From 65500843d487dc53e3785b66c27b236b420e43ec Mon Sep 17 00:00:00 2001 From: Charles Gauthier <85585494+charlesgauthier-udm@users.noreply.github.com> Date: Fri, 1 Sep 2023 14:42:05 -0400 Subject: [PATCH] Warn of SpatialAverager error over large region and densify polygons (#293) * Added warning for error over large regions and a function to densify polygons to prevent errors * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * Remove unnecessary import * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * Remove unnecessary import again, made a mistake last time * Moved densify_poly func to util.py and simplified it using shapely.segmentize. Simplified the length check for the warning in the SpatialAverager * Updated requirements to shapely>=2 so that shapely.segmentize can be used * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * Remove unnecessary imports since change to shapely.segmentize * add note about segmentize in the nb * Added details in CHANGES.rst * Moved the poly length check to its own function and modified the warning * Added test to verify that a warning is raised for large polys * Removed densify_polys function * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * Removed unnecessary test * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * Suggested changes to xesmf/frontend.py Co-authored-by: Pascal Bourgault * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * Update CHANGES.rst * Removed shapely>=2 dependencies --------- Co-authored-by: charlesgauthier-udm Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com> Co-authored-by: David Huard Co-authored-by: Pascal Bourgault --- CHANGES.rst | 1 + binder/environment.yml | 1 + ci/environment-upstream-dev.yml | 1 + doc/notebooks/Spatial_Averaging.ipynb | 6 ++++++ xesmf/frontend.py | 21 +++++++++++++++++++++ xesmf/tests/test_frontend.py | 7 +++++++ 6 files changed, 37 insertions(+) diff --git a/CHANGES.rst b/CHANGES.rst index 76aca736..68fea715 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -6,6 +6,7 @@ What's new New features ~~~~~~~~~~~~ +* Added a check in SpatialAverager that warns user if they are using polygons with long segments that could cause errors (:pull:`293`). By `Charles Gauthier `_ * Added the ability to generate regridding weights in parallel. By `Charles Gauthier `_ * Added the ability to chunk over horizontal/core dimensions. Added a `output_chunks` argument to the `Regridder` that allows the user to specify the chunking of the output data. By `Charles Gauthier `_ diff --git a/binder/environment.yml b/binder/environment.yml index 18e3a674..bcb6ee6e 100644 --- a/binder/environment.yml +++ b/binder/environment.yml @@ -7,6 +7,7 @@ dependencies: - dask - numpy - scipy + - shapely - matplotlib - cartopy - cf_xarray>=0.3.1 diff --git a/ci/environment-upstream-dev.yml b/ci/environment-upstream-dev.yml index 8f1aeabe..659aee2a 100644 --- a/ci/environment-upstream-dev.yml +++ b/ci/environment-upstream-dev.yml @@ -13,6 +13,7 @@ dependencies: - pydap - pytest - pytest-cov + - shapely - sparse>=0.8.0 - pip: - git+https://github.com/pydata/xarray.git diff --git a/doc/notebooks/Spatial_Averaging.ipynb b/doc/notebooks/Spatial_Averaging.ipynb index dc67110e..ad094cd6 100644 --- a/doc/notebooks/Spatial_Averaging.ipynb +++ b/doc/notebooks/Spatial_Averaging.ipynb @@ -28,6 +28,12 @@ "[regionmask](https://regionmask.readthedocs.io/) or\n", "[clisops](https://clisops.readthedocs.io).\n", "\n", + "Also, note that low resolution polygons such as large boxes might get distorted\n", + "when mapped on a sphere. Make sure polygon segments are at sufficiently high\n", + "resolution, on the order of 1°. The `shapely` package (v2) has a\n", + "[`segmentize` function](https://shapely.readthedocs.io/en/latest/reference/shapely.segmentize.html)\n", + "to add vertices to polygon segments.\n", + "\n", "The following example shows just how simple it is to compute the average over\n", "different countries. The notebook used `geopandas`, a simple and efficient\n", "container for geometries, and `descartes` for plotting maps. Make sure both\n", diff --git a/xesmf/frontend.py b/xesmf/frontend.py index 4b086583..1d14f796 100644 --- a/xesmf/frontend.py +++ b/xesmf/frontend.py @@ -8,6 +8,7 @@ import numpy as np import sparse as sps import xarray as xr +from shapely import LineString from xarray import DataArray, Dataset from .backend import Grid, LocStream, Mesh, add_corner, esmf_regrid_build, esmf_regrid_finalize @@ -1180,6 +1181,9 @@ def __init__( self._lon_out_name = 'lon' self._lat_out_name = 'lat' + # Check length of polys segments + self._check_polys_length(polys) + poly_centers = [poly.centroid.xy for poly in polys] self._lon_out = np.asarray([c[0][0] for c in poly_centers]) self._lat_out = np.asarray([c[1][0] for c in poly_centers]) @@ -1202,6 +1206,23 @@ def __init__( unmapped_to_nan=False, ) + @staticmethod + def _check_polys_length(polys, threshold=1): + # Check length of polys segments, issue warning if too long + check_polys, check_holes, _, _ = split_polygons_and_holes(polys) + check_polys.extend(check_holes) + poly_segments = [] + for check_poly in check_polys: + b = check_poly.boundary.coords + # Length of each segment + poly_segments.extend([LineString(b[k : k + 2]).length for k in range(len(b) - 1)]) + if np.any(np.array(poly_segments) > threshold): + warnings.warn( + f'`polys` contains large (> {threshold}°) segments. This could lead to errors over large regions. For a more accurate average, segmentize (densify) your shapes with `shapely.segmentize(polys, {threshold})`', + UserWarning, + stacklevel=2, + ) + def _compute_weights_and_area(self, mesh_out): """Return the weights and the area of the destination mesh cells.""" diff --git a/xesmf/tests/test_frontend.py b/xesmf/tests/test_frontend.py index 9937177a..1eff1371 100644 --- a/xesmf/tests/test_frontend.py +++ b/xesmf/tests/test_frontend.py @@ -952,3 +952,10 @@ def test_spatial_averager_mask(): savg = xe.SpatialAverager(dsm, [poly], geom_dim_name='my_geom') out = savg(dsm.abc) assert_allclose(out, 2, rtol=1e-3) + + +def test_densify_polys(): + # Check that using a large poly raises a warning + poly = Polygon([(-80, -40), (80, -40), (80, 40), (-80, 40)]) # Large poly + with pytest.warns(UserWarning): + xe.SpatialAverager(ds_in, [poly])