Merge pull request #604 from DHI/dfs2_slice

Dfs2 subsetting with physical coordinates

.github/workflows/build_docs.yml

@@ -1,8 +1,7 @@
 name: Build documentation (don't publish)
 
 on:
-  push:
-    pull_request:
+  pull_request:
       branches: [ main ]
 
 jobs:
@@ -11,9 +10,9 @@ jobs:
     runs-on: ubuntu-latest
 
     steps:
-    - uses: actions/checkout@v2
+    - uses: actions/checkout@v3
     - name: Set up Python 3.10
-      uses: actions/setup-python@v2
+      uses: actions/setup-python@v4
       with:
         python-version: "3.10"
 
@@ -30,4 +29,4 @@ jobs:
       uses: actions/upload-artifact@v2
       with:
         name: html
-        path: docs/_build/html
+        path: docs/_build/html

Makefile

@@ -1,25 +1,33 @@
 
-build: test
-	#python setup.py sdist bdist_wheel
+LIB = mikeio
+
+check: lint typecheck test
+
+build: typecheck test
 	python -m build
 
+lint:
+	ruff .
+
 test:
 	pytest --disable-warnings
 
+typecheck:
+	mypy $(LIB)/ --config-file pyproject.toml
+
+coverage: 
+	pytest --cov-report html --cov=$(LIB) tests/
+
 doctest:
 	pytest mikeio/dfs/*.py mikeio/dfsu/*.py mikeio/eum/*.py mikeio/pfs/*.py mikeio/spatial/_grid_geometry.py --doctest-modules
 	rm -f *.dfs* # remove temporary files, created from doctests
 
 perftest:
 	pytest tests/performance/ --durations=0
 
-typecheck:
-	mypy mikeio/
-
-coverage: 
-	pytest --cov-report html --cov=mikeio tests/
-
 docs: FORCE
 	cd docs; make html ;cd -
 
 FORCE:
+
+

docs/dfs2.md

@@ -16,6 +16,35 @@ items:
   0:  Elevation <Total Water Depth> (meter)
 ```
 
+## Subset in space
+
+The most convenient way to subset in space is to use the `sel` method, which returns a new (smaller) dataset, which can be further processed or written to disk using the `to_dfs` method. 
+
+```python
+>>> ds.geometry
+<mikeio.Grid2D>
+x: [12.2, 12.21, ..., 13.1] (nx=216, dx=0.004167)
+y: [55.2, 55.21, ..., 56.3] (ny=264, dy=0.004167)
+projection: LONG/LAT
+>>> ds_aoi = ds.sel(x=slice(12.5, 13.0), y=slice(55.5, 56.0))
+>>> ds_aoi.geometry
+<mikeio.Grid2D>
+x: [12.5, 12.5, ..., 12.99] (nx=120, dx=0.004167)
+y: [55.5, 55.5, ..., 55.99] (ny=120, dy=0.004167)
+projection: LONG/LAT
+```
+
+In order to specify an open-ended subset (i.e. where the end of the subset is the end of the domain), use `None` as the end of the slice. 
+
+```python
+>>> ds.sel(x=slice(None,13.0))
+<mikeio.Dataset>
+dims: (time:1, y:264, x:191)
+time: 2020-05-15 11:04:52 (time-invariant)
+geometry: Grid2D (ny=264, nx=191)
+items:
+  0:  Elevation <Total Water Depth> (meter)
+```
 
 ## Grid2D
 

mikeio/dataset/_dataarray.py

@@ -3,7 +3,7 @@
 from copy import deepcopy
 from datetime import datetime
 from functools import cached_property
-from typing import Iterable, Optional, Sequence, Tuple
+from typing import Iterable, Optional, Sequence, Tuple, Mapping
 
 
 import numpy as np
@@ -839,10 +839,14 @@ def sel(
         time: 1997-09-15 21:00:00 - 1997-09-16 03:00:00 (3 records)
         geometry: Dfsu2D (3700 elements, 2090 nodes)
         """
+        if any([isinstance(v, slice) for v in kwargs.values()]):
+            return self._sel_with_slice(kwargs)
+
         da = self
 
         # select in space
         if len(kwargs) > 0:
+
             idx = self.geometry.find_index(**kwargs)
             if isinstance(idx, tuple):
                 # TODO: support for dfs3
@@ -866,6 +870,27 @@ def sel(
             da = da[time]  # __getitem__ is 🚀
 
         return da
+
+    def _sel_with_slice(self, kwargs: Mapping[str,slice]) -> "DataArray":
+        for k, v in kwargs.items():
+            if isinstance(v, slice):
+                idx_start = self.geometry.find_index(**{k:v.start})
+                idx_stop = self.geometry.find_index(**{k:v.stop})
+                pos = 0
+                if isinstance(idx_start, tuple):
+                    if k == "x":
+                        pos = 0
+                    if k == "y":
+                        pos = 1
+
+                start = idx_start[pos][0] if idx_start is not None else None
+                stop = idx_stop[pos][0] if idx_stop is not None else None
+
+                idx = slice(start, stop)
+
+                self = self.isel(idx, axis=k)
+
+        return self
 
     def interp(
         # TODO find out optimal syntax to allow interpolation to single point, new time, grid, mesh...

mikeio/spatial/_FM_geometry.py

@@ -2,7 +2,7 @@
 import warnings
 from collections import namedtuple
 from functools import cached_property
-from typing import Collection, Sequence, Optional, List
+from typing import Collection, Optional, List
 
 import numpy as np
 from mikecore.DfsuFile import DfsuFileType  # type: ignore
@@ -830,7 +830,7 @@ def boundary_polylines(self) -> BoundaryPolylines:
         """Lists of closed polylines defining domain outline"""
         return self._get_boundary_polylines()
 
-    def contains(self, points) -> Sequence[bool]:
+    def contains(self, points):
         """test if a list of points are contained by mesh
 
         Parameters
@@ -1037,7 +1037,7 @@ def find_index(self, x=None, y=None, coords=None, area=None) -> np.ndarray:
             raise ValueError("Provide either coordinates or area")
 
     @staticmethod
-    def _inside_polygon(polygon, xy) -> bool:
+    def _inside_polygon(polygon, xy):
         import matplotlib.path as mp
 
         if polygon.ndim == 1:

mikeio/spatial/_FM_utils.py

@@ -97,7 +97,7 @@ def _plot_map(
     """
 
     import matplotlib.pyplot as plt
-    import matplotlib.cm as cm
+    import matplotlib
 
     VALID_PLOT_TYPES = (
         "mesh_only",
@@ -111,7 +111,7 @@ def _plot_map(
         ok_list = ", ".join(VALID_PLOT_TYPES)
         raise Exception(f"plot_type {plot_type} unknown! ({ok_list})")
 
-    cmap = cmap or cm.viridis
+    cmap = cmap or matplotlib.colormaps["viridis"]
 
     nc = node_coordinates
     ec = element_coordinates

mikeio/spatial/_grid_geometry.py

@@ -744,8 +744,12 @@ def find_index(
                 raise ValueError("x,y and coords cannot be given at the same time!")
             coords = np.column_stack([np.atleast_1d(x), np.atleast_1d(y)])
         elif x is not None:
+            if x < self.x[0] or x > self.x[-1]:
+                raise OutsideModelDomainError(x=x, y=None)
             return np.atleast_1d(np.argmin(np.abs(self.x - x))), None
         elif y is not None:
+            if y < self.y[0] or y > self.y[-1]:
+                raise OutsideModelDomainError(x=None, y=y)
             return None, np.atleast_1d(np.argmin(np.abs(self.y - y)))
 
         if coords is not None:

tests/test_dataarray.py

@@ -851,6 +851,25 @@ def test_modify_values_1d(da1):
     da1.isel([0, 4, 7]).values[1] = 10.0
     assert da1.values[4] != 10.0
 
+def test_get_2d_slice_with_sel(da_grid2d):
+    assert da_grid2d.shape == (10, 14, 7)
+    da3 = da_grid2d.sel(x=slice(10.0, 10.3))
+    assert da3.shape == (10, 14,3)
+    da4 = da_grid2d.sel(y=slice(-5.0, 0.0))
+    assert da4.shape == (10, 5, 7)
+
+    da5 = da_grid2d.sel(x=slice(10.0, 10.3), y=slice(-5.0,0.0))
+    assert da5.shape == (10,5,3)
+
+    da6 = da_grid2d.sel(x=slice(None, 10.3), y=slice(-4.0, None))
+    assert da6.shape == (10, 8, 3)
+
+def test_get_2d_outside_domain_raises_error(da_grid2d):
+    with pytest.raises(OutsideModelDomainError):
+        da_grid2d.sel(x=0.0)
+
+    with pytest.raises(OutsideModelDomainError):
+        da_grid2d.sel(x=slice(0.0,1.0))
 
 def test_modify_values_2d_all(da2):
     assert da2.shape == (10, 7)