Merge pull request #779 from DHI/test_dfsu_read_x_y

Add missing dfsu read x,y test

Author: ecomodeller

mikeio/dfsu/_dfsu.py

@@ -399,8 +399,8 @@ def read(
         time: int | str | slice | None = None,
         elements: Sequence[int] | np.ndarray | None = None,
         area: tuple[float, float, float, float] | None = None,
-        x: float | None = None,
-        y: float | None = None,
+        x: float | Sequence[float] | None = None,
+        y: float | Sequence[float] | None = None,
         keepdims: bool = False,
         dtype: Any = np.float32,
         error_bad_data: bool = True,
@@ -421,7 +421,7 @@ def read(
             Read only data inside (horizontal) area given as a
             bounding box (tuple with left, lower, right, upper)
             or as list of coordinates for a polygon, by default None
-        x, y: float, optional
+        x, y: float or list[float], optional
             Read only data for elements containing the (x,y) points(s),
             by default None
         elements: list[int], optional
@@ -447,8 +447,11 @@ def read(
         single_time_selected, time_steps = _valid_timesteps(dfs, time)
 
         _validate_elements_and_geometry_sel(elements, area=area, x=x, y=y)
-        if elements is None:
-            elements = self._parse_geometry_sel(area=area, x=x, y=y)
+        if area is not None:
+            elements = self.geometry._elements_in_area(area)
+
+        if (x is not None) or (y is not None):
+            elements = self.geometry.find_index(x=x, y=y)
 
         if elements is None:
             geometry = self.geometry
@@ -459,44 +462,34 @@ def read(
             geometry = self.geometry.elements_to_geometry(elements)
 
         item_numbers = _valid_item_numbers(dfs.ItemInfo, items)
-        items = _get_item_info(dfs.ItemInfo, item_numbers)
         n_items = len(item_numbers)
 
-        deletevalue = self.deletevalue
-
-        data_list = []
-
-        shape: tuple[int, ...]
-
         t_rel = np.zeros(len(time_steps))
 
         n_steps = len(time_steps)
-        shape = (
+        shape: tuple[int, ...] = (
             (n_elems,)
             if (single_time_selected and not keepdims)
             else (n_steps, n_elems)
         )
-        for item in range(n_items):
-            # Initialize an empty data block
-            data: np.ndarray = np.ndarray(shape=shape, dtype=dtype)
-            data_list.append(data)
+        data_list: list[np.ndarray] = [
+            np.ndarray(shape=shape, dtype=dtype) for _ in range(n_items)
+        ]
 
         for i in trange(n_steps, disable=not self.show_progress):
-            it = time_steps[i]
             for item in range(n_items):
-                dfs, d, t = _read_item_time_step(
+                dfs, d, t_rel[i] = _read_item_time_step(
                     dfs=dfs,
                     filename=self._filename,
                     time=time,
                     item_numbers=item_numbers,
-                    deletevalue=deletevalue,
+                    deletevalue=self.deletevalue,
                     shape=shape,
                     item=item,
-                    it=it,
+                    it=time_steps[i],
                     error_bad_data=error_bad_data,
                     fill_bad_data_value=fill_bad_data_value,
                 )
-                t_rel[i] = t
 
                 if elements is not None:
                     d = d[elements]
@@ -506,13 +499,9 @@ def read(
                 else:
                     data_list[item][i] = d
 
-        time = pd.to_datetime(t_rel, unit="s", origin=self.start_time)
-
         dfs.Close()
 
-        dims: tuple[str, ...]
-
-        dims = ("time", "element")
+        dims: tuple[str, ...] = ("time", "element")
 
         if single_time_selected and not keepdims:
             dims = ("element",)
@@ -522,10 +511,13 @@ def read(
             dims = tuple([d for d in dims if d != "element"])
             data_list = [np.squeeze(d, axis=-1) for d in data_list]
 
+        time = pd.to_datetime(t_rel, unit="s", origin=self.start_time)
+        item_infos = _get_item_info(dfs.ItemInfo, item_numbers)
+
         return Dataset(
             data_list,
             time,
-            items,
+            item_infos,
             geometry=geometry,
             dims=dims,
             validate=False,
@@ -558,53 +550,6 @@ def append(self, ds: Dataset, validate: bool = True) -> None:
         info = _get_dfsu_info(self._filename)
         self._n_timesteps = info.n_timesteps
 
-    def _parse_geometry_sel(
-        self,
-        area: tuple[float, float, float, float] | None,
-        x: float | None,
-        y: float | None,
-    ) -> np.ndarray | None:
-        """Parse geometry selection.
-
-        Parameters
-        ----------
-        area : list[float], optional
-            Read only data inside (horizontal) area given as a
-            bounding box (tuple with left, lower, right, upper)
-            or as list of coordinates for a polygon, by default None
-        x : float, optional
-            Read only data for elements containing the (x,y) points(s),
-            by default None
-        y : float, optional
-            Read only data for elements containing the (x,y) points(s),
-            by default None
-
-        Returns
-        -------
-        list[int]
-            List of element ids
-
-        Raises
-        ------
-        ValueError
-            If no elements are found in selection
-
-        """
-        elements = None
-
-        if area is not None:
-            elements = self.geometry._elements_in_area(area)
-
-        if (x is not None) or (y is not None):
-            elements = self.geometry.find_index(x=x, y=y)
-
-        if (x is not None) or (y is not None) or (area is not None):
-            # selection was attempted
-            if (elements is None) or len(elements) == 0:
-                raise ValueError("No elements in selection!")
-
-        return elements
-
     def get_overset_grid(
         self,
         dx: float | None = None,

mikeio/spatial/_FM_geometry.py

@@ -1077,14 +1077,16 @@ def _elements_in_area(
             xc = self.element_coordinates[:, 0]
             yc = self.element_coordinates[:, 1]
             mask = (xc >= x0) & (xc <= x1) & (yc >= y0) & (yc <= y1)
-            return np.where(mask)[0]
         elif self._area_is_polygon(area):
             polygon = np.array(area)
             xy = self.element_coordinates[:, :2]
             mask = self._inside_polygon(polygon, xy)
-            return np.where(mask)[0]
         else:
             raise ValueError("'area' must be bbox [x0,y0,x1,y1] or polygon")
+        elements = np.where(mask)[0]
+        if len(elements) == 0:
+            raise ValueError("No elements in selection!")
+        return elements
 
     def elements_to_geometry(
         self, elements: int | Sequence[int], keepdims: bool = False

tests/test_dfsu.py

@@ -259,6 +259,35 @@ def test_read_elements():
     assert ds2.Wind_speed.to_numpy()[0, 1] == pytest.approx(9.530759811401367)
 
 
+def test_read_x_y() -> None:
+    x = [1.49318531, 3.69276145]
+    y = [53.97088571, 54.08928194]
+    dfs = mikeio.open("tests/testdata/wind_north_sea.dfsu")
+    assert isinstance(dfs, mikeio.Dfsu2DH)
+    ds = dfs.read(x=x, y=y, time=0, keepdims=True)
+    assert isinstance(ds.geometry, mikeio.GeometryFM2D)
+    assert ds.geometry.element_coordinates[0][0] == pytest.approx(1.4931853081272184)
+    assert ds.Wind_speed.to_numpy()[0, 0] == pytest.approx(9.530759811401367)
+
+    x = x[::-1]
+    y = y[::-1]
+    ds2 = mikeio.read(
+        filename="tests/testdata/wind_north_sea.dfsu", x=x, y=y, time=0, keepdims=True
+    )
+    assert isinstance(ds2.geometry, mikeio.GeometryFM2D)
+    assert ds2.geometry.element_coordinates[1][0] == pytest.approx(1.4931853081272184)
+    assert ds2.Wind_speed.to_numpy()[0, 1] == pytest.approx(9.530759811401367)
+
+    x = 1.49318531
+    y = 53.97088571
+    ds3 = mikeio.read(
+        filename="tests/testdata/wind_north_sea.dfsu", x=x, y=y, time=0, keepdims=True
+    )
+    assert isinstance(ds3.geometry, mikeio.spatial.GeometryPoint2D)
+    assert ds3.geometry.x == pytest.approx(1.4931853081272184)
+    assert ds3.geometry.y == pytest.approx(53.97088571)
+
+
 def test_find_index_on_island():
     filename = "tests/testdata/FakeLake.dfsu"
     dfs = mikeio.open(filename)