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 <[email protected]>

* [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 <[email protected]>
Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Co-authored-by: David Huard <[email protected]>
Co-authored-by: Pascal Bourgault <[email protected]>

Author: 5 people

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 <https://github.com/charlesgauthier-udm>`_
 * Added the ability to generate regridding weights in parallel. By `Charles Gauthier <https://github.com/charlesgauthier-udm>`_
 * 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 <https://github.com/charlesgauthier-udm>`_

binder/environment.yml

@@ -7,6 +7,7 @@ dependencies:
   - dask
   - numpy
   - scipy
+  - shapely
   - matplotlib
   - cartopy
   - cf_xarray>=0.3.1

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

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",

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."""
 

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])