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 diff --git a/mikeio/spatial/grid_geometry.py b/mikeio/spatial/grid_geometry.py index 7021e1dc7..3c741e4d2 100644 --- a/mikeio/spatial/grid_geometry.py +++ b/mikeio/spatial/grid_geometry.py @@ -2,8 +2,7 @@ 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, GeometryPoint2D, @@ -11,7 +10,6 @@ GeometryUndefined, BoundingBox, ) -from ..eum import EUMType, EUMUnit def _check_equidistant(x: np.ndarray) -> None: @@ -652,58 +650,75 @@ def get_node_coordinates(self): gn = Grid2D(x=xn, y=yn) return gn.xy - def to_mesh(self, outfilename, projection=None, z=None): - """export grid to mesh file + def to_geometryFM(self, *, z=None, west=2, east=4, north=5, south=3): + """convert Grid2D to GeometryFM Parameters ---------- - outfilename : str - path of new mesh file - projection : str, optional - WKT projection string, by default 'LONG/LAT' - z : float or array(float), optional + z : float, optional bathymetry values for each node, by default 0 - if array: must have length=(nx+1)*(ny+1) + 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 """ - if projection is None: - projection = "LONG/LAT" + 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) + # node coordinates 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: + + 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 + 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, zn]) + 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, z=None): + """export grid to mesh file + + Parameters + ---------- + outfilename : str + path of new mesh file + z : float or array(float), optional + bathymetry values for each node, by default 0 + if array: must have length=(nx+1)*(ny+1) + """ + g = self.to_geometryFM() + + 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 == 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 - - builder = MeshBuilder() - builder.SetNodes(x, y, z, codes) - - elem_table = gn._to_element_table(index_base=1) - builder.SetElements(elem_table) - - builder.SetProjection(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 diff --git a/tests/test_geometry_grid.py b/tests/test_geometry_grid.py index 98cf42735..28b637091 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,42 @@ 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] + + 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_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"