From 2fe8ab853aae059abec115ce95e375e1169d4274 Mon Sep 17 00:00:00 2001 From: Bas des Tombe Date: Sun, 22 Dec 2024 08:47:28 +0100 Subject: [PATCH 01/11] Speed up gdf_to_bool_da and gdf_to_bool_ds functions. Added centroid checks and minimum area fraction filtering to docs. --- nlmod/dims/grid.py | 124 ++++++++++++++++++++++++++++++++++++++------- 1 file changed, 107 insertions(+), 17 deletions(-) diff --git a/nlmod/dims/grid.py b/nlmod/dims/grid.py index b85c9b66..33bf5be9 100644 --- a/nlmod/dims/grid.py +++ b/nlmod/dims/grid.py @@ -25,6 +25,7 @@ from scipy.interpolate import griddata from shapely.affinity import affine_transform from shapely.geometry import Point, Polygon +from shapely.ops import unary_union from tqdm import tqdm from .. import cache, util @@ -1610,10 +1611,23 @@ def aggregate_vector_per_cell(gdf, fields_methods, modelgrid=None): return celldata -def gdf_to_bool_da(gdf, ds, ix=None, buffer=0.0, **kwargs): - """Convert a GeoDataFrame with polygon geometries into a data array corresponding to - the modelgrid in which each cell is 1 (True) if one or more geometries are (partly) - in that cell. +def gdf_to_bool_da( + gdf, + ds, + ix=None, + buffer=0.0, + contains_centroid=False, + min_area_fraction=None, + **kwargs + ): + """Return True if grid cell is partly in polygons, False otherwise. + + This function returns True for grid cells that are partly in the polygons. If + contains_centroid is True, only cells are returned where the centroid is in the + polygon. If min_area_fraction is set, only cells are returned where the area of the + intersection is larger than min_area_fraction * cell_area. + + ix can be provided to speed up the process. If not provided it is computed from ds. Parameters ---------- @@ -1622,30 +1636,91 @@ def gdf_to_bool_da(gdf, ds, ix=None, buffer=0.0, **kwargs): ds : xr.Dataset Dataset with model data. ix : flopy.utils.GridIntersect, optional - If not provided it is computed from ds. + If not provided it is computed from ds. Speeds up the the function buffer : float, optional buffer around the geometries. The default is 0. + contains_centroid : bool, optional + if True, only store intersection result if cell centroid is + contained within intersection shape, only used if shape type is + "polygon" + min_area_fraction : float, optional + float defining minimum intersection area threshold, if intersection + area is smaller than min_frac_area * cell_area, do not store + intersection result. **kwargs : keyword arguments keyword arguments are passed to the intersect_*-methods. Returns ------- da : xr.DataArray - True if polygon is in cell, False otherwise. Grid dimensions according to ds. + True if polygon is in cell, False otherwise. """ - da = gdf_to_count_da(gdf, ds, ix=ix, buffer=buffer, **kwargs) > 0 + if isinstance(gdf, gpd.GeoDataFrame): + multipolygon = unary_union(gdf.geometry) + elif isinstance(gdf, shapely.geometry.base.BaseGeometry): + multipolygon = gdf + else: + msg = "gdf_to_bool_da() only support GeoDataFrame or shapely" + raise TypeError(msg) + + if buffer > 0.0: + multipolygon = multipolygon.buffer(buffer) + + if ix is None: + modelgrid = modelgrid_from_ds(ds) + ix = GridIntersect(modelgrid) + + if kwargs or contains_centroid or min_area_fraction is not None: + r = ix.intersect( + multipolygon, + contains_centroid=contains_centroid, + min_area_fraction=min_area_fraction, + **kwargs + ) + else: + r = ix.intersects(multipolygon) + + if ds.gridtype == "structured": + da = util.get_da_from_da_ds(ds, dims=("y", "x"), data=False) + elif ds.gridtype == "vertex": + da = util.get_da_from_da_ds(ds, dims=("icell2d",), data=False) + else: + msg = "gdf_to_bool_da() only support structured or vertex gridtypes" + raise ValueError(msg) + + if ds.gridtype == "structured": + ixs, iys = zip(*r["cellids"], strict=True) + da.values[ixs, iys] = True + elif ds.gridtype == "vertex": + da[r["cellids"].astype(int)] = True + return da -def gdf_to_bool_ds(gdf, ds, da_name, keep_coords=None, ix=None, buffer=0.0, **kwargs): - """Convert a GeoDataFrame with polygon geometries into a model dataset with a - data_array named 'da_name' in which each cell is 1 (True) if one or more geometries - are (partly) in that cell. +def gdf_to_bool_ds( + gdf, + ds, + da_name, + keep_coords=None, + ix=None, + buffer=0.0, + contains_centroid=False, + min_area_fraction=None, + **kwargs +): + """Return True if grid cell is partly in polygons, False otherwise. + + This function returns True for grid cells that are partly in the polygons. If + contains_centroid is True, only cells are returned where the centroid is in the + polygon. If min_area_fraction is set, only cells are returned where the area of the + intersection is larger than min_area_fraction * cell_area. + + ix can be provided to speed up the process. If not provided it is computed from ds. Parameters ---------- - gdf : geopandas.GeoDataFrame - polygon shapes with surface water. + gdf : geopandas.GeoDataFrame or shapely.geometry + shapes that will be rasterised. ds : xr.Dataset Dataset with model data. da_name : str @@ -1654,20 +1729,35 @@ def gdf_to_bool_ds(gdf, ds, da_name, keep_coords=None, ix=None, buffer=0.0, **kw the coordinates in ds the you want keep in your empty ds. If None all coordinates are kept from original ds. The default is None. ix : flopy.utils.GridIntersect, optional - If not provided it is computed from ds. + If not provided it is computed from ds. Speeds up the the function buffer : float, optional buffer around the geometries. The default is 0. + contains_centroid : bool, optional + if True, only store intersection result if cell centroid is + contained within intersection shape, only used if shape type is + "polygon" + min_area_fraction : float, optional + float defining minimum intersection area threshold, if intersection + area is smaller than min_frac_area * cell_area, do not store + intersection result. **kwargs : keyword arguments keyword arguments are passed to the intersect_*-methods. Returns ------- ds_out : xr.Dataset - Dataset with a single DataArray, this DataArray is 1 if polygon is in - cell, 0 otherwise. Grid dimensions according to ds and mfgrid. + True if polygon is in cell, False otherwise. """ ds_out = util.get_ds_empty(ds, keep_coords=keep_coords) - ds_out[da_name] = gdf_to_bool_da(gdf, ds, ix=ix, buffer=buffer, **kwargs) + ds_out[da_name] = gdf_to_bool_da( + gdf, + ds, + ix=ix, + buffer=buffer, + contains_centroid=contains_centroid, + min_area_fraction=min_area_fraction, + **kwargs, + ) return ds_out From 7585fb2bde8c1142a609e25ea7fd00b36e08504f Mon Sep 17 00:00:00 2001 From: Bas des Tombe Date: Sun, 22 Dec 2024 08:55:30 +0100 Subject: [PATCH 02/11] Remove trainling white spaces --- nlmod/dims/grid.py | 28 ++++++++++++++-------------- 1 file changed, 14 insertions(+), 14 deletions(-) diff --git a/nlmod/dims/grid.py b/nlmod/dims/grid.py index 33bf5be9..03c0c166 100644 --- a/nlmod/dims/grid.py +++ b/nlmod/dims/grid.py @@ -1612,12 +1612,12 @@ def aggregate_vector_per_cell(gdf, fields_methods, modelgrid=None): def gdf_to_bool_da( - gdf, - ds, - ix=None, - buffer=0.0, - contains_centroid=False, - min_area_fraction=None, + gdf, + ds, + ix=None, + buffer=0.0, + contains_centroid=False, + min_area_fraction=None, **kwargs ): """Return True if grid cell is partly in polygons, False otherwise. @@ -1698,14 +1698,14 @@ def gdf_to_bool_da( def gdf_to_bool_ds( - gdf, - ds, - da_name, - keep_coords=None, - ix=None, - buffer=0.0, - contains_centroid=False, - min_area_fraction=None, + gdf, + ds, + da_name, + keep_coords=None, + ix=None, + buffer=0.0, + contains_centroid=False, + min_area_fraction=None, **kwargs ): """Return True if grid cell is partly in polygons, False otherwise. From 198706ba2cd96dd6a7e2d2fbd5014e5f07e4210e Mon Sep 17 00:00:00 2001 From: Bas des Tombe Date: Sun, 22 Dec 2024 09:02:17 +0100 Subject: [PATCH 03/11] Remove more trailing spaces --- nlmod/dims/grid.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/nlmod/dims/grid.py b/nlmod/dims/grid.py index 03c0c166..e17639a4 100644 --- a/nlmod/dims/grid.py +++ b/nlmod/dims/grid.py @@ -1672,9 +1672,9 @@ def gdf_to_bool_da( if kwargs or contains_centroid or min_area_fraction is not None: r = ix.intersect( - multipolygon, - contains_centroid=contains_centroid, - min_area_fraction=min_area_fraction, + multipolygon, + contains_centroid=contains_centroid, + min_area_fraction=min_area_fraction, **kwargs ) else: From 60b57fd3c0e76b0de49ba31244a2fc602db69625 Mon Sep 17 00:00:00 2001 From: Bas des Tombe Date: Sun, 22 Dec 2024 10:54:47 +0100 Subject: [PATCH 04/11] GridIntersect initialization method in gdf_to_bool requires explicit method='vertex' --- nlmod/dims/grid.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/nlmod/dims/grid.py b/nlmod/dims/grid.py index e17639a4..f3708b3c 100644 --- a/nlmod/dims/grid.py +++ b/nlmod/dims/grid.py @@ -1668,7 +1668,7 @@ def gdf_to_bool_da( if ix is None: modelgrid = modelgrid_from_ds(ds) - ix = GridIntersect(modelgrid) + ix = GridIntersect(modelgrid, method="vertex") if kwargs or contains_centroid or min_area_fraction is not None: r = ix.intersect( From 415bd26673498cfc101cc1c79a9f12ae6d71c4d2 Mon Sep 17 00:00:00 2001 From: Bas des Tombe Date: Sun, 22 Dec 2024 11:03:01 +0100 Subject: [PATCH 05/11] Enhance buffer parameter documentation in gdf_to_bool_da and gdf_to_bool_ds functions to clarify its effect on geometry dilation and erosion. --- nlmod/dims/grid.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/nlmod/dims/grid.py b/nlmod/dims/grid.py index f3708b3c..2390a0ff 100644 --- a/nlmod/dims/grid.py +++ b/nlmod/dims/grid.py @@ -1638,7 +1638,8 @@ def gdf_to_bool_da( ix : flopy.utils.GridIntersect, optional If not provided it is computed from ds. Speeds up the the function buffer : float, optional - buffer around the geometries. The default is 0. + The distance to buffer around the geometries in meters. A positive distance + produces a dilation, a negative distance an erosion. The default is 0. contains_centroid : bool, optional if True, only store intersection result if cell centroid is contained within intersection shape, only used if shape type is @@ -1663,7 +1664,7 @@ def gdf_to_bool_da( msg = "gdf_to_bool_da() only support GeoDataFrame or shapely" raise TypeError(msg) - if buffer > 0.0: + if buffer != 0.0: multipolygon = multipolygon.buffer(buffer) if ix is None: @@ -1731,7 +1732,8 @@ def gdf_to_bool_ds( ix : flopy.utils.GridIntersect, optional If not provided it is computed from ds. Speeds up the the function buffer : float, optional - buffer around the geometries. The default is 0. + The distance to buffer around the geometries in meters. A positive distance + produces a dilation, a negative distance an erosion. The default is 0. contains_centroid : bool, optional if True, only store intersection result if cell centroid is contained within intersection shape, only used if shape type is From ebeba06e34ad9efee2fbd448d6fa478098785588 Mon Sep 17 00:00:00 2001 From: Bas des Tombe Date: Mon, 23 Dec 2024 08:12:21 +0100 Subject: [PATCH 06/11] empty geometry support in gdf_to_bool_da function --- nlmod/dims/grid.py | 20 ++++++++++++-------- 1 file changed, 12 insertions(+), 8 deletions(-) diff --git a/nlmod/dims/grid.py b/nlmod/dims/grid.py index 2390a0ff..dec3694b 100644 --- a/nlmod/dims/grid.py +++ b/nlmod/dims/grid.py @@ -1656,6 +1656,18 @@ def gdf_to_bool_da( da : xr.DataArray True if polygon is in cell, False otherwise. """ + if ds.gridtype == "structured": + da = util.get_da_from_da_ds(ds, dims=("y", "x"), data=False) + elif ds.gridtype == "vertex": + da = util.get_da_from_da_ds(ds, dims=("icell2d",), data=False) + else: + msg = "gdf_to_bool_da() only support structured or vertex gridtypes" + raise ValueError(msg) + + if len(gdf) == 0: + logger.warning("gdf passed to gdf_to_bool_da() is empty") + return da + if isinstance(gdf, gpd.GeoDataFrame): multipolygon = unary_union(gdf.geometry) elif isinstance(gdf, shapely.geometry.base.BaseGeometry): @@ -1681,14 +1693,6 @@ def gdf_to_bool_da( else: r = ix.intersects(multipolygon) - if ds.gridtype == "structured": - da = util.get_da_from_da_ds(ds, dims=("y", "x"), data=False) - elif ds.gridtype == "vertex": - da = util.get_da_from_da_ds(ds, dims=("icell2d",), data=False) - else: - msg = "gdf_to_bool_da() only support structured or vertex gridtypes" - raise ValueError(msg) - if ds.gridtype == "structured": ixs, iys = zip(*r["cellids"], strict=True) da.values[ixs, iys] = True From 7d0ae652bb1e2d7c5c1a77286e3b06b8ed78c736 Mon Sep 17 00:00:00 2001 From: Bas des Tombe Date: Mon, 23 Dec 2024 08:17:32 +0100 Subject: [PATCH 07/11] Removed whitespace --- nlmod/dims/grid.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/nlmod/dims/grid.py b/nlmod/dims/grid.py index dec3694b..e88f7150 100644 --- a/nlmod/dims/grid.py +++ b/nlmod/dims/grid.py @@ -1663,7 +1663,7 @@ def gdf_to_bool_da( else: msg = "gdf_to_bool_da() only support structured or vertex gridtypes" raise ValueError(msg) - + if len(gdf) == 0: logger.warning("gdf passed to gdf_to_bool_da() is empty") return da From ad2c9e6d2fc0cd22765ed1dada295a25d0d734d5 Mon Sep 17 00:00:00 2001 From: Bas des Tombe Date: Mon, 23 Dec 2024 09:56:32 +0100 Subject: [PATCH 08/11] Removed logging empty gdf and no union for singleton --- nlmod/dims/grid.py | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/nlmod/dims/grid.py b/nlmod/dims/grid.py index e88f7150..f4bb25d8 100644 --- a/nlmod/dims/grid.py +++ b/nlmod/dims/grid.py @@ -1664,12 +1664,13 @@ def gdf_to_bool_da( msg = "gdf_to_bool_da() only support structured or vertex gridtypes" raise ValueError(msg) - if len(gdf) == 0: - logger.warning("gdf passed to gdf_to_bool_da() is empty") - return da - if isinstance(gdf, gpd.GeoDataFrame): - multipolygon = unary_union(gdf.geometry) + if len(gdf) == 0: + return da + elif len(gdf) == 1: + multipolygon = gdf.geometry.values[0] + else: + multipolygon = unary_union(gdf.geometry) elif isinstance(gdf, shapely.geometry.base.BaseGeometry): multipolygon = gdf else: From e8d47cccdba7640c485c21d769125e260a05d4ad Mon Sep 17 00:00:00 2001 From: Bas des Tombe Date: Mon, 23 Dec 2024 11:35:18 +0100 Subject: [PATCH 09/11] Rotate the polygon instead of the modelgrid --- nlmod/dims/grid.py | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/nlmod/dims/grid.py b/nlmod/dims/grid.py index f4bb25d8..31acad10 100644 --- a/nlmod/dims/grid.py +++ b/nlmod/dims/grid.py @@ -1680,9 +1680,17 @@ def gdf_to_bool_da( if buffer != 0.0: multipolygon = multipolygon.buffer(buffer) + # Rotate the polygon instead of the modelgrid if ix is None: - modelgrid = modelgrid_from_ds(ds) + modelgrid = modelgrid_from_ds(ds, rotated=False) ix = GridIntersect(modelgrid, method="vertex") + + grid_rotation = ds.attrs.get("angrot", 0.0) + ix_rotation = ix.mfgrid.angrot + + if grid_rotation != 0.0 and ix_rotation == 0.0: + affine = get_affine_mod_to_world(ds) + multipolygon = affine_transform_gdf(multipolygon, affine) if kwargs or contains_centroid or min_area_fraction is not None: r = ix.intersect( From e7994d3c55809735e316492037388e4a29add397 Mon Sep 17 00:00:00 2001 From: Bas des Tombe Date: Mon, 23 Dec 2024 11:49:03 +0100 Subject: [PATCH 10/11] Removing whitespace --- nlmod/dims/grid.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/nlmod/dims/grid.py b/nlmod/dims/grid.py index 31acad10..75f90c9c 100644 --- a/nlmod/dims/grid.py +++ b/nlmod/dims/grid.py @@ -1684,7 +1684,7 @@ def gdf_to_bool_da( if ix is None: modelgrid = modelgrid_from_ds(ds, rotated=False) ix = GridIntersect(modelgrid, method="vertex") - + grid_rotation = ds.attrs.get("angrot", 0.0) ix_rotation = ix.mfgrid.angrot From 1a02343be3d54cb86a8f8e60a4a6dc761ffd7237 Mon Sep 17 00:00:00 2001 From: Bas des Tombe Date: Mon, 23 Dec 2024 14:02:15 +0100 Subject: [PATCH 11/11] Rotate grid in the correct direction --- nlmod/dims/grid.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/nlmod/dims/grid.py b/nlmod/dims/grid.py index 75f90c9c..52d34d03 100644 --- a/nlmod/dims/grid.py +++ b/nlmod/dims/grid.py @@ -1649,7 +1649,7 @@ def gdf_to_bool_da( area is smaller than min_frac_area * cell_area, do not store intersection result. **kwargs : keyword arguments - keyword arguments are passed to the intersect_*-methods. + keyword arguments are passed to the intersect-method. Returns ------- @@ -1689,8 +1689,8 @@ def gdf_to_bool_da( ix_rotation = ix.mfgrid.angrot if grid_rotation != 0.0 and ix_rotation == 0.0: - affine = get_affine_mod_to_world(ds) - multipolygon = affine_transform_gdf(multipolygon, affine) + affine = get_affine_world_to_mod(ds).to_shapely() + multipolygon = affine_transform(multipolygon, affine) if kwargs or contains_centroid or min_area_fraction is not None: r = ix.intersect( @@ -1702,10 +1702,10 @@ def gdf_to_bool_da( else: r = ix.intersects(multipolygon) - if ds.gridtype == "structured": + if r.size > 0 and ds.gridtype == "structured": ixs, iys = zip(*r["cellids"], strict=True) da.values[ixs, iys] = True - elif ds.gridtype == "vertex": + elif r.size > 0 and ds.gridtype == "vertex": da[r["cellids"].astype(int)] = True return da