Merge remote-tracking branch 'origin/v4-dev' into 2031-uxgrid-vertical-search

Author: fluidnumericsJoe

docs/examples/example_globcurrent.py

@@ -30,28 +30,6 @@ def set_globcurrent_fieldset(
     )
 
 
-@pytest.mark.v4remove(
-    reason="indices keryword is not supported in v4. Subsetting should be done on xarray level."
-)
-@pytest.mark.parametrize(
-    "use_xarray", [True, pytest.param(False, marks=pytest.mark.xfail)]
-)
-def test_globcurrent_fieldset(use_xarray):
-    fieldset = set_globcurrent_fieldset()
-    assert fieldset.U.lon.size == 81
-    assert fieldset.U.lat.size == 41
-    assert fieldset.V.lon.size == 81
-    assert fieldset.V.lat.size == 41
-
-    if not use_xarray:
-        indices = {"lon": [5], "lat": range(20, 30)}
-        fieldsetsub = set_globcurrent_fieldset(indices=indices)
-        assert np.allclose(fieldsetsub.U.lon, fieldset.U.lon[indices["lon"]])
-        assert np.allclose(fieldsetsub.U.lat, fieldset.U.lat[indices["lat"]])
-        assert np.allclose(fieldsetsub.V.lon, fieldset.V.lon[indices["lon"]])
-        assert np.allclose(fieldsetsub.V.lat, fieldset.V.lat[indices["lat"]])
-
-
 @pytest.mark.parametrize(
     "dt, lonstart, latstart", [(3600.0, 25, -35), (-3600.0, 20, -39)]
 )
@@ -103,55 +81,6 @@ def test_globcurrent_particles():
     assert abs(pset[0].lat - -35.3) < 1
 
 
-@pytest.mark.v4remove
-@pytest.mark.xfail(
-    reason="Can't patch metadata without using xarray. v4 will natively use xarray anyway. GH1919."
-)
-@pytest.mark.parametrize("dt", [-300, 300])
-def test_globcurrent_xarray_vs_netcdf(dt):
-    fieldsetNetcdf = set_globcurrent_fieldset(use_xarray=False)
-    fieldsetxarray = set_globcurrent_fieldset(use_xarray=True)
-    lonstart, latstart, runtime = (25, -35, timedelta(days=7))
-
-    psetN = parcels.ParticleSet(
-        fieldsetNetcdf, pclass=parcels.Particle, lon=lonstart, lat=latstart
-    )
-    psetN.execute(parcels.AdvectionRK4, runtime=runtime, dt=dt)
-
-    psetX = parcels.ParticleSet(
-        fieldsetxarray, pclass=parcels.Particle, lon=lonstart, lat=latstart
-    )
-    psetX.execute(parcels.AdvectionRK4, runtime=runtime, dt=dt)
-
-    assert np.allclose(psetN[0].lon, psetX[0].lon)
-    assert np.allclose(psetN[0].lat, psetX[0].lat)
-
-
-@pytest.mark.v4remove
-@pytest.mark.xfail(
-    reason="Timeslices will be removed in v4, as users will be able to use xarray directly."
-)
-@pytest.mark.parametrize("dt", [-300, 300])
-def test_globcurrent_netcdf_timestamps(dt):
-    fieldsetNetcdf = set_globcurrent_fieldset()
-    timestamps = fieldsetNetcdf.U.grid.timeslices
-    fieldsetTimestamps = set_globcurrent_fieldset(timestamps=timestamps)
-    lonstart, latstart, runtime = (25, -35, timedelta(days=7))
-
-    psetN = parcels.ParticleSet(
-        fieldsetNetcdf, pclass=parcels.Particle, lon=lonstart, lat=latstart
-    )
-    psetN.execute(parcels.AdvectionRK4, runtime=runtime, dt=dt)
-
-    psetT = parcels.ParticleSet(
-        fieldsetTimestamps, pclass=parcels.Particle, lon=lonstart, lat=latstart
-    )
-    psetT.execute(parcels.AdvectionRK4, runtime=runtime, dt=dt)
-
-    assert np.allclose(psetN.lon[0], psetT.lon[0])
-    assert np.allclose(psetN.lat[0], psetT.lat[0])
-
-
 def test__particles_init_time():
     fieldset = set_globcurrent_fieldset()
 

parcels/_datasets/structured/generic.py

@@ -7,6 +7,8 @@
 
 __all__ = ["T", "X", "Y", "Z", "datasets"]
 
+TIME = xr.date_range("2000", "2001", T)
+
 
 def _rotated_curvilinear_grid():
     XG = np.arange(X)
@@ -44,7 +46,7 @@ def _rotated_curvilinear_grid():
                 {"axis": "Z"},
             ),
             "depth": (["ZG"], np.arange(Z), {"axis": "Z"}),
-            "time": (["time"], xr.date_range("2000", "2001", T), {"axis": "T"}),
+            "time": (["time"], TIME, {"axis": "T"}),
             "lon": (
                 ["YG", "XG"],
                 LON,
@@ -119,7 +121,7 @@ def _unrolled_cone_curvilinear_grid():
                 {"axis": "Z"},
             ),
             "depth": (["ZG"], np.arange(Z), {"axis": "Z"}),
-            "time": (["time"], xr.date_range("2000", "2001", T), {"axis": "T"}),
+            "time": (["time"], TIME, {"axis": "T"}),
             "lon": (
                 ["YG", "XG"],
                 LON,
@@ -175,7 +177,7 @@ def _unrolled_cone_curvilinear_grid():
             "lon": (["XG"], 2 * np.pi / X * np.arange(0, X)),
             "lat": (["YG"], 2 * np.pi / (Y) * np.arange(0, Y)),
             "depth": (["ZG"], np.arange(Z)),
-            "time": (["time"], xr.date_range("2000", "2001", T), {"axis": "T"}),
+            "time": (["time"], TIME, {"axis": "T"}),
         },
     ),
     "2d_left_unrolled_cone": _unrolled_cone_curvilinear_grid(),

parcels/_datasets/unstructured/generic.py

@@ -1,15 +1,16 @@
 import math
-from datetime import datetime, timedelta
 
 import numpy as np
-import pandas as pd
 import uxarray as ux
+import xarray as xr
 
 __all__ = ["Nx", "datasets"]
 
+T = 13
 Nx = 20
 vmax = 1.0
 delta = 0.1
+TIME = xr.date_range("2000", "2001", T)
 
 
 def _stommel_gyre_delaunay():
@@ -57,7 +58,7 @@ def _stommel_gyre_delaunay():
         uxgrid=uxgrid,
         dims=["time", "nz1", "n_node"],
         coords=dict(
-            time=(["time"], pd.to_datetime(["2000-01-01"])),
+            time=(["time"], [TIME[0]]),
             nz1=(["nz1"], [0]),
         ),
         attrs=dict(
@@ -70,7 +71,7 @@ def _stommel_gyre_delaunay():
         uxgrid=uxgrid,
         dims=["time", "nz1", "n_node"],
         coords=dict(
-            time=(["time"], pd.to_datetime(["2000-01-01"])),
+            time=(["time"], [TIME[0]]),
             nz1=(["nz1"], [0]),
         ),
         attrs=dict(
@@ -83,7 +84,7 @@ def _stommel_gyre_delaunay():
         uxgrid=uxgrid,
         dims=["time", "nz1", "n_node"],
         coords=dict(
-            time=(["time"], pd.to_datetime(["2000-01-01"])),
+            time=(["time"], [TIME[0]]),
             nz1=(["nz1"], [0]),
         ),
         attrs=dict(description="pressure", units="N/m^2", location="node", mesh="delaunay", Conventions="UGRID-1.0"),
@@ -108,8 +109,6 @@ 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 = (
@@ -127,25 +126,23 @@ def _fesom2_square_delaunay_uniform_z_coordinate():
 
     # Define arrays U (zonal), V (meridional) and P (sea surface height)
     U = np.ones(
-        (num_days, nz1, uxgrid.n_face), dtype=np.float64
+        (T, nz1, uxgrid.n_face), dtype=np.float64
     )  # Lateral velocity is on the element centers and face centers
     V = np.ones(
-        (num_days, nz1, uxgrid.n_face), dtype=np.float64
+        (T, nz1, uxgrid.n_face), dtype=np.float64
     )  # Lateral velocity is on the element centers and face centers
     W = np.zeros(
-        (num_days, nz, uxgrid.n_node), dtype=np.float64
+        (T, nz, uxgrid.n_node), dtype=np.float64
     )  # Vertical velocity is on the element faces and face vertices
-    P = np.ones(
-        (num_days, nz1, uxgrid.n_node), dtype=np.float64
-    )  # Pressure is on the element centers and face vertices
+    P = np.ones((T, 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"], date_array),
+            time=(["time"], TIME),
             nz1=(["nz1"], zc),
         ),
         attrs=dict(
@@ -158,7 +155,7 @@ def _fesom2_square_delaunay_uniform_z_coordinate():
         uxgrid=uxgrid,
         dims=["time", "nz1", "n_face"],
         coords=dict(
-            time=(["time"], date_array),
+            time=(["time"], TIME),
             nz1=(["nz1"], zc),
         ),
         attrs=dict(
@@ -171,7 +168,7 @@ def _fesom2_square_delaunay_uniform_z_coordinate():
         uxgrid=uxgrid,
         dims=["time", "nz", "n_node"],
         coords=dict(
-            time=(["time"], date_array),
+            time=(["time"], TIME),
             nz=(["nz"], zf),
         ),
         attrs=dict(
@@ -184,7 +181,7 @@ def _fesom2_square_delaunay_uniform_z_coordinate():
         uxgrid=uxgrid,
         dims=["time", "nz1", "n_node"],
         coords=dict(
-            time=(["time"], date_array),
+            time=(["time"], TIME),
             nz1=(["nz1"], zc),
         ),
         attrs=dict(description="pressure", units="N/m^2", location="node", mesh="delaunay", Conventions="UGRID-1.0"),

parcels/_reprs.py

@@ -9,7 +9,7 @@
     from parcels import Field, FieldSet, ParticleSet
 
 
-def field_repr(field: Field) -> str:
+def field_repr(field: Field) -> str:  # TODO v4: Rework or remove entirely
     """Return a pretty repr for Field"""
     out = f"""<{type(field).__name__}>
     name            : {field.name!r}
@@ -63,7 +63,7 @@ def particleset_repr(pset: ParticleSet) -> str:
     return textwrap.dedent(out).strip()
 
 
-def fieldset_repr(fieldset: FieldSet) -> str:
+def fieldset_repr(fieldset: FieldSet) -> str:  # TODO v4: Rework or remove entirely
     """Return a pretty repr for FieldSet"""
     fields_repr = "\n".join([repr(f) for f in fieldset.get_fields()])
 

parcels/field.py

@@ -12,7 +12,7 @@
 
 from parcels._core.utils.time import TimeInterval
 from parcels._core.utils.unstructured import get_vertical_location_from_dims
-from parcels._reprs import default_repr, field_repr
+from parcels._reprs import default_repr
 from parcels._typing import (
     Mesh,
     VectorType,
@@ -212,9 +212,6 @@ def __init__(
         if "time" not in self.data.dims:
             raise ValueError("Field is missing a 'time' dimension. ")
 
-    def __repr__(self):
-        return field_repr(self)
-
     @property
     def units(self):
         return self._units

parcels/fieldset.py

@@ -9,7 +9,6 @@
 
 from parcels._core.utils.time import get_datetime_type_calendar
 from parcels._core.utils.time import is_compatible as datetime_is_compatible
-from parcels._reprs import fieldset_repr
 from parcels._typing import Mesh
 from parcels.field import Field, VectorField
 from parcels.v4.grid import Grid
@@ -79,9 +78,6 @@ def time_interval(self):
         time_intervals = (t for t in time_intervals if t is not None)
         return functools.reduce(lambda x, y: x.intersection(y), time_intervals)
 
-    def __repr__(self):
-        return fieldset_repr(self)
-
     def dimrange(self, dim):
         """Returns maximum value of a dimension (lon, lat, depth or time)
         on 'left' side and minimum value on 'right' side for all grids

parcels/v4/gridadapter.py

@@ -70,13 +70,17 @@ def depth(self):
         return self.grid._ds["depth"].values
 
     @property
-    def time(self):
+    def _datetimes(self):
         try:
             axis = self.grid.axes["T"]
         except KeyError:
             return np.zeros(1)
         return get_time(axis)
 
+    @property
+    def time(self):
+        return self._datetimes.astype(np.float64) / 1e9
+
     @property
     def xdim(self):
         return get_dimensionality(self.grid.axes.get("X"))
@@ -95,7 +99,7 @@ def tdim(self):
 
     @property
     def time_origin(self):
-        return TimeConverter(self.time[0])
+        return TimeConverter(self._datetimes[0])
 
     @property
     def _z4d(self) -> Literal[0, 1]: