Small refactors of VectorField and pragma comments

.gitignore

@@ -1,7 +1,6 @@
 build/*
 docs/_build/*
 docs/_downloads
-gh-pages/*
 lib/
 bin/
 parcels_examples
@@ -19,15 +18,13 @@ out-*
 .DS_store
 
 .vscode
-.idea/*
 .env
 Profile.prof
 **/.ipynb_checkpoints/*
 .cache/*
-.pytest_cache/*
 .eggs/*
 parcels.egg-info/*
 dist/parcels*.egg
-parcels/examples/particle*.png
 parcels/_version_setup.py
-/.pytest_cache/
+.pytest_cache
+.coverage

docs/examples/example_dask_chunk_OCMs.py

@@ -17,7 +17,7 @@ def compute_nemo_particle_advection(fieldset, mode):
     lonp = 2.5 * np.ones(npart)
     latp = [i for i in 52.0 + (-1e-3 + np.random.rand(npart) * 2.0 * 1e-3)]
 
-    def periodicBC(particle, fieldSet, time):
+    def periodicBC(particle, fieldSet, time):  # pragma: no cover
         if particle.lon > 15.0:
             particle_dlon -= 15.0  # noqa
         if particle.lon < 0:
@@ -652,7 +652,7 @@ class MyParticle(ptype[mode]):
 
     pset = parcels.ParticleSet(fieldset, pclass=MyParticle, lon=2.5, lat=52)
 
-    def Sample2D(particle, fieldset, time):
+    def Sample2D(particle, fieldset, time):  # pragma: no cover
         particle.sample_var_curvilinear += fieldset.nav_lon[
             time, particle.depth, particle.lat, particle.lon
         ]

docs/examples/example_globcurrent.py

@@ -129,7 +129,7 @@ def test_globcurrent_time_periodic(mode, rundays):
             fieldset, pclass=MyParticle, lon=25, lat=-35, time=fieldset.U.grid.time[0]
         )
 
-        def SampleU(particle, fieldset, time):
+        def SampleU(particle, fieldset, time):  # pragma: no cover
             u, v = fieldset.UV[time, particle.depth, particle.lat, particle.lon]
             particle.sample_var += u
 
@@ -255,7 +255,7 @@ def test_globcurrent_startparticles_between_time_arrays(mode, dt, with_starttime
 
     MyParticle = ptype[mode].add_variable("sample_var", initial=0.0)
 
-    def SampleP(particle, fieldset, time):
+    def SampleP(particle, fieldset, time):  # pragma: no cover
         particle.sample_var += fieldset.P[
             time, particle.depth, particle.lat, particle.lon
         ]
@@ -288,7 +288,7 @@ def test_globcurrent_particle_independence(mode, rundays=5):
     fieldset = set_globcurrent_fieldset()
     time0 = fieldset.U.grid.time[0]
 
-    def DeleteP0(particle, fieldset, time):
+    def DeleteP0(particle, fieldset, time):  # pragma: no cover
         if particle.id == 0:
             particle.delete()
 
@@ -353,7 +353,7 @@ def test_error_outputdt_not_multiple_dt(mode, tmpdir):
     pset = parcels.ParticleSet(fieldset, pclass=ptype[mode], lon=[0], lat=[0])
     ofile = pset.ParticleFile(name=filepath, outputdt=timedelta(days=1))
 
-    def DoNothing(particle, fieldset, time):
+    def DoNothing(particle, fieldset, time):  # pragma: no cover
         pass
 
     with pytest.raises(ValueError):

docs/examples/example_mitgcm.py

@@ -29,7 +29,7 @@ def run_mitgcm_zonally_reentrant(mode: Literal["scipy", "jit"], path: Path):
     fieldset.add_periodic_halo(zonal=True)
     fieldset.add_constant("domain_width", 1000000)
 
-    def periodicBC(particle, fieldset, time):
+    def periodicBC(particle, fieldset, time):  # pragma: no cover
         if particle.lon < 0:
             particle_dlon += fieldset.domain_width  # noqa
         elif particle.lon > fieldset.domain_width:

docs/examples/example_moving_eddies.py

@@ -269,7 +269,7 @@ def test_periodic_and_computeTimeChunk_eddies(mode):
         fieldset=fieldset, pclass=ptype[mode], lon=[3.3, 3.3], lat=[46.0, 47.8]
     )
 
-    def periodicBC(particle, fieldset, time):
+    def periodicBC(particle, fieldset, time):  # pragma: no cover
         if particle.lon < fieldset.halo_west:
             particle_dlon += fieldset.halo_east - fieldset.halo_west  # noqa
         elif particle.lon > fieldset.halo_east:
@@ -279,7 +279,7 @@ def periodicBC(particle, fieldset, time):
         elif particle.lat > fieldset.halo_north:
             particle_dlat -= fieldset.halo_north - fieldset.halo_south
 
-    def slowlySouthWestward(particle, fieldset, time):
+    def slowlySouthWestward(particle, fieldset, time):  # pragma: no cover
         particle_dlon -= 5 * particle.dt / 1e5  # noqa
         particle_dlat -= 3 * particle.dt / 1e5  # noqa
 

docs/examples/example_nemo_curvilinear.py

@@ -47,7 +47,7 @@ def run_nemo_curvilinear(mode, outfile, advtype="RK4"):
         latp = np.linspace(-70, 70, npart)
         runtime = timedelta(days=15)
 
-    def periodicBC(particle, fieldSet, time):
+    def periodicBC(particle, fieldSet, time):  # pragma: no cover
         if particle.lon > 180:
             particle_dlon -= 360  # noqa
 

docs/examples/example_peninsula.py

@@ -91,7 +91,7 @@ def peninsula_fieldset(xdim, ydim, mesh="flat", grid_type="A"):
     return fieldset
 
 
-def UpdateP(particle, fieldset, time):
+def UpdateP(particle, fieldset, time):  # pragma: no cover
     if time == 0:
         particle.p_start = fieldset.P[time, particle.depth, particle.lat, particle.lon]
     particle.p = fieldset.P[time, particle.depth, particle.lat, particle.lon]

docs/examples/example_stommel.py

@@ -77,13 +77,13 @@ def stommel_fieldset(xdim=200, ydim=200, grid_type="A"):
     return fieldset
 
 
-def UpdateP(particle, fieldset, time):
+def UpdateP(particle, fieldset, time):  # pragma: no cover
     if time == 0:
         particle.p_start = fieldset.P[time, particle.depth, particle.lat, particle.lon]
     particle.p = fieldset.P[time, particle.depth, particle.lat, particle.lon]
 
 
-def AgeP(particle, fieldset, time):
+def AgeP(particle, fieldset, time):  # pragma: no cover
     particle.age += particle.dt
     if particle.age > fieldset.maxage:
         particle.delete()

parcels/application_kernels/EOSseawaterproperties.py

@@ -5,7 +5,7 @@
 __all__ = ["AdiabticTemperatureGradient", "PressureFromLatDepth", "PtempFromTemp", "TempFromPtemp", "UNESCODensity"]
 
 
-def PressureFromLatDepth(particle, fieldset, time):
+def PressureFromLatDepth(particle, fieldset, time):  # pragma: no cover
     """
     Calculates pressure in dbars from depth in meters and latitude.
 
@@ -28,7 +28,7 @@
     particle.pressure = ((1 - C1) - math.pow(((math.pow((1 - C1), 2)) - (8.84e-6 * particle.depth)), 0.5)) / 4.42e-6
 
 
-def AdiabticTemperatureGradient(particle, fieldset, time):
+def AdiabticTemperatureGradient(particle, fieldset, time):  # pragma: no cover
     """Calculates adiabatic temperature gradient as per UNESCO 1983 routines.
 
 
@@ -81,7 +81,7 @@
     )
 
 
-def PtempFromTemp(particle, fieldset, time):
+def PtempFromTemp(particle, fieldset, time):  # pragma: no cover
     """
     Calculates potential temperature as per UNESCO 1983 report.
 
@@ -188,7 +188,7 @@
     particle.potemp = (th + (del_th - 2 * q) / 6) / 1.00024
 
 
-def TempFromPtemp(particle, fieldset, time):
+def TempFromPtemp(particle, fieldset, time):  # pragma: no cover
     """
     Calculates temperature from potential temperature at the reference
     pressure PR and in situ pressure P.
@@ -295,7 +295,7 @@
     particle.temp = (th + (del_th - 2 * q) / 6) / 1.00024
 
 
-def UNESCODensity(particle, fieldset, time):
+def UNESCODensity(particle, fieldset, time):  # pragma: no cover
     # This is a kernel which calculates the UNESCO density
     # (https://link.springer.com/content/pdf/bbm%3A978-3-319-18908-6%2F1.pdf),
     # from pressure, temperature and salinity.

parcels/application_kernels/TEOSseawaterdensity.py

@@ -5,7 +5,7 @@
 __all__ = ["PolyTEOS10_bsq"]
 
 
-def PolyTEOS10_bsq(particle, fieldset, time):
+def PolyTEOS10_bsq(particle, fieldset, time):  # pragma: no cover
     """Calculates density based on the polyTEOS10-bsq algorithm from Appendix A.2 of
     https://www.sciencedirect.com/science/article/pii/S1463500315000566
     requires fieldset.abs_salinity and fieldset.cons_temperature Fields in the fieldset

parcels/application_kernels/advection.py

@@ -14,7 +14,7 @@
 ]
 
 
-def AdvectionRK4(particle, fieldset, time):
+def AdvectionRK4(particle, fieldset, time):  # pragma: no cover
     """Advection of particles using fourth-order Runge-Kutta integration."""
     (u1, v1) = fieldset.UV[particle]
     lon1, lat1 = (particle.lon + u1 * 0.5 * particle.dt, particle.lat + v1 * 0.5 * particle.dt)
@@ -27,7 +27,7 @@
     particle_dlat += (v1 + 2 * v2 + 2 * v3 + v4) / 6.0 * particle.dt  # noqa
 
 
-def AdvectionRK4_3D(particle, fieldset, time):
+def AdvectionRK4_3D(particle, fieldset, time):  # pragma: no cover
     """Advection of particles using fourth-order Runge-Kutta integration including vertical velocity."""
     (u1, v1, w1) = fieldset.UVW[particle]
     lon1 = particle.lon + u1 * 0.5 * particle.dt
@@ -47,7 +47,7 @@
     particle_ddepth += (w1 + 2 * w2 + 2 * w3 + w4) / 6 * particle.dt  # noqa
 
 
-def AdvectionRK4_3D_CROCO(particle, fieldset, time):
+def AdvectionRK4_3D_CROCO(particle, fieldset, time):  # pragma: no cover
     """Advection of particles using fourth-order Runge-Kutta integration including vertical velocity.
     This kernel assumes the vertical velocity is the 'w' field from CROCO output and works on sigma-layers.
     """
@@ -92,14 +92,14 @@
     ) / 6
 
 
-def AdvectionEE(particle, fieldset, time):
+def AdvectionEE(particle, fieldset, time):  # pragma: no cover
     """Advection of particles using Explicit Euler (aka Euler Forward) integration."""
     (u1, v1) = fieldset.UV[particle]
     particle_dlon += u1 * particle.dt  # noqa
     particle_dlat += v1 * particle.dt  # noqa
 
 
-def AdvectionRK45(particle, fieldset, time):
+def AdvectionRK45(particle, fieldset, time):  # pragma: no cover
     """Advection of particles using adaptive Runge-Kutta 4/5 integration.
 
     Note that this kernel requires a Particle Class that has an extra Variable 'next_dt'
@@ -161,7 +161,7 @@
         return StatusCode.Repeat
 
 
-def AdvectionAnalytical(particle, fieldset, time):
+def AdvectionAnalytical(particle, fieldset, time):  # pragma: no cover
     """Advection of particles using 'analytical advection' integration.
 
     Based on Ariane/TRACMASS algorithm, as detailed in e.g. Doos et al (https://doi.org/10.5194/gmd-10-1733-2017).
@@ -232,14 +232,14 @@
     else:
         dz = 1.0
 
-    c1 = fieldset.UV.dist(py[0], py[1], px[0], px[1], grid.mesh, np.dot(i_u.phi2D_lin(0.0, xsi), py))
-    c2 = fieldset.UV.dist(py[1], py[2], px[1], px[2], grid.mesh, np.dot(i_u.phi2D_lin(eta, 1.0), py))
-    c3 = fieldset.UV.dist(py[2], py[3], px[2], px[3], grid.mesh, np.dot(i_u.phi2D_lin(1.0, xsi), py))
-    c4 = fieldset.UV.dist(py[3], py[0], px[3], px[0], grid.mesh, np.dot(i_u.phi2D_lin(eta, 0.0), py))
+    c1 = i_u.geodesic_distance(py[0], py[1], px[0], px[1], grid.mesh, np.dot(i_u.phi2D_lin(0.0, xsi), py))
+    c2 = i_u.geodesic_distance(py[1], py[2], px[1], px[2], grid.mesh, np.dot(i_u.phi2D_lin(eta, 1.0), py))
+    c3 = i_u.geodesic_distance(py[2], py[3], px[2], px[3], grid.mesh, np.dot(i_u.phi2D_lin(1.0, xsi), py))
+    c4 = i_u.geodesic_distance(py[3], py[0], px[3], px[0], grid.mesh, np.dot(i_u.phi2D_lin(eta, 0.0), py))
     rad = np.pi / 180.0
     deg2m = 1852 * 60.0
     meshJac = (deg2m * deg2m * math.cos(rad * particle.lat)) if grid.mesh == "spherical" else 1
-    dxdy = fieldset.UV.jacobian(py, px, eta, xsi) * meshJac
+    dxdy = i_u.compute_jacobian_determinant(py, px, eta, xsi) * meshJac
 
     if withW:
         U0 = direction * fieldset.U.data[ti, zi + 1, yi + 1, xi] * c4 * dz

parcels/application_kernels/advectiondiffusion.py

@@ -10,7 +10,7 @@
 __all__ = ["AdvectionDiffusionEM", "AdvectionDiffusionM1", "DiffusionUniformKh"]
 
 
-def AdvectionDiffusionM1(particle, fieldset, time):
+def AdvectionDiffusionM1(particle, fieldset, time):  # pragma: no cover
     """Kernel for 2D advection-diffusion, solved using the Milstein scheme at first order (M1).
 
     Assumes that fieldset has fields `Kh_zonal` and `Kh_meridional`
@@ -47,7 +47,7 @@
     particle_dlat += v * particle.dt + 0.5 * dKdy * (dWy**2 + particle.dt) + by * dWy  # noqa
 
 
-def AdvectionDiffusionEM(particle, fieldset, time):
+def AdvectionDiffusionEM(particle, fieldset, time):  # pragma: no cover
     """Kernel for 2D advection-diffusion, solved using the Euler-Maruyama scheme (EM).
 
     Assumes that fieldset has fields `Kh_zonal` and `Kh_meridional`
@@ -83,7 +83,7 @@
     particle_dlat += ay * particle.dt + by * dWy  # noqa
 
 
-def DiffusionUniformKh(particle, fieldset, time):
+def DiffusionUniformKh(particle, fieldset, time):  # pragma: no cover
     """Kernel for simple 2D diffusion where diffusivity (Kh) is assumed uniform.
 
     Assumes that fieldset has constant fields `Kh_zonal` and `Kh_meridional`.

parcels/field.py

@@ -1977,28 +1977,17 @@ def _check_grid_dimensions(grid1, grid2):
             and np.allclose(grid1.time_full, grid2.time_full)
         )
 
-    def dist(self, lat1: float, lat2: float, lon1: float, lon2: float, mesh: Mesh, lat: float):
-        if mesh == "spherical":
-            rad = np.pi / 180.0
-            deg2m = 1852 * 60.0
-            return np.sqrt(((lon2 - lon1) * deg2m * math.cos(rad * lat)) ** 2 + ((lat2 - lat1) * deg2m) ** 2)
-        else:
-            return np.sqrt((lon2 - lon1) ** 2 + (lat2 - lat1) ** 2)
-
-    def jacobian(self, py: np.ndarray, px: np.ndarray, eta: float, xsi: float):
-        dphidxsi = [eta - 1, 1 - eta, eta, -eta]
-        dphideta = [xsi - 1, -xsi, xsi, 1 - xsi]
+    @deprecated_made_private  # TODO: Remove 6 months after v3.2.0
+    def dist(self, *args, **kwargs):
+        raise NotImplementedError
 
-        dxdxsi = np.dot(px, dphidxsi)
-        dxdeta = np.dot(px, dphideta)
-        dydxsi = np.dot(py, dphidxsi)
-        dydeta = np.dot(py, dphideta)
-        jac = dxdxsi * dydeta - dxdeta * dydxsi
-        return jac
+    @deprecated_made_private  # TODO: Remove 6 months after v3.2.0
+    def jacobian(self, *args, **kwargs):
+        raise NotImplementedError
 
     def spatial_c_grid_interpolation2D(self, ti, z, y, x, time, particle=None, applyConversion=True):
         grid = self.U.grid
-        (zeta, eta, xsi, zi, yi, xi) = self.U._search_indices(time, z, y, x, ti, particle=particle)
+        (_, eta, xsi, zi, yi, xi) = self.U._search_indices(time, z, y, x, ti, particle=particle)
 
         if grid._gtype in [GridType.RectilinearSGrid, GridType.RectilinearZGrid]:
             px = np.array([grid.lon[xi], grid.lon[xi + 1], grid.lon[xi + 1], grid.lon[xi]])
@@ -2014,10 +2003,10 @@ def spatial_c_grid_interpolation2D(self, ti, z, y, x, time, particle=None, apply
             px[1:] = np.where(-px[1:] + px[0] > 180, px[1:] + 360, px[1:])
         xx = (1 - xsi) * (1 - eta) * px[0] + xsi * (1 - eta) * px[1] + xsi * eta * px[2] + (1 - xsi) * eta * px[3]
         assert abs(xx - x) < 1e-4
-        c1 = self.dist(py[0], py[1], px[0], px[1], grid.mesh, np.dot(i_u.phi2D_lin(0.0, xsi), py))
-        c2 = self.dist(py[1], py[2], px[1], px[2], grid.mesh, np.dot(i_u.phi2D_lin(eta, 1.0), py))
-        c3 = self.dist(py[2], py[3], px[2], px[3], grid.mesh, np.dot(i_u.phi2D_lin(1.0, xsi), py))
-        c4 = self.dist(py[3], py[0], px[3], px[0], grid.mesh, np.dot(i_u.phi2D_lin(eta, 0.0), py))
+        c1 = i_u.geodesic_distance(py[0], py[1], px[0], px[1], grid.mesh, np.dot(i_u.phi2D_lin(0.0, xsi), py))
+        c2 = i_u.geodesic_distance(py[1], py[2], px[1], px[2], grid.mesh, np.dot(i_u.phi2D_lin(eta, 1.0), py))
+        c3 = i_u.geodesic_distance(py[2], py[3], px[2], px[3], grid.mesh, np.dot(i_u.phi2D_lin(1.0, xsi), py))
+        c4 = i_u.geodesic_distance(py[3], py[0], px[3], px[0], grid.mesh, np.dot(i_u.phi2D_lin(eta, 0.0), py))
         if grid.zdim == 1:
             if self.gridindexingtype == "nemo":
                 U0 = self.U.data[ti, yi + 1, xi] * c4
@@ -2049,7 +2038,7 @@ def spatial_c_grid_interpolation2D(self, ti, z, y, x, time, particle=None, apply
         else:
             meshJac = deg2m if grid.mesh == "spherical" else 1
 
-        jac = self.jacobian(py, px, eta, xsi) * meshJac
+        jac = i_u.compute_jacobian_determinant(py, px, eta, xsi) * meshJac
 
         u = (
             (-(1 - eta) * U - (1 - xsi) * V) * px[0]

parcels/kernel.py

@@ -300,7 +300,7 @@
 
     def add_scipy_positionupdate_kernels(self):
         # Adding kernels that set and update the coordinate changes
-        def Setcoords(particle, fieldset, time):
+        def Setcoords(particle, fieldset, time):  # pragma: no cover
             particle_dlon = 0  # noqa
             particle_dlat = 0  # noqa
             particle_ddepth = 0  # noqa
@@ -309,7 +309,7 @@
             particle.depth = particle.depth_nextloop
             particle.time = particle.time_nextloop
 
-        def Updatecoords(particle, fieldset, time):
+        def Updatecoords(particle, fieldset, time):  # pragma: no cover
             particle.lon_nextloop = particle.lon + particle_dlon  # type: ignore[name-defined] # noqa
             particle.lat_nextloop = particle.lat + particle_dlat  # type: ignore[name-defined] # noqa
             particle.depth_nextloop = particle.depth + particle_ddepth  # type: ignore[name-defined] # noqa