From 35b8cf0b768fad7b49006582b1b76469cefdd433 Mon Sep 17 00:00:00 2001 From: Henrik Andersson Date: Mon, 13 May 2024 14:34:00 +0200 Subject: [PATCH 1/4] Convert dfsu2d to UGRID netcdf --- mikeio/dataset/_dataarray.py | 34 ++++++++++++++++++++++++++++++++++ mikeio/dataset/_dataset.py | 9 +++++++++ mypy.ini | 3 +++ pyproject.toml | 2 +- tests/test_dfsu.py | 20 ++++++++++++++++++++ 5 files changed, 67 insertions(+), 1 deletion(-) diff --git a/mikeio/dataset/_dataarray.py b/mikeio/dataset/_dataarray.py index e15a9b1bc..82e2a9cb3 100644 --- a/mikeio/dataset/_dataarray.py +++ b/mikeio/dataset/_dataarray.py @@ -1873,6 +1873,40 @@ def to_xarray(self) -> "xarray.DataArray": ) return xr_da + def to_uxarray(self): + import xugrid as xu + import xarray as xr + + assert isinstance( + self.geometry, GeometryFM2D + ), "Only flexible mesh 2D data is supported" + # face_code_connectivity is slightly different than the element_table + conn = np.zeros(shape=(self.geometry.n_elements, 4), dtype=np.int32) + conn[:] = -1 + + g = self.geometry + + for i in range(g.n_elements): + if len(g.element_table[i]) == 3: + conn[i, 0:3] = g.element_table[i] + else: + conn[i] = g.element_table[i] + + xn = g.element_coordinates[:, 0] + yn = g.element_coordinates[:, 1] + + grid = xu.Ugrid2d(xn, yn, fill_value=-1, face_node_connectivity=conn) + + da = xr.DataArray( + name=self.name, + data=self.to_numpy(), + coords={"time": self.time}, + dims=["time", grid.face_dimension], + ) + + uda = xu.UgridDataArray(da, grid) + return uda + # =============================================== def __repr__(self) -> str: diff --git a/mikeio/dataset/_dataset.py b/mikeio/dataset/_dataset.py index 6da0ff971..03319a4ff 100644 --- a/mikeio/dataset/_dataset.py +++ b/mikeio/dataset/_dataset.py @@ -1898,6 +1898,15 @@ def to_xarray(self) -> "xarray.Dataset": data = {da.name: da.to_xarray() for da in self} return xarray.Dataset(data) + def to_xugrid(self): + import xugrid as xu + + data = [da.to_uxarray() for da in self] + ds = xu.UgridDataset(grids=data[0].grid) + for da in data: + ds[da.name] = da + return ds + # =============================================== def __repr__(self) -> str: diff --git a/mypy.ini b/mypy.ini index f91b72461..12ccc3ab6 100644 --- a/mypy.ini +++ b/mypy.ini @@ -33,4 +33,7 @@ ignore_missing_imports = True [mypy-yaml.*] ignore_missing_imports = True +[mypy-xugrid.*] +ignore_missing_imports = True + diff --git a/pyproject.toml b/pyproject.toml index 1539327c4..836e15ba3 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -58,7 +58,7 @@ dev = ["pytest", "mypy==1.6.1", ] -test = ["pytest", "pytest-cov", "xarray","mypy==1.6.1","shapely","pyproj"] +test = ["pytest", "pytest-cov", "xarray","mypy==1.6.1","shapely","pyproj","uxgrid"] notebooks= [ "nbformat", diff --git a/tests/test_dfsu.py b/tests/test_dfsu.py index 678f82af3..ef2762af7 100644 --- a/tests/test_dfsu.py +++ b/tests/test_dfsu.py @@ -960,3 +960,23 @@ def test_writing_non_equdistant_dfsu_is_not_possible(tmp_path): with pytest.raises(ValueError, match="equidistant"): fp = tmp_path / "not_gonna_work.dfsu" dss.to_dfs(fp) + + +def test_convert_dfsu2d_to_xugrid(tmp_path): + + import xugrid + + ds = mikeio.read("tests/testdata/NorthSea_HD_and_windspeed.dfsu") + ds.sel(y=55, x=0).isel(time=-1).to_numpy() == pytest.approx(0.25443718) + xu_ds = ds.to_xugrid() + + xu_ds.ugrid.sel(y=55, x=0)["Surface elevation"].isel( + time=-1 + ).to_numpy() == pytest.approx(0.25443718) + assert xu_ds.time[-1] == ds.time[-1] + + fp = tmp_path / "north_sea.nc" + xu_ds.ugrid.to_netcdf(fp) + + xu_ds2 = xugrid.open_dataset(fp) + assert xu_ds2.time[-1] == ds.time[-1] From d2159a7ae2e9b1cffa8e83f5875605746969f91c Mon Sep 17 00:00:00 2001 From: Henrik Andersson Date: Mon, 13 May 2024 14:35:58 +0200 Subject: [PATCH 2/4] it is called xugrid --- mikeio/dataset/_dataarray.py | 2 ++ pyproject.toml | 2 +- 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/mikeio/dataset/_dataarray.py b/mikeio/dataset/_dataarray.py index 82e2a9cb3..b7deaf4fa 100644 --- a/mikeio/dataset/_dataarray.py +++ b/mikeio/dataset/_dataarray.py @@ -1897,6 +1897,8 @@ def to_uxarray(self): grid = xu.Ugrid2d(xn, yn, fill_value=-1, face_node_connectivity=conn) + # TODO support data with time axis + da = xr.DataArray( name=self.name, data=self.to_numpy(), diff --git a/pyproject.toml b/pyproject.toml index 836e15ba3..4572c85d5 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -58,7 +58,7 @@ dev = ["pytest", "mypy==1.6.1", ] -test = ["pytest", "pytest-cov", "xarray","mypy==1.6.1","shapely","pyproj","uxgrid"] +test = ["pytest", "pytest-cov", "xarray","mypy==1.6.1","shapely","pyproj","xugrid"] notebooks= [ "nbformat", From 846155e4861bc6c612f9cbd9879408bfc52c8248 Mon Sep 17 00:00:00 2001 From: Henrik Andersson Date: Mon, 13 May 2024 14:39:31 +0200 Subject: [PATCH 3/4] Try 3.9 --- .github/workflows/full_test.yml | 2 +- .github/workflows/legacy_test.yml | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/full_test.yml b/.github/workflows/full_test.yml index 26a261d3a..256d06d7e 100644 --- a/.github/workflows/full_test.yml +++ b/.github/workflows/full_test.yml @@ -17,7 +17,7 @@ jobs: strategy: matrix: os: [ubuntu-latest, windows-latest] - python-version: [3.8, "3.12"] + python-version: [3.9, "3.12"] steps: - uses: actions/checkout@v3 diff --git a/.github/workflows/legacy_test.yml b/.github/workflows/legacy_test.yml index 63f746dce..77ac50f87 100644 --- a/.github/workflows/legacy_test.yml +++ b/.github/workflows/legacy_test.yml @@ -17,7 +17,7 @@ jobs: - name: Set up Python uses: actions/setup-python@v4 with: - python-version: 3.8 + python-version: 3.9 - name: Install older dependencies run: | From f33eb25c0dd642c7a532e09391ccf1b3d8776fe6 Mon Sep 17 00:00:00 2001 From: Henrik Andersson Date: Mon, 13 May 2024 15:30:17 +0200 Subject: [PATCH 4/4] Use node coords --- mikeio/dataset/_dataarray.py | 46 +++++++++++++++++++++++------------- mikeio/dataset/_dataset.py | 2 +- tests/test_dfsu.py | 6 +++++ 3 files changed, 36 insertions(+), 18 deletions(-) diff --git a/mikeio/dataset/_dataarray.py b/mikeio/dataset/_dataarray.py index b7deaf4fa..a03486c0b 100644 --- a/mikeio/dataset/_dataarray.py +++ b/mikeio/dataset/_dataarray.py @@ -1873,7 +1873,7 @@ def to_xarray(self) -> "xarray.DataArray": ) return xr_da - def to_uxarray(self): + def to_xugrid(self): import xugrid as xu import xarray as xr @@ -1881,30 +1881,42 @@ def to_uxarray(self): self.geometry, GeometryFM2D ), "Only flexible mesh 2D data is supported" # face_code_connectivity is slightly different than the element_table - conn = np.zeros(shape=(self.geometry.n_elements, 4), dtype=np.int32) - conn[:] = -1 g = self.geometry - for i in range(g.n_elements): - if len(g.element_table[i]) == 3: + if g.is_tri_only: + conn = np.zeros(shape=(self.geometry.n_elements, 3), dtype=np.int32) + for i in range(g.n_elements): conn[i, 0:3] = g.element_table[i] - else: - conn[i] = g.element_table[i] + else: + conn = np.zeros(shape=(self.geometry.n_elements, 4), dtype=np.int32) + conn[:] = -1 - xn = g.element_coordinates[:, 0] - yn = g.element_coordinates[:, 1] + for i in range(g.n_elements): + if len(g.element_table[i]) == 3: + conn[i, 0:3] = g.element_table[i] + else: + conn[i] = g.element_table[i] - grid = xu.Ugrid2d(xn, yn, fill_value=-1, face_node_connectivity=conn) + xn = g.node_coordinates[:, 0] + yn = g.node_coordinates[:, 1] - # TODO support data with time axis + grid = xu.Ugrid2d(xn, yn, fill_value=-1, face_node_connectivity=conn) - da = xr.DataArray( - name=self.name, - data=self.to_numpy(), - coords={"time": self.time}, - dims=["time", grid.face_dimension], - ) + if "time" not in self.dims: + da = xr.DataArray( + name=self.name, + data=self.to_numpy(), + coords={}, + dims=[grid.face_dimension], + ) + else: + da = xr.DataArray( + name=self.name, + data=self.to_numpy(), + coords={"time": self.time}, + dims=["time", grid.face_dimension], + ) uda = xu.UgridDataArray(da, grid) return uda diff --git a/mikeio/dataset/_dataset.py b/mikeio/dataset/_dataset.py index 03319a4ff..7283ff8b1 100644 --- a/mikeio/dataset/_dataset.py +++ b/mikeio/dataset/_dataset.py @@ -1901,7 +1901,7 @@ def to_xarray(self) -> "xarray.Dataset": def to_xugrid(self): import xugrid as xu - data = [da.to_uxarray() for da in self] + data = [da.to_xugrid() for da in self] ds = xu.UgridDataset(grids=data[0].grid) for da in data: ds[da.name] = da diff --git a/tests/test_dfsu.py b/tests/test_dfsu.py index ef2762af7..9d9826018 100644 --- a/tests/test_dfsu.py +++ b/tests/test_dfsu.py @@ -980,3 +980,9 @@ def test_convert_dfsu2d_to_xugrid(tmp_path): xu_ds2 = xugrid.open_dataset(fp) assert xu_ds2.time[-1] == ds.time[-1] + + +def test_mixed_mesh_time_invariant_to_xugrid(): + ds = mikeio.read("tests/testdata/FakeLake.dfsu", time=0) + xu_ds = ds.to_xugrid() + assert len(xu_ds.mesh2d_nFaces) == 1011