Merge pull request #668 from DHI/non-utm

Dfs{2,3} local coordinates

mikeio/dfs/_dfs2.py

@@ -113,12 +113,13 @@ class Dfs2(_Dfs123):
     def __init__(
         self,
         filename: str | Path,
-        type: Literal["horizontal", "spectral"] = "horizontal",
+        type: Literal["horizontal", "spectral", "vertical"] = "horizontal",
     ):
         filename = str(filename)
         super().__init__(filename)
 
         is_spectral = type == "spectral"
+        is_vertical = type == "vertical"
         dfs = DfsFileFactory.Dfs2FileOpen(str(filename))
 
         x0 = dfs.SpatialAxis.X0 if is_spectral else 0.0
@@ -137,6 +138,7 @@ def __init__(
             orientation=orientation,
             origin=origin,
             is_spectral=is_spectral,
+            is_vertical=is_vertical,
         )
         dfs.Close()
         self._validate_no_orientation_in_geo()

mikeio/spatial/_grid_geometry.py

@@ -426,11 +426,12 @@ def __init__(
         dy: float | None = None,
         ny: int | None = None,
         bbox: Tuple[float, float, float, float] | None = None,
-        projection: str = "NON-UTM",
+        projection: str = "LONG/LAT",
         origin: Tuple[float, float] | None = None,
         orientation: float = 0.0,
         axis_names: Tuple[str, str] = ("x", "y"),
         is_spectral: bool = False,
+        is_vertical: bool = False,
     ):
         """Create equidistant 2D spatial geometry
 
@@ -464,6 +465,8 @@ def __init__(
             names of x and y axes, by default ("x", "y")
         is_spectral : bool, optional
             if True, the grid is spectral, by default False
+        is_vertical : bool, optional
+            if True, the grid is vertical, by default False
 
         Examples
         --------
@@ -491,7 +494,13 @@ def __init__(
             dy = self._dx if dy is None else dy
             self._y0, self._dy, self._ny = _parse_grid_axis("y", y, y0, dy, ny)
 
+        if self.is_local_coordinates and not (is_spectral or is_vertical):
+            self._x0 = self._x0 + self._dx / 2
+            self._y0 = self._y0 + self._dy / 2
+            self._shift_origin_on_write = False
+
         self.is_spectral = is_spectral
+        self.is_vertical = is_vertical
 
         self.plot = _Grid2DPlotter(self)
 
@@ -945,6 +954,7 @@ def _index_to_Grid2D(
                 projection=self.projection,
                 orientation=self._orientation,
                 is_spectral=self.is_spectral,
+                is_vertical=self.is_vertical,
                 origin=origin,
             )
 
@@ -1120,6 +1130,10 @@ def __init__(
         self._origin = origin
         self._orientation = orientation
 
+        if self.is_local_coordinates:
+            self._x0 = self._x0 + self._dx / 2
+            self._y0 = self._y0 + self._dy / 2
+
     @property
     def default_dims(self) -> Tuple[str, ...]:
         return ("z", "y", "x")

tests/test_dfs2.py

@@ -40,7 +40,7 @@ def dfs2_pt_spectrum_linearf():
 @pytest.fixture
 def dfs2_vertical_nonutm():
     filepath = Path("tests/testdata/hd_vertical_slice.dfs2")
-    return mikeio.open(filepath)
+    return mikeio.open(filepath, type="vertical")
 
 
 @pytest.fixture
@@ -50,7 +50,7 @@ def dfs2_gebco():
 
 
 def test_get_time_without_reading_data():
-    dfs = mikeio.open("tests/testdata/hd_vertical_slice.dfs2")
+    dfs = mikeio.open("tests/testdata/hd_vertical_slice.dfs2", type="vertical")
 
     assert isinstance(dfs.time, pd.DatetimeIndex)
     assert len(dfs.time) == 13
@@ -281,8 +281,10 @@ def test_properties_vertical_nonutm(dfs2_vertical_nonutm):
 
 def test_isel_vertical_nonutm(dfs2_vertical_nonutm):
     ds = dfs2_vertical_nonutm.read()
+    assert ds.geometry.is_vertical
     dssel = ds.isel(y=slice(45, None))
     g = dssel.geometry
+    assert g.is_vertical
     assert g._x0 == 0
     assert g._y0 == 0  # TODO: should this be 45?
     assert g.x[0] == 0
@@ -346,7 +348,9 @@ def test_properties_pt_spectrum_linearf(dfs2_pt_spectrum_linearf):
 
 
 def test_dir_wave_spectra_relative_time_axis():
-    ds = mikeio.read("tests/testdata/dir_wave_analysis_spectra.dfs2")
+    ds = mikeio.open(
+        "tests/testdata/dir_wave_analysis_spectra.dfs2", type="spectral"
+    ).read()
     assert ds.n_items == 1
     assert ds.geometry.nx == 128
     assert ds.geometry.ny == 37
@@ -819,3 +823,44 @@ def test_read_dfs2_static_dt_zero():
 
     assert ds2.shape == (2, 2)
     assert "time" not in ds2.dims
+
+
+def test_write_read_local_coordinates(tmp_path):
+    da = mikeio.DataArray(
+        np.array([[1, 2, 3], [4, 5, 6]]),
+        geometry=mikeio.Grid2D(nx=3, ny=2, dx=0.5, projection="NON-UTM"),
+    )
+    fp = tmp_path / "local_coordinates.dfs2"
+    da.to_dfs(fp)
+
+    ds = mikeio.read(fp)
+    assert da.geometry == ds.geometry
+
+
+def test_to_xarray():
+    # data is not important here
+    data = np.array([[1, 2, 3], [4, 5, 6]])
+
+    # geographic coordinates
+    dag = mikeio.DataArray(
+        data=data,
+        geometry=mikeio.Grid2D(nx=3, ny=2, dx=0.5, projection="LONG/LAT"),
+    )
+    assert dag.geometry.x[0] == pytest.approx(0.0)
+    assert dag.geometry.y[0] == pytest.approx(0.0)
+    xr_dag = dag.to_xarray()
+    assert xr_dag.x[0] == pytest.approx(0.0)
+    assert xr_dag.y[0] == pytest.approx(0.0)
+
+    # local coordinates
+    da = mikeio.DataArray(
+        data=data,
+        geometry=mikeio.Grid2D(nx=3, ny=2, dx=0.5, projection="NON-UTM"),
+    )
+    # local coordinates (=NON-UTM) have a different convention, geometry.x still refers to element centers
+    assert da.geometry.x[0] == pytest.approx(0.25)
+    assert da.geometry.y[0] == pytest.approx(0.25)
+
+    xr_da = da.to_xarray()
+    assert xr_da.x[0] == pytest.approx(0.25)
+    assert xr_da.y[0] == pytest.approx(0.25)

tests/test_dfs3.py

@@ -201,3 +201,50 @@ def test_MIKE_SHE_dfs3_output():
     assert g2.x[0] == g.x[0] + 30 * g.dx
     assert g2.y[0] == g.y[0]  # + 35 * g.dy
     assert g2.origin == pytest.approx((g2.x[0], g2.y[0]))
+
+
+def test_write_read_local_coordinates(tmp_path):
+    da = mikeio.DataArray(
+        np.array([[[1, 2, 3], [4, 5, 6]], [[7, 8, 9], [10, 11, 12]]]),
+        geometry=mikeio.Grid3D(
+            nx=3, ny=2, nz=2, dx=0.5, dy=1, dz=1, projection="NON-UTM"
+        ),
+    )
+    fp = tmp_path / "local_coordinates.dfs3"
+    da.to_dfs(fp)
+
+    ds = mikeio.read(fp)
+    assert da.geometry == ds.geometry
+
+
+def test_to_xarray():
+    # data is not important here
+    data = np.array([[[1, 2, 3], [4, 5, 6]], [[7, 8, 9], [10, 11, 12]]])
+
+    # geographic coordinates
+    dag = mikeio.DataArray(
+        data=data,
+        geometry=mikeio.Grid3D(
+            nx=3, ny=2, nz=2, dy=0.5, dz=1, dx=0.5, projection="LONG/LAT"
+        ),
+    )
+    assert dag.geometry.x[0] == pytest.approx(0.0)
+    assert dag.geometry.y[0] == pytest.approx(0.0)
+    xr_dag = dag.to_xarray()
+    assert xr_dag.x[0] == pytest.approx(0.0)
+    assert xr_dag.y[0] == pytest.approx(0.0)
+
+    # local coordinates
+    da = mikeio.DataArray(
+        data=data,
+        geometry=mikeio.Grid3D(
+            nx=3, ny=2, nz=2, dy=0.5, dz=1, dx=0.5, projection="NON-UTM"
+        ),
+    )
+    # local coordinates (=NON-UTM) have a different convention, geometry.x still refers to element centers
+    assert da.geometry.x[0] == pytest.approx(0.25)
+    assert da.geometry.y[0] == pytest.approx(0.25)
+
+    xr_da = da.to_xarray()
+    assert xr_da.x[0] == pytest.approx(0.25)
+    assert xr_da.y[0] == pytest.approx(0.25)

tests/test_geometry_grid.py

@@ -306,7 +306,7 @@ def test_to_geometryFM():
     assert isinstance(g, GeometryFM2D)
     assert g.n_elements == nx * ny
     assert g.n_nodes == (nx + 1) * (ny + 1)
-    assert g.projection_string == "NON-UTM"
+    assert g.projection_string == "LONG/LAT"
 
     xe = g.element_coordinates[:, 0]
     ye = g.element_coordinates[:, 1]
@@ -398,7 +398,7 @@ def test_grid2d_equality():
 
     assert g1 == g2
 
-    g3 = Grid2D(dx=0.1, nx=2, dy=0.2, ny=4, projection="LONG/LAT")
+    g3 = Grid2D(dx=0.1, nx=2, dy=0.2, ny=4, projection="NON-UTM")
     g4 = Grid2D(dx=0.1, nx=2, dy=0.2, ny=4)
 
     assert g3 != g4