diff --git a/.github/workflows/build_docs.yml b/.github/workflows/build_docs.yml index d8589330f..38a48fb69 100644 --- a/.github/workflows/build_docs.yml +++ b/.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 \ No newline at end of file + path: docs/_build/html diff --git a/Makefile b/Makefile index b7823e3b6..f17987469 100644 --- a/Makefile +++ b/Makefile @@ -1,11 +1,23 @@ -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 @@ -13,13 +25,9 @@ doctest: 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: + + diff --git a/docs/dfs2.md b/docs/dfs2.md index 8ed411dbf..3c1530816 100644 --- a/docs/dfs2.md +++ b/docs/dfs2.md @@ -16,6 +16,35 @@ items: 0: Elevation (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 + +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 + +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)) + +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 (meter) +``` ## Grid2D diff --git a/mikeio/dataset/_dataarray.py b/mikeio/dataset/_dataarray.py index b3c35a047..5433b3d22 100644 --- a/mikeio/dataset/_dataarray.py +++ b/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... diff --git a/mikeio/spatial/_FM_geometry.py b/mikeio/spatial/_FM_geometry.py index daca0f11b..ba5b0eafb 100644 --- a/mikeio/spatial/_FM_geometry.py +++ b/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: diff --git a/mikeio/spatial/_FM_utils.py b/mikeio/spatial/_FM_utils.py index 65603bd27..498ba716d 100644 --- a/mikeio/spatial/_FM_utils.py +++ b/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 diff --git a/mikeio/spatial/_grid_geometry.py b/mikeio/spatial/_grid_geometry.py index d18a13248..7bb20bd21 100644 --- a/mikeio/spatial/_grid_geometry.py +++ b/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: diff --git a/tests/test_dataarray.py b/tests/test_dataarray.py index 8f8803560..1f571bc4f 100644 --- a/tests/test_dataarray.py +++ b/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) diff --git a/tests/testdata/consistency/oresundHD.dfs2 b/tests/testdata/consistency/oresundHD.dfs2 index cb141ee58..8107abfab 100644 Binary files a/tests/testdata/consistency/oresundHD.dfs2 and b/tests/testdata/consistency/oresundHD.dfs2 differ diff --git a/tests/testdata/consistency/oresundHD2.dfs2 b/tests/testdata/consistency/oresundHD2.dfs2 deleted file mode 100644 index 8107abfab..000000000 Binary files a/tests/testdata/consistency/oresundHD2.dfs2 and /dev/null differ