Repository: Ouranosinc/xclim

Issue:
Create land *domain* to store streamflow indices
### Description
<!--Describe what you were trying to get done.
Tell us what happened, what went wrong, and what you expected to happen.-->

We have only *atmos* at the moment. 

Also, the documentation does not include the streamflow related indices or the frequency analysis code. 

Discussion:

Patch:
diff --git a/HISTORY.rst b/HISTORY.rst
index a42ac318..65cfc757 100644
--- a/HISTORY.rst
+++ b/HISTORY.rst
@@ -2,19 +2,24 @@
 History
 =======
 
+0.12.x-beta (2019-11-18)
+------------------------
+* 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.
+* Created land module for standardized access to streamflow indices.
+
 0.11.x-beta (2019-10-17)
 ------------------------
 * Added type hinting to call signatures of many functions for more explicit type-checking.
-* Now supporting explicit builds for Windows OS via Travis CI.
 * Added Kmeans clustering ensemble reduction algorithms.
 * Added utilities for converting between wind velocity (sfcWind) and wind components (uas, vas) arrays.
-* Added type hinting to call signatures of many functions for more explicit type-checking
-* Now supporting explicit builds for Windows OS via Travis CI
-* Fix failing test with Python 3.7
-* Fixed bug in subset.subset_bbox that could add unwanted coordinates/dims to some variables when applied to an entire dataset
-* Reformatted packaging configuration to pure Py3 wheel that ignore tests and test data
-* Now officially supporting Python3.8
-* Enhancement to precip_accumulation() to allow estimated amounts solid (or liquid) phase precipitation
+* Added type hinting to call signatures of many functions for more explicit type-checking.
+* Now supporting explicit builds for Windows OS via Travis CI.
+* Fix failing test with Python 3.7.
+* Fixed bug in subset.subset_bbox that could add unwanted coordinates/dims to some variables when applied to an entire dataset.
+* Reformatted packaging configuration to pure Py3 wheel that ignore tests and test data.
+* Now officially supporting Python3.8!
+* Enhancement to precip_accumulation() to allow estimated amounts solid (or liquid) phase precipitation.
 * Bugfix for frequency analysis choking on time series with NaNs only.
 
 0.10.x-beta (2019-06-18)
@@ -36,9 +41,9 @@ History
 * Now employing PEP8 + Black compatible autoformatting.
 * Added Windows and macOS images to Travis CI build ensemble.
 * Added variable thresholds for tasmax and tasmin in daily_freezethaw_events.
-* Updated subset.py to use date formatted strings ("%Y", "%Y%m" etc.) in temporal subsetting
+* Updated subset.py to use date formatted strings ("%Y", "%Y%m" etc.) in temporal subsetting.
 * Clean-up of day-of-year resampling. Precipitation percentile threshold will work without a doy index.
-* Addressed deprecations for xarray 0.13.0
+* Addressed deprecations for xarray 0.13.0.
 * Added a decorator function that verifies validity and reformats subset calls using start_date or end_date signatures.
 * Fixed a bug where 'lon' or 'lon_bounds' would return false values if either signatures were set to 0.
 
diff --git a/docs/index.rst b/docs/index.rst
index 45822ea0..a5d73478 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -10,6 +10,7 @@ xclim Official Documentation
    usage
    notebooks/example
    indices
+   modules
    icclim
    utilities
    contributing
diff --git a/docs/modules.rst b/docs/modules.rst
new file mode 100644
index 00000000..4e7a0bd0
--- /dev/null
+++ b/docs/modules.rst
@@ -0,0 +1,35 @@
+Modules
+=======
+
+Atmos
+-----
+
+.. warning::
+    This is a work in progress!
+
+.. automodule:: xclim.atmos._precip
+   :members:
+   :private-members:
+   :inherited-members:
+   :show-inheritance:
+   :undoc-members:
+
+.. automodule:: xclim.atmos._temperature
+   :members:
+   :private-members:
+   :inherited-members:
+   :show-inheritance:
+   :undoc-members:
+
+Land
+----
+
+.. warning::
+    This is a work in progress!
+
+.. automodule:: xclim.land._streamflow
+   :members:
+   :private-members:
+   :inherited-members:
+   :show-inheritance:
+   :undoc-members:
diff --git a/setup.cfg b/setup.cfg
index 8792cf74..00ce2ce1 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -1,5 +1,5 @@
 [bumpversion]
-current_version = 0.11.4
+current_version = 0.12.0
 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 c134f0a7..3bd00a04 100644
--- a/setup.py
+++ b/setup.py
@@ -10,7 +10,7 @@ URL = "https://github.com/Ouranosinc/xclim"
 AUTHOR = "Travis Logan"
 AUTHOR_EMAIL = "logan.travis@ouranos.ca"
 REQUIRES_PYTHON = ">=3.5.0"
-VERSION = "0.11.4"
+VERSION = "0.12.0"
 LICENSE = "Apache Software License 2.0"
 
 with open("README.rst") as readme_file:
diff --git a/xclim/__init__.py b/xclim/__init__.py
index 78824dde..c623efc7 100644
--- a/xclim/__init__.py
+++ b/xclim/__init__.py
@@ -9,7 +9,7 @@ from xclim import indices
 
 __author__ = """Travis Logan"""
 __email__ = "logan.travis@ouranos.ca"
-__version__ = "0.11.4"
+__version__ = "0.12.0"
 
 
 def build_module(name, objs, doc="", source=None, mode="ignore"):
diff --git a/xclim/land/__init__.py b/xclim/land/__init__.py
new file mode 100644
index 00000000..8ae9353f
--- /dev/null
+++ b/xclim/land/__init__.py
@@ -0,0 +1,1 @@
+from ._streamflow import *
diff --git a/xclim/land/_streamflow.py b/xclim/land/_streamflow.py
new file mode 100644
index 00000000..93f37890
--- /dev/null
+++ b/xclim/land/_streamflow.py
@@ -0,0 +1,112 @@
+# -*- coding: utf-8 -*-
+import numpy as np
+
+from xclim import checks
+from xclim import generic
+from xclim.indices import base_flow_index
+from xclim.utils import Indicator
+from xclim.utils import wrapped_partial
+
+# from boltons.funcutils import FunctionBuilder
+# import calendar
+
+__all__ = ["base_flow_index", "freq_analysis", "stats", "fit", "doy_qmax", "doy_qmin"]
+
+
+class Streamflow(Indicator):
+    units = "m^3 s-1"
+    context = "hydro"
+    standard_name = "discharge"
+
+    @staticmethod
+    def compute(*args, **kwds):
+        pass
+
+    def cfprobe(self, q):
+        checks.check_valid(q, "standard_name", "streamflow")
+
+
+class Stats(Streamflow):
+    def missing(self, *args, **kwds):
+        """Return whether an output is considered missing or not."""
+        from functools import reduce
+
+        indexer = kwds["indexer"]
+        freq = kwds["freq"] or generic.default_freq(**indexer)
+
+        miss = (checks.missing_any(da, freq, **indexer) for da in args)
+        return reduce(np.logical_or, miss)
+
+
+# Disable the missing value check because the output here is not a time series.
+class FA(Streamflow):
+    def missing(self, *args, **kwds):
+        """Return whether an output is considered missing or not."""
+        return False
+
+
+# Disable the daily checks because the inputs are period extremas.
+class Fit(FA):
+    def validate(self, da):
+        pass
+
+
+base_flow_index = Streamflow(
+    identifier="base_flow_index",
+    units="",
+    long_name="Base flow index",
+    compute=base_flow_index,
+)
+
+
+freq_analysis = FA(
+    identifier="freq_analysis",
+    var_name="q{window}{mode}{indexer}",
+    long_name="N-year return period {mode} {indexer} {window}-day flow",
+    description="Streamflow frequency analysis for the {mode} {indexer} {window}-day flow "
+    "estimated using the {dist} distribution.",
+    compute=generic.frequency_analysis,
+)
+
+
+stats = Stats(
+    identifier="stats",
+    var_name="q{indexer}{op}",
+    long_name="{freq} {op} of {indexer} daily flow ",
+    description="{freq} {op} of {indexer} daily flow",
+    compute=generic.select_resample_op,
+)
+
+
+fit = Fit(
+    identifier="fit",
+    var_name="params",
+    units="",
+    standard_name="{dist} parameters",
+    long_name="{dist} distribution parameters",
+    description="Parameters of the {dist} distribution fitted ",
+    cell_methods="time: fit",
+    compute=generic.fit,
+)
+
+
+doy_qmax = Streamflow(
+    identifier="doy_qmax",
+    var_name="q{indexer}_doy_qmax",
+    long_name="Day of the year of the maximum over {indexer}",
+    description="Day of the year of the maximum over {indexer}",
+    units="",
+    _partial=True,
+    compute=wrapped_partial(generic.select_resample_op, op=generic.doymax),
+)
+
+
+doy_qmin = Streamflow(
+    identifier="doy_qmin",
+    var_name="q{indexer}_doy_qmin",
+    long_name="Day of the year of the minimum over {indexer}",
+    description="Day of the year of the minimum over {indexer}",
+    units="",
+    _partial=True,
+    compute=wrapped_partial(generic.select_resample_op, op=generic.doymin),
+)
diff --git a/xclim/streamflow.py b/xclim/streamflow.py
index 89e0ca98..8f6f1258 100644
--- a/xclim/streamflow.py
+++ b/xclim/streamflow.py
@@ -1,110 +1,16 @@
-# -*- coding: utf-8 -*-
-import numpy as np
-
-from xclim import checks
-from xclim import generic
-from xclim.indices import base_flow_index
-from xclim.utils import Indicator
-from xclim.utils import wrapped_partial
-
-# from boltons.funcutils import FunctionBuilder
-# import calendar
-
-
-class Streamflow(Indicator):
-    units = "m^3 s-1"
-    context = "hydro"
-    standard_name = "discharge"
-
-    @staticmethod
-    def compute(*args, **kwds):
-        pass
-
-    def cfprobe(self, q):
-        checks.check_valid(q, "standard_name", "streamflow")
-
-
-class Stats(Streamflow):
-    def missing(self, *args, **kwds):
-        """Return whether an output is considered missing or not."""
-        from functools import reduce
-
-        indexer = kwds["indexer"]
-        freq = kwds["freq"] or generic.default_freq(**indexer)
-
-        miss = (checks.missing_any(da, freq, **indexer) for da in args)
-        return reduce(np.logical_or, miss)
-
-
-# Disable the missing value check because the output here is not a time series.
-class FA(Streamflow):
-    def missing(self, *args, **kwds):
-        """Return whether an output is considered missing or not."""
-        return False
-
-
-# Disable the daily checks because the inputs are period extremas.
-class Fit(FA):
-    def validate(self, da):
-        pass
-
-
-base_flow_index = Streamflow(
-    identifier="base_flow_index",
-    units="",
-    long_name="Base flow index",
-    compute=base_flow_index,
-)
-
-
-freq_analysis = FA(
-    identifier="freq_analysis",
-    var_name="q{window}{mode}{indexer}",
-    long_name="N-year return period {mode} {indexer} {window}-day flow",
-    description="Streamflow frequency analysis for the {mode} {indexer} {window}-day flow "
-    "estimated using the {dist} distribution.",
-    compute=generic.frequency_analysis,
-)
-
-
-stats = Stats(
-    identifier="stats",
-    var_name="q{indexer}{op}",
-    long_name="{freq} {op} of {indexer} daily flow ",
-    description="{freq} {op} of {indexer} daily flow",
-    compute=generic.select_resample_op,
-)
-
-
-fit = Fit(
-    identifier="fit",
-    var_name="params",
-    units="",
-    standard_name="{dist} parameters",
-    long_name="{dist} distribution parameters",
-    description="Parameters of the {dist} distribution fitted ",
-    cell_methods="time: fit",
-    compute=generic.fit,
-)
-
-
-doy_qmax = Streamflow(
-    identifier="doy_qmax",
-    var_name="q{indexer}_doy_qmax",
-    long_name="Day of the year of the maximum over {indexer}",
-    description="Day of the year of the maximum over {indexer}",
-    units="",
-    _partial=True,
-    compute=wrapped_partial(generic.select_resample_op, op=generic.doymax),
-)
-
-
-doy_qmin = Streamflow(
-    identifier="doy_qmin",
-    var_name="q{indexer}_doy_qmin",
-    long_name="Day of the year of the minimum over {indexer}",
-    description="Day of the year of the minimum over {indexer}",
-    units="",
-    _partial=True,
-    compute=wrapped_partial(generic.select_resample_op, op=generic.doymin),
+import warnings
+
+from .land import base_flow_index
+from .land import doy_qmax
+from .land import doy_qmin
+from .land import fit
+from .land import freq_analysis
+from .land import stats
+
+warnings.warn(
+    "{} will be deprecated in xclim v0.13.x. Please begin using the 'xclim.land' module.".format(
+        __name__,
+    ),
+    DeprecationWarning,
+    stacklevel=2,
 )
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()