Merge pull request #2050 from OceanParcels/patch-no-time-interval-index-search

Patch no time interval index search

Author: fluidnumericsJoe

parcels/_datasets/unstructured/generic.py

@@ -44,12 +44,12 @@ def _stommel_gyre_delaunay():
     uxgrid.attrs["Conventions"] = "UGRID-1.0"
 
     # Define arrays U (zonal), V (meridional) and P (sea surface height)
-    U = np.zeros((1, nz1, lat.size), dtype=np.float64)
-    V = np.zeros((1, nz1, lat.size), dtype=np.float64)
+    U = np.zeros((1, nz1, uxgrid.n_face), dtype=np.float64)
+    V = np.zeros((1, nz1, uxgrid.n_face), dtype=np.float64)
     W = np.zeros((1, nz, lat.size), dtype=np.float64)
-    P = np.zeros((1, nz1, lat.size), dtype=np.float64)
+    P = np.zeros((1, nz1, uxgrid.n_face), dtype=np.float64)
 
-    for i, (x, y) in enumerate(zip(lon_flat, lat_flat, strict=False)):
+    for i, (x, y) in enumerate(zip(uxgrid.face_lon, uxgrid.face_lat, strict=False)):
         xi = x / 60.0
         yi = y / 60.0
 
@@ -61,10 +61,10 @@ def _stommel_gyre_delaunay():
         data=U,
         name="U",
         uxgrid=uxgrid,
-        dims=["time", "nz1", "n_node"],
+        dims=["time", "nz1", "n_face"],
         coords=dict(
             time=(["time"], [TIME[0]]),
-            nz1=(["nz1"], [0]),
+            nz1=(["nz1"], zc),
         ),
         attrs=dict(
             description="zonal velocity", units="m/s", location="node", mesh="delaunay", Conventions="UGRID-1.0"
@@ -74,7 +74,7 @@ def _stommel_gyre_delaunay():
         data=V,
         name="V",
         uxgrid=uxgrid,
-        dims=["time", "nz1", "n_node"],
+        dims=["time", "nz1", "n_face"],
         coords=dict(
             time=(["time"], [TIME[0]]),
             nz1=(["nz1"], zc),
@@ -100,7 +100,7 @@ def _stommel_gyre_delaunay():
         data=P,
         name="p",
         uxgrid=uxgrid,
-        dims=["time", "nz1", "n_node"],
+        dims=["time", "nz1", "n_face"],
         coords=dict(
             time=(["time"], [TIME[0]]),
             nz1=(["nz1"], zc),

parcels/_index_search.py

@@ -38,7 +38,7 @@ def _search_time_index(field: Field, time: datetime):
     if the sampled value is outside the time value range.
     """
     if field.time_interval is None:
-        return 0
+        return 0, 0
 
     if time not in field.time_interval:
         _raise_time_extrapolation_error(time, field=None)

tests/v4/test_field.py

@@ -147,6 +147,20 @@ def test_field_unstructured_z_linear():
     assert np.isclose(W.eval(time=ds.time[0].values, z=900.0, y=30.0, x=30.0, applyConversion=False), 900.0)
 
 
+def test_field_constant_in_time():
+    """Tests field evaluation for a field with no time interval (i.e., constant in time)."""
+    ds = datasets_unstructured["stommel_gyre_delaunay"]
+    grid = UxGrid(ds.uxgrid, z=ds.coords["nz"])
+    # Note that the vertical coordinate is required to be the position of the layer interfaces ("nz"), not the mid-layers ("nz1")
+    P = Field(name="p", data=ds.p, grid=grid, interp_method=UXPiecewiseConstantFace)
+
+    # Assert that the field can be evaluated at any time, and returns the same value
+    time = np.datetime64("2000-01-01T00:00:00")
+    P1 = P.eval(time=time, z=10.0, y=30.0, x=30.0, applyConversion=False)
+    P2 = P.eval(time=time + np.timedelta64(1, "D"), z=10.0, y=30.0, x=30.0, applyConversion=False)
+    assert np.isclose(P1, P2)
+
+
 def test_field_unstructured_grid_creation(): ...