Repository: Ouranosinc/xclim

Issue:
add optional threshold in utils.subset_gridpoint
* xclim version:
* Python version:
* Operating System:

### Description

Describe what you were trying to get done.
Tell us what happened, what went wrong, and what you expected to happen.

### What I Did

```
Paste the command(s) you ran and the output.
If there was a crash, please include the traceback here.
```


Discussion:

Patch:
diff --git a/HISTORY.rst b/HISTORY.rst
index ced91a46..0a2c589e 100644
--- a/HISTORY.rst
+++ b/HISTORY.rst
@@ -2,6 +2,11 @@
 History
 =======
 
+0.12.x-beta (2019-11-15)
+------------------------
+* Added a distance function computing the geodesic distance to a point.
+* Added a `tolerance` argument to `subset_gridpoint` raising an error if distance to closest point is larger than tolerance.
+
 0.11.x-beta (2019-10-17)
 ------------------------
 * Added type hinting to call signatures of many functions for more explicit type-checking.
diff --git a/xclim/subset.py b/xclim/subset.py
index c6d8f935..a5f50462 100644
--- a/xclim/subset.py
+++ b/xclim/subset.py
@@ -313,6 +313,7 @@ def subset_gridpoint(
     lat: Optional[float] = None,
     start_date: Optional[str] = None,
     end_date: Optional[str] = None,
+    tolerance: Optional[float] = None,
 ) -> Union[xarray.DataArray, xarray.Dataset]:
     """Extract a nearest gridpoint from datarray based on lat lon coordinate.
 
@@ -342,6 +343,8 @@ def subset_gridpoint(
     end_yr : int
       Deprecated
         Last year of the subset. Defaults to last year of input data-array.
+    tolerance : Optional[float]
+      Raise error if the distance to the nearest gridpoint is larger than tolerance in meters.
 
     Returns
     -------
@@ -374,23 +377,20 @@ def subset_gridpoint(
             # if 'lon' and 'lat' are present as data dimensions use the .sel method.
             if "lat" in dims and "lon" in dims:
                 da = da.sel(lat=lat, lon=lon, method="nearest")
+
+                if tolerance is not None:
+                    # Calculate the geodesic distance between grid points and the point of interest.
+                    dist = distance(da, lon, lat)
+
             else:
-                g = Geod(ellps="WGS84")  # WGS84 ellipsoid - decent globaly
-                lon1 = da.lon.values
-                lat1 = da.lat.values
-                shp_orig = lon1.shape
-                lon1 = np.reshape(lon1, lon1.size)
-                lat1 = np.reshape(lat1, lat1.size)
-                # calculate geodesic distance between grid points and point of interest
-                az12, az21, dist = g.inv(
-                    lon1,
-                    lat1,
-                    np.broadcast_to(lon, lon1.shape),
-                    np.broadcast_to(lat, lat1.shape),
-                )
-                dist = dist.reshape(shp_orig)
-                iy, ix = np.unravel_index(np.argmin(dist, axis=None), dist.shape)
-                xydims = [x for x in da.lon.dims]
+                # Calculate the geodesic distance between grid points and the point of interest.
+                dist = distance(da, lon, lat)
+
+                # Find the indices for the closest point
+                iy, ix = np.unravel_index(dist.argmin(), dist.shape)
+
+                # Select data from closest point
+                xydims = [x for x in dist.dims]
                 args = dict()
                 args[xydims[0]] = iy
                 args[xydims[1]] = ix
@@ -402,6 +402,14 @@ def subset_gridpoint(
                 )
             )
 
+    if tolerance is not None:
+        if dist.min() > tolerance:
+            raise ValueError(
+                "Distance to closest point ({}) is larger than tolerance ({})".format(
+                    dist, tolerance
+                )
+            )
+
     if start_date or end_date:
         da = subset_time(da, start_date=start_date, end_date=end_date)
 
@@ -458,3 +466,38 @@ def subset_time(
     """
 
     return da.sel(time=slice(start_date, end_date))
+
+
+def distance(da, lon, lat):
+    """Return distance to point in meters.
+
+    Parameters
+    ----------
+    da : Union[xarray.DataArray, xarray.Dataset]
+      Input data.
+    lon : Optional[float]
+      Longitude coordinate.
+    lat : Optional[float]
+      Latitude coordinate.
+
+    Returns
+    -------
+    xarray.DataArray
+      Distance in meters to point.
+
+    Note
+    ----
+    To get the indices from closest point, use
+    >>> d = distance(da)
+    >>> k = d.argmin()
+    >>> i, j = numpy.unravel_index(k, d.shape)
+
+    """
+    g = Geod(ellps="WGS84")  # WGS84 ellipsoid - decent globaly
+
+    def func(lons, lats):
+        return g.inv(*np.broadcast_arrays(lons, lats, lon, lat))[2]
+
+    out = xarray.apply_ufunc(func, da.lon, da.lat)
+    out.attrs["units"] = "m"
+    return out


Test patch:
diff --git a/tests/test_subset.py b/tests/test_subset.py
index e1851b03..8a47743a 100644
--- a/tests/test_subset.py
+++ b/tests/test_subset.py
@@ -245,6 +245,15 @@ class TestSubsetGridPoint:
         with pytest.raises(Exception):
             subset.subset_gridpoint(da, lon=-72.4, lat=46.1)
 
+    def test_tolerance(self):
+        da = xr.open_dataset(self.nc_poslons).tas
+        lon = -72.5
+        lat = 46.2
+        with pytest.raises(ValueError):
+            subset.subset_gridpoint(da, lon=lon, lat=lat, tolerance=1)
+
+        subset.subset_gridpoint(da, lon=lon, lat=lat, tolerance=1e5)
+
 
 class TestSubsetBbox:
     nc_poslons = os.path.join(
@@ -506,3 +515,37 @@ class TestSubsetBbox:
             ' support "start_date" and "end_date" (type: str) using formats of "%Y", "%Y-%m" or "%Y-%m-%d".'
             not in [q.message for q in record]
         )
+
+
+class TestDistance:
+    def test_values(self):
+        # Check values are OK. Values taken from pyproj test.
+        boston_lat = 42.0 + (15.0 / 60.0)
+        boston_lon = -71.0 - (7.0 / 60.0)
+        portland_lat = 45.0 + (31.0 / 60.0)
+        portland_lon = -123.0 - (41.0 / 60.0)
+
+        da = xr.DataArray(
+            0, coords={"lon": [boston_lon], "lat": [boston_lat]}, dims=["lon", "lat"]
+        )
+        d = subset.distance(da, lon=portland_lon, lat=portland_lat)
+        np.testing.assert_almost_equal(d, 4164074.239, decimal=3)
+
+    def test_broadcasting(self):
+        # Check output dimensions match lons and lats.
+        lon = np.linspace(-180, 180, 20)
+        lat = np.linspace(-90, 90, 30)
+        da = xr.Dataset(data_vars=None, coords={"lon": lon, "lat": lat})
+        da["data"] = xr.DataArray(
+            np.random.rand(lon.size, lat.size), dims=["lon", "lat"]
+        )
+
+        d = subset.distance(da, -34, 56)
+        assert d.dims == da.data.dims
+        assert d.shape == da.data.shape
+        assert d.units == "m"
+
+        # Example of how to get the actual 2D indices.
+        k = d.argmin()
+        i, j = np.unravel_index(k, da.data.shape)
+        assert d[i, j] == d.min()