From 38a13b493936e28aabaeddb8cced1de2812981e8 Mon Sep 17 00:00:00 2001 From: Vecko <36369090+VeckoTheGecko@users.noreply.github.com> Date: Mon, 10 Mar 2025 14:04:29 +0100 Subject: [PATCH] Remove NestedField --- docs/documentation/index.rst | 1 - docs/examples/tutorial_NestedFields.ipynb | 242 ---------------------- parcels/field.py | 79 +------ parcels/fieldset.py | 43 +--- parcels/interaction/interactionkernel.py | 4 +- parcels/kernel.py | 8 +- parcels/particleset.py | 11 +- tests/test_fieldset_sampling.py | 5 +- 8 files changed, 20 insertions(+), 373 deletions(-) delete mode 100644 docs/examples/tutorial_NestedFields.ipynb diff --git a/docs/documentation/index.rst b/docs/documentation/index.rst index 3e3fc9b967..78baeea3a5 100644 --- a/docs/documentation/index.rst +++ b/docs/documentation/index.rst @@ -20,7 +20,6 @@ Parcels has several documentation and tutorial Jupyter notebooks and scripts whi ../examples/tutorial_nemo_curvilinear.ipynb ../examples/tutorial_nemo_3D.ipynb ../examples/tutorial_croco_3D.ipynb - ../examples/tutorial_NestedFields.ipynb ../examples/tutorial_timevaryingdepthdimensions.ipynb ../examples/tutorial_periodic_boundaries.ipynb ../examples/tutorial_interpolation.ipynb diff --git a/docs/examples/tutorial_NestedFields.ipynb b/docs/examples/tutorial_NestedFields.ipynb deleted file mode 100644 index b399588fc4..0000000000 --- a/docs/examples/tutorial_NestedFields.ipynb +++ /dev/null @@ -1,242 +0,0 @@ -{ - "cells": [ - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# NestedFields\n" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "In some applications, you may have access to different fields that each cover only part of the region of interest. Then, you would like to combine them all together. You may also have a field covering the entire region and another one only covering part of it, but with a higher resolution. The set of those fields form what we call nested fields.\n", - "\n", - "It is possible to combine all those fields with kernels, either with different if/else statements depending on particle position.\n", - "\n", - "However, an easier way to work with nested fields in Parcels is to combine all those fields into one `NestedField` object. The Parcels code will then try to successively interpolate the different fields.\n", - "\n", - "For each Particle, the algorithm is the following:\n", - "\n", - "1. Interpolate the particle onto the first `Field` in the `NestedFields` list.\n", - "\n", - "2. If the interpolation succeeds or if an error other than `ErrorOutOfBounds` is thrown, the function is stopped.\n", - "\n", - "3. If an `ErrorOutOfBounds` is thrown, try step 1) again with the next `Field` in the `NestedFields` list\n", - "\n", - "4. If interpolation on the last `Field` in the `NestedFields` list also returns an `ErrorOutOfBounds`, then the Particle is flagged as OutOfBounds.\n", - "\n", - "This algorithm means that **the order of the fields in the 'NestedField' matters**. In particular, the smallest/finest resolution fields have to be listed _before_ the larger/coarser resolution fields.\n", - "\n", - "Note also that `pset.execute()` can be _much_ slower on `NestedField` objects than on normal `Fields`. This is because the search for the correct `Field` will always start with the first `Field` in the list. This means that especially in cases where many particles are in the outermost `Field`, simulations can become very slow.\n", - "\n", - "This tutorial shows how to use these `NestedField` with a very idealised example.\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "import matplotlib.pyplot as plt\n", - "import numpy as np\n", - "import xarray as xr\n", - "\n", - "import parcels" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "First define a zonal and meridional velocity field defined on a high resolution (dx = 100m) 2kmx2km grid with a flat mesh. The zonal velocity is uniform and 1 m/s, and the meridional velocity is equal to 0.5 _ cos(lon / 200 _ pi / 2) m/s.\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "dim = 21\n", - "lon = np.linspace(0.0, 2e3, dim, dtype=np.float32)\n", - "lat = np.linspace(0.0, 2e3, dim, dtype=np.float32)\n", - "lon_g, lat_g = np.meshgrid(lon, lat)\n", - "V1_data = np.cos(lon_g / 200 * np.pi / 2)\n", - "U1 = parcels.Field(\"U1\", np.ones((dim, dim), dtype=np.float32), lon=lon, lat=lat)\n", - "V1 = parcels.Field(\"V1\", V1_data, grid=U1.grid)" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Now define the same velocity field on a low resolution (dx = 2km) 20kmx4km grid.\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "xdim = 11\n", - "ydim = 3\n", - "lon = np.linspace(-2e3, 18e3, xdim, dtype=np.float32)\n", - "lat = np.linspace(-1e3, 3e3, ydim, dtype=np.float32)\n", - "lon_g, lat_g = np.meshgrid(lon, lat)\n", - "V2_data = np.cos(lon_g / 200 * np.pi / 2)\n", - "U2 = parcels.Field(\"U2\", np.ones((ydim, xdim), dtype=np.float32), lon=lon, lat=lat)\n", - "V2 = parcels.Field(\"V2\", V2_data, grid=U2.grid)" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "We now combine those fields into a `NestedField` and create the fieldset\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "U = parcels.NestedField(\"U\", [U1, U2])\n", - "V = parcels.NestedField(\"V\", [V1, V2])\n", - "fieldset = parcels.FieldSet(U, V)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "pset = parcels.ParticleSet(fieldset, pclass=parcels.Particle, lon=[0], lat=[1000])\n", - "\n", - "output_file = pset.ParticleFile(\n", - " name=\"NestedFieldParticle.zarr\", outputdt=50, chunks=(1, 100)\n", - ")\n", - "\n", - "pset.execute(parcels.AdvectionRK4, runtime=14000, dt=10, output_file=output_file)\n", - "\n", - "ds = xr.open_zarr(\"NestedFieldParticle.zarr\")\n", - "plt.plot(ds.lon.T, ds.lat.T, \".-\")\n", - "plt.plot([0, 2e3, 2e3, 0, 0], [0, 0, 2e3, 2e3, 0], c=\"orange\")\n", - "plt.plot([-2e3, 18e3, 18e3, -2e3, -2e3], [-1e3, -1e3, 3e3, 3e3, -1e3], c=\"green\");" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "As we observe, there is a change of dynamic at lon=2000, which corresponds to the change of grid.\n", - "\n", - "The analytical solution to the problem:\n", - "\n", - "\\begin{align}\n", - "dx/dt &= 1;\\\\\n", - "dy/dt &= \\cos(x \\pi/400);\\\\\n", - "\\text{with } x(0) &= 0, y(0) = 1000\n", - "\\end{align}\n", - "\n", - "is\n", - "\n", - "\\begin{align}\n", - "x(t) &= t;\\\\\n", - "y(t) &= 1000 + 400/\\pi \\sin(t \\pi / 400)\n", - "\\end{align}\n", - "\n", - "which is captured by the High Resolution field (orange area) but not the Low Resolution one (green area).\n" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Keep track of the field interpolated\n", - "\n", - "For different reasons, you may want to keep track of the field you have interpolated. You can do that easily by creating another field that share the grid with original fields.\n", - "Watch out that this operation has a cost of a full interpolation operation.\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Need to redefine fieldset\n", - "fieldset = parcels.FieldSet(U, V)\n", - "\n", - "ones_array1 = np.ones((U1.grid.ydim, U1.grid.xdim), dtype=np.float32)\n", - "F1 = parcels.Field(\"F1\", ones_array1, grid=U1.grid)\n", - "\n", - "ones_array2 = np.ones((U2.grid.ydim, U2.grid.xdim), dtype=np.float32)\n", - "F2 = parcels.Field(\"F2\", 2 * ones_array2, grid=U2.grid)\n", - "\n", - "F = parcels.NestedField(\"F\", [F1, F2])\n", - "fieldset.add_field(F)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "def SampleNestedFieldIndex(particle, fieldset, time):\n", - " particle.f = fieldset.F[time, particle.depth, particle.lat, particle.lon]\n", - "\n", - "\n", - "SampleParticle = parcels.Particle.add_variable(\"f\", dtype=np.int32)\n", - "\n", - "pset = parcels.ParticleSet(fieldset, pclass=SampleParticle, lon=[1000], lat=[500])\n", - "pset.execute(SampleNestedFieldIndex, runtime=1)\n", - "print(\n", - " f\"Particle ({pset[0].lon:g}, {pset[0].lat:g}) \"\n", - " f\"interpolates Field #{int(pset[0].f)}\"\n", - ")\n", - "\n", - "pset[0].lon = 10000\n", - "pset.execute(SampleNestedFieldIndex, runtime=1)\n", - "print(\n", - " f\"Particle ({pset[0].lon:g}, {pset[0].lat:g}) \"\n", - " f\"interpolates Field #{int(pset[0].f)}\"\n", - ")" - ] - } - ], - "metadata": { - "kernelspec": { - "display_name": "parcels", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.12.3" - } - }, - "nbformat": 4, - "nbformat_minor": 2 -} diff --git a/parcels/field.py b/parcels/field.py index d7056d37bb..12ac0bde4f 100644 --- a/parcels/field.py +++ b/parcels/field.py @@ -53,7 +53,7 @@ from parcels.fieldset import FieldSet -__all__ = ["Field", "NestedField", "VectorField"] +__all__ = ["Field", "VectorField"] def _isParticle(key): @@ -140,11 +140,6 @@ class Field: Write the Field in NetCDF format at the same frequency as the ParticleFile outputdt, using a filenaming scheme based on the ParticleFile name - Examples - -------- - For usage examples see the following tutorials: - - * `Nested Fields <../examples/tutorial_NestedFields.ipynb>`__ """ allow_time_extrapolation: bool @@ -1352,75 +1347,3 @@ def __getitem__(self, key): return self.eval(*key) except tuple(AllParcelsErrorCodes.keys()) as error: return _deal_with_errors(error, key, vector_type=self.vector_type) - - -class NestedField(list): - """NestedField is a class that allows for interpolation of fields on different grids of potentially varying resolution. - - The NestedField class is a list of Fields where the first Field that contains the particle within the domain is then used for interpolation. - This induces that the order of the fields in the list matters. - Each one it its turn, a field is interpolated: if the interpolation succeeds or if an error other - than `ErrorOutOfBounds` is thrown, the function is stopped. Otherwise, next field is interpolated. - NestedField returns an `ErrorOutOfBounds` only if last field is as well out of boundaries. - NestedField is composed of either Fields or VectorFields. - - Parameters - ---------- - name : str - Name of the NestedField - F : list of Field - List of fields (order matters). F can be a scalar Field, a VectorField, or the zonal component (U) of the VectorField - V : list of Field - List of fields defining the meridional component of a VectorField, if F is the zonal component. (default: None) - W : list of Field - List of fields defining the vertical component of a VectorField, if F and V are the zonal and meridional components (default: None) - - - Examples - -------- - See `here <../examples/tutorial_NestedFields.ipynb>`__ - for a detailed tutorial - - """ - - def __init__(self, name: str, F, V=None, W=None): - if V is None: - if isinstance(F[0], VectorField): - vector_type = F[0].vector_type - for Fi in F: - assert isinstance(Fi, Field) or ( - isinstance(Fi, VectorField) and Fi.vector_type == vector_type - ), "Components of a NestedField must be Field or VectorField" - self.append(Fi) - elif W is None: - for i, Fi, Vi in zip(range(len(F)), F, V, strict=True): - assert isinstance(Fi, Field) and isinstance( - Vi, Field - ), "F, and V components of a NestedField must be Field" - self.append(VectorField(f"{name}_{i}", Fi, Vi)) - else: - for i, Fi, Vi, Wi in zip(range(len(F)), F, V, W, strict=True): - assert ( - isinstance(Fi, Field) and isinstance(Vi, Field) and isinstance(Wi, Field) - ), "F, V and W components of a NestedField must be Field" - self.append(VectorField(f"{name}_{i}", Fi, Vi, Wi)) - self.name = name - - def __getitem__(self, key): - if isinstance(key, int): - return list.__getitem__(self, key) - else: - for iField in range(len(self)): - try: - if _isParticle(key): - val = list.__getitem__(self, iField).eval(key.time, key.depth, key.lat, key.lon, particle=None) - else: - val = list.__getitem__(self, iField).eval(*key) - break - except tuple(AllParcelsErrorCodes.keys()) as error: - if iField == len(self) - 1: - vector_type = self[iField].vector_type if isinstance(self[iField], VectorField) else None - return _deal_with_errors(error, key, vector_type=vector_type) - else: - pass - return val diff --git a/parcels/fieldset.py b/parcels/fieldset.py index a7d2eb0fa4..cf47357e4d 100644 --- a/parcels/fieldset.py +++ b/parcels/fieldset.py @@ -9,7 +9,7 @@ from parcels._compat import MPI from parcels._typing import GridIndexingType, InterpMethodOption, Mesh -from parcels.field import Field, NestedField, VectorField +from parcels.field import Field, VectorField from parcels.grid import Grid from parcels.gridset import GridSet from parcels.particlefile import ParticleFile @@ -35,7 +35,7 @@ class FieldSet: in custom kernels. """ - def __init__(self, U: Field | NestedField | None, V: Field | NestedField | None, fields=None): + def __init__(self, U: Field | None, V: Field | None, fields=None): self.gridset = GridSet() self._completed: bool = False self._particlefile: ParticleFile | None = None @@ -150,7 +150,7 @@ def from_data( v = fields.pop("V", None) return cls(u, v, fields=fields) - def add_field(self, field: Field | NestedField, name: str | None = None): + def add_field(self, field: Field, name: str | None = None): """Add a :class:`parcels.field.Field` object to the FieldSet. Parameters @@ -166,8 +166,6 @@ def add_field(self, field: Field | NestedField, name: str | None = None): -------- For usage examples see the following tutorials: - * `Nested Fields <../examples/tutorial_NestedFields.ipynb>`__ - * `Unit converters <../examples/tutorial_unitconverters.ipynb>`__ (Default value = None) """ @@ -179,11 +177,6 @@ def add_field(self, field: Field | NestedField, name: str | None = None): if hasattr(self, name): # check if Field with same name already exists when adding new Field raise RuntimeError(f"FieldSet already has a Field with name '{name}'") - if isinstance(field, NestedField): - setattr(self, name, field) - for fld in field: - self.gridset.add_grid(fld) - fld.fieldset = self else: setattr(self, name, field) self.gridset.add_grid(field) @@ -222,21 +215,12 @@ def add_vector_field(self, vfield): if isinstance(v, Field) and (v not in self.get_fields()): self.add_field(v) vfield.fieldset = self - if isinstance(vfield, NestedField): - for f in vfield: - f.fieldset = self def _add_UVfield(self): if not hasattr(self, "UV") and hasattr(self, "U") and hasattr(self, "V"): - if isinstance(self.U, NestedField): - self.add_vector_field(NestedField("UV", self.U, self.V)) - else: - self.add_vector_field(VectorField("UV", self.U, self.V)) + self.add_vector_field(VectorField("UV", self.U, self.V)) if not hasattr(self, "UVW") and hasattr(self, "W"): - if isinstance(self.U, NestedField): - self.add_vector_field(NestedField("UVW", self.U, self.V, self.W)) - else: - self.add_vector_field(VectorField("UVW", self.U, self.V, self.W)) + self.add_vector_field(VectorField("UVW", self.U, self.V, self.W)) def _check_complete(self): assert self.U, 'FieldSet does not have a Field named "U"' @@ -268,13 +252,8 @@ def check_velocityfields(U, V, W): if V.gridindexingtype != U.gridindexingtype or (W and W.gridindexingtype != U.gridindexingtype): raise ValueError("Not all velocity Fields have the same gridindexingtype") - if isinstance(self.U, NestedField): - w = self.W if hasattr(self, "W") else [None] * len(self.U) - for U, V, W in zip(self.U, self.V, w, strict=True): - check_velocityfields(U, V, W) - else: - W = self.W if hasattr(self, "W") else None - check_velocityfields(self.U, self.V, W) + W = self.W if hasattr(self, "W") else None + check_velocityfields(self.U, self.V, W) for g in self.gridset.grids: g._check_zonal_periodic() @@ -288,7 +267,7 @@ def check_velocityfields(U, V, W): self._add_UVfield() for f in self.get_fields(): - if isinstance(f, (VectorField, NestedField)) or f._dataFiles is None: + if isinstance(f, VectorField) or f._dataFiles is None: continue self._completed = True @@ -1104,12 +1083,6 @@ def get_fields(self) -> list[Field | VectorField]: if type(v) in [Field, VectorField]: if v not in fields: fields.append(v) - elif isinstance(v, NestedField): - if v not in fields: - fields.append(v) - for v2 in v: - if v2 not in fields: - fields.append(v2) return fields def add_constant(self, name, value): diff --git a/parcels/interaction/interactionkernel.py b/parcels/interaction/interactionkernel.py index f2192c18f6..0f51c821df 100644 --- a/parcels/interaction/interactionkernel.py +++ b/parcels/interaction/interactionkernel.py @@ -5,7 +5,7 @@ import numpy as np from parcels._compat import MPI -from parcels.field import NestedField, VectorField +from parcels.field import VectorField from parcels.kernel import BaseKernel from parcels.tools.statuscodes import StatusCode @@ -121,7 +121,7 @@ def execute_python(self, pset, endtime, dt): """ if self.fieldset is not None: for f in self.fieldset.get_fields(): - if isinstance(f, (VectorField, NestedField)): + if isinstance(f, VectorField): continue f.data = np.array(f.data) diff --git a/parcels/kernel.py b/parcels/kernel.py index 2bfcd41dd3..300a9efc97 100644 --- a/parcels/kernel.py +++ b/parcels/kernel.py @@ -16,7 +16,7 @@ AdvectionRK4_3D_CROCO, AdvectionRK45, ) -from parcels.field import Field, NestedField, VectorField +from parcels.field import Field, VectorField from parcels.grid import GridType from parcels.tools.statuscodes import ( StatusCode, @@ -216,10 +216,6 @@ def check_fieldsets_in_kernels(self, pyfunc): # TODO v4: this can go into anoth and self._fieldset.W._scaling_factor > 0 ): warning = True - if isinstance(self._fieldset.W, NestedField): - for f in self._fieldset.W: - if f._creation_log != "from_nemo" and f._scaling_factor is not None and f._scaling_factor > 0: - warning = True if warning: warnings.warn( "Note that in AdvectionRK4_3D, vertical velocity is assumed positive towards increasing z. " @@ -337,7 +333,7 @@ def execute(self, pset, endtime, dt): if pset.fieldset is not None: for f in self.fieldset.get_fields(): - if isinstance(f, (VectorField, NestedField)): + if isinstance(f, VectorField): continue f.data = np.array(f.data) diff --git a/parcels/particleset.py b/parcels/particleset.py index e6946e2d8c..2655cbdf88 100644 --- a/parcels/particleset.py +++ b/parcels/particleset.py @@ -12,7 +12,7 @@ from parcels._compat import MPI from parcels.application_kernels.advection import AdvectionRK4 -from parcels.field import Field, NestedField +from parcels.field import Field from parcels.grid import CurvilinearGrid, GridType from parcels.interaction.interactionkernel import InteractionKernel from parcels.interaction.neighborsearch import ( @@ -305,13 +305,8 @@ def __getitem__(self, index): @staticmethod def lonlatdepth_dtype_from_field_interp_method(field): - if isinstance(field, NestedField): - for f in field: - if f.interp_method == "cgrid_velocity": - return np.float64 - else: - if field.interp_method == "cgrid_velocity": - return np.float64 + if field.interp_method == "cgrid_velocity": + return np.float64 return np.float32 @property diff --git a/tests/test_fieldset_sampling.py b/tests/test_fieldset_sampling.py index 1b9f674e5f..c8033a658d 100644 --- a/tests/test_fieldset_sampling.py +++ b/tests/test_fieldset_sampling.py @@ -12,7 +12,6 @@ Field, FieldSet, Geographic, - NestedField, Particle, ParticleSet, StatusCode, @@ -758,7 +757,11 @@ def test_multiple_grid_addlater_error(): assert fail +@pytest.mark.v4alpha +@pytest.mark.xfail(reason="Implementation of NestedFields is being reconsidered in v4.") def test_nestedfields(): + from parcels import NestedField + xdim = 10 ydim = 20