diff --git a/opendrift/readers/reader_shape.py b/opendrift/readers/reader_shape.py index 9fc3f1304..1f590cb0e 100644 --- a/opendrift/readers/reader_shape.py +++ b/opendrift/readers/reader_shape.py @@ -62,36 +62,61 @@ 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] @@ -99,11 +124,11 @@ def __init__(self, shapes, proj4_str = '+proj=lonlat +ellps=WGS84', invert=False 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): @@ -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 diff --git a/tests/readers/test_reader_shape.py b/tests/readers/test_reader_shape.py index 64f13e811..ab558f48b 100644 --- a/tests/readers/test_reader_shape.py +++ b/tests/readers/test_reader_shape.py @@ -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 @@ -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 + diff --git a/tests/test_data/shp/test_more_vars.cpg b/tests/test_data/shp/test_more_vars.cpg new file mode 100644 index 000000000..3ad133c04 --- /dev/null +++ b/tests/test_data/shp/test_more_vars.cpg @@ -0,0 +1 @@ +UTF-8 \ No newline at end of file diff --git a/tests/test_data/shp/test_more_vars.dbf b/tests/test_data/shp/test_more_vars.dbf new file mode 100644 index 000000000..6408526ce Binary files /dev/null and b/tests/test_data/shp/test_more_vars.dbf differ diff --git a/tests/test_data/shp/test_more_vars.prj b/tests/test_data/shp/test_more_vars.prj new file mode 100644 index 000000000..f45cbadf0 --- /dev/null +++ b/tests/test_data/shp/test_more_vars.prj @@ -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]] \ No newline at end of file diff --git a/tests/test_data/shp/test_more_vars.qix b/tests/test_data/shp/test_more_vars.qix new file mode 100644 index 000000000..4e60bdc77 Binary files /dev/null and b/tests/test_data/shp/test_more_vars.qix differ diff --git a/tests/test_data/shp/test_more_vars.qmd b/tests/test_data/shp/test_more_vars.qmd new file mode 100644 index 000000000..071726c19 --- /dev/null +++ b/tests/test_data/shp/test_more_vars.qmd @@ -0,0 +1,27 @@ + + + + + + + + + + + + + + + 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]] + +proj=longlat +datum=WGS84 +no_defs + 3452 + 4326 + EPSG:4326 + WGS 84 + longlat + EPSG:7030 + true + + + + diff --git a/tests/test_data/shp/test_more_vars.shp b/tests/test_data/shp/test_more_vars.shp new file mode 100644 index 000000000..240b419ed Binary files /dev/null and b/tests/test_data/shp/test_more_vars.shp differ diff --git a/tests/test_data/shp/test_more_vars.shx b/tests/test_data/shp/test_more_vars.shx new file mode 100644 index 000000000..2e6de532a Binary files /dev/null and b/tests/test_data/shp/test_more_vars.shx differ diff --git a/tests/test_data/shp/test_vars.cpg b/tests/test_data/shp/test_vars.cpg new file mode 100644 index 000000000..3ad133c04 --- /dev/null +++ b/tests/test_data/shp/test_vars.cpg @@ -0,0 +1 @@ +UTF-8 \ No newline at end of file diff --git a/tests/test_data/shp/test_vars.dbf b/tests/test_data/shp/test_vars.dbf new file mode 100644 index 000000000..6408526ce Binary files /dev/null and b/tests/test_data/shp/test_vars.dbf differ diff --git a/tests/test_data/shp/test_vars.prj b/tests/test_data/shp/test_vars.prj new file mode 100644 index 000000000..f45cbadf0 --- /dev/null +++ b/tests/test_data/shp/test_vars.prj @@ -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]] \ No newline at end of file diff --git a/tests/test_data/shp/test_vars.qix b/tests/test_data/shp/test_vars.qix new file mode 100644 index 000000000..5fc35560c Binary files /dev/null and b/tests/test_data/shp/test_vars.qix differ diff --git a/tests/test_data/shp/test_vars.qmd b/tests/test_data/shp/test_vars.qmd new file mode 100644 index 000000000..e5752fb61 --- /dev/null +++ b/tests/test_data/shp/test_vars.qmd @@ -0,0 +1,27 @@ + + + + + + + + + + + + + + + + +proj=longlat +datum=WGS84 +no_defs + 0 + 0 + + + + + false + + + + diff --git a/tests/test_data/shp/test_vars.shp b/tests/test_data/shp/test_vars.shp new file mode 100644 index 000000000..3fc831bf6 Binary files /dev/null and b/tests/test_data/shp/test_vars.shp differ diff --git a/tests/test_data/shp/test_vars.shx b/tests/test_data/shp/test_vars.shx new file mode 100644 index 000000000..953ac4bfc Binary files /dev/null and b/tests/test_data/shp/test_vars.shx differ diff --git a/tests/test_data/shp/test_vars_landmask.cpg b/tests/test_data/shp/test_vars_landmask.cpg new file mode 100644 index 000000000..3ad133c04 --- /dev/null +++ b/tests/test_data/shp/test_vars_landmask.cpg @@ -0,0 +1 @@ +UTF-8 \ No newline at end of file diff --git a/tests/test_data/shp/test_vars_landmask.dbf b/tests/test_data/shp/test_vars_landmask.dbf new file mode 100644 index 000000000..994d8c6a6 Binary files /dev/null and b/tests/test_data/shp/test_vars_landmask.dbf differ diff --git a/tests/test_data/shp/test_vars_landmask.prj b/tests/test_data/shp/test_vars_landmask.prj new file mode 100644 index 000000000..f45cbadf0 --- /dev/null +++ b/tests/test_data/shp/test_vars_landmask.prj @@ -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]] \ No newline at end of file diff --git a/tests/test_data/shp/test_vars_landmask.qix b/tests/test_data/shp/test_vars_landmask.qix new file mode 100644 index 000000000..f40e1b2da Binary files /dev/null and b/tests/test_data/shp/test_vars_landmask.qix differ diff --git a/tests/test_data/shp/test_vars_landmask.qmd b/tests/test_data/shp/test_vars_landmask.qmd new file mode 100644 index 000000000..071726c19 --- /dev/null +++ b/tests/test_data/shp/test_vars_landmask.qmd @@ -0,0 +1,27 @@ + + + + + + + + + + + + + + + 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]] + +proj=longlat +datum=WGS84 +no_defs + 3452 + 4326 + EPSG:4326 + WGS 84 + longlat + EPSG:7030 + true + + + + diff --git a/tests/test_data/shp/test_vars_landmask.shp b/tests/test_data/shp/test_vars_landmask.shp new file mode 100644 index 000000000..db561c8d4 Binary files /dev/null and b/tests/test_data/shp/test_vars_landmask.shp differ diff --git a/tests/test_data/shp/test_vars_landmask.shx b/tests/test_data/shp/test_vars_landmask.shx new file mode 100644 index 000000000..b44edd6ed Binary files /dev/null and b/tests/test_data/shp/test_vars_landmask.shx differ