Rename particle position attributes from "lon" and "lat" to "x" and "y" (#2719)

* Refactor lon lat to x y

Update particle definition and all related code

- particle.py particle variables
- ParticleSet args
- all internal mentions of .lon and .lat
- test suite code

* Update repr

* Update interpolators

* Update error messages

* Update local varnames in kernels

* Update migration guide

* Review feedback

* Restore v3 files

* Update docs/user_guide/v4-migration.md

Co-authored-by: Erik van Sebille <e.vansebille@uu.nl>

---------

Co-authored-by: Erik van Sebille <e.vansebille@uu.nl>

Author: VeckoTheGecko

docs/user_guide/v4-migration.md

@@ -18,6 +18,7 @@ Version 4 of Parcels is unreleased at the moment. The information in this migrat
 - Users need to explicitly use `convert_z_to_sigma_croco` in sampling kernels (such as the `AdvectionRK4_3D_CROCO` or `SampleOMegaCroco` kernels) when working with CROCO data, as the automatic conversion from depth to sigma grids under the hood has been removed.
 - We added a new AdvectionRK2 Kernel. The AdvectionRK4 kernel is still available, but RK2 is now the recommended default advection scheme as it is faster while the accuracy is comparable for most applications. See also the Choosing an integration method tutorial.
 - Functions shouldn't be converted to Kernels before adding to a pset.execute() call. Instead, simply pass the function(s) as a list to pset.execute().
+- Kernel variables `lat` and `lon` have been renamed to `y` and `x`, and `dlat` and `dlon` have been renamed to `dy` and `dx`. These changes are also reflected on the ParticleSet as well as the ParticleFile output.
 
 ## FieldSet
 

src/parcels/_compat.py

@@ -1,7 +1,7 @@
 """Import helpers for compatability between installations."""
 
 
-# for compat with v3 of parcels when users provide `initial=attrgetter("lon")` to a Variable
+# for compat with v3 of parcels when users provide `initial=attrgetter("...")` to a Variable
 # so that particle initial state matches another variable
 class _AttrgetterHelper:
     """

src/parcels/_core/field.py

@@ -184,7 +184,7 @@ def __getitem__(self, key):
         self._check_velocitysampling()
         try:
             if isinstance(key, ParticleSetView):
-                return self.eval(key.time, key.z, key.lat, key.lon, key)
+                return self.eval(key.time, key.z, key.y, key.x, key)
             else:
                 return self.eval(*key)
         except tuple(AllParcelsErrorCodes.keys()) as error:
@@ -293,7 +293,7 @@ def eval(self, time: datetime, z, y, x, particles=None):
     def __getitem__(self, key):
         try:
             if isinstance(key, ParticleSetView):
-                return self.eval(key.time, key.z, key.lat, key.lon, key)
+                return self.eval(key.time, key.z, key.y, key.x, key)
             else:
                 return self.eval(*key)
         except tuple(AllParcelsErrorCodes.keys()) as error:
@@ -389,7 +389,7 @@ def _assert_same_time_interval(fields: Sequence[Field]) -> None:
 
 def _get_positions(field: Field, time, z, y, x, particles, _ei) -> tuple[dict, dict]:
     """Initialize and populate particle_positions and grid_positions dictionaries"""
-    particle_positions = {"time": time, "z": z, "lat": y, "lon": x}
+    particle_positions = {"time": time, "z": z, "y": y, "x": x}
     grid_positions = {}
     grid_positions.update(_search_time_index(field, time))
     grid_positions.update(field.grid.search(z, y, x, ei=_ei))

src/parcels/_core/kernel.py

@@ -106,13 +106,13 @@ def remove_deleted(self, pset):
             pset.remove_indices(indices)
 
     def _position_update(self, particles, fieldset):
-        particles.lon += particles.dlon
-        particles.lat += particles.dlat
+        particles.x += particles.dx
+        particles.y += particles.dy
         particles.z += particles.dz
         particles.time += particles.dt
 
-        particles.dlon = 0
-        particles.dlat = 0
+        particles.dx = 0
+        particles.dy = 0
         particles.dz = 0
 
         if hasattr(self.fieldset, "RK45_tol"):
@@ -244,6 +244,6 @@ def execute(self, pset, endtime, dt):
                     if error_code == StatusCode.ErrorOutsideTimeInterval:
                         error_func(pset[inds].time)
                     else:
-                        error_func(pset[inds].z, pset[inds].lat, pset[inds].lon)
+                        error_func(pset[inds].z, pset[inds].y, pset[inds].x)
 
         return pset

src/parcels/_core/particle.py

@@ -141,18 +141,26 @@ def get_default_particle(spatial_dtype: type[np.float32] | type[np.float64]) ->
                 attrs={"standard_name": "vertical coordinate", "units": "m", "positive": "down"},
             ),
             Variable(
-                "lat",
+                "y",
                 dtype=spatial_dtype,
-                attrs={"standard_name": "latitude", "units": "degrees_north", "axis": "Y"},
+                attrs={
+                    "standard_name": "latitude",
+                    "units": "degrees_north",
+                    "axis": "Y",
+                },  # TODO v4: https://github.com/Parcels-code/Parcels/issues/2720
             ),
             Variable(
-                "lon",
+                "x",
                 dtype=spatial_dtype,
-                attrs={"standard_name": "longitude", "units": "degrees_east", "axis": "X"},
+                attrs={
+                    "standard_name": "longitude",
+                    "units": "degrees_east",
+                    "axis": "X",
+                },  # TODO v4: https://github.com/Parcels-code/Parcels/issues/2720
             ),
             Variable("dz", dtype=spatial_dtype, to_write=False),
-            Variable("dlat", dtype=spatial_dtype, to_write=False),
-            Variable("dlon", dtype=spatial_dtype, to_write=False),
+            Variable("dy", dtype=spatial_dtype, to_write=False),
+            Variable("dx", dtype=spatial_dtype, to_write=False),
             Variable(
                 "particle_id",
                 dtype=np.int64,

src/parcels/_core/particleset.py

@@ -60,8 +60,8 @@ def __init__(
         pclass=Particle,
         time=None,
         z=None,
-        lat=None,
-        lon=None,
+        y=None,
+        x=None,
         particle_ids=None,
         **kwargs,
     ):
@@ -70,21 +70,21 @@ def __init__(
 
         self.fieldset = fieldset
         time = np.empty(shape=0) if time is None else np.array(time).flatten()
-        lat = np.empty(shape=0) if lat is None else np.array(lat).flatten()
-        lon = np.empty(shape=0) if lon is None else np.array(lon).flatten()
+        y = np.empty(shape=0) if y is None else np.array(y).flatten()
+        x = np.empty(shape=0) if x is None else np.array(x).flatten()
 
         if particle_ids is None:
-            particle_ids = np.arange(lon.size)
+            particle_ids = np.arange(x.size)
 
         if z is None:
             minz = 0
             for field in self.fieldset.fields.values():
                 if field.grid.depth is not None:
                     minz = min(minz, field.grid.depth[0])
-            z = np.ones(lon.size) * minz
+            z = np.ones(x.size) * minz
         else:
             z = np.array(z).flatten()
-        assert lon.size == lat.size and lon.size == z.size, "lon, lat, z don't all have the same lenghts"
+        assert x.size == y.size and x.size == z.size, "lon, lat, z don't all have the same lenghts"
 
         if time is None or len(time) == 0:
             # do not set a time yet (because sign_dt not known)
@@ -95,26 +95,26 @@ def __init__(
             time = timedelta_to_float(time)
         else:
             raise TypeError("particle time must be a datetime, timedelta, or date object")
-        time = np.repeat(time, lon.size) if time.size == 1 else time
+        time = np.repeat(time, x.size) if time.size == 1 else time
 
-        assert lon.size == time.size, "time and positions (lon, lat, z) do not have the same lengths."
+        assert x.size == time.size, "time and positions (lon, lat, z) do not have the same lengths."
 
         if fieldset.time_interval:
             _warn_particle_times_outside_fieldset_time_bounds(time, fieldset.time_interval)
 
         for kwvar in kwargs:
             kwargs[kwvar] = np.array(kwargs[kwvar]).flatten()
-            assert lon.size == kwargs[kwvar].size, f"{kwvar} and positions (lon, lat, z) don't have the same lengths."
+            assert x.size == kwargs[kwvar].size, f"{kwvar} and positions (lon, lat, z) don't have the same lengths."
 
         self._data = create_particle_data(
             pclass=pclass,
-            nparticles=lon.size,
+            nparticles=x.size,
             ngrids=len(fieldset.gridset),
             initial=dict(
                 time=time,
                 z=z,
-                lat=lat,
-                lon=lon,
+                y=y,
+                x=x,
                 particle_id=particle_ids,
             ),
         )
@@ -245,7 +245,7 @@ def remove_indices(self, indices):
     def populate_indices(self):
         """Pre-populate guesses of particle ei (element id) indices"""
         for i, grid in enumerate(self.fieldset.gridset):
-            grid_positions = grid.search(self.z, self.lat, self.lon)
+            grid_positions = grid.search(self.z, self.y, self.x)
             self._data["ei"][:, i] = grid.ravel_index(
                 {
                     "X": grid_positions["X"]["index"],

src/parcels/_core/statuscodes.py

@@ -41,7 +41,7 @@ class FieldInterpolationError(RuntimeError):
 
 
 def _raise_field_interpolation_error(z, y, x):
-    raise FieldInterpolationError(f"Field interpolation returned NaN at (z={z}, lat={y}, lon={x})")
+    raise FieldInterpolationError(f"Field interpolation returned NaN at (z={z}, y={y}, x={x})")
 
 
 class FieldOutOfBoundError(RuntimeError):
@@ -51,7 +51,7 @@ class FieldOutOfBoundError(RuntimeError):
 
 
 def _raise_field_out_of_bound_error(z, y, x):
-    raise FieldOutOfBoundError(f"Field sampled out-of-bound, at (z={z}, lat={y}, lon={x})")
+    raise FieldOutOfBoundError(f"Field sampled out-of-bound, at (z={z}, y={y}, x={x})")
 
 
 class FieldOutOfBoundSurfaceError(RuntimeError):
@@ -65,7 +65,7 @@ def format_out(val):
         return "unknown" if val is None else val
 
     raise FieldOutOfBoundSurfaceError(
-        f"Field sampled out-of-bound at the surface, at (z={format_out(z)}, lat={format_out(y)}, lon={format_out(x)})"
+        f"Field sampled out-of-bound at the surface, at (z={format_out(z)}, y={format_out(y)}, x={format_out(x)})"
     )
 
 
@@ -82,7 +82,7 @@ class GridSearchingError(RuntimeError):
 
 
 def _raise_grid_searching_error(z, y, x):
-    raise GridSearchingError(f"Grid searching failed at (z={z}, lat={y}, lon={x})")
+    raise GridSearchingError(f"Grid searching failed at (z={z}, y={y}, x={x})")
 
 
 class GeneralError(RuntimeError):
@@ -92,7 +92,7 @@ class GeneralError(RuntimeError):
 
 
 def _raise_general_error(z, y, x):
-    raise GeneralError(f"General error occurred at (z={z}, lat={y}, lon={x})")
+    raise GeneralError(f"General error occurred at (z={z}, y={y}, x={x})")
 
 
 class OutsideTimeInterval(RuntimeError):

src/parcels/_reprs.py

@@ -97,8 +97,8 @@ def particleset_repr(pset: ParticleSet) -> str:
 def particlesetview_repr(pview: Any) -> str:
     """Return a pretty repr for ParticleSetView"""
     time_string = "not_yet_set" if pview.time is None or np.isnan(pview.time) else f"{pview.time:f}"
-    out = f"P[{pview.particle_id}]: time={time_string}, z={pview.z:f}, lat={pview.lat:f}, lon={pview.lon:f}"
-    vars = [v.name for v in pview._pclass.variables if v.to_write is True and v.name not in ["lon", "lat", "z", "time"]]
+    out = f"P[{pview.particle_id}]: time={time_string}, z={pview.z:f}, y={pview.y:f}, x={pview.x:f}"
+    vars = [v.name for v in pview._pclass.variables if v.to_write is True and v.name not in ["z", "y", "x", "time"]]
     for var in vars:
         out += f", {var}={getattr(pview, var):f}"
 

src/parcels/interpolators/_uxinterpolators.py

@@ -171,7 +171,7 @@ def interp(
         u = vectorfield.U.interp_method.interp(particle_positions, grid_positions, vectorfield.U)
         v = vectorfield.V.interp_method.interp(particle_positions, grid_positions, vectorfield.V)
         if vectorfield.grid._mesh == "spherical":
-            u /= 1852 * 60 * np.cos(np.deg2rad(particle_positions["lat"]))
+            u /= 1852 * 60 * np.cos(np.deg2rad(particle_positions["y"]))
             v /= 1852 * 60
 
         if "3D" in vectorfield.vector_type:

src/parcels/interpolators/_xinterpolators.py

@@ -132,7 +132,7 @@ def interp(
         field: Field,
     ):
         """Returning the single value of a Constant Field (with a size=(1,1,1,1) array)"""
-        return field.data[0, 0, 0, 0].values * np.ones_like(particle_positions["lon"])
+        return field.data[0, 0, 0, 0].values * np.ones_like(particle_positions["x"])
 
 
 class XLinear_Velocity(VectorInterpolator):  # noqa:  N801
@@ -149,7 +149,7 @@ def interp(
         u = _xlinear.interp(particle_positions, grid_positions, vectorfield.U)
         v = _xlinear.interp(particle_positions, grid_positions, vectorfield.V)
         if vectorfield.grid._mesh == "spherical":
-            u /= 1852 * 60 * np.cos(np.deg2rad(particle_positions["lat"]))
+            u /= 1852 * 60 * np.cos(np.deg2rad(particle_positions["y"]))
             v /= 1852 * 60
 
         if vectorfield.W:
@@ -304,7 +304,7 @@ def _compute_corner_data(data, selection_dict) -> np.ndarray:
             v = v.compute()
 
         if grid._mesh == "spherical":
-            conversion = 1852 * 60.0 * np.cos(np.deg2rad(particle_positions["lat"]))
+            conversion = 1852 * 60.0 * np.cos(np.deg2rad(particle_positions["y"]))
             u /= conversion
             v /= conversion