Repository: Ouranosinc/xclim

Issue:
Subset with individual features
* xclim version: 0.13.0

It would be great to have an option for us to explode the features of a shapefile and subset spatially according to their boundaries. This has been implemented in a few libraries but it wouldn't take much work to modify our current `subset_shape` and `create_mask` to have this behaviour.

Discussion:
Zeitsperre: I think you may have forgotten to add the multi-feature json to the PR. Once that's in, I'll take a look. Thanks!
coveralls: ## Pull Request Test Coverage Report for [Build 1480](https://coveralls.io/builds/28099692)

* **15** of **15**   **(100.0%)**  changed or added relevant lines in **1** file are covered.
* No unchanged relevant lines lost coverage.
* Overall coverage increased (+**0.02%**) to **90.836%**

---



|  Totals | [![Coverage Status](https://coveralls.io/builds/28099692/badge)](https://coveralls.io/builds/28099692) |
| :-- | --: |
| Change from base [Build 1473](https://coveralls.io/builds/28036063): |  0.02% |
| Covered Lines: | 1685 |
| Relevant Lines: | 1855 |

---
##### 💛  - [Coveralls](https://coveralls.io)


Patch:
diff --git a/HISTORY.rst b/HISTORY.rst
index 50dd6f92..ec7dbd31 100644
--- a/HISTORY.rst
+++ b/HISTORY.rst
@@ -2,6 +2,11 @@
 History
 =======
 
+0.14.x
+------
+* Fix in `subset.wrap_lons_and_split_at_greenwich` to preserve multi-region dataframes.
+
+
 0.13.x (2020-01-10)
 -------------------
 * Documentation improvements: list of indicators, RTD theme, notebook example.
diff --git a/setup.cfg b/setup.cfg
index 1bb16fe4..df1ca2e7 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -1,5 +1,5 @@
 [bumpversion]
-current_version = 0.13.0
+current_version = 0.13.1-beta
 commit = False
 tag = True
 parse = (?P<major>\d+)\.(?P<minor>\d+).(?P<patch>\d+)(\-(?P<release>[a-z]+))?
diff --git a/setup.py b/setup.py
index 38f7694b..b3d84518 100644
--- a/setup.py
+++ b/setup.py
@@ -12,7 +12,7 @@ URL = "https://github.com/Ouranosinc/xclim"
 AUTHOR = "Travis Logan"
 AUTHOR_EMAIL = "logan.travis@ouranos.ca"
 REQUIRES_PYTHON = ">=3.6.0"
-VERSION = "0.13.0"
+VERSION = "0.13.1-beta"
 LICENSE = "Apache Software License 2.0"
 
 with open("README.rst") as readme_file:
diff --git a/xclim/__init__.py b/xclim/__init__.py
index f8a0186b..f6e7e099 100644
--- a/xclim/__init__.py
+++ b/xclim/__init__.py
@@ -12,7 +12,7 @@ from xclim import seaIce
 
 __author__ = """Travis Logan"""
 __email__ = "logan.travis@ouranos.ca"
-__version__ = "0.13.0"
+__version__ = "0.13.1-beta"
 
 
 def build_module(name, objs, doc="", source=None, mode="ignore"):
diff --git a/xclim/subset.py b/xclim/subset.py
index ff58380e..a61b36f8 100644
--- a/xclim/subset.py
+++ b/xclim/subset.py
@@ -193,37 +193,38 @@ def wrap_lons_and_split_at_greenwich(func):
                     stacklevel=4,
                 )
             split_flag = False
-            if (poly.geometry.total_bounds[0] < 0) and (
-                poly.geometry.total_bounds[2] > 0
-            ):
-                split_flag = True
-                warnings.warn(
-                    "Geometry crosses the Greenwich Meridian. Proceeding to split polygon at Greenwich."
-                    " This feature is experimental. Output might not be accurate.",
-                    UserWarning,
-                    stacklevel=4,
-                )
-
-                # Create a meridian line at Greenwich, split polygons at this line and erase a buffer line
-                union = Polygon(cascaded_union(poly.geometry))
-                meridian = LineString([Point(0, 90), Point(0, -90)])
-                buffered = meridian.buffer(0.000000001)
-                split_polygons = split(union, meridian)
-                # TODO: This doesn't seem to be thread safe in Travis CI on macOS. Merits testing with a local machine.
-                buffered_split_polygons = [
-                    feat for feat in split_polygons.difference(buffered)
-                ]
-
-                # Load split features into a new GeoDataFrame with WGS84 CRS
-                split_gdf = gpd.GeoDataFrame(
-                    list(range(len(buffered_split_polygons))),
-                    geometry=buffered_split_polygons,
-                    crs={"init": "epsg:4326"},
-                )
-                # split_gdf.crs = CRS.from_epsg(4326)
-                split_gdf.columns = ["index", "geometry"]
-
-                poly = split_gdf
+            for (index, feature) in poly.iterrows():
+                if (feature.geometry.bounds[0] < 0) and (
+                    feature.geometry.bounds[2] > 0
+                ):
+                    split_flag = True
+                    warnings.warn(
+                        "Geometry crosses the Greenwich Meridian. Proceeding to split polygon at Greenwich."
+                        " This feature is experimental. Output might not be accurate.",
+                        UserWarning,
+                        stacklevel=4,
+                    )
+
+                    # Create a meridian line at Greenwich, split polygons at this line and erase a buffer line
+                    union = Polygon(cascaded_union(feature.geometry))
+                    meridian = LineString([Point(0, 90), Point(0, -90)])
+                    buffered = meridian.buffer(0.000000001)
+                    split_polygons = split(union, meridian)
+                    # TODO: This doesn't seem to be thread safe in Travis CI on macOS. Merits testing with a local machine.
+                    buffered_split_polygons = [
+                        feat for feat in split_polygons.difference(buffered)
+                    ]
+
+                    # Cannot assign iterable with `at` (pydata/pandas#26333) so a small hack:
+                    # Load split features into a new GeoDataFrame with WGS84 CRS
+                    split_gdf = gpd.GeoDataFrame(
+                        geometry=[cascaded_union(buffered_split_polygons)],
+                        crs={"init": "epsg:4326"},
+                    )
+                    poly.at[[index], "geometry"] = split_gdf.geometry.values
+                    # split_gdf.columns = ["index", "geometry"]
+
+                    # feature = split_gdf
 
             # Reproject features in WGS84 CSR to use 0 to 360 as longitudinal values
             poly = poly.to_crs(
@@ -238,6 +239,7 @@ def wrap_lons_and_split_at_greenwich(func):
                 )
                 poly = gpd.GeoDataFrame(poly.buffer(0.000000001), columns=["geometry"])
                 poly.crs = crs1
+
             kwargs["poly"] = poly
 
         return func(*args, **kwargs)
@@ -253,7 +255,10 @@ def create_mask(
     poly: gpd.GeoDataFrame = None,
     wrap_lons: bool = False,
 ):
-    """
+    """Creates a mask with values corresponding to the features in a GeoDataFrame.
+
+    The returned mask's points have the value of the first geometry of `poly` they fall in.
+
     Parameters
     ----------
     x_dim : xarray.DataArray
@@ -268,7 +273,33 @@ def create_mask(
     Returns
     -------
     xarray.DataArray
+
+    Examples
+    --------
+    >>> from xclim import subset
+    >>> import xarray as xr
+    >>> import geopandas as gpd
+    >>> ds = xr.open_dataset('example.nc')
+    >>> polys = gpd.read_file('regions.json')
+    Get a mask from all polygons in 'regions.json'
+    >>> mask = subset.create_mask(x_dim=ds.lon, y_dim=ds.lat, poly=polys)
+    >>> ds = ds.assign_coords(regions=mask)
+    Operations can be applied to each regions with  `groupby`. Ex:
+    >>> ds = ds.groupby('regions').mean()
+    Extra step to retrieve the names of those polygons stored in the "id" column
+    >>> region_names = xr.DataArray(polys.id, dims=('regions',)))
+    >>> ds = ds.assign_coords(regions_names=region_names)
     """
+    # Check for intersections
+    for i, (inda, pola) in enumerate(poly.iterrows()):
+        for (indb, polb) in poly.loc[i + 1 :].iterrows():
+            if pola.geometry.intersects(polb.geometry):
+                warnings.warn(
+                    f"List of shapes contains overlap between {inda} and {indb}. Only {inda} will be used.",
+                    UserWarning,
+                    stacklevel=4,
+                )
+
     if len(x_dim.shape) == 1 & len(y_dim.shape) == 1:
         # create a 2d grid of lon, lat values
         lon1, lat1 = np.meshgrid(


Test patch:
diff --git a/tests/test_subset.py b/tests/test_subset.py
index e3eacf7d..54462eea 100644
--- a/tests/test_subset.py
+++ b/tests/test_subset.py
@@ -1,5 +1,6 @@
 import os
 
+import geopandas as gpd
 import numpy as np
 import pytest
 import xarray as xr
@@ -526,6 +527,7 @@ class TestSubsetShape:
     eastern_canada_geojson = os.path.join(TESTS_DATA, "cmip5", "eastern_canada.json")
     southern_qc_geojson = os.path.join(TESTS_DATA, "cmip5", "southern_qc_geojson.json")
     small_geojson = os.path.join(TESTS_DATA, "cmip5", "small_geojson.json")
+    multi_regions_geojson = os.path.join(TESTS_DATA, "cmip5", "multi_regions.json")
 
     def compare_vals(self, ds, sub, vari, flag_2d=False):
         # check subsetted values against original
@@ -649,6 +651,17 @@ class TestSubsetShape:
         assert len(sub.lon.values) == 3
         assert len(sub.lat.values) == 3
 
+    def test_mask_multiregions(self):
+        ds = xr.open_dataset(self.nc_file)
+        regions = gpd.read_file(self.multi_regions_geojson)
+
+        mask = subset.create_mask(
+            x_dim=ds.lon, y_dim=ds.lat, poly=regions, wrap_lons=True
+        )
+        vals, counts = np.unique(mask.values[mask.notnull()], return_counts=True)
+        assert all(vals == [0, 1, 2])
+        assert all(counts == [58, 250, 22])
+
 
 class TestDistance:
     def test_values(self):
diff --git a/tests/testdata/cmip5/multi_regions.json b/tests/testdata/cmip5/multi_regions.json
new file mode 100644
index 00000000..8c0d9a7a
--- /dev/null
+++ b/tests/testdata/cmip5/multi_regions.json
@@ -0,0 +1,162 @@
+{
+  "type": "FeatureCollection",
+  "features": [
+    {
+      "type": "Feature",
+      "properties": {"id": "Québec"},
+      "geometry": {
+        "type": "Polygon",
+        "coordinates": [
+          [
+            [
+              -79.8046875,
+              63.39152174400882
+            ],
+            [
+              -81.9140625,
+              48.22467264956519
+            ],
+            [
+              -75.234375,
+              39.639537564366684
+            ],
+            [
+              -62.9296875,
+              48.69096039092549
+            ],
+            [
+              -48.515625,
+              54.57206165565852
+            ],
+            [
+              -64.6875,
+              63.23362741232569
+            ],
+            [
+              -78.046875,
+              63.39152174400882
+            ],
+            [
+              -79.8046875,
+              63.39152174400882
+            ]
+          ]
+        ]
+      }
+    },
+    {
+      "type": "Feature",
+      "properties": {"id": "Europe"},
+      "geometry": {
+        "type": "Polygon",
+        "coordinates": [
+          [
+            [
+              -18.28125,
+              68.52823492039876
+            ],
+            [
+              -27.0703125,
+              66.37275500247455
+            ],
+            [
+              -18.28125,
+              61.10078883158897
+            ],
+            [
+              -11.25,
+              55.57834467218206
+            ],
+            [
+              -13.0078125,
+              46.800059446787316
+            ],
+            [
+              -10.1953125,
+              32.84267363195431
+            ],
+            [
+              9.140625,
+              37.16031654673677
+            ],
+            [
+              29.53125,
+              33.137551192346145
+            ],
+            [
+              32.6953125,
+              42.8115217450979
+            ],
+            [
+              43.2421875,
+              48.45835188280866
+            ],
+            [
+              38.3203125,
+              65.6582745198266
+            ],
+            [
+              39.0234375,
+              70.8446726342528
+            ],
+            [
+              16.5234375,
+              73.42842364106816
+            ],
+            [
+              -10.546875,
+              71.07405646336098
+            ],
+            [
+              -18.28125,
+              68.00757101804004
+            ],
+            [
+              -18.28125,
+              68.52823492039876
+            ]
+          ]
+        ]
+      }
+    },
+    {
+      "type": "Feature",
+      "properties": {"id": "Newfoundland"},
+      "geometry": {
+        "type": "Polygon",
+        "coordinates": [
+          [
+            [
+              -46.40625,
+              57.70414723434193
+            ],
+            [
+              -57.65624999999999,
+              54.77534585936447
+            ],
+            [
+              -64.6875,
+              46.800059446787316
+            ],
+            [
+              -58.00781249999999,
+              40.713955826286046
+            ],
+            [
+              -47.4609375,
+              45.583289756006316
+            ],
+            [
+              -42.890625,
+              53.74871079689897
+            ],
+            [
+              -46.40625,
+              57.70414723434193
+            ]
+          ]
+        ]
+      }
+    }
+  ]
+}