From 6650461630eb260dc9d536c1835be3e5bf48c68f Mon Sep 17 00:00:00 2001
From: Kristen Thyng <kthyng@gmail.com>
Date: Thu, 19 Jan 2023 16:55:29 -0600
Subject: [PATCH 1/7] two main changes in sel2d/sel2dcf

1. a mask can be input to limit the lons/lats from the DataArray/Dataset that are used in searching for the nearest point, in case the nearest model point is on land but we still want a valid model point returned.
2. incorporating changes from xoak that optional return distance of the model point(s) from the requested point(s).
Changed tests, etc. Have not updated docs yet.
---
 extract_model/extract_model.py | 45 ++++++++++++++++---
 extract_model/utils.py         | 79 ++++++++++++++++++----------------
 tests/test_accessor.py         |  6 ++-
 tests/test_em.py               | 43 ++++++++++++++++--
 tests/test_sel2d.py            | 24 ++++++-----
 5 files changed, 138 insertions(+), 59 deletions(-)

diff --git a/extract_model/extract_model.py b/extract_model/extract_model.py
index 883c74e..71590d3 100644
--- a/extract_model/extract_model.py
+++ b/extract_model/extract_model.py
@@ -2,6 +2,7 @@
 Main file for this code. The main code is in `select`, and the rest is to help with variable name management.
 """
 
+from typing import Optional
 import warnings
 
 from numbers import Number
@@ -11,6 +12,8 @@
 import xarray as xr
 import xoak  # noqa: F401
 
+from xarray import DataArray
+
 
 try:
     import xesmf as xe
@@ -367,8 +370,8 @@ def select(
     return da.squeeze(), weights
 
 
-def sel2d(var, **kwargs):
-    """Find the value of the var at closest location to inputs.
+def sel2d(var, mask: Optional[DataArray] = None, distances_name: str="distance", **kwargs):
+    """Find the value of the var at closest location to inputs, optionally respecting mask.
 
     This is meant to mimic `xarray` `.sel()` in API and idea, except that the horizontal selection is done for 2D coordinates instead of 1D coordinates, since `xarray` cannot yet handle 2D coordinates. This wraps `xoak`.
 
@@ -380,6 +383,8 @@ def sel2d(var, **kwargs):
     Like in `xarray.sel()`, input values can be numbers, lists, or arrays, and arrays can be single or multidimensional. For longitude and latitude, however, input values cannot be slices.
 
     Can also pass through `xarray.sel()` information for other dimension selections.
+    
+    Optionally input mask so that if requested lon/lat is on land, the nearest valid model point will be returned. Otherwise nan's will be returned. If requested lon/lat is outside domain but not on land, the nearest model output will be returned regardless.
 
     Parameters
     ----------
@@ -389,10 +394,14 @@ def sel2d(var, **kwargs):
         >>> em.sel2d(da, ...)
         instead of `ds.variable` directly. Then subsequent calls will be faster. See `xoak` for more information.
         A Dataset will "remember" the index calculated for whichever grid coordinates were first requested and subsequently run faster for requests on that grid (and not run for other grids).
+    mask : DataArray, optional
+        If input, mask is applied to lon/lat so that if requested lon/lat is on land, the nearest valid model point will be returned. Otherwise nan's will be returned. If requested lon/lat is outside domain but not on land, the nearest model output will be returned regardless.
+    distances_name : str, optional
+        Provide a name in which to save the distances from xoak; there will be one per lon/lat location found. If None, distances won't be returned in object.
 
     Returns
     -------
-    An xarray object of the same type as input as var which is selected in horizontal coordinates to input locations and, in input, to time and vertical selections. If not selected, other dimensions are brought along.
+    An xarray object of the same type as input as var which is selected in horizontal coordinates to input locations and, in input, to time and vertical selections. If not selected, other dimensions are brought along. If distances_name is not None, Dataset is returned.
 
     Notes
     -----
@@ -445,23 +454,47 @@ def sel2d(var, **kwargs):
     # create Dataset
     ds_to_find = xr.Dataset({"lat_to_find": (dims, lats), "lon_to_find": (dims, lons)})
 
+    if mask is not None:
+        
+        # Assume mask is 2D — but not true for wetting/drying
+    
+        # find indices representing mask
+        eta, xi = np.where(mask.values)
+    
+        # make advanced indexer to flatten arrays
+        var_flat = var.cf.isel(X=xr.DataArray(xi, dims="loc"), Y=xr.DataArray(eta, dims="loc"))
+    
+        var = var_flat.copy()
+    
     if var.xoak.index is None:
         var.xoak.set_index([latname, lonname], "sklearn_geo_balltree")
     elif (latname, lonname) != var.xoak._index_coords:
         raise ValueError(
             f"Index has been built for grid with coords {var.xoak._index_coords} but coord names input are ({latname}, {lonname})."
         )
+    elif var.xoak.index is not None:
+        pass
+    else:
+        warnings.warn(
+            "Maybe a mask is not present or being properly identified in var. You could use `use_mask=False`.",
+            RuntimeWarning,
+        )
+
 
     # perform selection
     output = var.xoak.sel(
-        {latname: ds_to_find.lat_to_find, lonname: ds_to_find.lon_to_find}
+        {latname: ds_to_find.lat_to_find, lonname: ds_to_find.lon_to_find}, distances_name=distances_name
     )
+    
+    # distances between input points and nearest points
+    # distances = var.xoak._index.query(np.array([*zip(lats,lons)]))['distances'][:,0]
+    # import pdb; pdb.set_trace()
 
     with xr.set_options(keep_attrs=True):
         return output.sel(**kwargs)
 
 
-def sel2dcf(var, **kwargs):
+def sel2dcf(var, mask: Optional[DataArray] = None, distances_name: str="distance", **kwargs):
     """Find nearest value(s) on 2D horizontal grid using cf-xarray names.
 
     Use "longitude" and "latitude" for those coordinate names.
@@ -511,7 +544,7 @@ def sel2dcf(var, **kwargs):
 
     new_kwargs.update(kwargs)
 
-    return sel2d(var, **new_kwargs)
+    return sel2d(var, mask=mask, distances_name=distances_name, **new_kwargs)
 
 
 def selZ(var, depths):
diff --git a/extract_model/utils.py b/extract_model/utils.py
index e10d06b..5c64ad1 100644
--- a/extract_model/utils.py
+++ b/extract_model/utils.py
@@ -482,7 +482,7 @@ def order(da):
     )
 
 
-def preprocess_roms(ds):
+def preprocess_roms(ds, interp_vertical: bool = True):
     """Preprocess ROMS model output for use with cf-xarray.
 
     Also fixes any other known issues with model output.
@@ -490,12 +490,14 @@ def preprocess_roms(ds):
     Parameters
     ----------
     ds: xarray Dataset
+    
+    interp_vertical=True
 
     Returns
     -------
     Same Dataset but with some metadata added and/or altered.
     """
-
+    
     # add axes attributes for dimensions
     dims = [dim for dim in ds.dims if dim.startswith("s_")]
     for dim in dims:
@@ -546,39 +548,40 @@ def preprocess_roms(ds):
             ds["s_w"].attrs["standard_name"] = "ocean_s_coordinate_g2"
 
     # calculate vertical coord
-    name_dict = {}
-    if "s_rho" in ds.dims:
-        name_dict["s_rho"] = "z_rho"
-        if "positive" in ds.s_rho.attrs:
-            ds.s_rho.attrs.pop("positive")
-    if "s_w" in ds.dims:
-        name_dict["s_w"] = "z_w"
-        if "positive" in ds.s_w.attrs:
-            ds.s_w.attrs.pop("positive")
-    ds.cf.decode_vertical_coords(outnames=name_dict)
-
-    # fix attrs
-    for zname in ["z_rho", "z_w"]:  # name_dict.values():
-        if zname in ds:
-            ds[
-                zname
-            ].attrs = {}  # coord inherits from one of the vars going into calculation
-            ds[zname].attrs["positive"] = "up"
-            ds[zname].attrs["units"] = "m"
-            ds[zname] = order(ds[zname])
-
-    # replace s_rho with z_rho, etc, to make z_rho the vertical coord
-    for sname, zname in name_dict.items():
-        for var in ds.data_vars:
-            if ds[var].ndim == 4:
-                if "coordinates" in ds[var].encoding:
-                    coords = ds[var].encoding["coordinates"]
-                    if sname in coords:  # replace if present
-                        coords = coords.replace(sname, zname)
-                    else:  # still add z_rho or z_w
-                        if zname in ds.coords and ds[zname].shape == ds[var].shape:
-                            coords += f" {zname}"
-                    ds[var].encoding["coordinates"] = coords
+    if interp_vertical:
+        name_dict = {}
+        if "s_rho" in ds.dims:
+            name_dict["s_rho"] = "z_rho"
+            if "positive" in ds.s_rho.attrs:
+                ds.s_rho.attrs.pop("positive")
+        if "s_w" in ds.dims:
+            name_dict["s_w"] = "z_w"
+            if "positive" in ds.s_w.attrs:
+                ds.s_w.attrs.pop("positive")
+        ds.cf.decode_vertical_coords(outnames=name_dict)
+
+        # fix attrs
+        for zname in ["z_rho", "z_w"]:  # name_dict.values():
+            if zname in ds:
+                ds[
+                    zname
+                ].attrs = {}  # coord inherits from one of the vars going into calculation
+                ds[zname].attrs["positive"] = "up"
+                ds[zname].attrs["units"] = "m"
+                ds[zname] = order(ds[zname])
+
+        # replace s_rho with z_rho, etc, to make z_rho the vertical coord
+        for sname, zname in name_dict.items():
+            for var in ds.data_vars:
+                if ds[var].ndim == 4:
+                    if "coordinates" in ds[var].encoding:
+                        coords = ds[var].encoding["coordinates"]
+                        if sname in coords:  # replace if present
+                            coords = coords.replace(sname, zname)
+                        else:  # still add z_rho or z_w
+                            if zname in ds.coords and ds[zname].shape == ds[var].shape:
+                                coords += f" {zname}"
+                        ds[var].encoding["coordinates"] = coords
 
     # # easier to remove "coordinates" attribute from any variables than add it to all
     # for var in ds.data_vars:
@@ -706,12 +709,14 @@ def preprocess_rtofs(ds):
     raise NotImplementedError
 
 
-def preprocess(ds, model_type=None):
+def preprocess(ds, model_type=None, kwargs=None):
     """A preprocess function for reading in with xarray.
 
     This tries to address known model shortcomings in a generic way so that
     `cf-xarray` will work generally, including decoding vertical coordinates.
     """
+    
+    kwargs = kwargs or {}
 
     # This is an internal attribute used by netCDF which xarray doesn't know or care about, but can
     # be returned from THREDDS.
@@ -739,7 +744,7 @@ def preprocess(ds, model_type=None):
         model_type = guess_model_type(ds)
 
     if model_type in preprocess_map:
-        return preprocess_map[model_type](ds)
+        return preprocess_map[model_type](ds, **kwargs)
 
     return ds
 
diff --git a/tests/test_accessor.py b/tests/test_accessor.py
index 474e1f9..5f75fa8 100644
--- a/tests/test_accessor.py
+++ b/tests/test_accessor.py
@@ -19,6 +19,7 @@ def test_2dsel():
     model = models[3]
     da = model["da"]
     i, j = model["i"], model["j"]
+    varname = da.name
 
     if da.cf["longitude"].ndim == 1:
         longitude = float(da.cf["X"][i])
@@ -37,7 +38,8 @@ def test_2dsel():
     da_check = da.cf.isel(X=i, Y=j)
 
     # checks that the resultant model output is the same
-    assert np.allclose(da_sel2d.squeeze(), da_check)
+    assert np.allclose(da_sel2d[varname].squeeze(), da_check)
 
     da_test = da.em.sel2dcf(longitude=lon_comp, latitude=lat_comp)
-    assert np.allclose(da_sel2d, da_test)
+    assert np.allclose(da_sel2d[varname], da_test[varname])
+    assert np.allclose(da_sel2d["distance"], da_test["distance"])
diff --git a/tests/test_em.py b/tests/test_em.py
index d8b0534..f85c078 100644
--- a/tests/test_em.py
+++ b/tests/test_em.py
@@ -75,6 +75,7 @@ def test_sel2d(model):
 
     da = model["da"]
     i, j = model["i"], model["j"]
+    varname = da.name
 
     if da.cf["longitude"].ndim == 1:
         # sel2d is for 2D horizontal grids
@@ -87,8 +88,9 @@ def test_sel2d(model):
         latitude = float(da.cf["latitude"][j, i])
 
         # take a nearby point to test function
-        lon_comp = longitude - 0.001
-        lat_comp = latitude - 0.001
+        lon_comp = longitude
+        dlat = 0.001
+        lat_comp = latitude - dlat
 
         inputs = {
             da.cf["longitude"].name: lon_comp,
@@ -97,7 +99,9 @@ def test_sel2d(model):
         da_sel2d = em.sel2d(da, **inputs)
         da_check = da.cf.isel(X=i, Y=j)
 
-        assert np.allclose(da_sel2d.squeeze(), da_check)
+        assert np.allclose(da_sel2d[varname].squeeze(), da_check)
+        # 6371 is radius of earth in km
+        assert np.allclose(da_sel2d["distance"], np.deg2rad(dlat)*6371)
 
 
 @pytest.mark.parametrize("model", models, ids=lambda x: x["name"])
@@ -329,3 +333,36 @@ def test_preprocess(model):
     conds = [True if axis in da.cf.axes else axis for axis in axes]
 
     assert all(conds)
+
+
+def test_sel2d_simple_2D():
+
+    ds = xr.Dataset(
+        coords={
+            'lon': (('eta','xi'), np.array([[0, 1],[2, 3]]), {"standard_name": "longitude"}),
+            'lat': (('eta','xi'), np.array([[4, 5],[6, 7]]), {"standard_name": "latitude"}),
+            'eta': (('eta'), [0,1], {"axis": "Y"}),
+            'xi': (('xi'), [0,1], {"axis": "X"}),
+        }
+    )
+
+    # check distance when ran with exact grid point
+    ds_out = em.sel2d(ds, lon=0, lat=4)
+    assert np.allclose(ds_out["distance"], 0)
+
+    # check that alternative function call returns exact results
+    ds_outcf = em.sel2dcf(ds, longitude=0, latitude=4)
+    assert ds_out == ds_outcf
+
+    # use mask leaving one valid point to check behavior with mask
+    mask = (ds.cf["longitude"] == 3).astype(int)
+    ds_out = em.sel2d(ds, lon=0, lat=4, mask=mask)
+    assert ds_out.lon == 3
+    assert ds_out.lat == 7
+
+    ds_outcf = em.sel2dcf(ds, longitude=0, latitude=4, mask=mask)
+    assert ds_out == ds_outcf
+    
+    # if distance_name=None, no distance returned
+    ds_out = em.sel2d(ds, lon=0, lat=4, distances_name=None)
+    assert "distance" not in ds_out.variables
diff --git a/tests/test_sel2d.py b/tests/test_sel2d.py
index a9e6020..23560f6 100644
--- a/tests/test_sel2d.py
+++ b/tests/test_sel2d.py
@@ -27,9 +27,9 @@
 i, j = model["i"], model["j"]
 lon = float(da.cf["longitude"][j, i])
 lat = float(da.cf["latitude"][j, i])
+varname = da.name
 
 
-# import pdb; pdb.set_trace()
 def test_lon_lat_types():
     """Test sel2d with different types for lon/lat"""
 
@@ -37,15 +37,15 @@ def test_lon_lat_types():
 
     # Floats
     da_test = em.sel2d(da, lon_rho=lon, lat_rho=lat).squeeze()
-    assert np.allclose(da_check, da_test)
+    assert np.allclose(da_check, da_test[varname])
 
     # List
     da_test = em.sel2d(da, lon_rho=[lon], lat_rho=[lat]).squeeze()
-    assert np.allclose(da_check, da_test)
+    assert np.allclose(da_check, da_test[varname])
 
     # Array
     da_test = em.sel2d(da, lon_rho=np.array([lon]), lat_rho=np.array([lat])).squeeze()
-    assert np.allclose(da_check, da_test)
+    assert np.allclose(da_check, da_test[varname])
 
 
 def test_2D():
@@ -62,7 +62,7 @@ def test_2D():
 
     da_test = em.sel2d(da, lon_rho=Lon, lat_rho=Lat).squeeze()
 
-    assert np.allclose(da_check, da_test)
+    assert np.allclose(da_check, da_test[varname])
 
 
 def test_index_reuse():
@@ -95,7 +95,8 @@ def test_ll_name_reversal():
 
     da1 = em.sel2d(da, lon_rho=lon, lat_rho=lat).squeeze()
     da2 = em.sel2d(da, lat_rho=lat, lon_rho=lon).squeeze()
-    assert np.allclose(da1, da2)
+    assert np.allclose(da1[varname], da2[varname])
+    assert np.allclose(da1["distance"], da2["distance"])
 
 
 def test_sel_time():
@@ -106,7 +107,7 @@ def test_sel_time():
 
     da_test = em.sel2d(da, lon_rho=lon, lat_rho=lat, ocean_time=0)
 
-    assert np.allclose(da_check, da_test)
+    assert np.allclose(da_check, da_test[varname])
 
     # Won't run in different input order
     with pytest.raises(ValueError):
@@ -118,14 +119,15 @@ def test_cf_versions():
 
     da_check = em.sel2d(da, lon_rho=lon, lat_rho=lat)
     da_test = em.sel2dcf(da, longitude=lon, latitude=lat)
-    assert np.allclose(da_check, da_test)
+    assert np.allclose(da_check[varname], da_test[varname])
+    assert np.allclose(da_check["distance"], da_test["distance"])
 
     da_test = em.sel2dcf(da, latitude=lat, longitude=lon)
-    assert np.allclose(da_check, da_test)
+    assert np.allclose(da_check[varname], da_test[varname])
 
     da_check = em.sel2d(da, lon_rho=lon, lat_rho=lat, ocean_time=0)
     da_test = em.sel2dcf(da, latitude=lat, longitude=lon, T=0)
-    assert np.allclose(da_check, da_test)
+    assert np.allclose(da_check[varname], da_test[varname])
 
     da_test = em.sel2dcf(da, T=0, longitude=lon, latitude=lat)
-    assert np.allclose(da_check, da_test)
+    assert np.allclose(da_check[varname], da_test[varname])

From 4cd037699e3a9cd15c978b7c67e2d40284190565 Mon Sep 17 00:00:00 2001
From: Kristen Thyng <kthyng@gmail.com>
Date: Thu, 19 Jan 2023 16:55:54 -0600
Subject: [PATCH 2/7] precommit

---
 extract_model/extract_model.py | 32 +++++++++++++++++++-------------
 extract_model/utils.py         | 10 ++++++----
 tests/test_em.py               | 20 ++++++++++++++------
 3 files changed, 39 insertions(+), 23 deletions(-)

diff --git a/extract_model/extract_model.py b/extract_model/extract_model.py
index 71590d3..e9db785 100644
--- a/extract_model/extract_model.py
+++ b/extract_model/extract_model.py
@@ -2,10 +2,10 @@
 Main file for this code. The main code is in `select`, and the rest is to help with variable name management.
 """
 
-from typing import Optional
 import warnings
 
 from numbers import Number
+from typing import Optional
 
 import cf_xarray  # noqa: F401
 import numpy as np
@@ -370,7 +370,9 @@ def select(
     return da.squeeze(), weights
 
 
-def sel2d(var, mask: Optional[DataArray] = None, distances_name: str="distance", **kwargs):
+def sel2d(
+    var, mask: Optional[DataArray] = None, distances_name: str = "distance", **kwargs
+):
     """Find the value of the var at closest location to inputs, optionally respecting mask.
 
     This is meant to mimic `xarray` `.sel()` in API and idea, except that the horizontal selection is done for 2D coordinates instead of 1D coordinates, since `xarray` cannot yet handle 2D coordinates. This wraps `xoak`.
@@ -383,7 +385,7 @@ def sel2d(var, mask: Optional[DataArray] = None, distances_name: str="distance",
     Like in `xarray.sel()`, input values can be numbers, lists, or arrays, and arrays can be single or multidimensional. For longitude and latitude, however, input values cannot be slices.
 
     Can also pass through `xarray.sel()` information for other dimension selections.
-    
+
     Optionally input mask so that if requested lon/lat is on land, the nearest valid model point will be returned. Otherwise nan's will be returned. If requested lon/lat is outside domain but not on land, the nearest model output will be returned regardless.
 
     Parameters
@@ -455,17 +457,19 @@ def sel2d(var, mask: Optional[DataArray] = None, distances_name: str="distance",
     ds_to_find = xr.Dataset({"lat_to_find": (dims, lats), "lon_to_find": (dims, lons)})
 
     if mask is not None:
-        
+
         # Assume mask is 2D — but not true for wetting/drying
-    
+
         # find indices representing mask
         eta, xi = np.where(mask.values)
-    
+
         # make advanced indexer to flatten arrays
-        var_flat = var.cf.isel(X=xr.DataArray(xi, dims="loc"), Y=xr.DataArray(eta, dims="loc"))
-    
+        var_flat = var.cf.isel(
+            X=xr.DataArray(xi, dims="loc"), Y=xr.DataArray(eta, dims="loc")
+        )
+
         var = var_flat.copy()
-    
+
     if var.xoak.index is None:
         var.xoak.set_index([latname, lonname], "sklearn_geo_balltree")
     elif (latname, lonname) != var.xoak._index_coords:
@@ -480,12 +484,12 @@ def sel2d(var, mask: Optional[DataArray] = None, distances_name: str="distance",
             RuntimeWarning,
         )
 
-
     # perform selection
     output = var.xoak.sel(
-        {latname: ds_to_find.lat_to_find, lonname: ds_to_find.lon_to_find}, distances_name=distances_name
+        {latname: ds_to_find.lat_to_find, lonname: ds_to_find.lon_to_find},
+        distances_name=distances_name,
     )
-    
+
     # distances between input points and nearest points
     # distances = var.xoak._index.query(np.array([*zip(lats,lons)]))['distances'][:,0]
     # import pdb; pdb.set_trace()
@@ -494,7 +498,9 @@ def sel2d(var, mask: Optional[DataArray] = None, distances_name: str="distance",
         return output.sel(**kwargs)
 
 
-def sel2dcf(var, mask: Optional[DataArray] = None, distances_name: str="distance", **kwargs):
+def sel2dcf(
+    var, mask: Optional[DataArray] = None, distances_name: str = "distance", **kwargs
+):
     """Find nearest value(s) on 2D horizontal grid using cf-xarray names.
 
     Use "longitude" and "latitude" for those coordinate names.
diff --git a/extract_model/utils.py b/extract_model/utils.py
index 5c64ad1..1a56885 100644
--- a/extract_model/utils.py
+++ b/extract_model/utils.py
@@ -490,14 +490,14 @@ def preprocess_roms(ds, interp_vertical: bool = True):
     Parameters
     ----------
     ds: xarray Dataset
-    
+
     interp_vertical=True
 
     Returns
     -------
     Same Dataset but with some metadata added and/or altered.
     """
-    
+
     # add axes attributes for dimensions
     dims = [dim for dim in ds.dims if dim.startswith("s_")]
     for dim in dims:
@@ -565,7 +565,9 @@ def preprocess_roms(ds, interp_vertical: bool = True):
             if zname in ds:
                 ds[
                     zname
-                ].attrs = {}  # coord inherits from one of the vars going into calculation
+                ].attrs = (
+                    {}
+                )  # coord inherits from one of the vars going into calculation
                 ds[zname].attrs["positive"] = "up"
                 ds[zname].attrs["units"] = "m"
                 ds[zname] = order(ds[zname])
@@ -715,7 +717,7 @@ def preprocess(ds, model_type=None, kwargs=None):
     This tries to address known model shortcomings in a generic way so that
     `cf-xarray` will work generally, including decoding vertical coordinates.
     """
-    
+
     kwargs = kwargs or {}
 
     # This is an internal attribute used by netCDF which xarray doesn't know or care about, but can
diff --git a/tests/test_em.py b/tests/test_em.py
index f85c078..4ba11ef 100644
--- a/tests/test_em.py
+++ b/tests/test_em.py
@@ -101,7 +101,7 @@ def test_sel2d(model):
 
         assert np.allclose(da_sel2d[varname].squeeze(), da_check)
         # 6371 is radius of earth in km
-        assert np.allclose(da_sel2d["distance"], np.deg2rad(dlat)*6371)
+        assert np.allclose(da_sel2d["distance"], np.deg2rad(dlat) * 6371)
 
 
 @pytest.mark.parametrize("model", models, ids=lambda x: x["name"])
@@ -339,10 +339,18 @@ def test_sel2d_simple_2D():
 
     ds = xr.Dataset(
         coords={
-            'lon': (('eta','xi'), np.array([[0, 1],[2, 3]]), {"standard_name": "longitude"}),
-            'lat': (('eta','xi'), np.array([[4, 5],[6, 7]]), {"standard_name": "latitude"}),
-            'eta': (('eta'), [0,1], {"axis": "Y"}),
-            'xi': (('xi'), [0,1], {"axis": "X"}),
+            "lon": (
+                ("eta", "xi"),
+                np.array([[0, 1], [2, 3]]),
+                {"standard_name": "longitude"},
+            ),
+            "lat": (
+                ("eta", "xi"),
+                np.array([[4, 5], [6, 7]]),
+                {"standard_name": "latitude"},
+            ),
+            "eta": (("eta"), [0, 1], {"axis": "Y"}),
+            "xi": (("xi"), [0, 1], {"axis": "X"}),
         }
     )
 
@@ -362,7 +370,7 @@ def test_sel2d_simple_2D():
 
     ds_outcf = em.sel2dcf(ds, longitude=0, latitude=4, mask=mask)
     assert ds_out == ds_outcf
-    
+
     # if distance_name=None, no distance returned
     ds_out = em.sel2d(ds, lon=0, lat=4, distances_name=None)
     assert "distance" not in ds_out.variables

From 14bfa4b8d980d667b570a02675833d636f26a6ca Mon Sep 17 00:00:00 2001
From: Kristen Thyng <kthyng@gmail.com>
Date: Thu, 19 Jan 2023 16:58:04 -0600
Subject: [PATCH 3/7] updated whats new

---
 docs/whats_new.rst | 8 ++++++++
 1 file changed, 8 insertions(+)

diff --git a/docs/whats_new.rst b/docs/whats_new.rst
index 72d6ed5..ba11787 100644
--- a/docs/whats_new.rst
+++ b/docs/whats_new.rst
@@ -1,6 +1,14 @@
 :mod:`What's New`
 ----------------------------
 
+v1.1.0 (January 19, 2023)
+=========================
+
+Two main changes in `sel2d` / `sel2dcf`:
+
+* a mask can be input to limit the lons/lats from the DataArray/Dataset that are used in searching for the nearest point, in case the nearest model point is on land but we still want a valid model point returned.
+* incorporating changes from xoak that optional return distance of the model point(s) from the requested point(s).
+
 v1.0.0 (December 9, 2022)
 =========================
 * Simplified dependencies

From 497e256e23da31b8b9a7d680c55fdf97d104ca6b Mon Sep 17 00:00:00 2001
From: Kristen Thyng <kthyng@gmail.com>
Date: Thu, 19 Jan 2023 17:03:33 -0600
Subject: [PATCH 4/7] fix buld problem

https://github.com/conda-incubator/setup-miniconda/issues/274#issuecomment-1384764160
---
 .github/workflows/test.yaml | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/.github/workflows/test.yaml b/.github/workflows/test.yaml
index 09512dc..d1cbc7a 100644
--- a/.github/workflows/test.yaml
+++ b/.github/workflows/test.yaml
@@ -40,7 +40,8 @@ jobs:
         if: ${{ runner.os != 'Windows' }}
         uses: conda-incubator/setup-miniconda@v2
         with:
-          mamba-version: "*" # activate this to build with mamba.
+          miniforge-variant: Mambaforge
+          python-version: ${{ matrix.python-version }}
           channels: conda-forge, defaults # These need to be specified to use mamba
           channel-priority: true
           environment-file: ci/environment-py${{ matrix.python-version }}.yml

From d50f10447bece3ca878d66bab4206428278b74f8 Mon Sep 17 00:00:00 2001
From: Kristen Thyng <kthyng@gmail.com>
Date: Thu, 19 Jan 2023 17:57:14 -0600
Subject: [PATCH 5/7] updating builds

---
 .github/workflows/test.yaml                                  | 5 +++--
 ci/{environment-py3.7-win.yml => environment-py3.10-win.yml} | 2 +-
 ci/{environment-py3.7.yml => environment-py3.10.yml}         | 2 +-
 3 files changed, 5 insertions(+), 4 deletions(-)
 rename ci/{environment-py3.7-win.yml => environment-py3.10-win.yml} (94%)
 rename ci/{environment-py3.7.yml => environment-py3.10.yml} (96%)

diff --git a/.github/workflows/test.yaml b/.github/workflows/test.yaml
index d1cbc7a..2831950 100644
--- a/.github/workflows/test.yaml
+++ b/.github/workflows/test.yaml
@@ -9,7 +9,7 @@ jobs:
       fail-fast: false
       matrix:
         os: ["macos-latest", "ubuntu-latest", "windows-latest"]
-        python-version: ["3.7", "3.8", "3.9"]
+        python-version: ["3.8", "3.9", "3.10"]
     steps:
       - name: Checkout source
         uses: actions/checkout@v3
@@ -53,7 +53,8 @@ jobs:
         if: ${{ runner.os == 'Windows' }}
         uses: conda-incubator/setup-miniconda@v2
         with:
-          mamba-version: "*" # activate this to build with mamba.
+          miniforge-variant: Mambaforge
+          python-version: ${{ matrix.python-version }}
           channels: conda-forge, defaults # These need to be specified to use mamba
           channel-priority: true
           environment-file: ci/environment-py${{ matrix.python-version }}-win.yml
diff --git a/ci/environment-py3.7-win.yml b/ci/environment-py3.10-win.yml
similarity index 94%
rename from ci/environment-py3.7-win.yml
rename to ci/environment-py3.10-win.yml
index 5785c5d..324ec49 100644
--- a/ci/environment-py3.7-win.yml
+++ b/ci/environment-py3.10-win.yml
@@ -2,7 +2,7 @@ name: test-env-win
 channels:
   - conda-forge
 dependencies:
-  - python=3.7
+  - python=3.10
   - cf_xarray>=0.6
   - dask
   - netcdf4
diff --git a/ci/environment-py3.7.yml b/ci/environment-py3.10.yml
similarity index 96%
rename from ci/environment-py3.7.yml
rename to ci/environment-py3.10.yml
index f7eaba6..f65387b 100644
--- a/ci/environment-py3.7.yml
+++ b/ci/environment-py3.10.yml
@@ -2,7 +2,7 @@ name: test-env-mac-unix
 channels:
   - conda-forge
 dependencies:
-  - python=3.7
+  - python=3.10
   - cf_xarray>=0.6
   - dask <=2022.05.0 # for xESMF, https://github.com/axiom-data-science/extract_model/issues/49
   - netcdf4

From 8179d7d8b7efa49a8b5cb2404f517cd199d736ef Mon Sep 17 00:00:00 2001
From: Kristen Thyng <kthyng@gmail.com>
Date: Thu, 19 Jan 2023 18:04:41 -0600
Subject: [PATCH 6/7] xoak needs to be dev version

---
 ci/environment-py3.10-win.yml | 3 ++-
 ci/environment-py3.10.yml     | 3 ++-
 ci/environment-py3.8-win.yml  | 3 ++-
 ci/environment-py3.8.yml      | 3 ++-
 ci/environment-py3.9-win.yml  | 3 ++-
 ci/environment-py3.9.yml      | 3 ++-
 6 files changed, 12 insertions(+), 6 deletions(-)

diff --git a/ci/environment-py3.10-win.yml b/ci/environment-py3.10-win.yml
index 324ec49..31006e6 100644
--- a/ci/environment-py3.10-win.yml
+++ b/ci/environment-py3.10-win.yml
@@ -11,10 +11,11 @@ dependencies:
   - requests
   - scikit-learn  # used by xoak for tree
   - xarray
-  - xoak
+  # - xoak
   - pytest
   - pytest-benchmark
   - pip:
     - codecov
     - pytest-cov
     - coverage[toml]
+    - git+https://github.com/xarray-contrib/xoak
diff --git a/ci/environment-py3.10.yml b/ci/environment-py3.10.yml
index f65387b..3c5d013 100644
--- a/ci/environment-py3.10.yml
+++ b/ci/environment-py3.10.yml
@@ -12,10 +12,11 @@ dependencies:
   - scikit-learn  # used by xoak for tree
   - xarray
   - xesmf
-  - xoak
+  # - xoak
   - pytest
   - pytest-benchmark
   - pip:
     - codecov
     - pytest-cov
     - coverage[toml]
+    - git+https://github.com/xarray-contrib/xoak
diff --git a/ci/environment-py3.8-win.yml b/ci/environment-py3.8-win.yml
index 4b581fe..de2f319 100644
--- a/ci/environment-py3.8-win.yml
+++ b/ci/environment-py3.8-win.yml
@@ -11,10 +11,11 @@ dependencies:
   - requests
   - scikit-learn  # used by xoak for tree
   - xarray
-  - xoak
+  # - xoak
   - pytest
   - pytest-benchmark
   - pip:
     - codecov
     - pytest-cov
     - coverage[toml]
+    - git+https://github.com/xarray-contrib/xoak
diff --git a/ci/environment-py3.8.yml b/ci/environment-py3.8.yml
index c2e6026..7487908 100644
--- a/ci/environment-py3.8.yml
+++ b/ci/environment-py3.8.yml
@@ -12,10 +12,11 @@ dependencies:
   - scikit-learn  # used by xoak for tree
   - xarray
   - xesmf
-  - xoak
+  # - xoak
   - pytest
   - pytest-benchmark
   - pip:
     - codecov
     - pytest-cov
     - coverage[toml]
+    - git+https://github.com/xarray-contrib/xoak
diff --git a/ci/environment-py3.9-win.yml b/ci/environment-py3.9-win.yml
index f683edc..0b13f34 100644
--- a/ci/environment-py3.9-win.yml
+++ b/ci/environment-py3.9-win.yml
@@ -11,10 +11,11 @@ dependencies:
   - requests
   - scikit-learn  # used by xoak for tree
   - xarray
-  - xoak
+  # - xoak
   - pytest
   - pytest-benchmark
   - pip:
     - codecov
     - pytest-cov
     - coverage[toml]
+    - git+https://github.com/xarray-contrib/xoak
diff --git a/ci/environment-py3.9.yml b/ci/environment-py3.9.yml
index f7b6721..c118a46 100644
--- a/ci/environment-py3.9.yml
+++ b/ci/environment-py3.9.yml
@@ -12,10 +12,11 @@ dependencies:
   - scikit-learn  # used by xoak for tree
   - xarray
   - xesmf
-  - xoak
+  # - xoak
   - pytest
   - pytest-benchmark
   - pip:
     - codecov
     - pytest-cov
     - coverage[toml]
+    - git+https://github.com/xarray-contrib/xoak

From c89824a6d6767704256d6e4aeed7d065d973e677 Mon Sep 17 00:00:00 2001
From: Kristen Thyng <kthyng@gmail.com>
Date: Thu, 19 Jan 2023 18:07:44 -0600
Subject: [PATCH 7/7] need my branch of xoak for now

---
 ci/environment-py3.10-win.yml | 2 +-
 ci/environment-py3.10.yml     | 2 +-
 ci/environment-py3.8-win.yml  | 2 +-
 ci/environment-py3.8.yml      | 2 +-
 ci/environment-py3.9-win.yml  | 2 +-
 ci/environment-py3.9.yml      | 2 +-
 6 files changed, 6 insertions(+), 6 deletions(-)

diff --git a/ci/environment-py3.10-win.yml b/ci/environment-py3.10-win.yml
index 31006e6..a58d4aa 100644
--- a/ci/environment-py3.10-win.yml
+++ b/ci/environment-py3.10-win.yml
@@ -18,4 +18,4 @@ dependencies:
     - codecov
     - pytest-cov
     - coverage[toml]
-    - git+https://github.com/xarray-contrib/xoak
+    - git+https://github.com/kthyng/xoak@include_distances
diff --git a/ci/environment-py3.10.yml b/ci/environment-py3.10.yml
index 3c5d013..3d24938 100644
--- a/ci/environment-py3.10.yml
+++ b/ci/environment-py3.10.yml
@@ -19,4 +19,4 @@ dependencies:
     - codecov
     - pytest-cov
     - coverage[toml]
-    - git+https://github.com/xarray-contrib/xoak
+    - git+https://github.com/kthyng/xoak@include_distances
diff --git a/ci/environment-py3.8-win.yml b/ci/environment-py3.8-win.yml
index de2f319..ebd3c14 100644
--- a/ci/environment-py3.8-win.yml
+++ b/ci/environment-py3.8-win.yml
@@ -18,4 +18,4 @@ dependencies:
     - codecov
     - pytest-cov
     - coverage[toml]
-    - git+https://github.com/xarray-contrib/xoak
+    - git+https://github.com/kthyng/xoak@include_distances
diff --git a/ci/environment-py3.8.yml b/ci/environment-py3.8.yml
index 7487908..ee6f036 100644
--- a/ci/environment-py3.8.yml
+++ b/ci/environment-py3.8.yml
@@ -19,4 +19,4 @@ dependencies:
     - codecov
     - pytest-cov
     - coverage[toml]
-    - git+https://github.com/xarray-contrib/xoak
+    - git+https://github.com/kthyng/xoak@include_distances
diff --git a/ci/environment-py3.9-win.yml b/ci/environment-py3.9-win.yml
index 0b13f34..c3dcfde 100644
--- a/ci/environment-py3.9-win.yml
+++ b/ci/environment-py3.9-win.yml
@@ -18,4 +18,4 @@ dependencies:
     - codecov
     - pytest-cov
     - coverage[toml]
-    - git+https://github.com/xarray-contrib/xoak
+    - git+https://github.com/kthyng/xoak@include_distances
diff --git a/ci/environment-py3.9.yml b/ci/environment-py3.9.yml
index c118a46..f799cb0 100644
--- a/ci/environment-py3.9.yml
+++ b/ci/environment-py3.9.yml
@@ -19,4 +19,4 @@ dependencies:
     - codecov
     - pytest-cov
     - coverage[toml]
-    - git+https://github.com/xarray-contrib/xoak
+    - git+https://github.com/kthyng/xoak@include_distances