Repository: DHI/mikeio
Pull request #483 (2022-11-24 06:01:48)

Xarray point add missing coordinates for GeometryPoint{2,3}D


Patch:
diff --git a/mikeio/dataarray.py b/mikeio/dataarray.py
index b514d7d4..3aa73b45 100644
--- a/mikeio/dataarray.py
+++ b/mikeio/dataarray.py
@@ -1488,7 +1488,7 @@ class DataArray(DataUtilsMixin, TimeSeries):
             if isinstance(self.geometry, Grid2D):  # TODO DIY bilinear interpolation
                 xr_da = self.to_xarray()
                 dai = xr_da.interp(x=x, y=y).values
-                geometry = GeometryPoint2D(x=x, y=y)
+                geometry = GeometryPoint2D(x=x, y=y, projection=self.geometry.projection)
             elif isinstance(self.geometry, Grid1D):
                 if interpolant is None:
                     interpolant = self.geometry.get_spatial_interpolant(coords)
@@ -1501,9 +1501,9 @@ class DataArray(DataUtilsMixin, TimeSeries):
                     )
                 dai = self.geometry.interp2d(self, *interpolant).flatten()
                 if z is None:
-                    geometry = GeometryPoint2D(x=x, y=y)
+                    geometry = GeometryPoint2D(x=x, y=y, projection=self.geometry.projection)
                 else:
-                    geometry = GeometryPoint3D(x=x, y=y, z=z)
+                    geometry = GeometryPoint3D(x=x, y=y, z=z, projection=self.geometry.projection)
 
             da = DataArray(
                 data=dai, time=self.time, geometry=geometry, item=deepcopy(self.item)
@@ -2245,6 +2245,13 @@ class DataArray(DataUtilsMixin, TimeSeries):
             coords["element"] = xr.DataArray(
                 data=self.geometry.element_ids, dims="element"
             )
+        elif isinstance(self.geometry, GeometryPoint2D):
+            coords["x"] = self.geometry.x
+            coords["y"] = self.geometry.y
+        elif isinstance(self.geometry, GeometryPoint3D):
+            coords["x"] = self.geometry.x
+            coords["y"] = self.geometry.y
+            coords["z"] = self.geometry.z
 
         xr_da = xr.DataArray(
             data=self.values,
diff --git a/mikeio/spatial/FM_geometry.py b/mikeio/spatial/FM_geometry.py
index f884b307..ca3c24af 100644
--- a/mikeio/spatial/FM_geometry.py
+++ b/mikeio/spatial/FM_geometry.py
@@ -1004,9 +1004,9 @@ class GeometryFM(_Geometry):
             coords = self.element_coordinates[idx].flatten()
 
             if self.is_layered:
-                return GeometryPoint3D(*coords)
+                return GeometryPoint3D(*coords, projection=self.projection)
             else:
-                return GeometryPoint2D(coords[0], coords[1])
+                return GeometryPoint2D(coords[0], coords[1], projection=self.projection)
         else:
             if self._type == DfsuFileType.DfsuSpectral1D:
                 return self._nodes_to_geometry(nodes=idx)
diff --git a/mikeio/spatial/geometry.py b/mikeio/spatial/geometry.py
index af124983..237300eb 100644
--- a/mikeio/spatial/geometry.py
+++ b/mikeio/spatial/geometry.py
@@ -95,7 +95,8 @@ class GeometryUndefined(_Geometry):
 
 
 class GeometryPoint2D(_Geometry):
-    def __init__(self, x: float, y: float):
+    def __init__(self, x: float, y: float, projection = None):
+        self._projstr = projection
         self.x = x
         self.y = y
 
@@ -108,7 +109,8 @@ class GeometryPoint2D(_Geometry):
 
 
 class GeometryPoint3D(_Geometry):
-    def __init__(self, x: float, y: float, z: float):
+    def __init__(self, x: float, y: float, z: float, projection = None):
+        self._projstr = projection
         self.x = x
         self.y = y
         self.z = z


Test patch:
diff --git a/tests/test_dfsu.py b/tests/test_dfsu.py
index 484deb89..609a1461 100644
--- a/tests/test_dfsu.py
+++ b/tests/test_dfsu.py
@@ -924,12 +924,29 @@ def test_dataset_interp():
     y = 6184000
 
     dai = da.interp(x=x, y=y)
-
     assert isinstance(dai, DataArray)
     assert dai.shape == (ds.n_timesteps,)
     assert dai.name == da.name
     assert dai.geometry.x == x
     assert dai.geometry.y == y
+    assert dai.geometry.projection == ds.geometry.projection
+
+def test_dataset_interp_to_xarray():
+    ds = mikeio.read("tests/testdata/oresundHD_run1.dfsu")
+
+    assert not ds.geometry.is_geo
+
+    x = 360000
+    y = 6184000
+
+    dsi = ds.interp(x=x, y=y)
+
+    xr_dsi = dsi.to_xarray()
+    assert float(xr_dsi.x) == pytest.approx(x)
+    assert float(xr_dsi.y) == pytest.approx(y)
+
+
+
 
 
 def test_interp_like_grid():
diff --git a/tests/test_dfsu_layered.py b/tests/test_dfsu_layered.py
index 11928c44..dc550e7e 100644
--- a/tests/test_dfsu_layered.py
+++ b/tests/test_dfsu_layered.py
@@ -225,6 +225,7 @@ def test_read_dfsu3d_xyz():
     ds = dfs.read()  # all data in file
     dspt1 = ds.sel(x=x, y=y, z=z)
     assert isinstance(dspt1.geometry, GeometryPoint3D)
+    assert dspt1.geometry.projection == ds.geometry.projection
 
     dspt2 = dfs.read(x=x, y=y, z=z)
     assert isinstance(dspt2.geometry, GeometryPoint3D)
@@ -237,6 +238,20 @@ def test_read_dfsu3d_xyz():
     assert dspt3[0].values == dspt1[0].values[-1]
     # 20.531237
 
+def test_read_dfsu3d_xyz_to_xarray():
+    filename = "tests/testdata/oresund_sigma_z.dfsu"
+    dfs = mikeio.open(filename)
+
+    (x, y, z) = (333934.1, 6158101.5, -5)
+
+    ds = dfs.read()  # all data in file
+    dspt1 = ds.sel(x=x, y=y, z=z)
+    
+    xr_ds = dspt1.to_xarray()
+    assert float(xr_ds.x) == pytest.approx(x)
+    assert float(xr_ds.y) == pytest.approx(y)
+    assert float(xr_ds.z) == pytest.approx(z)
+
 
 def test_read_column_select_single_time_plot():
     filename = "tests/testdata/oresund_sigma_z.dfsu"