Repository: Ouranosinc/xclim

Issue:
Error converting array to scalars in subset.subset_shape
* xclim version: 0.16.0
* Python version: PAVICS
* Operating System: PAVICS

### Description
I was subsetting a shapefile to be used to calculate ensemble percentiles. The shapefile has multiple watershed sub-boundaries.


### What I Did
```python 
import geopandas as gpd
shp = 'GIS-Data/Beaver.shp'
gdf = gpd.GeoDataFrame.from_file(shp)

# use more percentiles in this example
perc_val = [5, 10, 25, 50, 75, 90, 95]

fig = gdf.plot(facecolor='white')
with ProgressBar():
    dsEnsPoly = subset.subset_shape(dsEns, shape=shp)
    dsPercPoly = ensembles.ensemble_percentiles(dsEnsPoly, values=perc_val)
```
### What I Received
```python
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-30-5aaec4dd7e8c> in <module>
     10 fig = gdf.plot(facecolor='white')
     11 with ProgressBar():
---> 12     dsEnsPoly = subset.subset_shape(dsEns, shape=shp)
     13     dsPercPoly = ensembles.ensemble_percentiles(dsEnsPoly, values=perc_val)
     14     dsPercPoly.sel(time='2041').tg_mean_p50.plot()

/usr/local/envs/birdy/lib/python3.7/site-packages/xclim/subset.py in func_checker(*args, **kwargs)
    206             formatted_args.append(argument)
    207 
--> 208         final = func(*formatted_args, **kwargs)
    209 
    210         for k, v in conv.items():

/usr/local/envs/birdy/lib/python3.7/site-packages/xclim/subset.py in subset_shape(ds, shape, raster_crs, shape_crs, buffer, wrap_lons, start_date, end_date)
    477     # Get the shape's bounding box
    478     bounds = poly.bounds
--> 479     lon_bnds = (float(bounds.minx.values), float(bounds.maxx.values))
    480     lat_bnds = (float(bounds.miny.values), float(bounds.maxy.values))
    481 

TypeError: only size-1 arrays can be converted to Python scalars
```

### Further Investigation
```python
poly = gpd.GeoDataFrame.from_file('GIS-Data/Aguasabon_Project.shp')
poly
bounds = poly.bounds
bounds
#float(min(bounds.minx.values))
#float(max(bounds.maxx.values))
```
minx	miny	maxx	maxy
0	-87.464326	49.216421	-86.632188	49.816493
1	-87.441197	49.067714	-86.342397	50.158592
2	-87.347309	48.771573	-86.930992	49.280212

The last two comments could fix the issue if they meet the intent of the rest of the code.


Discussion:

Patch:
diff --git a/docs/notebooks/example.ipynb b/docs/notebooks/example.ipynb
index 344812d6..6bc3484a 100644
--- a/docs/notebooks/example.ipynb
+++ b/docs/notebooks/example.ipynb
@@ -51,40 +51,27 @@
    "source": [
     "## Setting up the Dask client: parallel processing\n",
     "\n",
-    "This step is not mandatory as dask will fall back to its \"single machine scheduler\" if a `Client` is not created (the case is this notebook). However, this default scheduler doesn't allow you to set the number of threads or a memory limit and doesn't start the dashboard, which can be quite useful to understand your task's progress.\n",
+    "First we create a pool of workers that will wait for jobs. The `xarray` library will automatically connect to these workers and and dispatch them jobs that can be run in parallel.  \n",
     "\n",
-    "The `xarray` library will automatically connect to the dask scheduler when needed.\n",
+    "The dashboard link lets you see in real time how busy those workers are.\n",
     "\n",
-    "The dashboard link lets you see in real time how busy those workers are. More info : [dask distributed documentation](https://distributed.dask.org/en/latest/).\n",
+    "* [dask distributed documentation](https://distributed.dask.org/en/latest/)\n",
     "\n",
-    "```python\n",
+    "This step is not mandatory as dask will fall back to its \"single machine scheduler\" if a `Client` is not created. However, this default scheduler doesn't allow you to set the number of threads or a memory limit and doesn't start the dashboard, which can be quite useful to understand your task's progress.\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
     "from distributed import Client\n",
     "\n",
     "# Depending on your workstation specifications, you may need to adjust these values.\n",
     "# On a single machine, n_workers=1 is usually better.\n",
     "client=Client(n_workers=1, threads_per_worker=4, memory_limit=\"4GB\") \n",
-    "client\n",
-    "```\n",
-    "\n",
-    "<table style=\"border: 2px solid white;\">\n",
-    "<tr>\n",
-    "<td style=\"vertical-align: top; border: 0px solid white\">\n",
-    "<h3 style=\"text-align: left;\">Client</h3>\n",
-    "<ul style=\"text-align: left; list-style: none; margin: 0; padding: 0;\">\n",
-    "  <li><b>Scheduler: </b>tcp://127.0.0.1:36783</li>\n",
-    "  <li><b>Dashboard: </b><a href='http://127.0.0.1:35847/status' target='_blank'>http://127.0.0.1:35847/status</a>\n",
-    "</ul>\n",
-    "</td>\n",
-    "<td style=\"vertical-align: top; border: 0px solid white\">\n",
-    "<h3 style=\"text-align: left;\">Cluster</h3>\n",
-    "<ul style=\"text-align: left; list-style:none; margin: 0; padding: 0;\">\n",
-    "  <li><b>Workers: </b>1</li>\n",
-    "  <li><b>Cores: </b>4</li>\n",
-    "  <li><b>Memory: </b>4.00 GB</li>\n",
-    "</ul>\n",
-    "</td>\n",
-    "</tr>\n",
-    "</table>"
+    "client"
    ]
   },
   {
@@ -120,8 +107,8 @@
    "outputs": [],
    "source": [
     "# Chunking in memory along the time dimension.\n",
-    "# Note that the data type is a 'dask.array'. xarray will automatically use the dask scheduler.\n",
-    "ds = xr.open_dataset(data_url, chunks={'time': 365, 'lat': 84, 'lon': 75}, drop_variables=['ts', 'time_vectors'])\n",
+    "# Note that the data type is a 'dask.array'. xarray will automatically use client workers.\n",
+    "ds = xr.open_dataset(data_url, chunks={'time': 365, 'lat': 168, 'lon': 150}, drop_variables=['ts', 'time_vectors'])\n",
     "ds"
    ]
   },
@@ -301,7 +288,7 @@
    "outputs": [],
    "source": [
     "%%time\n",
-    "output_file = output_folder / 'test_tx_days_above.nc'\n",
+    "output_file = output_folder / 'test_tx_max.nc'\n",
     "dsOut.to_netcdf(output_file)"
    ]
   },
@@ -333,7 +320,7 @@
    "outputs": [],
    "source": [
     "# rechunk data in memory for the entire grid  \n",
-    "ds2c = ds2.chunk(chunks={'time': 365})\n",
+    "ds2c = ds2.chunk(chunks={'time':4 * 365})\n",
     "print(ds2c.chunks['time'])"
    ]
   },
@@ -344,7 +331,7 @@
    "outputs": [],
    "source": [
     "%%time\n",
-    "out = xc.atmos.tx_days_above(ds2c.tasmax, thresh='30 C', freq='YS')\n",
+    "out = xc.atmos.tx_max(ds2c.tasmax, freq='YS')\n",
     "dsOut = xr.Dataset(data_vars=None, coords=out.coords, attrs=ds.attrs)\n",
     "dsOut[out.name] = out\n",
     "\n",
@@ -366,7 +353,93 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "# ds4 = ds3.load()"
+    "ds4 = ds3.load()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Unit handling in `xclim`\n",
+    "\n",
+    "A lot of effort has been placed into automatic handling of input data units. `xclim` will automatically detect the input variable(s) units (e.g. °C versus °K or mm/s versus mm/day etc.) and adjust on-the-fly in order to calculate indices in the consistent manner. This comes with the obvious caveat that input data requires metadata attribute for units\n",
+    "\n",
+    "In the example below, we compute weekly total precipitation in mm using inputs of mm/s and mm/d. As you see, the output is identical."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Compute with the original mm s-1 data\n",
+    "out1 = xc.atmos.precip_accumulation(ds4.pr, freq='MS')\n",
+    "# Create a copy of the data converted to mm d-1\n",
+    "pr_mmd = ds4.pr * 3600 * 24\n",
+    "pr_mmd.attrs[\"units\"] = \"mm d-1\"\n",
+    "out2 = xc.atmos.precip_accumulation(pr_mmd, freq='MS')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# import plotting stuff\n",
+    "import matplotlib.pyplot as plt\n",
+    "%matplotlib inline\n",
+    "plt.style.use('seaborn')\n",
+    "plt.rcParams['figure.figsize'] = (11, 5)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "plt.figure()\n",
+    "out1.plot(label='From mm s-1', linestyle='-')\n",
+    "out2.plot(label='From mm d-1', linestyle='none', marker='o')\n",
+    "plt.legend()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Threshold indices\n",
+    "\n",
+    "`xclim` unit handling also applies to threshold indicators. Users can provide threshold in units of choice and `xclim` will adjust automatically. For example determining the number of days with tasmax > 20°C users can define a threshold input of '20 C' or '20 degC' even if input data is in Kelvin. Alernatively users can even provide a threshold in Kelvin '293.15 K' (if they really wanted to)."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Create a copy of the data converted to C\n",
+    "tasmax_C = ds4.tasmax - 273.15\n",
+    "tasmax_C.attrs['units'] = 'C'\n",
+    "\n",
+    "# Using Kelvin data, threshold in Celsius\n",
+    "out1 = xc.atmos.tx_days_above(ds4.tasmax, thresh='20 C', freq='MS')\n",
+    "\n",
+    "# Using Celsius data\n",
+    "out2 = xc.atmos.tx_days_above(tasmax_C, thresh='20 C', freq='MS')\n",
+    "\n",
+    "# Using Celsius but with threshold in Kelvin\n",
+    "out3 = xc.atmos.tx_days_above(tasmax_C, thresh='293.15 K', freq='MS')\n",
+    "\n",
+    "# Plot and see that it's all identical:\n",
+    "plt.figure()\n",
+    "out1.plot(label='K and degC', linestyle='-')\n",
+    "out2.plot(label='degC and degC', marker='s', markersize=10, linestyle='none')\n",
+    "out3.plot(label='degC and K', marker='o', linestyle='none')\n",
+    "plt.legend()"
    ]
   },
   {
@@ -378,7 +451,6 @@
   }
  ],
  "metadata": {
-  "celltoolbar": "Raw Cell Format",
   "kernelspec": {
    "display_name": "Python 3",
    "language": "python",
@@ -394,7 +466,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.8.1"
+   "version": "3.8.0"
   }
  },
  "nbformat": 4,
diff --git a/xclim/indicators/atmos/_conversion.py b/xclim/indicators/atmos/_conversion.py
index 771b04de..b61ef4fc 100644
--- a/xclim/indicators/atmos/_conversion.py
+++ b/xclim/indicators/atmos/_conversion.py
@@ -50,7 +50,7 @@ saturation_vapor_pressure = Converter(
     units="Pa",
     long_name="Saturation vapor pressure",
     description=lambda **kws: (
-        "The saturation vapor pressure was calculated from the temperature "
+        "The saturation vapor pressure was calculated from a temperature "
         "according to the {method} method."
     )
     + (
@@ -68,13 +68,18 @@ relative_humidity_from_dewpoint = Converter(
     units="%",
     long_name="Relative Humidity",
     standard_name="relative_humidity",
-    description="Computed from temperature and dew point temperature.",
+    description=lambda **kws: (
+        "Computed from temperature, and dew point temperature through the "
+        "saturation vapor pressures, which were calculated "
+        "according to the {method} method."
+    )
+    + (
+        " The computation was done in reference to ice for temperatures below {ice_thresh}."
+        if kws["ice_thresh"] is not None
+        else ""
+    ),
     compute=wrapped_partial(
-        indices.relative_humidity,
-        huss=None,
-        ps=None,
-        method="dewpoint",
-        invalid_values="mask",
+        indices.relative_humidity, huss=None, ps=None, invalid_values="mask",
     ),
 )
 
diff --git a/xclim/indices/_conversion.py b/xclim/indices/_conversion.py
index 8c70b548..e023d789 100644
--- a/xclim/indices/_conversion.py
+++ b/xclim/indices/_conversion.py
@@ -181,7 +181,7 @@ def saturation_vapor_pressure(
 
     - "goffgratch46" or "GG46", based on [goffgratch46]_, values and equation taken from [voemel]_.
     - "sonntag90" or "SO90", taken from [sonntag90]_.
-    - "tetens30" or "TE30", based on [tetens30], values and equation taken from [voemel]_.
+    - "tetens30" or "TE30", based on [tetens30]_, values and equation taken from [voemel]_.
     - "wmo08" or "WMO08", taken from [wmo08]_.
 
 
@@ -191,7 +191,7 @@ def saturation_vapor_pressure(
     .. [sonntag90] Sonntag, D. (1990). Important new values of the physical constants of 1986, vapour pressure formulations based on the ITS-90, and psychrometer formulae. Zeitschrift für Meteorologie, 40(5), 340-344.
     .. [tetens30] Tetens, O. 1930. Über einige meteorologische Begriffe. Z. Geophys 6: 207-309.
     .. [voemel] http://cires1.colorado.edu/~voemel/vp.html
-    .. [wmo2008] World Meteorological Organization. (2008). Guide to meteorological instruments and methods of observation. Geneva, Switzerland: World Meteorological Organization. https://www.weather.gov/media/epz/mesonet/CWOP-WMO8.pdf
+    .. [wmo08] World Meteorological Organization. (2008). Guide to meteorological instruments and methods of observation. Geneva, Switzerland: World Meteorological Organization. https://www.weather.gov/media/epz/mesonet/CWOP-WMO8.pdf
     """
     if ice_thresh is not None:
         thresh = convert_units_to(ice_thresh, "degK")
@@ -279,71 +279,84 @@ def relative_humidity(
     method: str = "sonntag90",
     invalid_values: str = "clip",
 ) -> xr.DataArray:
-    """Compute relative humidity from temperature and either dewpoint temperature
-    or from specific humidity and pressure (through the saturation vapor pressure)
+    r"""
+    Compute relative humidity from temperature and either dewpoint temperature or specific humidity and pressure through the saturation vapor pressure.
 
     Parameters
     ----------
     tas : xr.DataArray
         Temperature array
     dtas : xr.DataArray
-        Dewpoint temperature, if specified, "method" must be set to "dewpoint".
+        Dewpoint temperature, if specified, overrides huss and ps.
     huss : xr.DataArray
         Specific Humidity
     ps : xr.DataArray
         Air Pressure
     ice_thresh : str
         Threshold temperature under which to switch to equations in reference to ice instead of water.
-        If None (default) everything is computed with reference to water. Does nothing if 'method' is "dewpoint,"
-    method : {"dewpoint", "goffgratch46", "sonntag90", "tetens30", "wmo08"}
+        If None (default) everything is computed with reference to water. Does nothing if 'method' is "bohren98".
+    method : {"bohren98", "goffgratch46", "sonntag90", "tetens30", "wmo08"}
         Which method to use, see notes of this function and of `saturation_vapor_pressure`.
     invalid_values : {"clip", "mask", None}
-        What to do with values outside the 0-100 range.
-        If "clip" (default), clips everything to 0 - 100,
-        if "mask", replaces values outside the range by np.nan,
-        if None, does nothing.
+        What to do with values outside the 0-100 range. If "clip" (default), clips everything to 0 - 100,
+        if "mask", replaces values outside the range by np.nan, and if `None` , does nothing.
 
     Notes
     -----
-
     In the following, let :math:`T`, :math:`T_d`, :math:`q` and :math:`p` be the temperature,
     the dew point temperature, the specific humidity and the air pressure.
 
-    **For the "dewpoint" method** : With :math:`L` the Enthalpy of vaporization of water
-    and :math:`R_w` the gas constant for water vapor, the relative humidity is computed as:
+    **For the "bohren98" method** : This method does not use the saturation vapor pressure directly,
+    but rather uses an approximation of the ratio of :math:`\frac{e_{sat}(T_d)}{e_{sat}(T)}`.
+    With :math:`L` the enthalpy of vaporization of water and :math:`R_w` the gas constant for water vapor,
+    the relative humidity is computed as:
 
     .. math::
 
-        RH = e^{\\frac{-L (T - T_d)}{R_wTT_d}}
+        RH = e^{\frac{-L (T - T_d)}{R_wTT_d}}
 
-    Formula taken from [Lawrence_2005]_.
+    From [BohrenAlbrecht1998]_, formula taken from [Lawrence2005]_. :math:`L = 2.5\times 10^{-6}` J kg-1, exact for :math:`T = 273.15` K, is used.
 
     **Other methods**: With :math:`w`, :math:`w_{sat}`, :math:`e_{sat}` the mixing ratio,
-    the saturation mixing ratio and the saturation vapor pressure, relative humidity is computed as:
+    the saturation mixing ratio and the saturation vapor pressure.
+    If the dewpoint temperature is given, relative humidity is computed as:
 
-        ... math::
+    .. math::
 
-            RH = 100\\frac{w}{w_{sat}}
-            w = \\frac{q}{1-q}
-            w_{sat} = 0.622\\frac{e_{sat}}{P - e_{sat}}
+        RH = 100\frac{e_{sat}(T_d)}{e_{sat}(T)}
 
-    The methods differ by how :math:`e_{sat}` is computed. See the doc of `xclim.core.utils.saturation_vapor_pressure`.
+    Otherwise, the specific humidity and the air pressure must be given so relative humidity can be computed as:
+
+    .. math::
+
+        RH = 100\frac{w}{w_{sat}}
+        w = \frac{q}{1-q}
+        w_{sat} = 0.622\frac{e_{sat}}{P - e_{sat}}
+
+    The methods differ by how :math:`e_{sat}` is computed. See the doc of :py:meth:`xclim.core.utils.saturation_vapor_pressure`.
 
     References
     ----------
-    .. [Lawrence_2005] Lawrence, M.G. (2005). The Relationship between Relative Humidity and the Dewpoint Temperature in Moist Air: A Simple Conversion and Applications. Bull. Amer. Meteor. Soc., 86, 225–234, https://doi.org/10.1175/BAMS-86-2-225
+    .. [Lawrence2005] Lawrence, M.G. (2005). The Relationship between Relative Humidity and the Dewpoint Temperature in Moist Air: A Simple Conversion and Applications. Bull. Amer. Meteor. Soc., 86, 225–234, https://doi.org/10.1175/BAMS-86-2-225
+    .. [BohrenAlbrecht1998] Craig F. Bohren, Bruce A. Albrecht. Atmospheric Thermodynamics. Oxford University Press, 1998.
     """
-    if dtas is not None and method != "dewpoint":
-        raise ValueError(
-            "If the dewpoint temperature (dtas) is passed, method must be set to 'dewpoint'"
-        )
 
-    if method == "dewpoint":
+    if method in ("bohren98", "BA90"):
+        if dtas is None:
+            raise ValueError("To use method 'bohren98' (BA98), dewpoint must be given.")
         dtas = convert_units_to(dtas, "degK")
         tas = convert_units_to(tas, "degK")
         L = 2.501e6
         Rw = (461.5,)
         rh = 100 * np.exp(-L * (tas - dtas) / (Rw * tas * dtas))
+    elif dtas is not None:
+        e_sat_dt = saturation_vapor_pressure(
+            tas=dtas, ice_thresh=ice_thresh, method=method
+        )
+        e_sat_t = saturation_vapor_pressure(
+            tas=tas, ice_thresh=ice_thresh, method=method
+        )
+        rh = 100 * e_sat_dt / e_sat_t
     else:
         ps = convert_units_to(ps, "Pa")
         huss = convert_units_to(huss, "")
@@ -386,7 +399,7 @@ def specific_humidity(
     ice_thresh : str
         Threshold temperature under which to switch to equations in reference to ice instead of water.
         If None (default) everything is computed with reference to water.
-    method : {"dewpoint", "goffgratch46", "sonntag90", "tetens30", "wmo08"}
+    method : {"goffgratch46", "sonntag90", "tetens30", "wmo08"}
         Which method to use, see notes of this function and of `saturation_vapor_pressure`.
     invalid_values : {"clip", "mask", None}
         What to do with values larger than the saturation specific humidity and lower than 0.
diff --git a/xclim/subset.py b/xclim/subset.py
index 281e6977..cb764868 100644
--- a/xclim/subset.py
+++ b/xclim/subset.py
@@ -555,9 +555,9 @@ def subset_shape(
         poly.geometry = poly.buffer(buffer)
 
     # Get the shape's bounding box.
-    bounds = poly.bounds
-    lon_bnds = (float(bounds.minx.values), float(bounds.maxx.values))
-    lat_bnds = (float(bounds.miny.values), float(bounds.maxy.values))
+    minx, miny, maxx, maxy = poly.total_bounds
+    lon_bnds = (minx, maxx)
+    lat_bnds = (miny, maxy)
 
     # If polygon doesn't cross prime meridian, subset bbox first to reduce processing time
     # Only case not implemented is when lon_bnds cross the 0 deg meridian but dataset grid has all positive lons


Test patch:
diff --git a/tests/test_atmos.py b/tests/test_atmos.py
index 5d0f57a1..670e16ed 100644
--- a/tests/test_atmos.py
+++ b/tests/test_atmos.py
@@ -28,7 +28,7 @@ def test_relative_humidity_dewpoint(tas_series, rh_series):
             dtas=tas_series(np.array([-15, -10, -2, 5, 10, 20, 29, 20, 30]) + K2C),
         ),
         rh_series([np.nan, 100, 93, 71, 52, 73, 94, 31, 20]),
-        rtol=0.01,
+        rtol=0.02,
         atol=1,
     )
 
diff --git a/tests/test_indices.py b/tests/test_indices.py
index 21e7a7ce..10e7afbb 100644
--- a/tests/test_indices.py
+++ b/tests/test_indices.py
@@ -1150,20 +1150,25 @@ class TestWindConversion:
         )
 
 
+@pytest.mark.parametrize(
+    "method", ["bohren98", "tetens30", "sonntag90", "goffgratch46", "wmo08"]
+)
 @pytest.mark.parametrize(
     "invalid_values,exp0", [("clip", 100), ("mask", np.nan), (None, 151)]
 )
-def test_relative_humidity_dewpoint(tas_series, rh_series, invalid_values, exp0):
+def test_relative_humidity_dewpoint(
+    tas_series, rh_series, method, invalid_values, exp0
+):
     np.testing.assert_allclose(
         xci.relative_humidity(
             tas=tas_series(np.array([-20, -10, -1, 10, 20, 25, 30, 40, 60]) + K2C),
             dtas=tas_series(np.array([-15, -10, -2, 5, 10, 20, 29, 20, 30]) + K2C),
-            method="dewpoint",
+            method=method,
             invalid_values=invalid_values,
         ),
         # Expected values obtained by hand calculation
         rh_series([exp0, 100, 93, 71, 52, 73, 94, 31, 20]),
-        rtol=0.01,
+        rtol=0.02,
         atol=1,
     )
 
diff --git a/tests/test_subset.py b/tests/test_subset.py
index 25c73870..4f0e2c95 100644
--- a/tests/test_subset.py
+++ b/tests/test_subset.py
@@ -769,6 +769,13 @@ class TestSubsetShape:
         assert all(vals == [0, 1, 2])
         assert all(counts == [58, 250, 22])
 
+    def test_subset_multiregions(self):
+        ds = xr.open_dataset(self.nc_file)
+        regions = gpd.read_file(self.multi_regions_geojson)
+        regions.set_index("id")
+        ds_sub = subset.subset_shape(ds, shape=regions)
+        assert ds_sub.notnull().sum() == 58 + 250 + 22
+
 
 class TestDistance:
     def test_values(self):