Removing transpose option in Field creation

As not needed anymore when using xarray under the hood in Parcels v4

Author: erikvansebille

docs/examples/example_moving_eddies.py

@@ -65,9 +65,9 @@ def cosd(x):
     dy = (lat[1] - lat[0]) * 1852 * 60 if mesh == "spherical" else lat[1] - lat[0]
 
     # Define arrays U (zonal), V (meridional), and P (sea surface height) on A-grid
-    U = np.zeros((lon.size, lat.size, time.size), dtype=np.float32)
-    V = np.zeros((lon.size, lat.size, time.size), dtype=np.float32)
-    P = np.zeros((lon.size, lat.size, time.size), dtype=np.float32)
+    U = np.zeros((time.size, lat.size, lon.size), dtype=np.float32)
+    V = np.zeros((time.size, lat.size, lon.size), dtype=np.float32)
+    P = np.zeros((time.size, lat.size, lon.size), dtype=np.float32)
 
     # Some constants
     corio_0 = 1.0e-4  # Coriolis parameter
@@ -78,32 +78,32 @@ def cosd(x):
     dX = eddyspeed * 86400 / dx  # Grid cell movement of eddy max each day
     dY = eddyspeed * 86400 / dy  # Grid cell movement of eddy max each day
 
-    [x, y] = np.mgrid[: lon.size, : lat.size]
+    [y, x] = np.mgrid[: lat.size, : lon.size]
     for t in range(time.size):
         hymax_1 = lat.size / 7.0
         hxmax_1 = 0.75 * lon.size - dX * t
         hymax_2 = 3.0 * lat.size / 7.0 + dY * t
         hxmax_2 = 0.75 * lon.size - dX * t
 
-        P[:, :, t] = h0 * np.exp(
+        P[t, :, :] = h0 * np.exp(
             -((x - hxmax_1) ** 2) / (sig * lon.size / 4.0) ** 2
             - (y - hymax_1) ** 2 / (sig * lat.size / 7.0) ** 2
         )
-        P[:, :, t] += h0 * np.exp(
+        P[t, :, :] += h0 * np.exp(
             -((x - hxmax_2) ** 2) / (sig * lon.size / 4.0) ** 2
             - (y - hymax_2) ** 2 / (sig * lat.size / 7.0) ** 2
         )
 
-        V[:-1, :, t] = -np.diff(P[:, :, t], axis=0) / dx / corio_0 * g
-        V[-1, :, t] = V[-2, :, t]  # Fill in the last column
+        V[t, :, :-1] = -np.diff(P[t, :, :], axis=1) / dx / corio_0 * g
+        V[t, :, -1] = V[t, :, -2]  # Fill in the last column
 
-        U[:, :-1, t] = np.diff(P[:, :, t], axis=1) / dy / corio_0 * g
-        U[:, -1, t] = U[:, -2, t]  # Fill in the last row
+        U[t, :-1, :] = np.diff(P[t, :, :], axis=0) / dy / corio_0 * g
+        U[t, -1, :] = U[t, -2, :]  # Fill in the last row
 
     data = {"U": U, "V": V, "P": P}
     dimensions = {"lon": lon, "lat": lat, "time": time}
 
-    fieldset = parcels.FieldSet.from_data(data, dimensions, transpose=True, mesh=mesh)
+    fieldset = parcels.FieldSet.from_data(data, dimensions, mesh=mesh)
 
     # setting some constants for AdvectionRK45 kernel
     fieldset.RK45_min_dt = 1e-3

parcels/field.py

@@ -108,13 +108,7 @@ class Field:
     name : str
         Name of the field
     data : np.ndarray
-        2D, 3D or 4D numpy array of field data.
-
-        1. If data shape is [xdim, ydim], [xdim, ydim, zdim], [xdim, ydim, tdim] or [xdim, ydim, zdim, tdim],
-           whichever is relevant for the dataset, use the flag transpose=True
-        2. If data shape is [ydim, xdim], [zdim, ydim, xdim], [tdim, ydim, xdim] or [tdim, zdim, ydim, xdim],
-           use the flag transpose=False
-        3. If data has any other shape, you first need to reorder it
+        2D, 3D or 4D numpy array of field data with shape [ydim, xdim], [zdim, ydim, xdim], [tdim, ydim, xdim] or [tdim, zdim, ydim, xdim],
     lon : np.ndarray or list
         Longitude coordinates (numpy vector or array) of the field (only if grid is None)
     lat : np.ndarray or list
@@ -138,8 +132,6 @@ class Field:
         mesh and time_origin information. Can be constructed from any of the Grid objects
     fieldtype : str
         Type of Field to be used for UnitConverter (either 'U', 'V', 'Kh_zonal', 'Kh_meridional' or None)
-    transpose : bool
-        Transpose data to required (lon, lat) layout
     vmin : float
         Minimum allowed value on the field. Data below this value are set to zero
     vmax : float
@@ -183,7 +175,6 @@ def __init__(
         mesh: Mesh = "flat",
         timestamps=None,
         fieldtype=None,
-        transpose: bool = False,
         vmin: float | None = None,
         vmax: float | None = None,
         cast_data_dtype: type[np.float32] | type[np.float64] | Literal["float32", "float64"] = "float32",
@@ -271,7 +262,7 @@ def __init__(
             )
 
         if not self.grid.defer_load:
-            self.data = self._reshape(self.data, transpose)
+            self.data = self._reshape(self.data)
             self._loaded_time_indices = range(self.grid.tdim)
 
             # Hack around the fact that NaN and ridiculously large values
@@ -668,7 +659,7 @@ def from_netcdf(
                         errormessage = (
                             f"Field {filebuffer.name} expecting a data shape of [tdim={grid.tdim}, zdim={grid.zdim}, "
                             f"ydim={grid.ydim - 2 * grid.meridional_halo}, xdim={grid.xdim - 2 * grid.zonal_halo}] "
-                            f"but got shape {buffer_data.shape}. Flag transpose=True could help to reorder the data."
+                            f"but got shape {buffer_data.shape}."
                         )
                         assert buffer_data.shape[0] == grid.tdim, errormessage
                         assert buffer_data.shape[2] == grid.ydim - 2 * grid.meridional_halo, errormessage
@@ -768,7 +759,7 @@ def from_xarray(
             **kwargs,
         )
 
-    def _reshape(self, data, transpose=False):
+    def _reshape(self, data):
         # Ensure that field data is the right data type
         if not isinstance(data, (np.ndarray, da.core.Array)):
             data = np.array(data)
@@ -777,8 +768,6 @@ def _reshape(self, data, transpose=False):
         elif (self.cast_data_dtype == np.float64) and (data.dtype != np.float64):
             data = data.astype(np.float64)
         lib = np if isinstance(data, np.ndarray) else da
-        if transpose:
-            data = lib.transpose(data)
         if self.grid._lat_flipped:
             data = lib.flip(data, axis=-2)
 
@@ -803,10 +792,7 @@ def _reshape(self, data, transpose=False):
             if len(data.shape) == 4:
                 data = data.reshape(sum(((data.shape[0],), data.shape[2:]), ()))
         if len(data.shape) == 4:
-            errormessage = (
-                f"Field {self.name} expecting a data shape of [tdim, zdim, ydim, xdim]. "
-                "Flag transpose=True could help to reorder the data."
-            )
+            errormessage = f"Field {self.name} expecting a data shape of [tdim, zdim, ydim, xdim]. "
             assert data.shape[0] == self.grid.tdim, errormessage
             assert data.shape[2] == self.grid.ydim - 2 * self.grid.meridional_halo, errormessage
             assert data.shape[3] == self.grid.xdim - 2 * self.grid.zonal_halo, errormessage
@@ -819,10 +805,7 @@ def _reshape(self, data, transpose=False):
                 self.grid.tdim,
                 self.grid.ydim - 2 * self.grid.meridional_halo,
                 self.grid.xdim - 2 * self.grid.zonal_halo,
-            ), (
-                f"Field {self.name} expecting a data shape of [tdim, ydim, xdim]. "
-                "Flag transpose=True could help to reorder the data."
-            )
+            ), f"Field {self.name} expecting a data shape of [tdim, ydim, xdim]. "
         if self.grid.meridional_halo > 0 or self.grid.zonal_halo > 0:
             data = self.add_periodic_halo(
                 zonal=self.grid.zonal_halo > 0,

parcels/fieldset.py

@@ -74,7 +74,6 @@ def from_data(
         cls,
         data,
         dimensions,
-        transpose=False,
         mesh: Mesh = "spherical",
         allow_time_extrapolation: bool | None = None,
         **kwargs,
@@ -86,21 +85,14 @@ def from_data(
         data :
             Dictionary mapping field names to numpy arrays.
             Note that at least a 'U' and 'V' numpy array need to be given, and that
-            the built-in Advection kernels assume that U and V are in m/s
-
-            1. If data shape is [xdim, ydim], [xdim, ydim, zdim], [xdim, ydim, tdim] or [xdim, ydim, zdim, tdim],
-               whichever is relevant for the dataset, use the flag transpose=True
-            2. If data shape is [ydim, xdim], [zdim, ydim, xdim], [tdim, ydim, xdim] or [tdim, zdim, ydim, xdim],
-               use the flag transpose=False (default value)
-            3. If data has any other shape, you first need to reorder it
+            the built-in Advection kernels assume that U and V are in m/s.
+            Data shape is either [ydim, xdim], [zdim, ydim, xdim], [tdim, ydim, xdim] or [tdim, zdim, ydim, xdim],
         dimensions : dict
             Dictionary mapping field dimensions (lon,
             lat, depth, time) to numpy arrays.
             Note that dimensions can also be a dictionary of dictionaries if
             dimension names are different for each variable
             (e.g. dimensions['U'], dimensions['V'], etc).
-        transpose : bool
-            Whether to transpose data on read-in (Default value = False)
         mesh : str
             String indicating the type of mesh coordinates and
             units used during velocity interpolation, see also `this tutorial <../examples/tutorial_unitconverters.ipynb>`__:
@@ -154,7 +146,6 @@ def from_data(
                 name,
                 datafld,
                 grid=grid,
-                transpose=transpose,
                 allow_time_extrapolation=allow_time_extrapolation,
                 **kwargs,
             )

tests/test_advection.py

@@ -61,22 +61,22 @@ def test_advection_zonal(lon, lat, depth):
     """Particles at high latitude move geographically faster due to the pole correction in `GeographicPolar`."""
     npart = 10
     data2D = {
-        "U": np.ones((lon.size, lat.size), dtype=np.float32),
-        "V": np.zeros((lon.size, lat.size), dtype=np.float32),
+        "U": np.ones((lat.size, lon.size), dtype=np.float32),
+        "V": np.zeros((lat.size, lon.size), dtype=np.float32),
     }
     data3D = {
-        "U": np.ones((lon.size, lat.size, depth.size), dtype=np.float32),
-        "V": np.zeros((lon.size, lat.size, depth.size), dtype=np.float32),
+        "U": np.ones((depth.size, lat.size, lon.size), dtype=np.float32),
+        "V": np.zeros((depth.size, lat.size, lon.size), dtype=np.float32),
     }
     dimensions = {"lon": lon, "lat": lat}
-    fieldset2D = FieldSet.from_data(data2D, dimensions, mesh="spherical", transpose=True)
+    fieldset2D = FieldSet.from_data(data2D, dimensions, mesh="spherical")
 
     pset2D = ParticleSet(fieldset2D, pclass=Particle, lon=np.zeros(npart) + 20.0, lat=np.linspace(0, 80, npart))
     pset2D.execute(AdvectionRK4, runtime=timedelta(hours=2), dt=timedelta(seconds=30))
     assert (np.diff(pset2D.lon) > 1.0e-4).all()
 
     dimensions["depth"] = depth
-    fieldset3D = FieldSet.from_data(data3D, dimensions, mesh="spherical", transpose=True)
+    fieldset3D = FieldSet.from_data(data3D, dimensions, mesh="spherical")
     pset3D = ParticleSet(
         fieldset3D,
         pclass=Particle,
@@ -91,9 +91,9 @@ def test_advection_zonal(lon, lat, depth):
 def test_advection_meridional(lon, lat):
     """Particles at high latitude move geographically faster due to the pole correction in `GeographicPolar`."""
     npart = 10
-    data = {"U": np.zeros((lon.size, lat.size), dtype=np.float32), "V": np.ones((lon.size, lat.size), dtype=np.float32)}
+    data = {"U": np.zeros((lat.size, lon.size), dtype=np.float32), "V": np.ones((lat.size, lon.size), dtype=np.float32)}
     dimensions = {"lon": lon, "lat": lat}
-    fieldset = FieldSet.from_data(data, dimensions, mesh="spherical", transpose=True)
+    fieldset = FieldSet.from_data(data, dimensions, mesh="spherical")
 
     pset = ParticleSet(fieldset, pclass=Particle, lon=np.linspace(-60, 60, npart), lat=np.linspace(0, 30, npart))
     delta_lat = np.diff(pset.lat)
@@ -110,9 +110,9 @@ def test_advection_3D():
         "lat": np.linspace(0.0, 1e4, ydim, dtype=np.float32),
         "depth": np.linspace(0.0, 1.0, zdim, dtype=np.float32),
     }
-    data = {"U": np.ones((xdim, ydim, zdim), dtype=np.float32), "V": np.zeros((xdim, ydim, zdim), dtype=np.float32)}
-    data["U"][:, :, 0] = 0.0
-    fieldset = FieldSet.from_data(data, dimensions, mesh="flat", transpose=True)
+    data = {"U": np.ones((zdim, ydim, xdim), dtype=np.float32), "V": np.zeros((zdim, ydim, xdim), dtype=np.float32)}
+    data["U"][0, :, :] = 0.0
+    fieldset = FieldSet.from_data(data, dimensions, mesh="flat")
 
     pset = ParticleSet(
         fieldset, pclass=Particle, lon=np.zeros(npart), lat=np.zeros(npart) + 1e2, depth=np.linspace(0, 1, npart)
@@ -171,11 +171,11 @@ def SubmergeParticle(particle, fieldset, time):  # pragma: no cover
 def test_advection_RK45(lon, lat, rk45_tol):
     npart = 10
     data2D = {
-        "U": np.ones((lon.size, lat.size), dtype=np.float32),
-        "V": np.zeros((lon.size, lat.size), dtype=np.float32),
+        "U": np.ones((lat.size, lon.size), dtype=np.float32),
+        "V": np.zeros((lat.size, lon.size), dtype=np.float32),
     }
     dimensions = {"lon": lon, "lat": lat}
-    fieldset = FieldSet.from_data(data2D, dimensions, mesh="spherical", transpose=True)
+    fieldset = FieldSet.from_data(data2D, dimensions, mesh="spherical")
     fieldset.add_constant("RK45_tol", rk45_tol)
 
     dt = timedelta(seconds=30).total_seconds()
@@ -266,8 +266,8 @@ def create_periodic_fieldset(xdim, ydim, uvel, vvel):
         "lat": np.linspace(0.0, 1.0, ydim + 1, dtype=np.float32)[1:],
     }
 
-    data = {"U": uvel * np.ones((xdim, ydim), dtype=np.float32), "V": vvel * np.ones((xdim, ydim), dtype=np.float32)}
-    return FieldSet.from_data(data, dimensions, mesh="spherical", transpose=True)
+    data = {"U": uvel * np.ones((ydim, xdim), dtype=np.float32), "V": vvel * np.ones((ydim, xdim), dtype=np.float32)}
+    return FieldSet.from_data(data, dimensions, mesh="spherical")
 
 
 def periodicBC(particle, fieldset, time):  # pragma: no cover
@@ -373,10 +373,10 @@ def create_fieldset_stationary(xdim=100, ydim=100, maxtime=timedelta(hours=6)):
         "time": time,
     }
     data = {
-        "U": np.ones((xdim, ydim, 1), dtype=np.float32) * u_0 * np.cos(f * time),
-        "V": np.ones((xdim, ydim, 1), dtype=np.float32) * -u_0 * np.sin(f * time),
+        "U": np.transpose(np.ones((xdim, ydim, 1), dtype=np.float32) * u_0 * np.cos(f * time)),
+        "V": np.transpose(np.ones((xdim, ydim, 1), dtype=np.float32) * -u_0 * np.sin(f * time)),
     }
-    fieldset = FieldSet.from_data(data, dimensions, mesh="flat", transpose=True)
+    fieldset = FieldSet.from_data(data, dimensions, mesh="flat")
     # setting some constants for AdvectionRK45 kernel
     fieldset.RK45_min_dt = 1e-3
     fieldset.RK45_max_dt = 1e2
@@ -435,13 +435,13 @@ def test_stationary_eddy_vertical():
     lon_data = np.linspace(0, 25000, xdim, dtype=np.float32)
     lat_data = np.linspace(0, 25000, ydim, dtype=np.float32)
     time_data = np.arange(0.0, 6 * 3600 + 1e-5, 60.0, dtype=np.float64)
-    fld1 = np.ones((xdim, ydim, 1), dtype=np.float32) * u_0 * np.cos(f * time_data)
-    fld2 = np.ones((xdim, ydim, 1), dtype=np.float32) * -u_0 * np.sin(f * time_data)
-    fldzero = np.zeros((xdim, ydim, 1), dtype=np.float32) * time_data
+    fld1 = np.transpose(np.ones((xdim, ydim, 1), dtype=np.float32) * u_0 * np.cos(f * time_data))
+    fld2 = np.transpose(np.ones((xdim, ydim, 1), dtype=np.float32) * -u_0 * np.sin(f * time_data))
+    fldzero = np.transpose(np.zeros((xdim, ydim, 1), dtype=np.float32) * time_data)
 
     dimensions = {"lon": lon_data, "lat": lat_data, "time": time_data}
     data = {"U": fld1, "V": fldzero, "W": fld2}
-    fieldset = FieldSet.from_data(data, dimensions, mesh="flat", transpose=True)
+    fieldset = FieldSet.from_data(data, dimensions, mesh="flat")
 
     pset = ParticleSet(fieldset, pclass=Particle, lon=lon, lat=lat, depth=depth)
     pset.execute(AdvectionRK4_3D, dt=dt, endtime=endtime)
@@ -453,7 +453,7 @@ def test_stationary_eddy_vertical():
     assert np.allclose(pset.depth, exp_depth, rtol=1e-5)
 
     data = {"U": fldzero, "V": fld2, "W": fld1}
-    fieldset = FieldSet.from_data(data, dimensions, mesh="flat", transpose=True)
+    fieldset = FieldSet.from_data(data, dimensions, mesh="flat")
 
     pset = ParticleSet(fieldset, pclass=Particle, lon=lon, lat=lat, depth=depth)
     pset.execute(AdvectionRK4_3D, dt=dt, endtime=endtime)
@@ -483,10 +483,10 @@ def create_fieldset_moving(xdim=100, ydim=100, maxtime=timedelta(hours=6)):
         "time": time,
     }
     data = {
-        "U": np.ones((xdim, ydim, 1), dtype=np.float32) * u_g + (u_0 - u_g) * np.cos(f * time),
-        "V": np.ones((xdim, ydim, 1), dtype=np.float32) * -(u_0 - u_g) * np.sin(f * time),
+        "U": np.transpose(np.ones((xdim, ydim, 1), dtype=np.float32) * u_g + (u_0 - u_g) * np.cos(f * time)),
+        "V": np.transpose(np.ones((xdim, ydim, 1), dtype=np.float32) * -(u_0 - u_g) * np.sin(f * time)),
     }
-    return FieldSet.from_data(data, dimensions, mesh="flat", transpose=True)
+    return FieldSet.from_data(data, dimensions, mesh="flat")
 
 
 @pytest.fixture
@@ -555,11 +555,15 @@ def create_fieldset_decaying(xdim=100, ydim=100, maxtime=timedelta(hours=6)):
         "time": time,
     }
     data = {
-        "U": np.ones((xdim, ydim, 1), dtype=np.float32) * u_g * np.exp(-gamma_g * time)
-        + (u_0 - u_g) * np.exp(-gamma * time) * np.cos(f * time),
-        "V": np.ones((xdim, ydim, 1), dtype=np.float32) * -(u_0 - u_g) * np.exp(-gamma * time) * np.sin(f * time),
+        "U": np.transpose(
+            np.ones((xdim, ydim, 1), dtype=np.float32) * u_g * np.exp(-gamma_g * time)
+            + (u_0 - u_g) * np.exp(-gamma * time) * np.cos(f * time)
+        ),
+        "V": np.transpose(
+            np.ones((xdim, ydim, 1), dtype=np.float32) * -(u_0 - u_g) * np.exp(-gamma * time) * np.sin(f * time)
+        ),
     }
-    return FieldSet.from_data(data, dimensions, mesh="flat", transpose=True)
+    return FieldSet.from_data(data, dimensions, mesh="flat")
 
 
 @pytest.fixture

tests/test_data/create_testfields.py

@@ -29,12 +29,12 @@ def generate_testfieldset(xdim, ydim, zdim, tdim):
     lat = np.linspace(0.0, 1.0, ydim, dtype=np.float32)
     depth = np.linspace(0.0, 0.5, zdim, dtype=np.float32)
     time = np.linspace(0.0, tdim, tdim, dtype=np.float64)
-    U = np.ones((xdim, ydim, zdim, tdim), dtype=np.float32)
-    V = np.zeros((xdim, ydim, zdim, tdim), dtype=np.float32)
-    P = 2.0 * np.ones((xdim, ydim, zdim, tdim), dtype=np.float32)
+    U = np.ones((tdim, zdim, ydim, xdim), dtype=np.float32)
+    V = np.zeros((tdim, zdim, ydim, xdim), dtype=np.float32)
+    P = 2.0 * np.ones((tdim, zdim, ydim, xdim), dtype=np.float32)
     data = {"U": U, "V": V, "P": P}
     dimensions = {"lon": lon, "lat": lat, "depth": depth, "time": time}
-    fieldset = FieldSet.from_data(data, dimensions, mesh="flat", transpose=True)
+    fieldset = FieldSet.from_data(data, dimensions, mesh="flat")
     fieldset.write("testfields")
 
 
@@ -66,7 +66,7 @@ def generate_perlin_testfield():
     if asizeof is not None:
         print(f"Perlin U-field requires {U.size * U.itemsize} bytes of memory.")
         print(f"Perlin V-field requires {V.size * V.itemsize} bytes of memory.")
-    fieldset = FieldSet.from_data(data, dimensions, mesh="spherical", transpose=False)
+    fieldset = FieldSet.from_data(data, dimensions, mesh="spherical")
     # fieldset.write("perlinfields")  # can also be used, but then has a ghost depth dimension
     write_simple_2Dt(fieldset.U, os.path.join(os.path.dirname(__file__), "perlinfields"), varname="vozocrtx")
     write_simple_2Dt(fieldset.V, os.path.join(os.path.dirname(__file__), "perlinfields"), varname="vomecrty")