Merge pull request #374 from DHI/grid2d-to-geometryfm

Grid2D to GeometryFM

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
 

mikeio/spatial/grid_geometry.py

@@ -2,16 +2,14 @@
 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,
     GeometryPoint3D,
     GeometryUndefined,
     BoundingBox,
 )
-from ..eum import EUMType, EUMUnit
 
 
 def _check_equidistant(x: np.ndarray) -> None:
@@ -749,58 +747,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

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"