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

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

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

pip install -e .[dev]
parse_log_pytest
requirements.txt
pytest pytest-cov coveralls
apt-get update apt-get install -y gcc
3.6
requirements_dev.txt
pytest --no-header -rA --tb=line --color=no -p no:cacheprovider -W ignore::DeprecationWarning
attrs==22.2.0
bleach==4.1.0
bokeh==2.3.3
boltons==23.1.1
Bottleneck==1.2.1
bump2version==1.0.1
bumpversion==0.6.0
certifi==2021.5.30
cffi==1.15.1
cfgv==3.3.1
cftime==1.0.3.4
charset-normalizer==2.0.12
click==8.0.4
cloudpickle==2.2.1
colorama==0.4.5
contextvars==2.4
coverage==6.2
coveralls==3.3.1
cryptography==40.0.2
dask==2021.3.0
distlib==0.3.9
distributed==2021.3.0
docopt==0.6.2
docutils==0.18.1
filelock==3.4.1
flake8==5.0.4
fsspec==2022.1.0
HeapDict==1.0.1
identify==2.4.4
idna==3.10
immutables==0.19
importlib-metadata==4.2.0
importlib-resources==5.2.3
iniconfig==1.1.1
jeepney==0.7.1
Jinja2==3.0.3
joblib==1.1.1
keyring==23.4.1
locket==1.0.0
MarkupSafe==2.0.1
mccabe==0.7.0
msgpack==1.0.5
netCDF4==1.6.1
nodeenv==1.6.0
numpy==1.19.5
packaging==21.3
pandas==1.1.5
partd==1.2.0
Pillow==8.4.0
Pint==0.17
pkginfo==1.10.0
platformdirs==2.4.0
pluggy==1.0.0
pre-commit==2.17.0
psutil==7.0.0
py==1.11.0
pycodestyle==2.9.1
pycparser==2.21
pyflakes==2.5.0
Pygments==2.14.0
pyparsing==3.1.4
pyproj==2.4.1
pytest==7.0.1
pytest-cov==4.0.0
pytest-runner==5.3.2
python-dateutil==2.9.0.post0
pytz==2025.2
PyYAML==6.0.1
readme-renderer==34.0
requests==2.27.1
requests-toolbelt==1.0.0
rfc3986==1.5.0
scikit-learn==0.24.2
scipy==1.5.4
SecretStorage==3.3.3
six==1.17.0
sortedcontainers==2.4.0
tblib==1.7.0
threadpoolctl==3.1.0
toml==0.10.2
tomli==1.2.3
toolz==0.12.0
tornado==6.1
tox==3.28.0
tqdm==4.64.1
twine==3.8.0
typing_extensions==4.1.1
urllib3==1.26.20
virtualenv==20.16.2
watchdog==2.3.1
webencodings==0.5.1
xarray==0.13.0
-e git+https://github.com/Ouranosinc/xclim.git@608bac4c9976df778fb9fa442ef725c45dde9f36#egg=xclim
zict==2.1.0
zipp==3.6.0

name: xclim
channels:
  - defaults
  - https://repo.anaconda.com/pkgs/main
  - https://repo.anaconda.com/pkgs/r
  - conda-forge
dependencies:
  - _libgcc_mutex=0.1=main
  - _openmp_mutex=5.1=1_gnu
  - ca-certificates=2025.2.25=h06a4308_0
  - certifi=2021.5.30=py36h06a4308_0
  - ld_impl_linux-64=2.40=h12ee557_0
  - libffi=3.3=he6710b0_2
  - libgcc-ng=11.2.0=h1234567_1
  - libgomp=11.2.0=h1234567_1
  - libstdcxx-ng=11.2.0=h1234567_1
  - ncurses=6.4=h6a678d5_0
  - openssl=1.1.1w=h7f8727e_0
  - pip=21.2.2=py36h06a4308_0
  - python=3.6.13=h12debd9_1
  - readline=8.2=h5eee18b_0
  - setuptools=58.0.4=py36h06a4308_0
  - sqlite=3.45.3=h5eee18b_0
  - tk=8.6.14=h39e8969_0
  - wheel=0.37.1=pyhd3eb1b0_0
  - xz=5.6.4=h5eee18b_1
  - zlib=1.2.13=h5eee18b_1
  - pip:
      - attrs==22.2.0
      - bleach==4.1.0
      - bokeh==2.3.3
      - boltons==23.1.1
      - bottleneck==1.2.1
      - bump2version==1.0.1
      - bumpversion==0.6.0
      - cffi==1.15.1
      - cfgv==3.3.1
      - cftime==1.0.3.4
      - charset-normalizer==2.0.12
      - click==8.0.4
      - cloudpickle==2.2.1
      - colorama==0.4.5
      - contextvars==2.4
      - coverage==6.2
      - coveralls==3.3.1
      - cryptography==40.0.2
      - dask==2021.3.0
      - distlib==0.3.9
      - distributed==2021.3.0
      - docopt==0.6.2
      - docutils==0.18.1
      - filelock==3.4.1
      - flake8==5.0.4
      - fsspec==2022.1.0
      - heapdict==1.0.1
      - identify==2.4.4
      - idna==3.10
      - immutables==0.19
      - importlib-metadata==4.2.0
      - importlib-resources==5.2.3
      - iniconfig==1.1.1
      - jeepney==0.7.1
      - jinja2==3.0.3
      - joblib==1.1.1
      - keyring==23.4.1
      - locket==1.0.0
      - markupsafe==2.0.1
      - mccabe==0.7.0
      - msgpack==1.0.5
      - netcdf4==1.6.1
      - nodeenv==1.6.0
      - numpy==1.19.5
      - packaging==21.3
      - pandas==1.1.5
      - partd==1.2.0
      - pillow==8.4.0
      - pint==0.17
      - pkginfo==1.10.0
      - platformdirs==2.4.0
      - pluggy==1.0.0
      - pre-commit==2.17.0
      - psutil==7.0.0
      - py==1.11.0
      - pycodestyle==2.9.1
      - pycparser==2.21
      - pyflakes==2.5.0
      - pygments==2.14.0
      - pyparsing==3.1.4
      - pyproj==2.4.1
      - pytest==7.0.1
      - pytest-cov==4.0.0
      - pytest-runner==5.3.2
      - python-dateutil==2.9.0.post0
      - pytz==2025.2
      - pyyaml==6.0.1
      - readme-renderer==34.0
      - requests==2.27.1
      - requests-toolbelt==1.0.0
      - rfc3986==1.5.0
      - scikit-learn==0.24.2
      - scipy==1.5.4
      - secretstorage==3.3.3
      - six==1.17.0
      - sortedcontainers==2.4.0
      - tblib==1.7.0
      - threadpoolctl==3.1.0
      - toml==0.10.2
      - tomli==1.2.3
      - toolz==0.12.0
      - tornado==6.1
      - tox==3.28.0
      - tqdm==4.64.1
      - twine==3.8.0
      - typing-extensions==4.1.1
      - urllib3==1.26.20
      - virtualenv==20.16.2
      - watchdog==2.3.1
      - webencodings==0.5.1
      - xarray==0.13.0
      - zict==2.1.0
      - zipp==3.6.0
prefix: /opt/conda/envs/xclim

tests/test_subset.py::TestDistance::test_values tests/test_subset.py::TestDistance::test_broadcasting
tests/test_subset.py::TestSubsetGridPoint::test_simple tests/test_subset.py::TestSubsetGridPoint::test_time_simple tests/test_subset.py::TestSubsetGridPoint::test_time_dates_outofbounds tests/test_subset.py::TestSubsetGridPoint::test_time_start_only tests/test_subset.py::TestSubsetGridPoint::test_time_end_only tests/test_subset.py::TestSubsetGridPoint::test_time_incomplete_years tests/test_subset.py::TestSubsetGridPoint::test_irregular tests/test_subset.py::TestSubsetGridPoint::test_positive_lons tests/test_subset.py::TestSubsetGridPoint::test_raise tests/test_subset.py::TestSubsetGridPoint::test_tolerance tests/test_subset.py::TestSubsetBbox::test_simple tests/test_subset.py::TestSubsetBbox::test_irregular tests/test_subset.py::TestSubsetBbox::test_irregular_datset tests/test_subset.py::TestSubsetBbox::test_single_bounds_curvilinear tests/test_subset.py::TestSubsetBbox::test_positive_lons tests/test_subset.py::TestSubsetBbox::test_time tests/test_subset.py::TestSubsetBbox::test_raise tests/test_subset.py::TestSubsetBbox::test_warnings
tests/test_subset.py::TestSubsetGridPoint::test_dataset tests/test_subset.py::TestSubsetBbox::test_dataset tests/test_subset.py::TestSubsetBbox::test_inverted_coords tests/test_subset.py::TestSubsetBbox::test_single_bounds_rectilinear
HISTORY.rst xclim/subset.py
HISTORY.rst xclim/subset.py