[#2010] Add Conventions to sample dataset; add associated test

Author: fluidnumericsJoe

parcels/_datasets/unstructured/generic.py

@@ -1,4 +1,5 @@
 import math
+from datetime import datetime, timedelta
 
 import numpy as np
 import pandas as pd
@@ -35,6 +36,7 @@ def _stommel_gyre_delaunay():
         method="regional_delaunay",
         boundary_points=boundary_points,
     )
+    uxgrid.attrs["Conventions"] = "UGRID-1.0"
 
     # Define arrays U (zonal), V (meridional) and P (sea surface height)
     U = np.zeros((1, 1, lat.size), dtype=np.float64)
@@ -106,6 +108,8 @@ def _fesom2_square_delaunay_uniform_z_coordinate():
     zc = 0.5 * (zf[:-1] + zf[1:])  # Vertical element centers
     nz = zf.size
     nz1 = zc.size
+    num_days = 5
+    date_array = [datetime(2000, 1, 1) + timedelta(days=i) for i in range(num_days)]
 
     # mask any point on one of the boundaries
     mask = (
@@ -119,26 +123,29 @@ def _fesom2_square_delaunay_uniform_z_coordinate():
         method="regional_delaunay",
         boundary_points=boundary_points,
     )
+    uxgrid.attrs["Conventions"] = "UGRID-1.0"
 
     # Define arrays U (zonal), V (meridional) and P (sea surface height)
     U = np.ones(
-        (1, nz1, uxgrid.n_face), dtype=np.float64
+        (num_days, nz1, uxgrid.n_face), dtype=np.float64
     )  # Lateral velocity is on the element centers and face centers
     V = np.ones(
-        (1, nz1, uxgrid.n_face), dtype=np.float64
+        (num_days, nz1, uxgrid.n_face), dtype=np.float64
     )  # Lateral velocity is on the element centers and face centers
     W = np.zeros(
-        (1, nz, uxgrid.n_node), dtype=np.float64
+        (num_days, nz, uxgrid.n_node), dtype=np.float64
     )  # Vertical velocity is on the element faces and face vertices
-    P = np.ones((1, nz1, uxgrid.n_node), dtype=np.float64)  # Pressure is on the element centers and face vertices
+    P = np.ones(
+        (num_days, nz1, uxgrid.n_node), dtype=np.float64
+    )  # Pressure is on the element centers and face vertices
 
     u = ux.UxDataArray(
         data=U,
         name="U",
         uxgrid=uxgrid,
         dims=["time", "nz1", "n_face"],
         coords=dict(
-            time=(["time"], pd.to_datetime(["2000-01-01"])),
+            time=(["time"], date_array),
             nz1=(["nz1"], zc),
         ),
         attrs=dict(
@@ -151,7 +158,7 @@ def _fesom2_square_delaunay_uniform_z_coordinate():
         uxgrid=uxgrid,
         dims=["time", "nz1", "n_face"],
         coords=dict(
-            time=(["time"], pd.to_datetime(["2000-01-01"])),
+            time=(["time"], date_array),
             nz1=(["nz1"], zc),
         ),
         attrs=dict(
@@ -164,7 +171,7 @@ def _fesom2_square_delaunay_uniform_z_coordinate():
         uxgrid=uxgrid,
         dims=["time", "nz", "n_node"],
         coords=dict(
-            time=(["time"], pd.to_datetime(["2000-01-01"])),
+            time=(["time"], date_array),
             nz1=(["nz"], zf),
         ),
         attrs=dict(
@@ -177,7 +184,7 @@ def _fesom2_square_delaunay_uniform_z_coordinate():
         uxgrid=uxgrid,
         dims=["time", "nz1", "n_node"],
         coords=dict(
-            time=(["time"], pd.to_datetime(["2000-01-01"])),
+            time=(["time"], date_array),
             nz1=(["nz1"], zc),
         ),
         attrs=dict(description="pressure", units="N/m^2", location="node", mesh="delaunay", Conventions="UGRID-1.0"),

tests/v4/test_uxarray_fieldset.py

@@ -13,6 +13,7 @@
     VectorField,
     download_example_dataset,
 )
+from parcels._datasets.unstructured.generic import datasets as datasets_unstructured
 
 
 @pytest.fixture
@@ -88,3 +89,26 @@ def test_fesom_channel(ds_fesom_channel, uvw_fesom_channel):
 
     pset = ParticleSet(fieldset, pclass=Particle)
     pset.execute(endtime=timedelta(days=1), dt=timedelta(hours=1))
+
+
+@pytest.mark.skip(reason="_index_search.search_time_indices needs major refactoring")
+def test_fesom2_square_delaunay_uniform_z_coordinate_eval():
+    """
+    Test the evaluation of a fieldset with a FESOM2 square Delaunay grid and uniform z-coordinate.
+    Ensures that the fieldset can be created and evaluated correctly.
+    Since the underlying data is constant, we can check that the values are as expected.
+    """
+    ds = datasets_unstructured["fesom2_square_delaunay_uniform_z_coordinate"]
+    UVW = VectorField(
+        name="UVW",
+        U=Field(name="U", data=ds.U, grid=ds.uxgrid, interp_method=UXPiecewiseConstantFace),
+        V=Field(name="V", data=ds.V, grid=ds.uxgrid, interp_method=UXPiecewiseConstantFace),
+        W=Field(name="W", data=ds.W, grid=ds.uxgrid, interp_method=UXPiecewiseLinearNode),
+    )
+    P = Field(name="p", data=ds.p, grid=ds.uxgrid, interp_method=UXPiecewiseConstantFace)
+    fieldset = FieldSet([UVW, P, UVW.U, UVW.V, UVW.W])
+
+    assert fieldset.U.eval(time=ds.time[0].values, z=1.0, y=30.0, x=30.0, applyConversion=False) == 1.0
+    assert fieldset.V.eval(time=ds.time[0].values, z=1.0, y=30.0, x=30.0, applyConversion=False) == 1.0
+    assert fieldset.W.eval(time=ds.time[0].values, z=1.0, y=30.0, x=30.0, applyConversion=False) == 0.0
+    assert fieldset.P.eval(time=ds.time[0].values, z=1.0, y=30.0, x=30.0, applyConversion=False) == 1.0