Skip to content

Shape reader: multiple variables alowed #1460

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 1 commit into
base: master
Choose a base branch
from
+185 −24
Open

Shape reader: multiple variables alowed #1460

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
71 changes: 49 additions & 22 deletions opendrift/readers/reader_shape.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -62,48 +62,73 @@ def from_shpfiles(shpfiles, proj4_str = '+proj=lonlat +ellps=WGS84', invert=Fals
shpfiles = shpfiles
else:
shpfiles = [ shpfiles ]

# reading shapefiles
shp_iters = []
for shp in shpfiles:
logger.debug("Reading shapefile: %s" % shp)
from cartopy import io
reader = io.shapereader.Reader(shp)
shp_iters.append(reader.geometries())

shp_iters.append(reader.records())
return Reader(itertools.chain(*shp_iters), proj4_str, invert=invert)

def __init__(self, shapes, proj4_str = '+proj=lonlat +ellps=WGS84', invert=False):
def _ingest_landmask_only(self):
assert len(self.polys) > 0, "no geometries loaded"
logger.info("Pre-processing %d geometries" % len(self.polys))
self.land = {
self.variables[0]: shapely.ops.unary_union(self.polys)
}

def _ingest_multivariable(self):
self.land = {}
for var in self.variables:
features = filter(lambda x: x.attributes['variable'] == var,
self.records)
self.land[var] = shapely.ops.unary_union([x.geometry
for x in features])

def __init__(self, shapes, proj4_str = '+proj=lonlat +ellps=WGS84', invert=False):
self.invert = invert # True if polygons are lakes and not land areas
self.proj4 = proj4_str
self.crs = pyproj.CRS(self.proj4)

super (Reader, self).__init__ ()

# Depth
self.z = None

# loading and pre-processing shapes
self.polys = list(shapes) # reads geometries if Shape Readers
assert len(self.polys) > 0, "no geometries loaded"

logger.info("Pre-processing %d geometries" % len(self.polys))
self.land = shapely.ops.unary_union(self.polys)

self.xmin, self.ymin, self.xmax, self.ymax = self.land.bounds
try:
# Test if shapes are records or geometries
self.records = list(shapes)
self.polys = [x.geometry for x in self.records]
except AttributeError:
# This is a classic land mask file
self.polys = list(shapes) # reads geometries if Shape Readers
self.records = None
self._ingest_landmask_only()
else:
try:
# When records are passed, check for the `variable` attribute
add_vars = set([x.attributes['variable'] for x in self.records])
except KeyError:
# It's a plain land mask file without attributes
self._ingest_landmask_only()
else:
# It's a multi-variable shapefile
self.variables = list(add_vars)
self._ingest_multivariable()
# Regardless of the number of variables, the extent is always the sum
# of the polygons
whole = shapely.ops.unary_union(self.polys)
self.xmin, self.ymin, self.xmax, self.ymax = whole.bounds
self.xmin, self.ymin = self.lonlat2xy(self.xmin, self.ymin)
self.xmax, self.ymax = self.lonlat2xy(self.xmax, self.ymax)
self.xmin = self.xmin[0]
self.xmax = self.xmax[0]
self.ymin = self.ymin[0]
self.ymax = self.ymax[0]

def __on_land__(self, x, y):
def __on_land__(self, x, y, var):
if self.invert is False:
return shapely.vectorized.contains(self.land, x, y)
return shapely.vectorized.contains(self.land[var], x, y)
else: # Inverse if polygons are lakes and not land areas
return 1 - shapely.vectorized.contains(self.land, x, y)
return 1 - shapely.vectorized.contains(self.land[var], x, y)

def get_variables(self, requestedVariables, time = None,
x = None, y = None, z = None):
Expand All @@ -121,8 +146,10 @@ def get_variables(self, requestedVariables, time = None,
Binary mask of point x, y on land.

"""

self.check_arguments(requestedVariables, time, x, y, z)
return { 'x' : x, 'y' : y, 'land_binary_mask': self.__on_land__(x,y) }

requestedVariables = \
self.check_arguments(requestedVariables, time, x, y, z)[0]
variables = { 'x' : x, 'y' : y}
for var in requestedVariables:
variables[var] = self.__on_land__(x, y, var)
return variables

51 changes: 49 additions & 2 deletions tests/readers/test_reader_shape.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -27,8 +27,8 @@ def test_on_land():
name='admin_0_countries')
r = reader_shape.Reader.from_shpfiles(shpfilename)

assert r.__on_land__(np.array([10]), np.array([60])) == [True]
assert r.__on_land__(np.array([5]), np.array([60])) == [False]
assert r.__on_land__(np.array([10]), np.array([60]), r.variables[0]) == [True]
assert r.__on_land__(np.array([5]), np.array([60]), r.variables[0]) == [False]


@need_ne_shapes
Expand Down Expand Up @@ -58,3 +58,50 @@ def test_global_array(test_data):
print(oc.env.readers)
assert len(
oc.env.readers) == 2 # make sure opendrift doesn't add default basemap

def test_get_shape_variable():
oc = OceanDrift(loglevel=30)
shpfilename = [
oc.test_data_folder() + 'shp/test_vars.shp',
oc.test_data_folder() + 'shp/test_more_vars.shp',
]
reader = reader_shape.Reader.from_shpfiles(shpfilename)
test = reader.get_variables(['var0', 'var1'],
x=14.267009, y=68.528457)
assert test['var0'] == False and test['var1'] == True, \
"Failed on reading test_vars.shp"
test = reader.get_variables(['var0', 'var1'],
x=13.979286, y=68.164405)
assert test['var0'] == True and test['var1'] == False, \
"Failed on reading test_more_vars.shp"

def test_global_var_array():
oc = OceanDrift(loglevel=30)
shpfilename = [
oc.test_data_folder() + 'shp/test_vars.shp',
oc.test_data_folder() + 'shp/test_vars_landmask.shp',
oc.test_data_folder() + 'shp/test_more_vars.shp',
]
reader_landmask = reader_shape.Reader.from_shpfiles(shpfilename)

reader_nordic = reader_ROMS_native.Reader(
oc.test_data_folder() + \
'2Feb2016_Nordic_sigma_3d/Nordic-4km_SLEVELS_avg_00_subset2Feb2016.nc')

lon = np.array([15., 5.])
lat = np.array([65.6, 65.6])

# global
oc = OceanDrift(loglevel=00)
oc.set_config('general:use_auto_landmask', False)
oc.add_reader([reader_nordic, reader_landmask])
oc.env.finalize()
en, en_prof, missing = oc.env.get_environment(['land_binary_mask'],
reader_nordic.start_time, lon,
lat, np.array([0, 0]), None)

np.testing.assert_array_equal(en.land_binary_mask, np.array([True, False]))
print(oc.env.readers)
assert len(
oc.env.readers) == 2 # make sure opendrift doesn't add default basemap

1 change: 1 addition & 0 deletions tests/test_data/shp/test_more_vars.cpg
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1 @@
UTF-8
Binary file added tests/test_data/shp/test_more_vars.dbf
View file
Open in desktop
Binary file not shown.
1 change: 1 addition & 0 deletions tests/test_data/shp/test_more_vars.prj
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1 @@
GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]
Binary file added tests/test_data/shp/test_more_vars.qix
View file
Open in desktop
Binary file not shown.
27 changes: 27 additions & 0 deletions tests/test_data/shp/test_more_vars.qmd
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,27 @@
<!DOCTYPE qgis PUBLIC 'http://mrcc.com/qgis.dtd' 'SYSTEM'>
<qgis version="3.34.11-Prizren">
<identifier></identifier>
<parentidentifier></parentidentifier>
<language></language>
<type></type>
<title></title>
<abstract></abstract>
<links/>
<dates/>
<fees></fees>
<encoding></encoding>
<crs>
<spatialrefsys nativeFormat="Wkt">
<wkt>GEOGCRS["WGS 84",ENSEMBLE["World Geodetic System 1984 ensemble",MEMBER["World Geodetic System 1984 (Transit)"],MEMBER["World Geodetic System 1984 (G730)"],MEMBER["World Geodetic System 1984 (G873)"],MEMBER["World Geodetic System 1984 (G1150)"],MEMBER["World Geodetic System 1984 (G1674)"],MEMBER["World Geodetic System 1984 (G1762)"],MEMBER["World Geodetic System 1984 (G2139)"],ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]],ENSEMBLEACCURACY[2.0]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],CS[ellipsoidal,2],AXIS["geodetic latitude (Lat)",north,ORDER[1],ANGLEUNIT["degree",0.0174532925199433]],AXIS["geodetic longitude (Lon)",east,ORDER[2],ANGLEUNIT["degree",0.0174532925199433]],USAGE[SCOPE["Horizontal component of 3D system."],AREA["World."],BBOX[-90,-180,90,180]],ID["EPSG",4326]]</wkt>
<proj4>+proj=longlat +datum=WGS84 +no_defs</proj4>
<srsid>3452</srsid>
<srid>4326</srid>
<authid>EPSG:4326</authid>
<description>WGS 84</description>
<projectionacronym>longlat</projectionacronym>
<ellipsoidacronym>EPSG:7030</ellipsoidacronym>
<geographicflag>true</geographicflag>
</spatialrefsys>
</crs>
<extent/>
</qgis>
Binary file added tests/test_data/shp/test_more_vars.shp
View file
Open in desktop
Binary file not shown.
Binary file added tests/test_data/shp/test_more_vars.shx
View file
Open in desktop
Binary file not shown.
1 change: 1 addition & 0 deletions tests/test_data/shp/test_vars.cpg
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1 @@
UTF-8
Binary file added tests/test_data/shp/test_vars.dbf
View file
Open in desktop
Binary file not shown.
1 change: 1 addition & 0 deletions tests/test_data/shp/test_vars.prj
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1 @@
GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]
Binary file added tests/test_data/shp/test_vars.qix
View file
Open in desktop
Binary file not shown.
27 changes: 27 additions & 0 deletions tests/test_data/shp/test_vars.qmd
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,27 @@
<!DOCTYPE qgis PUBLIC 'http://mrcc.com/qgis.dtd' 'SYSTEM'>
<qgis version="3.34.11-Prizren">
<identifier></identifier>
<parentidentifier></parentidentifier>
<language></language>
<type></type>
<title></title>
<abstract></abstract>
<links/>
<dates/>
<fees></fees>
<encoding></encoding>
<crs>
<spatialrefsys nativeFormat="Wkt">
<wkt></wkt>
<proj4>+proj=longlat +datum=WGS84 +no_defs</proj4>
<srsid>0</srsid>
<srid>0</srid>
<authid></authid>
<description></description>
<projectionacronym></projectionacronym>
<ellipsoidacronym></ellipsoidacronym>
<geographicflag>false</geographicflag>
</spatialrefsys>
</crs>
<extent/>
</qgis>
Binary file added tests/test_data/shp/test_vars.shp
View file
Open in desktop
Binary file not shown.
Binary file added tests/test_data/shp/test_vars.shx
View file
Open in desktop
Binary file not shown.
1 change: 1 addition & 0 deletions tests/test_data/shp/test_vars_landmask.cpg
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1 @@
UTF-8
Binary file added tests/test_data/shp/test_vars_landmask.dbf
View file
Open in desktop
Binary file not shown.
1 change: 1 addition & 0 deletions tests/test_data/shp/test_vars_landmask.prj
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1 @@
GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]
Binary file added tests/test_data/shp/test_vars_landmask.qix
View file
Open in desktop
Binary file not shown.
27 changes: 27 additions & 0 deletions tests/test_data/shp/test_vars_landmask.qmd
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,27 @@
<!DOCTYPE qgis PUBLIC 'http://mrcc.com/qgis.dtd' 'SYSTEM'>
<qgis version="3.34.11-Prizren">
<identifier></identifier>
<parentidentifier></parentidentifier>
<language></language>
<type></type>
<title></title>
<abstract></abstract>
<links/>
<dates/>
<fees></fees>
<encoding></encoding>
<crs>
<spatialrefsys nativeFormat="Wkt">
<wkt>GEOGCRS["WGS 84",ENSEMBLE["World Geodetic System 1984 ensemble",MEMBER["World Geodetic System 1984 (Transit)"],MEMBER["World Geodetic System 1984 (G730)"],MEMBER["World Geodetic System 1984 (G873)"],MEMBER["World Geodetic System 1984 (G1150)"],MEMBER["World Geodetic System 1984 (G1674)"],MEMBER["World Geodetic System 1984 (G1762)"],MEMBER["World Geodetic System 1984 (G2139)"],ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]],ENSEMBLEACCURACY[2.0]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],CS[ellipsoidal,2],AXIS["geodetic latitude (Lat)",north,ORDER[1],ANGLEUNIT["degree",0.0174532925199433]],AXIS["geodetic longitude (Lon)",east,ORDER[2],ANGLEUNIT["degree",0.0174532925199433]],USAGE[SCOPE["Horizontal component of 3D system."],AREA["World."],BBOX[-90,-180,90,180]],ID["EPSG",4326]]</wkt>
<proj4>+proj=longlat +datum=WGS84 +no_defs</proj4>
<srsid>3452</srsid>
<srid>4326</srid>
<authid>EPSG:4326</authid>
<description>WGS 84</description>
<projectionacronym>longlat</projectionacronym>
<ellipsoidacronym>EPSG:7030</ellipsoidacronym>
<geographicflag>true</geographicflag>
</spatialrefsys>
</crs>
<extent/>
</qgis>
Binary file added tests/test_data/shp/test_vars_landmask.shp
View file
Open in desktop
Binary file not shown.
Binary file added tests/test_data/shp/test_vars_landmask.shx
View file
Open in desktop
Binary file not shown.