diff --git a/docs/user_guide/examples/explanation_kernelloop.md b/docs/user_guide/examples/explanation_kernelloop.md index c326d9906..143c4aab8 100644 --- a/docs/user_guide/examples/explanation_kernelloop.md +++ b/docs/user_guide/examples/explanation_kernelloop.md @@ -43,6 +43,7 @@ Besides having commutable Kernels, the main advantage of this implementation is Below is a simple example of some particles at the surface of the ocean. We create an idealised zonal wind flow that will "push" a particle that is already affected by the surface currents. The Kernel loop ensures that these two forces act at the same time and location. ```{code-cell} +:tags: [hide-output] import matplotlib.pyplot as plt import numpy as np import xarray as xr @@ -69,21 +70,18 @@ ds_fields["VWind"] = xr.DataArray( fieldset = parcels.FieldSet.from_copernicusmarine(ds_fields) -# Set unit converters for custom wind fields -fieldset.UWind.units = parcels.GeographicPolar() -fieldset.VWind.units = parcels.Geographic() +# Create a vecorfield for the wind +windvector = parcels.VectorField("Wind", fieldset.UWind, fieldset.VWind) +fieldset.add_field(windvector) ``` Now we define a wind kernel that uses a forward Euler method to apply the wind forcing. Note that we update the `particles.dlon` and `particles.dlat` variables, rather than `particles.lon` and `particles.lat` directly. ```{code-cell} def wind_kernel(particles, fieldset): - particles.dlon += ( - fieldset.UWind[particles] * particles.dt - ) - particles.dlat += ( - fieldset.VWind[particles] * particles.dt - ) + uwind, vwind = fieldset.Wind[particles] + particles.dlon += uwind * particles.dt + particles.dlat += vwind * particles.dt ``` First run a simulation where we apply kernels as `[AdvectionRK4, wind_kernel]` diff --git a/docs/user_guide/examples/tutorial_unitconverters.ipynb b/docs/user_guide/examples/tutorial_unitconverters.ipynb index a844833ea..a780c8550 100644 --- a/docs/user_guide/examples/tutorial_unitconverters.ipynb +++ b/docs/user_guide/examples/tutorial_unitconverters.ipynb @@ -5,7 +5,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "# 🖥️ Spherical grids and unit converters" + "# 🖥️ Spherical grids and velocity conversion" ] }, { @@ -23,7 +23,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Let's first import the relevant modules, and generate a simple dataset on a 2D spherical mesh, with `U`, `V` and `temperature` data arrays, with the velocities 1 m/s and the temperature 20C.\n" + "Let's first import the relevant modules, and generate a simple dataset on a 2D spherical mesh, with `U`, `V` and `temperature` data arrays, with the velocities 1 m s-1 and the temperature 20°C.\n" ] }, { @@ -67,7 +67,7 @@ "When using a `FieldSet` method for a specific dataset, such as `from_copernicusmarine()`, the grid information is known and parsed by Parcels, so we do not have to add the `mesh` argument.\n", "```\n", "\n", - "Plotting the `U` field indeed shows a uniform 1 m/s eastward flow.\n" + "Plotting the `U` field indeed shows a uniform 1 m s-1 eastward flow.\n" ] }, { @@ -104,7 +104,11 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "However, printing the velocites directly shows something perhaps surprising. Here, we use the square-bracket field-interpolation notation to print the field value at (5W, 40N, 0m depth) at time 0. _Note that sampling a velocity in Parcels is done by calling the `fieldset.UV` VectorField; see the [Field Sampling tutorial](https://docs.oceanparcels.org/en/latest/examples/tutorial_sampling.html#Sampling-velocity-fields) for more information._\n" + "However, printing the velocites directly shows something perhaps surprising. Here, we use the square-bracket field-interpolation notation to print the field value at (5W, 40N, 0m depth) at time 0. \n", + "\n", + "```{note}\n", + "Sampling a velocity in Parcels is done by calling the `fieldset.UV` VectorField; see also the section \"Sampling U and V separately\" below.\n", + "```\n" ] }, { @@ -117,8 +121,11 @@ "z = np.array([0])\n", "lat = np.array([40])\n", "lon = np.array([-5])\n", - "print(fieldset.UV[time, z, lat, lon])\n", - "print(fieldset.temperature[time, z, lat, lon])" + "u, v = fieldset.UV[time, z, lat, lon]\n", + "temp = fieldset.temperature[time, z, lat, lon]\n", + "\n", + "print(f\"(u, v) = ({u}, {v})\")\n", + "print(f\"temperature = {temp}\")" ] }, { @@ -126,19 +133,9 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "While the temperature field indeed is 20C, as we defined, these printed velocities are much smaller.\n", + "While the temperature field indeed is 20°C, as we defined, these printed velocities are much smaller.\n", "\n", - "This is because Parcels converts under the hood from m/s to degrees/s. This conversion is done with a `parcels.converters` object, which is stored in the `.units` attribute of each Field. Below, we print these\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "for fld in [fieldset.U, fieldset.V, fieldset.temperature]:\n", - " print(f\"{fld.name}: {fld.units}\")" + "This is because Parcels converts under the hood from m s-1 to degrees s-1." ] }, { @@ -146,9 +143,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "So the U field has a `GeographicPolar` UnitConverter object, the V field has a `Geographic` UnitConverter and the `temp` field has a `Unity` object.\n", - "\n", - "Indeed, if we multiply the value of the V field with 1852 \\* 60 (the number of meters in 1 degree of latitude), we get the expected 1 m/s.\n" + "Indeed, if we multiply the value of the U field with 1852 \\* 60 \\* cos(lat) (the number of meters in 1 degree of longitude), we get the expected 1 m s-1.\n" ] }, { @@ -158,7 +153,22 @@ "outputs": [], "source": [ "u, v = fieldset.UV[time, z, lat, lon]\n", - "print(v * 1852 * 60)" + "\n", + "u = u * 1852 * 60 * np.cos(np.deg2rad(lat))\n", + "v = v * 1852 * 60 * np.cos(np.deg2rad(lat))\n", + "print(f\"(u, v) = ({u}, {v})\")\n", + "\n", + "assert np.isclose(u, 1.0)\n", + "assert np.isclose(v, 1.0)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "```{note}\n", + "It may be surprising that the conversion factor depends on latitude also for the V velocity component. This is because Parcels assumes a flux-based interpolation, so that the velocity components are defined on the faces of the grid cells, which vary in size with latitude.\n", + "```" ] }, { @@ -166,7 +176,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Note that you can also interpolate the Field without a unit conversion, by using the `eval()` method and setting `applyConversion=False`, as below\n" + "You can also interpolate the Field to these values directly by using the `eval()` method and setting `applyConversion=False`, as below\n" ] }, { @@ -187,134 +197,17 @@ ] }, { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## UnitConverters for `mesh='flat'`\n" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "If longitudes and latitudes are given in meters, rather than degrees, simply add `mesh='flat'` when creating the XGrid object.\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "ds_flat = simple_UV_dataset(dims=(1, 1, nlat, nlon), mesh=\"flat\").isel(time=0, depth=0)\n", - "ds_flat[\"temperature\"] = ds_flat[\"U\"] + 20 # add temperature field of 20 deg\n", - "ds_flat[\"U\"].data[:] = 1.0 # set U to 1 m/s\n", - "ds_flat[\"V\"].data[:] = 1.0 # set V to 1 m/s\n", - "grid = parcels.XGrid.from_dataset(ds_flat, mesh=\"flat\")\n", - "U = parcels.Field(\"U\", ds_flat[\"U\"], grid, interp_method=parcels.interpolators.XLinear)\n", - "V = parcels.Field(\"V\", ds_flat[\"V\"], grid, interp_method=parcels.interpolators.XLinear)\n", - "UV = parcels.VectorField(\"UV\", U, V)\n", - "temperature = parcels.Field(\n", - " \"temperature\",\n", - " ds_flat[\"temperature\"],\n", - " grid,\n", - " interp_method=parcels.interpolators.XLinear,\n", - ")\n", - "fieldset_flat = parcels.FieldSet([U, V, UV, temperature])\n", - "\n", - "plt.pcolormesh(\n", - " fieldset_flat.U.grid.lon,\n", - " fieldset_flat.U.grid.lat,\n", - " fieldset_flat.U.data[0, 0, :, :],\n", - " vmin=0,\n", - " vmax=1,\n", - " shading=\"gouraud\",\n", - ")\n", - "plt.colorbar()\n", - "plt.show()\n", - "\n", - "print(\n", - " \"Velocities:\",\n", - " fieldset_flat.UV[time, z, lat, lon],\n", - ")\n", - "for fld in [fieldset_flat.U, fieldset_flat.V, fieldset_flat.temperature]:\n", - " print(f\"{fld.name}: {fld.units}\")" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Indeed, in this case all Fields have the same default `Unity` object.\n" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## UnitConverters for Diffusion fields\n" - ] - }, - { - "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ - "The units for Brownian diffusion are in $m^2/s$. If (and only if!) the diffusion fields are called \"Kh_zonal\" and \"Kh_meridional\", Parcels will automatically assign the correct Unitconverter objects to these fields.\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "kh_zonal = 100 # in m^2/s\n", - "kh_meridional = 100 # in m^2/s\n", - "\n", - "ds[\"Kh_zonal\"] = xr.DataArray(\n", - " data=kh_zonal * np.ones((nlat, nlon), dtype=np.float32), dims=[\"YG\", \"XG\"]\n", - ")\n", - "\n", - "kh_zonal_field = parcels.Field(\n", - " \"Kh_zonal\",\n", - " ds[\"Kh_zonal\"],\n", - " grid=fieldset.U.grid,\n", - " interp_method=parcels.interpolators.XLinear,\n", - ")\n", - "\n", - "fieldset.add_field(kh_zonal_field)\n", - "\n", - "ds[\"Kh_meridional\"] = xr.DataArray(\n", - " data=kh_meridional * np.ones((nlat, nlon), dtype=np.float32), dims=[\"YG\", \"XG\"]\n", - ")\n", - "\n", - "kh_meridional_field = parcels.Field(\n", - " \"Kh_meridional\",\n", - " ds[\"Kh_meridional\"],\n", - " grid=fieldset.U.grid,\n", - " interp_method=parcels.interpolators.XLinear,\n", - ")\n", - "\n", - "fieldset.add_field(kh_meridional_field)\n", - "\n", - "for fld in [fieldset.Kh_zonal, fieldset.Kh_meridional]:\n", - " val = fld[time, z, lat, lon]\n", - " print(f\"{fld.name}: {val} {fld.units}\")" + "## Don't sample U and V separately" ] }, { - "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ - "Here, the unitconverters are `GeographicPolarSquare` and `GeographicSquare`, respectively.\n", - "\n", - "Indeed, multiplying with $(1852\\cdot60)^2$ returns the original value\n" + "Sampling `U` and `V` separately will _not_ convert to degrees s-1, so these velocities cannot be used directly for advection on spherical coordinates. This is one of the main reasons to always use the `UV` VectorField for velocity sampling in Parcels.\n" ] }, { @@ -323,76 +216,21 @@ "metadata": {}, "outputs": [], "source": [ - "deg_to_m = 1852 * 60\n", - "print(fieldset.Kh_meridional[time, z, lat, lon] * deg_to_m**2)" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Adding a UnitConverter object to a Field\n" + "for fld in [fieldset.U, fieldset.V]:\n", + " print(f\"{fld.name}: {fld.eval(time, z, lat, lon)}\")" ] }, { - "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ - "So, to summarise, here is a table with all the conversions\n", - "\n", - "| Field name | Converter object | Conversion for `mesh='spherical'` | Conversion for `mesh='flat'` |\n", - "| ---------------- | ----------------------- | --------------------------------------------------------- | ---------------------------- |\n", - "| `\"U\"` | `GeographicPolar` | $1852 \\cdot 60 \\cdot \\cos(lat \\cdot \\frac{\\pi}{180})$ | 1 |\n", - "| `\"V\"` | `Geographic` | $1852 \\cdot 60$ | 1 |\n", - "| `\"Kh_zonal\"` | `GeographicPolarSquare` | $(1852 \\cdot 60 \\cdot \\cos(lat \\cdot \\frac{\\pi}{180}))^2$ | 1 |\n", - "| `\"Kh_meridional\"` | `GeographicSquare` | $(1852 \\cdot 60)^2$ | 1 |\n", - "| All other fields | `Unity` | 1 | 1 |\n", + "## Unit conversion for other fields such as diffusivity\n", "\n", - "Only four Field names are recognised and assigned an automatic UnitConverter object. This means that things might go very wrong when e.g. a velocity field is not called `U` or `V`.\n", + "For other fields such as diffusivity, Parcels does not apply similar unit conversions when using a `spherical` mesh. For example, if we define a diffusivity field with value 10 m2 s-1, Parcels will not convert this to degrees2 s-1 under the hood. \n", "\n", - "Fortunately, you can always add a UnitConverter later, as explained below:\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "ds[\"Ustokes\"] = xr.DataArray(\n", - " data=np.ones((nlat, nlon), dtype=np.float32), dims=[\"YG\", \"XG\"]\n", - ")\n", + "If you want to work with diffusivity in degrees2 s-1 (for example to move particles using a random walk), you will have to convert this yourself in your kernel. \n", "\n", - "fieldset.add_field(\n", - " parcels.Field(\n", - " \"Ustokes\",\n", - " ds[\"Ustokes\"],\n", - " grid=fieldset.U.grid,\n", - " interp_method=parcels.interpolators.XLinear,\n", - " )\n", - ")\n", - "print(fieldset.Ustokes[time, z, lat, lon])" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "This value for `Ustokes` of course is not as expected, since the mesh is spherical and hence this would mean 1 degree/s velocity. Assigning the correct `GeographicPolar` Unitconverter gives\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "fieldset.Ustokes.units = parcels.GeographicPolar()\n", - "print(fieldset.Ustokes[time, z, lat, lon])\n", - "print(fieldset.Ustokes[time, z, lat, lon] * 1852 * 60 * np.cos(40 * np.pi / 180))" + "Note that for the built-in `DiffusionUniformKh`, `AdvectionDiffusionM1` and `AdvectionDiffusionEM`, the conversion is done automatically." ] } ], diff --git a/docs/user_guide/examples_v3/example_moving_eddies.py b/docs/user_guide/examples_v3/example_moving_eddies.py index 108d6b079..348e21b32 100644 --- a/docs/user_guide/examples_v3/example_moving_eddies.py +++ b/docs/user_guide/examples_v3/example_moving_eddies.py @@ -25,8 +25,7 @@ def moving_eddies_fieldset(xdim=200, ydim=350, mesh="flat"): xdim : Vertical dimension of the generated fieldset (Default value = 200) mesh : str - String indicating the type of mesh coordinates and - units used during velocity interpolation: + String indicating the type of mesh coordinates used during velocity interpolation: 1. spherical: Lat and lon in degree, with a correction for zonal velocity U near the poles. diff --git a/docs/user_guide/examples_v3/example_peninsula.py b/docs/user_guide/examples_v3/example_peninsula.py index 718a53598..cccfdd7ce 100644 --- a/docs/user_guide/examples_v3/example_peninsula.py +++ b/docs/user_guide/examples_v3/example_peninsula.py @@ -24,8 +24,7 @@ def peninsula_fieldset(xdim, ydim, mesh="flat", grid_type="A"): xdim : Vertical dimension of the generated fieldset mesh : str - String indicating the type of mesh coordinates and - units used during velocity interpolation: + String indicating the type of mesh coordinates used during velocity interpolation: 1. spherical: Lat and lon in degree, with a correction for zonal velocity U near the poles. diff --git a/docs/user_guide/v4-migration.md b/docs/user_guide/v4-migration.md index a5685e58e..d548949b6 100644 --- a/docs/user_guide/v4-migration.md +++ b/docs/user_guide/v4-migration.md @@ -39,6 +39,7 @@ Version 4 of Parcels is unreleased at the moment. The information in this migrat - `Field.eval()` returns an array of floats instead of a single float (related to the vectorization) - `Field.eval()` does not throw OutOfBounds or other errors +- `applyConversion` has been renamed to `apply_conversion` and only works for VectorFields. Conversion of units should be handled in Kernels. ## GridSet @@ -46,4 +47,4 @@ Version 4 of Parcels is unreleased at the moment. The information in this migrat ## UnitConverters -- The default `UnitConverter` is now called `Unity()` +- UnitConverters have been removed. Instead, Interpolation functions should handle unit conversion internally, based on the value of `grid._mesh` ("spherical" or "flat"). diff --git a/src/parcels/__init__.py b/src/parcels/__init__.py index 177eceb7b..4647e98dd 100644 --- a/src/parcels/__init__.py +++ b/src/parcels/__init__.py @@ -24,14 +24,6 @@ from parcels._core.uxgrid import UxGrid from parcels._core.xgrid import XGrid -from parcels._core.converters import ( - Geographic, - GeographicPolar, - GeographicPolarSquare, - GeographicSquare, - Unity, -) - from parcels._core.statuscodes import ( AllParcelsErrorCodes, FieldInterpolationError, @@ -64,12 +56,6 @@ "BaseGrid", "UxGrid", "XGrid", - # Converters - "Geographic", - "GeographicPolar", - "GeographicPolarSquare", - "GeographicSquare", - "Unity", # Status codes and errors "AllParcelsErrorCodes", "FieldInterpolationError", diff --git a/src/parcels/_core/converters.py b/src/parcels/_core/converters.py deleted file mode 100644 index a3ecd965e..000000000 --- a/src/parcels/_core/converters.py +++ /dev/null @@ -1,117 +0,0 @@ -from __future__ import annotations - -from abc import ABC, abstractmethod -from math import pi - -import numpy as np -import numpy.typing as npt - -__all__ = [ - "Geographic", - "GeographicPolar", - "GeographicPolarSquare", - "GeographicSquare", - "UnitConverter", - "Unity", - "_convert_to_flat_array", - "_unitconverters_map", -] - - -def _convert_to_flat_array(var: npt.ArrayLike) -> npt.NDArray: - """Convert lists and single integers/floats to one-dimensional numpy arrays - - Parameters - ---------- - var : Array - list or numeric to convert to a one-dimensional numpy array - """ - return np.array(var).flatten() - - -class UnitConverter(ABC): - source_unit: str | None = None - target_unit: str | None = None - - @abstractmethod - def to_target(self, value, z, y, x): ... - - @abstractmethod - def to_source(self, value, z, y, x): ... - - -class Unity(UnitConverter): - """Interface class for spatial unit conversion during field sampling that performs no conversion.""" - - source_unit: None - target_unit: None - - def to_target(self, value, z, y, x): - return value - - def to_source(self, value, z, y, x): - return value - - -class Geographic(UnitConverter): - """Unit converter from geometric to geographic coordinates (m to degree)""" - - source_unit = "m" - target_unit = "degree" - - def to_target(self, value, z, y, x): - return value / 1000.0 / 1.852 / 60.0 - - def to_source(self, value, z, y, x): - return value * 1000.0 * 1.852 * 60.0 - - -class GeographicPolar(UnitConverter): - """Unit converter from geometric to geographic coordinates (m to degree) - with a correction to account for narrower grid cells closer to the poles. - """ - - source_unit = "m" - target_unit = "degree" - - def to_target(self, value, z, y, x): - return value / 1000.0 / 1.852 / 60.0 / np.cos(y * pi / 180) - - def to_source(self, value, z, y, x): - return value * 1000.0 * 1.852 * 60.0 * np.cos(y * pi / 180) - - -class GeographicSquare(UnitConverter): - """Square distance converter from geometric to geographic coordinates (m2 to degree2)""" - - source_unit = "m2" - target_unit = "degree2" - - def to_target(self, value, z, y, x): - return value / pow(1000.0 * 1.852 * 60.0, 2) - - def to_source(self, value, z, y, x): - return value * pow(1000.0 * 1.852 * 60.0, 2) - - -class GeographicPolarSquare(UnitConverter): - """Square distance converter from geometric to geographic coordinates (m2 to degree2) - with a correction to account for narrower grid cells closer to the poles. - """ - - source_unit = "m2" - target_unit = "degree2" - - def to_target(self, value, z, y, x): - return value / pow(1000.0 * 1.852 * 60.0 * np.cos(y * pi / 180), 2) - - def to_source(self, value, z, y, x): - return value * pow(1000.0 * 1.852 * 60.0 * np.cos(y * pi / 180), 2) - - -_unitconverters_map = { - "U": GeographicPolar(), - "V": Geographic(), - "Kh_zonal": GeographicPolarSquare(), - "Kh_meridional": GeographicSquare(), -} diff --git a/src/parcels/_core/field.py b/src/parcels/_core/field.py index a510209b3..94c011766 100644 --- a/src/parcels/_core/field.py +++ b/src/parcels/_core/field.py @@ -8,11 +8,6 @@ import uxarray as ux import xarray as xr -from parcels._core.converters import ( - UnitConverter, - Unity, - _unitconverters_map, -) from parcels._core.index_search import GRID_SEARCH_ERROR, LEFT_OUT_OF_BOUNDS, RIGHT_OUT_OF_BOUNDS, _search_time_index from parcels._core.particle import KernelParticle from parcels._core.statuscodes import ( @@ -135,25 +130,10 @@ def __init__( self.igrid = -1 # Default the grid index to -1 - if self.grid._mesh == "flat" or (self.name not in _unitconverters_map.keys()): - self.units = Unity() - elif self.grid._mesh == "spherical": - self.units = _unitconverters_map[self.name] - if self.data.shape[0] > 1: if "time" not in self.data.coords: raise ValueError("Field data is missing a 'time' coordinate.") - @property - def units(self): - return self._units - - @units.setter - def units(self, value): - if not isinstance(value, UnitConverter): - raise ValueError(f"Units must be a UnitConverter object, got {type(value)}") - self._units = value - @property def xdim(self): if type(self.data) is xr.DataArray: @@ -197,12 +177,30 @@ def _check_velocitysampling(self): stacklevel=2, ) - def eval(self, time: datetime, z, y, x, particles=None, applyConversion=True): + def eval(self, time: datetime, z, y, x, particles=None): """Interpolate field values in space and time. - We interpolate linearly in time and apply implicit unit - conversion to the result. Note that we defer to - scipy.interpolate to perform spatial interpolation. + Parameters + ---------- + time : float or array-like + Time(s) at which to sample the field. + z, y, x : scalar or array-like + Vertical (z), latitudinal (y) and longitudinal (x) positions to sample. + Inputs are promoted to 1-D arrays internally. + particles : ParticleSet, optional + If provided, used to associate results with particle indices and to + update particle state and element indices. Defaults to None. + + Returns + ------- + (value) : float or array-like + The interpolated value as a numpy.ndarray (or scalar) with the same + broadcasted shape as the input coordinates. + + Notes + ----- + - Particle states are updated for out-of-bounds, search errors and NaN + interpolation values. """ if particles is None: _ei = None @@ -218,8 +216,6 @@ def eval(self, time: datetime, z, y, x, particles=None, applyConversion=True): _update_particle_states_interp_value(particles, value) - if applyConversion: - value = self.units.to_target(value, z, y, x) return value def __getitem__(self, key): @@ -284,11 +280,34 @@ def vector_interp_method(self, method: Callable): self._vector_interp_method = method def eval(self, time: datetime, z, y, x, particles=None, applyConversion=True): - """Interpolate field values in space and time. - - We interpolate linearly in time and apply implicit unit - conversion to the result. Note that we defer to - scipy.interpolate to perform spatial interpolation. + """Interpolate vectorfield values in space and time. + + Parameters + ---------- + time : float or array-like + Time(s) at which to sample the field. + z, y, x : scalar or array-like + Vertical (z), latitudinal (y) and longitudinal (x) positions to sample. + Inputs are promoted to 1-D arrays internally. + particles : ParticleSet, optional + If provided, used to associate results with particle indices and to + update particle state and element indices. Defaults to None. + applyConversion : bool, default True + If True and the underlying grid is spherical, the horizontal velocity + components (U, V) are converted from m s^-1 to degrees s^-1. + If False, raw interpolated values are returned. + + Returns + ------- + (u, v, (w,)) : tuple or array-like + The interpolated velocity components: (u, v) for 2D vectors or (u, v, w) + for 3D vectors. Each element is a numpy.ndarray (or scalar) with the same + broadcasted shape as the input coordinates. + + Notes + ----- + - Particle states are updated for out-of-bounds, search errors and NaN + interpolation values. """ if particles is None: _ei = None @@ -315,8 +334,6 @@ def eval(self, time: datetime, z, y, x, particles=None, applyConversion=True): meshJac = 1852 * 60.0 * np.cos(np.deg2rad(y)) u = u / meshJac v = v / meshJac - if "3D" in self.vector_type: - w = self.W.units.to_target(w, z, y, x) for vel in (u, v, w): _update_particle_states_interp_value(particles, vel) diff --git a/src/parcels/_core/fieldset.py b/src/parcels/_core/fieldset.py index a5b8a1293..4e4a1a6bb 100644 --- a/src/parcels/_core/fieldset.py +++ b/src/parcels/_core/fieldset.py @@ -105,14 +105,6 @@ def add_field(self, field: Field, name: str | None = None): name : str Name of the :class:`parcels.field.Field` object to be added. Defaults to name in Field object. - - - Examples - -------- - For usage examples see the following tutorials: - - * `Unit converters <../examples/tutorial_unitconverters.ipynb>`__ (Default value = None) - """ if not isinstance(field, (Field, VectorField)): raise ValueError(f"Expected `field` to be a Field or VectorField object. Got {type(field)}") @@ -136,8 +128,8 @@ def add_constant_field(self, name: str, value, mesh: Mesh = "spherical"): value : Value of the constant field mesh : str - String indicating the type of mesh coordinates and - units used during velocity interpolation, see also `this tutorial <../examples/tutorial_unitconverters.ipynb>`__: + String indicating the type of mesh coordinates, + see also `this tutorial <../examples/tutorial_unitconverters.ipynb>`__: 1. spherical (default): Lat and lon in degree, with a correction for zonal velocity U near the poles. @@ -366,7 +358,7 @@ def from_sgrid_conventions( ds : xarray.Dataset xarray.Dataset with SGRID convention metadata. mesh : str - String indicating the type of mesh coordinates and units used during + String indicating the type of mesh coordinates used during velocity interpolation. Options are "spherical" or "flat". Returns diff --git a/src/parcels/_core/particleset.py b/src/parcels/_core/particleset.py index 4895725e5..deb224b76 100644 --- a/src/parcels/_core/particleset.py +++ b/src/parcels/_core/particleset.py @@ -9,7 +9,6 @@ from tqdm import tqdm from zarr.storage import DirectoryStore -from parcels._core.converters import _convert_to_flat_array from parcels._core.kernel import Kernel from parcels._core.particle import KernelParticle, Particle, create_particle_data from parcels._core.statuscodes import StatusCode @@ -73,9 +72,9 @@ def __init__( self._kernel = None self.fieldset = fieldset - lon = np.empty(shape=0) if lon is None else _convert_to_flat_array(lon) - lat = np.empty(shape=0) if lat is None else _convert_to_flat_array(lat) - time = np.empty(shape=0) if time is None else _convert_to_flat_array(time) + lon = np.empty(shape=0) if lon is None else np.array(lon).flatten() + lat = np.empty(shape=0) if lat is None else np.array(lat).flatten() + time = np.empty(shape=0) if time is None else np.array(time).flatten() if trajectory_ids is None: trajectory_ids = np.arange(lon.size) @@ -87,7 +86,7 @@ def __init__( minz = min(minz, field.grid.depth[0]) z = np.ones(lon.size) * minz else: - z = _convert_to_flat_array(z) + 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" if time is None or len(time) == 0: @@ -108,7 +107,7 @@ def __init__( for kwvar in kwargs: if kwvar not in ["partition_function"]: - kwargs[kwvar] = _convert_to_flat_array(kwargs[kwvar]) + 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." ) diff --git a/src/parcels/kernels/advectiondiffusion.py b/src/parcels/kernels/advectiondiffusion.py index f0466dd1a..51f6d2666 100644 --- a/src/parcels/kernels/advectiondiffusion.py +++ b/src/parcels/kernels/advectiondiffusion.py @@ -8,6 +8,16 @@ __all__ = ["AdvectionDiffusionEM", "AdvectionDiffusionM1", "DiffusionUniformKh"] +def meters_to_degrees_zonal(deg, lat): # pragma: no cover + """Convert square meters to square degrees longitude at a given latitude.""" + return deg / pow(1852 * 60.0 * np.cos(lat * np.pi / 180), 2) + + +def meters_to_degrees_meridional(deg): # pragma: no cover + """Convert square meters to square degrees latitude.""" + return deg / pow(1852 * 60.0, 2) + + def AdvectionDiffusionM1(particles, fieldset): # pragma: no cover """Kernel for 2D advection-diffusion, solved using the Milstein scheme at first order (M1). @@ -29,16 +39,28 @@ def AdvectionDiffusionM1(particles, fieldset): # pragma: no cover Kxp1 = fieldset.Kh_zonal[particles.time, particles.z, particles.lat, particles.lon + fieldset.dres, particles] Kxm1 = fieldset.Kh_zonal[particles.time, particles.z, particles.lat, particles.lon - fieldset.dres, particles] + if fieldset.Kh_zonal.grid._mesh == "spherical": + Kxp1 = meters_to_degrees_zonal(Kxp1, particles.lat) + Kxm1 = meters_to_degrees_zonal(Kxm1, particles.lat) dKdx = (Kxp1 - Kxm1) / (2 * fieldset.dres) u, v = fieldset.UV[particles.time, particles.z, particles.lat, particles.lon, particles] - bx = np.sqrt(2 * fieldset.Kh_zonal[particles.time, particles.z, particles.lat, particles.lon, particles]) + kh_zonal = fieldset.Kh_zonal[particles.time, particles.z, particles.lat, particles.lon, particles] + if fieldset.Kh_zonal.grid._mesh == "spherical": + kh_zonal = meters_to_degrees_zonal(kh_zonal, particles.lat) + bx = np.sqrt(2 * kh_zonal) Kyp1 = fieldset.Kh_meridional[particles.time, particles.z, particles.lat + fieldset.dres, particles.lon, particles] Kym1 = fieldset.Kh_meridional[particles.time, particles.z, particles.lat - fieldset.dres, particles.lon, particles] + if fieldset.Kh_meridional.grid._mesh == "spherical": + Kyp1 = meters_to_degrees_meridional(Kyp1) + Kym1 = meters_to_degrees_meridional(Kym1) dKdy = (Kyp1 - Kym1) / (2 * fieldset.dres) - by = np.sqrt(2 * fieldset.Kh_meridional[particles.time, particles.z, particles.lat, particles.lon, particles]) + kh_meridional = fieldset.Kh_meridional[particles.time, particles.z, particles.lat, particles.lon, particles] + if fieldset.Kh_meridional.grid._mesh == "spherical": + kh_meridional = meters_to_degrees_meridional(kh_meridional) + by = np.sqrt(2 * kh_meridional) # Particle positions are updated only after evaluating all terms. particles.dlon += u * particles.dt + 0.5 * dKdx * (dWx**2 + particles.dt) + bx * dWx @@ -66,15 +88,29 @@ def AdvectionDiffusionEM(particles, fieldset): # pragma: no cover Kxp1 = fieldset.Kh_zonal[particles.time, particles.z, particles.lat, particles.lon + fieldset.dres, particles] Kxm1 = fieldset.Kh_zonal[particles.time, particles.z, particles.lat, particles.lon - fieldset.dres, particles] + if fieldset.Kh_zonal.grid._mesh == "spherical": + Kxp1 = meters_to_degrees_zonal(Kxp1, particles.lat) + Kxm1 = meters_to_degrees_zonal(Kxm1, particles.lat) dKdx = (Kxp1 - Kxm1) / (2 * fieldset.dres) ax = u + dKdx - bx = np.sqrt(2 * fieldset.Kh_zonal[particles.time, particles.z, particles.lat, particles.lon, particles]) + + kh_zonal = fieldset.Kh_zonal[particles.time, particles.z, particles.lat, particles.lon, particles] + if fieldset.Kh_zonal.grid._mesh == "spherical": + kh_zonal = meters_to_degrees_zonal(kh_zonal, particles.lat) + bx = np.sqrt(2 * kh_zonal) Kyp1 = fieldset.Kh_meridional[particles.time, particles.z, particles.lat + fieldset.dres, particles.lon, particles] Kym1 = fieldset.Kh_meridional[particles.time, particles.z, particles.lat - fieldset.dres, particles.lon, particles] + if fieldset.Kh_meridional.grid._mesh == "spherical": + Kyp1 = meters_to_degrees_meridional(Kyp1) + Kym1 = meters_to_degrees_meridional(Kym1) dKdy = (Kyp1 - Kym1) / (2 * fieldset.dres) ay = v + dKdy - by = np.sqrt(2 * fieldset.Kh_meridional[particles.time, particles.z, particles.lat, particles.lon, particles]) + + kh_meridional = fieldset.Kh_meridional[particles.time, particles.z, particles.lat, particles.lon, particles] + if fieldset.Kh_meridional.grid._mesh == "spherical": + kh_meridional = meters_to_degrees_meridional(kh_meridional) + by = np.sqrt(2 * kh_meridional) # Particle positions are updated only after evaluating all terms. particles.dlon += ax * particles.dt + bx * dWx @@ -103,8 +139,15 @@ def DiffusionUniformKh(particles, fieldset): # pragma: no cover dWx = np.random.normal(0, np.sqrt(np.fabs(particles.dt))) dWy = np.random.normal(0, np.sqrt(np.fabs(particles.dt))) - bx = np.sqrt(2 * fieldset.Kh_zonal[particles]) - by = np.sqrt(2 * fieldset.Kh_meridional[particles]) + kh_zonal = fieldset.Kh_zonal[particles] + kh_meridional = fieldset.Kh_meridional[particles] + + if fieldset.Kh_zonal.grid._mesh == "spherical": + kh_zonal = meters_to_degrees_zonal(kh_zonal, particles.lat) + kh_meridional = meters_to_degrees_meridional(kh_meridional) + + bx = np.sqrt(2 * kh_zonal) + by = np.sqrt(2 * kh_meridional) particles.dlon += bx * dWx particles.dlat += by * dWy diff --git a/tests/test_diffusion.py b/tests/test_diffusion.py index f06ceba27..4f27373ab 100644 --- a/tests/test_diffusion.py +++ b/tests/test_diffusion.py @@ -7,11 +7,8 @@ from parcels import ( Field, FieldSet, - GeographicPolarSquare, - GeographicSquare, Particle, ParticleSet, - Unity, Variable, VectorField, XGrid, @@ -40,13 +37,6 @@ def test_fieldKh_Brownian(mesh): fieldset.add_constant_field("Kh_zonal", kh_zonal, mesh=mesh) fieldset.add_constant_field("Kh_meridional", kh_meridional, mesh=mesh) - if mesh == "spherical": - assert isinstance(fieldset.Kh_zonal.units, GeographicPolarSquare) - assert isinstance(fieldset.Kh_meridional.units, GeographicSquare) - else: - assert isinstance(fieldset.Kh_zonal.units, Unity) - assert isinstance(fieldset.Kh_meridional.units, Unity) - npart = 100 runtime = np.timedelta64(2, "h") @@ -58,10 +48,10 @@ def test_fieldKh_Brownian(mesh): expected_std_lat = np.sqrt(2 * kh_meridional * mesh_conversion**2 * timedelta_to_float(runtime)) tol = 500 * mesh_conversion # effectively 500 m errors - assert np.allclose(np.std(pset.lat), expected_std_lat, atol=tol) - assert np.allclose(np.std(pset.lon), expected_std_lon, atol=tol) - assert np.allclose(np.mean(pset.lon), 0, atol=tol) - assert np.allclose(np.mean(pset.lat), 0, atol=tol) + np.testing.assert_allclose(np.std(pset.lat), expected_std_lat, atol=tol) + np.testing.assert_allclose(np.std(pset.lon), expected_std_lon, atol=tol) + np.testing.assert_allclose(np.mean(pset.lon), 0, atol=tol) + np.testing.assert_allclose(np.mean(pset.lat), 0, atol=tol) @pytest.mark.parametrize("mesh", ["spherical", "flat"]) diff --git a/tests/test_field.py b/tests/test_field.py index 7795c953d..23609219a 100644 --- a/tests/test_field.py +++ b/tests/test_field.py @@ -201,31 +201,31 @@ def test_field_unstructured_z_linear(): # Test above first cell center - for piecewise constant, should return the depth of the first cell center assert np.isclose( - P.eval(time=[0], z=[10.0], y=[30.0], x=[30.0], applyConversion=False), + P.eval(time=[0], z=[10.0], y=[30.0], x=[30.0]), 55.555557, ) # Test below first cell center, but in the first layer - for piecewise constant, should return the depth of the first cell center assert np.isclose( - P.eval(time=[0], z=[65.0], y=[30.0], x=[30.0], applyConversion=False), + P.eval(time=[0], z=[65.0], y=[30.0], x=[30.0]), 55.555557, ) # Test bottom layer - for piecewise constant, should return the depth of the of the bottom layer cell center assert np.isclose( - P.eval(time=[0], z=[900.0], y=[30.0], x=[30.0], applyConversion=False), + P.eval(time=[0], z=[900.0], y=[30.0], x=[30.0]), 944.44445801, ) W = Field(name="W", data=ds.W, grid=grid, interp_method=UxPiecewiseLinearNode) assert np.isclose( - W.eval(time=[0], z=[10.0], y=[30.0], x=[30.0], applyConversion=False), + W.eval(time=[0], z=[10.0], y=[30.0], x=[30.0]), 10.0, ) assert np.isclose( - W.eval(time=[0], z=[65.0], y=[30.0], x=[30.0], applyConversion=False), + W.eval(time=[0], z=[65.0], y=[30.0], x=[30.0]), 65.0, ) assert np.isclose( - W.eval(time=[0], z=[900.0], y=[30.0], x=[30.0], applyConversion=False), + W.eval(time=[0], z=[900.0], y=[30.0], x=[30.0]), 900.0, ) @@ -239,13 +239,12 @@ def test_field_constant_in_time(): # Assert that the field can be evaluated at any time, and returns the same value time = np.datetime64("2000-01-01T00:00:00") - P1 = P.eval(time=time, z=[10.0], y=[30.0], x=[30.0], applyConversion=False) + P1 = P.eval(time=time, z=[10.0], y=[30.0], x=[30.0]) P2 = P.eval( time=time + np.timedelta64(1, "D"), z=[10.0], y=[30.0], x=[30.0], - applyConversion=False, ) assert np.isclose(P1, P2) diff --git a/tests/test_fieldset.py b/tests/test_fieldset.py index 131826647..1a2515bcb 100644 --- a/tests/test_fieldset.py +++ b/tests/test_fieldset.py @@ -6,7 +6,7 @@ import pytest import xarray as xr -from parcels import Field, Geographic, GeographicPolar, ParticleFile, ParticleSet, VectorField, XGrid +from parcels import Field, ParticleFile, ParticleSet, VectorField, XGrid from parcels._core.fieldset import CalendarError, FieldSet, _datetime_to_msg from parcels._datasets.structured.circulation_models import datasets as datasets_circulation_models from parcels._datasets.structured.generic import T as T_structured @@ -260,14 +260,6 @@ def test_fieldset_from_copernicusmarine(ds, caplog): assert "renamed it to 'V'" in caplog.text -@pytest.mark.parametrize("ds", _COPERNICUS_DATASETS) -def test_grid_mesh_units_from_copernicusmarine(ds): - fieldset = FieldSet.from_copernicusmarine(ds) - assert fieldset.U.grid._mesh == "spherical" - assert isinstance(fieldset.U.units, GeographicPolar) - assert isinstance(fieldset.V.units, Geographic) - - @pytest.mark.parametrize("ds", [datasets_circulation_models["ds_copernicusmarine"].copy()]) def test_fieldset_from_copernicusmarine_missing_axis(ds, caplog): del ds["latitude"].attrs["axis"] diff --git a/tests/test_interpolation.py b/tests/test_interpolation.py index c282fe7f0..26f1ad0cb 100644 --- a/tests/test_interpolation.py +++ b/tests/test_interpolation.py @@ -170,20 +170,19 @@ def test_interpolation_mesh_type(mesh, npart=10): U = Field("U", ds["U"], grid, interp_method=XLinear) V = Field("V", ds["V"], grid, interp_method=XLinear) UV = VectorField("UV", U, V) + fieldset = FieldSet([U, V, UV]) lat = 30.0 time = 0.0 u_expected = 1.0 if mesh == "flat" else 1.0 / (1852 * 60 * np.cos(np.radians(lat))) - assert np.isclose(U[time, 0, lat, 0], u_expected, atol=1e-7) - assert V[time, 0, lat, 0] == 0.0 + assert fieldset.U.eval(time, 0, lat, 0) == 1.0 + assert fieldset.V[time, 0, lat, 0] == 0.0 - u, v = UV[time, 0, lat, 0] + u, v = fieldset.UV[time, 0, lat, 0] assert np.isclose(u, u_expected, atol=1e-7) assert v == 0.0 - assert U.eval(time, 0, lat, 0, applyConversion=False) == 1 - interp_methods = { "linear": XLinear, diff --git a/tests/test_uxarray_fieldset.py b/tests/test_uxarray_fieldset.py index 05aa3c5b2..174a17eea 100644 --- a/tests/test_uxarray_fieldset.py +++ b/tests/test_uxarray_fieldset.py @@ -122,26 +122,29 @@ def test_fesom2_square_delaunay_uniform_z_coordinate_eval(): P = Field(name="p", data=ds.p, grid=grid, interp_method=UxPiecewiseLinearNode) fieldset = FieldSet([UVW, P, UVW.U, UVW.V, UVW.W]) + (u, v, w) = fieldset.UVW.eval(time=[0.0], z=[1.0], y=[30.0], x=[30.0], applyConversion=False) + assert np.allclose([u.item(), v.item(), w.item()], [1.0, 1.0, 0.0], rtol=1e-3, atol=1e-6) + assert np.isclose( - fieldset.U.eval(time=[0.0], z=[1.0], y=[30.0], x=[30.0], applyConversion=False), + fieldset.U.eval(time=[0.0], z=[1.0], y=[30.0], x=[30.0]), 1.0, rtol=1e-3, atol=1e-6, ) assert np.isclose( - fieldset.V.eval(time=[0.0], z=[1.0], y=[30.0], x=[30.0], applyConversion=False), + fieldset.V.eval(time=[0.0], z=[1.0], y=[30.0], x=[30.0]), 1.0, rtol=1e-3, atol=1e-6, ) assert np.isclose( - fieldset.W.eval(time=[0.0], z=[1.0], y=[30.0], x=[30.0], applyConversion=False), + fieldset.W.eval(time=[0.0], z=[1.0], y=[30.0], x=[30.0]), 0.0, rtol=1e-3, atol=1e-6, ) assert np.isclose( - fieldset.p.eval(time=[0.0], z=[1.0], y=[30.0], x=[30.0], applyConversion=False), + fieldset.p.eval(time=[0.0], z=[1.0], y=[30.0], x=[30.0]), 1.0, rtol=1e-3, atol=1e-6, @@ -163,7 +166,7 @@ def test_fesom2_square_delaunay_antimeridian_eval(): ) fieldset = FieldSet([P]) - assert np.isclose(fieldset.p.eval(time=[0], z=[1.0], y=[30.0], x=[-170.0], applyConversion=False), 1.0) - assert np.isclose(fieldset.p.eval(time=[0], z=[1.0], y=[30.0], x=[-180.0], applyConversion=False), 1.0) - assert np.isclose(fieldset.p.eval(time=[0], z=[1.0], y=[30.0], x=[180.0], applyConversion=False), 1.0) - assert np.isclose(fieldset.p.eval(time=[0], z=[1.0], y=[30.0], x=[170.0], applyConversion=False), 1.0) + assert np.isclose(fieldset.p.eval(time=[0], z=[1.0], y=[30.0], x=[-170.0]), 1.0) + assert np.isclose(fieldset.p.eval(time=[0], z=[1.0], y=[30.0], x=[-180.0]), 1.0) + assert np.isclose(fieldset.p.eval(time=[0], z=[1.0], y=[30.0], x=[180.0]), 1.0) + assert np.isclose(fieldset.p.eval(time=[0], z=[1.0], y=[30.0], x=[170.0]), 1.0)