From c4bc22b8d3d51182a6955bc7198256618384ec3f Mon Sep 17 00:00:00 2001 From: Jesper Sandvig Mariegaard <34088801+jsmariegaard@users.noreply.github.com> Date: Fri, 24 Jun 2022 10:19:36 +0000 Subject: [PATCH 1/9] fix boundary codes, round node_coordinates, default projection --- mikeio/spatial/grid_geometry.py | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/mikeio/spatial/grid_geometry.py b/mikeio/spatial/grid_geometry.py index 7021e1dc7..d7d2cec34 100644 --- a/mikeio/spatial/grid_geometry.py +++ b/mikeio/spatial/grid_geometry.py @@ -4,6 +4,7 @@ import numpy as np from mikecore.eum import eumQuantity from mikecore.MeshBuilder import MeshBuilder + from .geometry import ( _Geometry, GeometryPoint2D, @@ -666,7 +667,7 @@ def to_mesh(self, outfilename, projection=None, z=None): if array: must have length=(nx+1)*(ny+1) """ if projection is None: - projection = "LONG/LAT" + projection = self.projection # get node based grid xn = self._centers_to_nodes(self.x) @@ -687,14 +688,14 @@ def to_mesh(self, outfilename, projection=None, z=None): "z must either be scalar or have length of nodes ((nx+1)*(ny+1))" ) codes = np.zeros(n, dtype=int) - codes[y == gn.bbox.top] = 5 # north - codes[x == gn.bbox.right] = 4 # east - codes[y == gn.bbox.bottom] = 3 # south - codes[x == gn.bbox.left] = 2 # west - codes[(y == gn.bbox.top) & (x == gn.bbox.left)] = 5 # corner->north + codes[y == y[-1]] = 5 # north + codes[x == x[-1]] = 4 # east + codes[y == y[0]] = 3 # south + codes[x == x[0]] = 2 # west + codes[(y == y[-1]) & (x == x[0])] = 5 # corner->north builder = MeshBuilder() - builder.SetNodes(x, y, z, codes) + builder.SetNodes(np.round(x, 8), np.round(y, 8), z, codes) elem_table = gn._to_element_table(index_base=1) builder.SetElements(elem_table) From 9f452f884666078c59a4c48bbae91474f8626d3f Mon Sep 17 00:00:00 2001 From: Jesper Sandvig Mariegaard <34088801+jsmariegaard@users.noreply.github.com> Date: Fri, 24 Jun 2022 10:23:23 +0000 Subject: [PATCH 2/9] to_geometryFM --- mikeio/spatial/grid_geometry.py | 30 ++++++++++++++++++++++++++++++ 1 file changed, 30 insertions(+) diff --git a/mikeio/spatial/grid_geometry.py b/mikeio/spatial/grid_geometry.py index d7d2cec34..076dc3a63 100644 --- a/mikeio/spatial/grid_geometry.py +++ b/mikeio/spatial/grid_geometry.py @@ -653,6 +653,36 @@ def get_node_coordinates(self): gn = Grid2D(x=xn, y=yn) return gn.xy + def to_geometryFM(self): + """convert Grid2D to GeometryFM""" + from mikeio.spatial.FM_geometry import GeometryFM + + # get node based grid + xn = self._centers_to_nodes(self.x) + yn = self._centers_to_nodes(self.y) + gn = Grid2D(x=xn, y=yn) + + x = gn.xy[:, 0] + y = gn.xy[:, 1] + n = gn.nx * gn.ny + z = np.zeros_like(x) + + codes = np.zeros(n, dtype=int) + codes[y == y[-1]] = 5 # north + codes[x == x[-1]] = 4 # east + codes[y == y[0]] = 3 # south + codes[x == x[0]] = 2 # west + codes[(y == y[-1]) & (x == x[0])] = 5 # corner->north + + nc = np.column_stack([x, y, z]) + elem_table = gn._to_element_table(index_base=0) + return GeometryFM( + node_coordinates=nc, + element_table=elem_table, + codes=codes, + projection=self.projection, + ) + def to_mesh(self, outfilename, projection=None, z=None): """export grid to mesh file From ebf1ee5299f67d7e80122fef1590e1d60bf40f94 Mon Sep 17 00:00:00 2001 From: Jesper Sandvig Mariegaard <34088801+jsmariegaard@users.noreply.github.com> Date: Fri, 24 Jun 2022 14:50:59 +0000 Subject: [PATCH 3/9] remove duplicate import --- mikeio/dataarray.py | 1 - 1 file changed, 1 deletion(-) diff --git a/mikeio/dataarray.py b/mikeio/dataarray.py index 926388d58..7b475abd2 100644 --- a/mikeio/dataarray.py +++ b/mikeio/dataarray.py @@ -26,7 +26,6 @@ from mikecore.DfsuFile import DfsuFileType from .spatial.FM_utils import _plot_map, _plot_vertical_profile from mikecore.DfsuFile import DfsuFileType -from .spatial.FM_utils import _plot_map from .spectral_utils import plot_2dspectrum, calc_m0_from_spectrum from .data_utils import DataUtilsMixin From 4a8773cdafd4b4f0cf903ac5f20502c6163adad7 Mon Sep 17 00:00:00 2001 From: Jesper Sandvig Mariegaard <34088801+jsmariegaard@users.noreply.github.com> Date: Thu, 28 Jul 2022 08:06:31 +0200 Subject: [PATCH 4/9] remove projection argument --- mikeio/spatial/grid_geometry.py | 9 ++------- 1 file changed, 2 insertions(+), 7 deletions(-) diff --git a/mikeio/spatial/grid_geometry.py b/mikeio/spatial/grid_geometry.py index 076dc3a63..d5e5a574d 100644 --- a/mikeio/spatial/grid_geometry.py +++ b/mikeio/spatial/grid_geometry.py @@ -683,22 +683,17 @@ def to_geometryFM(self): projection=self.projection, ) - def to_mesh(self, outfilename, projection=None, z=None): + def to_mesh(self, outfilename, z=None): """export grid to mesh file Parameters ---------- outfilename : str path of new mesh file - projection : str, optional - WKT projection string, by default 'LONG/LAT' z : float or array(float), optional bathymetry values for each node, by default 0 if array: must have length=(nx+1)*(ny+1) """ - if projection is None: - projection = self.projection - # get node based grid xn = self._centers_to_nodes(self.x) yn = self._centers_to_nodes(self.y) @@ -730,7 +725,7 @@ def to_mesh(self, outfilename, projection=None, z=None): elem_table = gn._to_element_table(index_base=1) builder.SetElements(elem_table) - builder.SetProjection(projection) + builder.SetProjection(self.projection) quantity = eumQuantity.Create(EUMType.Bathymetry, EUMUnit.meter) builder.SetEumQuantity(quantity) newMesh = builder.CreateMesh() From 1dcba5a6bad315072e911d5b2e711e5dd06e8b66 Mon Sep 17 00:00:00 2001 From: Jesper Sandvig Mariegaard <34088801+jsmariegaard@users.noreply.github.com> Date: Thu, 28 Jul 2022 08:25:44 +0200 Subject: [PATCH 5/9] refactor to_mesh to use to_geometryFM --- mikeio/spatial/grid_geometry.py | 42 ++++++--------------------------- 1 file changed, 7 insertions(+), 35 deletions(-) diff --git a/mikeio/spatial/grid_geometry.py b/mikeio/spatial/grid_geometry.py index d5e5a574d..72beb78be 100644 --- a/mikeio/spatial/grid_geometry.py +++ b/mikeio/spatial/grid_geometry.py @@ -2,8 +2,6 @@ from typing import Sequence, Tuple, Union import warnings import numpy as np -from mikecore.eum import eumQuantity -from mikecore.MeshBuilder import MeshBuilder from .geometry import ( _Geometry, @@ -12,7 +10,6 @@ GeometryUndefined, BoundingBox, ) -from ..eum import EUMType, EUMUnit def _check_equidistant(x: np.ndarray) -> None: @@ -662,6 +659,7 @@ def to_geometryFM(self): yn = self._centers_to_nodes(self.y) gn = Grid2D(x=xn, y=yn) + # node coordinates x = gn.xy[:, 0] y = gn.xy[:, 1] n = gn.nx * gn.ny @@ -694,42 +692,16 @@ def to_mesh(self, outfilename, z=None): bathymetry values for each node, by default 0 if array: must have length=(nx+1)*(ny+1) """ - # get node based grid - xn = self._centers_to_nodes(self.x) - yn = self._centers_to_nodes(self.y) - gn = Grid2D(x=xn, y=yn) + g = self.to_geometryFM() - x = gn.xy[:, 0] - y = gn.xy[:, 1] - n = gn.nx * gn.ny - if z is None: - z = np.zeros(n) - else: - if np.isscalar(z): - z = z * np.ones(n) - else: - if len(z) != n: + if z is not None: + if not np.isscalar(z): + if len(z) != g.n_nodes: raise ValueError( "z must either be scalar or have length of nodes ((nx+1)*(ny+1))" ) - codes = np.zeros(n, dtype=int) - codes[y == y[-1]] = 5 # north - codes[x == x[-1]] = 4 # east - codes[y == y[0]] = 3 # south - codes[x == x[0]] = 2 # west - codes[(y == y[-1]) & (x == x[0])] = 5 # corner->north - - builder = MeshBuilder() - builder.SetNodes(np.round(x, 8), np.round(y, 8), z, codes) - - elem_table = gn._to_element_table(index_base=1) - builder.SetElements(elem_table) - - builder.SetProjection(self.projection) - quantity = eumQuantity.Create(EUMType.Bathymetry, EUMUnit.meter) - builder.SetEumQuantity(quantity) - newMesh = builder.CreateMesh() - newMesh.Write(outfilename) + g.node_coordinates[:, 2] = z + g.to_mesh(outfilename=outfilename) @dataclass From 80aa8200f4ed9e20cd184ac76f307ddc36304e77 Mon Sep 17 00:00:00 2001 From: Jesper Sandvig Mariegaard <34088801+jsmariegaard@users.noreply.github.com> Date: Thu, 28 Jul 2022 08:40:03 +0200 Subject: [PATCH 6/9] test_to_geometryFM --- tests/test_geometry_grid.py | 21 +++++++++++++++++++++ 1 file changed, 21 insertions(+) diff --git a/tests/test_geometry_grid.py b/tests/test_geometry_grid.py index 98cf42735..329184d73 100644 --- a/tests/test_geometry_grid.py +++ b/tests/test_geometry_grid.py @@ -5,6 +5,7 @@ import pytest from mikeio import Mesh from mikeio import Grid2D, Grid1D +from mikeio.spatial.FM_geometry import GeometryFM def test_create_nx_ny(): @@ -243,6 +244,26 @@ def test_find_index(): assert jj[2] == j2 +def test_to_geometryFM(): + nx = 5 + ny = 3 + grd = Grid2D(nx=nx, dx=1, ny=ny, dy=2) + g = grd.to_geometryFM() + assert isinstance(g, GeometryFM) + assert g.n_elements == nx * ny + assert g.n_nodes == (nx + 1) * (ny + 1) + assert g.projection_string == "NON-UTM" + + xe = g.element_coordinates[:, 0] + ye = g.element_coordinates[:, 1] + assert xe[0] == grd.x[0] + assert xe[ny - 1] == grd.x[0] + assert ye[0] == grd.y[0] + assert xe[-1] == grd.x[-1] + assert ye[ny - 1] == grd.y[-1] + assert ye[-1] == grd.y[-1] + + def test_to_mesh(tmp_path: Path): outfilename = tmp_path / "temp.mesh" From 2ef47baef9bb6b5173c6f0718ff1b19d1a4dea54 Mon Sep 17 00:00:00 2001 From: Jesper Sandvig Mariegaard <34088801+jsmariegaard@users.noreply.github.com> Date: Thu, 28 Jul 2022 08:47:21 +0200 Subject: [PATCH 7/9] test also codes --- tests/test_geometry_grid.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/tests/test_geometry_grid.py b/tests/test_geometry_grid.py index 329184d73..de317ebdd 100644 --- a/tests/test_geometry_grid.py +++ b/tests/test_geometry_grid.py @@ -263,6 +263,11 @@ def test_to_geometryFM(): assert ye[ny - 1] == grd.y[-1] assert ye[-1] == grd.y[-1] + assert g.codes[0] == 2 # west (lower left corner) + assert g.codes[1] == 3 # south + assert g.codes[-2] == 5 # north + assert g.codes[-1] == 4 # east (upper right corner) + def test_to_mesh(tmp_path: Path): outfilename = tmp_path / "temp.mesh" From 9e08774f91fe8b29b0782bdea522d176aea2fde8 Mon Sep 17 00:00:00 2001 From: Henrik Andersson Date: Wed, 3 Aug 2022 16:46:12 +0200 Subject: [PATCH 8/9] Option to specify bound codes and z value --- mikeio/spatial/grid_geometry.py | 33 +++++++++++++++++++++++++-------- tests/test_geometry_grid.py | 10 ++++++++++ 2 files changed, 35 insertions(+), 8 deletions(-) diff --git a/mikeio/spatial/grid_geometry.py b/mikeio/spatial/grid_geometry.py index 72beb78be..f3ece79da 100644 --- a/mikeio/spatial/grid_geometry.py +++ b/mikeio/spatial/grid_geometry.py @@ -650,8 +650,22 @@ def get_node_coordinates(self): gn = Grid2D(x=xn, y=yn) return gn.xy - def to_geometryFM(self): - """convert Grid2D to GeometryFM""" + def to_geometryFM(self,*,z=None, west=2, east=4, north=5,south=3): + """convert Grid2D to GeometryFM + + Parameters + ---------- + z : float, optional + bathymetry values for each node, by default 0 + west: int, optional + code value for west boundary + east: int, optional + code value for east boundary + north: int, optional + code value for north boundary + south: int, optional + code value for south boundary + """ from mikeio.spatial.FM_geometry import GeometryFM # get node based grid @@ -663,16 +677,19 @@ def to_geometryFM(self): x = gn.xy[:, 0] y = gn.xy[:, 1] n = gn.nx * gn.ny - z = np.zeros_like(x) + zn = np.zeros_like(x) + if z is not None: + zn[:] = z + codes = np.zeros(n, dtype=int) - codes[y == y[-1]] = 5 # north - codes[x == x[-1]] = 4 # east - codes[y == y[0]] = 3 # south - codes[x == x[0]] = 2 # west + codes[y == y[-1]] = north + codes[x == x[-1]] = east + codes[y == y[0]] = south + codes[x == x[0]] = west codes[(y == y[-1]) & (x == x[0])] = 5 # corner->north - nc = np.column_stack([x, y, z]) + nc = np.column_stack([x, y, zn]) elem_table = gn._to_element_table(index_base=0) return GeometryFM( node_coordinates=nc, diff --git a/tests/test_geometry_grid.py b/tests/test_geometry_grid.py index de317ebdd..d4c703c8f 100644 --- a/tests/test_geometry_grid.py +++ b/tests/test_geometry_grid.py @@ -268,6 +268,16 @@ def test_to_geometryFM(): assert g.codes[-2] == 5 # north assert g.codes[-1] == 4 # east (upper right corner) +def test_to_geometryFM_custom_z_custom_code(): + nx = 5 + ny = 3 + grd = Grid2D(nx=nx, dx=1, ny=ny, dy=2) + g = grd.to_geometryFM(z=-12.0, west=30) + assert isinstance(g, GeometryFM) + assert all(g.node_coordinates[:,2] == -12.0) + + assert g.codes[0] == 30 + def test_to_mesh(tmp_path: Path): outfilename = tmp_path / "temp.mesh" From 9f0a92231d61607ac8fcd2b6ed14b29293d70516 Mon Sep 17 00:00:00 2001 From: github-actions <${GITHUB_ACTOR}@users.noreply.github.com> Date: Wed, 3 Aug 2022 14:46:38 +0000 Subject: [PATCH 9/9] Formatted Python code with black --- mikeio/spatial/grid_geometry.py | 4 ++-- tests/test_geometry_grid.py | 3 ++- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/mikeio/spatial/grid_geometry.py b/mikeio/spatial/grid_geometry.py index f3ece79da..3c741e4d2 100644 --- a/mikeio/spatial/grid_geometry.py +++ b/mikeio/spatial/grid_geometry.py @@ -650,7 +650,7 @@ def get_node_coordinates(self): gn = Grid2D(x=xn, y=yn) return gn.xy - def to_geometryFM(self,*,z=None, west=2, east=4, north=5,south=3): + def to_geometryFM(self, *, z=None, west=2, east=4, north=5, south=3): """convert Grid2D to GeometryFM Parameters @@ -681,7 +681,7 @@ def to_geometryFM(self,*,z=None, west=2, east=4, north=5,south=3): zn = np.zeros_like(x) if z is not None: zn[:] = z - + codes = np.zeros(n, dtype=int) codes[y == y[-1]] = north codes[x == x[-1]] = east diff --git a/tests/test_geometry_grid.py b/tests/test_geometry_grid.py index d4c703c8f..28b637091 100644 --- a/tests/test_geometry_grid.py +++ b/tests/test_geometry_grid.py @@ -268,13 +268,14 @@ def test_to_geometryFM(): assert g.codes[-2] == 5 # north assert g.codes[-1] == 4 # east (upper right corner) + def test_to_geometryFM_custom_z_custom_code(): nx = 5 ny = 3 grd = Grid2D(nx=nx, dx=1, ny=ny, dy=2) g = grd.to_geometryFM(z=-12.0, west=30) assert isinstance(g, GeometryFM) - assert all(g.node_coordinates[:,2] == -12.0) + assert all(g.node_coordinates[:, 2] == -12.0) assert g.codes[0] == 30