Skip to content

Implementing FieldSet.from_croco() #2450

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
erikvansebille wants to merge 12 commits into
base: fieldset_from_nemo
Choose a base branch
from
Open

Implementing FieldSet.from_croco() #2450

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
12 commits
Select commit Hold shift + click to select a range
d3712b6
Fixing bug in index_search when curvilinear grid is flat
erikvansebille Jan 2, 2026
8e1ec53
Moving CROCO Sigma kernels to their own file
erikvansebille Jan 2, 2026
92f810c
Remove automatic replacement of CROCO kernel under the hood
erikvansebille Jan 2, 2026
cd55ab5
Moving CROCO tests to new test_sigmagrids file
erikvansebille Jan 2, 2026
b50c291
Updating croco 3D tutorial
erikvansebille Jan 2, 2026
a9672a5
Always attach a UV Field
erikvansebille Jan 2, 2026
19f43ba
Updating RK4 and RK2 CROCO
erikvansebille Jan 2, 2026
1d0a6bd
Adding first version of fieldset.from_croco()
erikvansebille Jan 2, 2026
a9f79d1
Merge branch 'fieldset_from_nemo' into fieldset_from_croco
erikvansebille Jan 5, 2026
8f8eaf8
Fixing CROCO interpolation to use XLinear for W
erikvansebille Jan 6, 2026
6b6cc66
Removing omega testing against v3
erikvansebille Jan 6, 2026
462e9b7
Merge branch 'fieldset_from_nemo' into fieldset_from_croco
erikvansebille Jan 6, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
163 changes: 79 additions & 84 deletions docs/user_guide/examples_v3/tutorial_croco_3D.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -34,16 +34,16 @@
"metadata": {},
"outputs": [],
"source": [
"import os\n",
"\n",
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"import xarray as xr\n",
"\n",
"import parcels\n",
"\n",
"example_dataset_folder = parcels.download_example_dataset(\"CROCOidealized_data\")\n",
"file = os.path.join(example_dataset_folder, \"CROCO_idealized.nc\")"
"data_folder = parcels.download_example_dataset(\"CROCOidealized_data\")\n",
"ds_fields = xr.open_dataset(data_folder / \"CROCO_idealized.nc\")\n",
"\n",
"ds_fields.load(); # Preload data to speed up access"
]
},
{
Expand All @@ -59,17 +59,16 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we create a FieldSet object using the `FieldSet.from_croco()` method. Note that CROCO is a C-grid (with similar indexing at MITgcm), so we need to provide the longitudes and latitudes of the $\\rho$-points of the grid (`lon_rho` and `lat_rho`). We also need to provide the sigma levels at the depth points (`s_w`). Finally, it is important to also provide the bathymetry field (`h`), which is needed to convert the depth levels of the particles to sigma-coordinates."
"Now we create a FieldSet object using the `FieldSet.from_croco()` method. Note that CROCO is a C-grid (with similar indexing at MITgcm)."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<div class=\"alert alert-info\">\n",
"\n",
"__Note__ that in the code below we use the `w` velocity field for vertical velocity. However, it is unclear whether this is always the right choice. CROCO (and ROMS) also output an `omega` field, which may be more appropriate to use. The idealised simulation below only works when using `w`, though. In other simulations, it is recommended to test whether `omega` provides more realistic results. See https://github.com/OceanParcels/Parcels/discussions/1728 for more information.\n",
"</div>"
"```{note}\n",
"In the code below, we use the `w` velocity field for vertical velocity. However, it is unclear whether this is always the right choice. CROCO (and ROMS) also output an `omega` field, which may be more appropriate to use. The idealised simulation below only works when using `w`, though. In other simulations, it is recommended to test whether `omega` provides more realistic results. See https://github.com/OceanParcels/Parcels/discussions/1728 for more information.\n",
"```"
]
},
{
Expand All @@ -78,29 +77,12 @@
"metadata": {},
"outputs": [],
"source": [
"variables = {\"U\": \"u\", \"V\": \"v\", \"W\": \"w\", \"H\": \"h\", \"Zeta\": \"zeta\", \"Cs_w\": \"Cs_w\"}\n",
"\n",
"lon_rho = \"x_rho\" # Note, this would be \"lon_rho\" for a dataset on a spherical grid\n",
"lat_rho = \"y_rho\" # Note ,this would be \"lat_rho\" for a dataset on a spherical grid\n",
"vars_to_use = [\"u\", \"v\", \"w\", \"h\", \"zeta\", \"Cs_w\"]\n",
"ds_fields = ds_fields[vars_to_use]\n",
"\n",
"dimensions = {\n",
" \"U\": {\"lon\": lon_rho, \"lat\": lat_rho, \"depth\": \"s_w\", \"time\": \"time\"},\n",
" \"V\": {\"lon\": lon_rho, \"lat\": lat_rho, \"depth\": \"s_w\", \"time\": \"time\"},\n",
" \"W\": {\"lon\": lon_rho, \"lat\": lat_rho, \"depth\": \"s_w\", \"time\": \"time\"},\n",
" \"H\": {\"lon\": lon_rho, \"lat\": lat_rho},\n",
" \"Zeta\": {\"lon\": lon_rho, \"lat\": lat_rho, \"time\": \"time\"},\n",
" \"Cs_w\": {\"depth\": \"s_w\"},\n",
"}\n",
"fieldset = parcels.FieldSet.from_croco(\n",
" file,\n",
" variables,\n",
" dimensions,\n",
" hc=xr.open_dataset(\n",
" file\n",
" ).hc.values, # Note this stretching parameter is needed for the vertical grid\n",
" allow_time_extrapolation=True, # Note, this is only needed for this specific example dataset, that has only one snapshot\n",
" mesh=\"flat\", # Note, this is only needed for this specific example dataset, that has been created on a 'flat' mesh (i.e. in km instead of in degrees)\n",
")"
"mesh = \"flat\" # Note, this is only needed for this specific example dataset, that has been created on a 'flat' mesh (i.e. in km instead of in degrees)\n",
"fieldset = parcels.FieldSet.from_croco(ds_fields, mesh=mesh)\n",
"fieldset.add_constant(\"hc\", ds_fields.hc)"
]
},
{
Expand All @@ -120,24 +102,27 @@
" [40e3, 80e3, 120e3],\n",
" [100, -10, -130, -250, -400, -850, -1400, -1550],\n",
")\n",
"Y = np.ones(X.size) * fieldset.U.grid.lat[25]\n",
"Y = np.ones(X.size) * 98000\n",
"\n",
"\n",
"def DeleteParticle(particle, fieldset, time):\n",
" if particle.state >= 50:\n",
" particle.delete()\n",
"def DeleteParticle(particles, fieldset):\n",
" any_error = particles.state >= 50\n",
" particles[any_error].state = parcels.StatusCode.Delete\n",
"\n",
"\n",
"pset = parcels.ParticleSet(\n",
" fieldset=fieldset, pclass=parcels.Particle, lon=X, lat=Y, depth=Z\n",
" fieldset=fieldset, pclass=parcels.Particle, lon=X, lat=Y, z=Z\n",
")\n",
"\n",
"outputfile = pset.ParticleFile(name=\"croco_particles3D.zarr\", outputdt=5000)\n",
"outputfile = parcels.ParticleFile(\n",
" store=\"croco_particles3D.zarr\",\n",
" outputdt=np.timedelta64(5000, \"s\"),\n",
")\n",
"\n",
"pset.execute(\n",
" [parcels.kernels.AdvectionRK4_3D, DeleteParticle],\n",
" runtime=5e4,\n",
" dt=100,\n",
" [parcels.kernels.AdvectionRK4_3D_CROCO, DeleteParticle],\n",
" runtime=np.timedelta64(50_000, \"s\"),\n",
" dt=np.timedelta64(100, \"s\"),\n",
" output_file=outputfile,\n",
")"
]
Expand All @@ -162,15 +147,20 @@
"fig, ax = plt.subplots(1, 1, figsize=(6, 4))\n",
"ds = xr.open_zarr(\"croco_particles3D.zarr\")\n",
"\n",
"ax.plot(X / 1e3, Z, \"k.\", label=\"Initial positions\")\n",
"ax.plot(ds.lon.T / 1e3, ds.z.T, \".-\")\n",
"\n",
"dsCROCO = xr.open_dataset(file)\n",
"for z in dsCROCO.s_w.values:\n",
" ax.plot(fieldset.H.lon / 1e3, fieldset.H.data[0, 25, :] * z, \"k\", linewidth=0.5)\n",
"for z in ds_fields.s_w.values:\n",
" ax.plot(\n",
" fieldset.h.grid.lon[0, :] / 1e3,\n",
" fieldset.h.data[0, 0, 25, :] * z,\n",
" \"k\",\n",
" linewidth=0.5,\n",
" )\n",
"ax.set_xlabel(\"X [km]\")\n",
"ax.set_xlim(30, 170)\n",
"ax.set_ylabel(\"Depth [m]\")\n",
"ax.set_title(\"Particles in idealized CROCO velocity field using 3D advection\")\n",
"ax.set_title(\"Particles in idealized CROCO velocity field using AdvectionRK4_3D_CROCO\")\n",
"plt.tight_layout()\n",
"plt.show()"
]
Expand All @@ -186,7 +176,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"It may be insightful to compare this 3D run with the `AdvectionRK4_3D` kernel with a run where the vertical velocity (`W`) is set to zero. In that case, the particles will not stay on their initial depth levels but instead follow sigma-layers."
"It may be insightful to compare this 3D run with the `AdvectionRK4_3D_CROCO` kernel with a run where the vertical velocity (`W`) is set to zero. In that case, the particles will not stay on their initial depth levels but instead follow sigma-layers."
]
},
{
Expand All @@ -195,32 +185,44 @@
"metadata": {},
"outputs": [],
"source": [
"import copy\n",
"\n",
"fieldset_noW = copy.copy(fieldset)\n",
"fieldset_noW = parcels.FieldSet.from_croco(ds_fields, mesh=mesh)\n",
"fieldset_noW.add_constant(\"hc\", ds_fields.hc)\n",
"fieldset_noW.W.data[:] = 0.0\n",
"\n",
"X, Z = np.meshgrid(\n",
" [40e3, 80e3, 120e3],\n",
" [100, -10, -130, -250, -400, -850, -1400, -1550],\n",
")\n",
"Y = np.ones(X.size) * 100\n",
"\n",
"pset_noW = parcels.ParticleSet(\n",
" fieldset=fieldset_noW, pclass=parcels.Particle, lon=X, lat=Y, depth=Z\n",
" fieldset=fieldset_noW, pclass=parcels.Particle, lon=X, lat=Y, z=Z\n",
")\n",
"\n",
"outputfile = pset.ParticleFile(name=\"croco_particles_noW.zarr\", outputdt=5000)\n",
"outputfile = parcels.ParticleFile(\n",
" store=\"croco_particles_noW.zarr\", outputdt=np.timedelta64(5000, \"s\")\n",
")\n",
"\n",
"pset_noW.execute(\n",
" [parcels.kernels.AdvectionRK4_3D, DeleteParticle],\n",
" runtime=5e4,\n",
" dt=100,\n",
" [parcels.kernels.AdvectionRK4_3D_CROCO, DeleteParticle],\n",
" runtime=np.timedelta64(50_000, \"s\"),\n",
" dt=np.timedelta64(100, \"s\"),\n",
" output_file=outputfile,\n",
")\n",
"\n",
"fig, ax = plt.subplots(1, 1, figsize=(6, 4))\n",
"ds = xr.open_zarr(\"croco_particles_noW.zarr\")\n",
"\n",
"ax.plot(X / 1e3, Z, \"k.\", label=\"Initial positions\")\n",
"ax.plot(ds.lon.T / 1e3, ds.z.T, \".-\")\n",
"\n",
"dsCROCO = xr.open_dataset(file)\n",
"for z in dsCROCO.s_w.values:\n",
" ax.plot(fieldset.H.lon / 1e3, fieldset.H.data[0, 25, :] * z, \"k\", linewidth=0.5)\n",
"for z in ds_fields.s_w.values:\n",
" ax.plot(\n",
" fieldset.h.grid.lon[0, :] / 1e3,\n",
" fieldset.h.data[0, 0, 25, :] * z,\n",
" \"k\",\n",
" linewidth=0.5,\n",
" )\n",
"ax.set_xlabel(\"X [km]\")\n",
"ax.set_xlim(30, 170)\n",
"ax.set_ylabel(\"Depth [m]\")\n",
Expand All @@ -240,56 +242,49 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"When using `FieldSet.from_croco()`, Parcels knows that depth needs to be converted to sigma-coordinates, before doing any interpolation. This is done under the hood, using code for interpolation (in this case a `T` Field) like\n",
"```python\n",
"# First calculate local sigma level of the particle, by linearly interpolating the scaling function that maps sigma to depth (using local ocean depth H, sea-surface Zeta and stretching parameters Cs_w and hc). See also https://croco-ocean.gitlabpages.inria.fr/croco_doc/model/model.grid.html#vertical-grid-parameters\n",
"h = fieldset.H[time, 0, particle.lat, particle.lon]\n",
"zeta = fieldset.H[time, 0, particle.lat, particle.lon]\n",
"sigma_levels = fieldset.U.grid.depth\n",
"z0 = fieldset.hc * sigma_levels + (h - fieldset.hc) * fieldset.Cs_w\n",
"zvec = z0 + zeta * (1 + (z0 / h))\n",
"zinds = zvec <= z\n",
"if z >= zvec[-1]:\n",
" zi = len(zvec) - 2\n",
"else:\n",
" zi = zinds.argmin() - 1 if z >= zvec[0] else 0\n",
"\n",
"sigma = sigma_levels[zi] + (z - zvec[zi]) * (sigma_levels[zi + 1] - sigma_levels[zi]) / (zvec[zi + 1] - zvec[zi])\n",
"When using `FieldSet.from_croco()`, Parcels needs to convert the particles.depth to sigma-coordinates, before doing any interpolation. This is done with the `convert_z_to_sigma_croco()` function. Interpolating onto a Field is then done like:\n",
"\n",
"# Now interpolate the field to the sigma level\n",
"temp = fieldset.T[time, sigma, particle.lat, particle.lon]\n",
"```python\n",
"def SampleTempCroco(particles, fieldset):\n",
" \"\"\"Sample temperature field on a CROCO sigma grid by first converting z to sigma levels.\"\"\"\n",
" sigma = convert_z_to_sigma_croco(fieldset, particles.time, particles.z, particles.lat, particles.lon, particles)\n",
" particles.temp = fieldset.T[particles.time, sigma, particles.lat, particles.lon, particles]\n",
"```"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"For the `AdvectionRK4_3D` kernel, Parcels will replace the kernel with `AdvectionRK4_3D_CROCO`, which works slightly different from the normal 3D advection kernel because it converts the vertical velocity in sigma-units. The conversion from depth to sigma is done at every time step, using the code shows above.\n",
"For Advection, you will need to use the `AdvectionRK4_3D_CROCO` kernel, which works slightly different from the normal 3D advection kernel because it converts the vertical velocity in sigma-units. The conversion from depth to sigma is done at every time step, using the `convert_z_to_sigma_croco()` function.\n",
"\n",
"In particular, the following algorithm is used (note that the RK4 version is slightly more complex than this Euler-Forward version, but the idea is identical)\n",
"In particular, the following algorithm is used (note that the RK4 version is slightly more complex than this Euler-Forward version, but the idea is identical). Also note that the vertical velocity is linearly interpolated here, which gives much better results than the default C-grid interpolation.\n",
"\n",
"```python\n",
"(u, v, w) = fieldset.UVW[time, particle.depth, particle.lat, particle.lon, particle] \n",
"sigma = particles.z / fieldset.h[particles.time, 0, particles.lat, particles.lon]\n",
"\n",
"sigma_levels = convert_z_to_sigma_croco(fieldset, particles.time, particles.z, particles.lat, particles.lon, particles)\n",
"(u, v) = fieldset.UV[time, sigma_levels, particle.lat, particle.lon, particle]\n",
"w = fieldset.W[time, sigma_levels, particle.lat, particle.lon, particle]\n",
"\n",
"# scaling the w with the sigma level of the particle\n",
"w_sigma = w * sigma / fieldset.H[time, particle.depth, particle.lat, particle.lon]\n",
"w_sigma = w * sigma / fieldset.h[particles.time, 0, particles.lat, particles.lon]\n",
"\n",
"lon_new = particle.lon + u*particle.dt\n",
"lat_new = particle.lat + v*particle.dt\n",
"lon_new = particles.lon + u*particles.dt\n",
"lat_new = particles.lat + v*particles.dt\n",
"\n",
"# calculating new sigma level\n",
"sigma_new = sigma + w_sigma*particle.dt \n",
"sigma_new = sigma + w_sigma*particles.dt \n",
"\n",
"# Converting back from sigma to depth, at _new_ location\n",
"depth_new = sigma_new * fieldset.H[time, particle.depth, lat_new, lon_new]\n",
"# converting back from sigma to depth, at _new_ location\n",
"depth_new = sigma_new * fieldset.h[particles.time, 0, lat_new, lon_new]\n",
"```"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "parcels",
"display_name": "docs",
"language": "python",
"name": "python3"
},
Expand All @@ -303,7 +298,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.12.3"
"version": "3.12.11"
}
},
"nbformat": 4,
Expand Down
2 changes: 1 addition & 1 deletion docs/user_guide/index.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -31,11 +31,11 @@ TODO: Add links to Reference API throughout
examples/explanation_grids.md
examples/tutorial_nemo_curvilinear.ipynb
examples/tutorial_nemo_3D.ipynb
examples/tutorial_croco_3D.ipynb
examples/tutorial_unitconverters.ipynb
```

<!-- examples/documentation_indexing.ipynb -->
<!-- examples/tutorial_croco_3D.ipynb -->
<!-- examples/tutorial_timevaryingdepthdimensions.ipynb -->

```{toctree}
Expand Down
1 change: 1 addition & 0 deletions docs/user_guide/v4-migration.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,7 @@ Version 4 of Parcels is unreleased at the moment. The information in this migrat
- `math` functions should be replaced with array compatible equivalents (e.g., `math.sin` -> `np.sin`). Instead of `ParcelsRandom` you should use numpy's random functions.
- `particle.depth` has been changed to `particles.z` to be consistent with the [CF conventions for trajectory data](https://cfconventions.org/cf-conventions/cf-conventions.html#trajectory-data), and to make Parcels also generalizable to atmospheric contexts.
- The `InteractionKernel` class has been removed. Since normal Kernels now have access to _all_ particles, particle-particle interaction can be performed within normal Kernels.
- Users need to explicitly use `convert_z_to_sigma_croco` in sampling kernels such as the `AdvectionRK4_3D_CROCO` kernel for CROCO fields, as the automatic conversion from depth to sigma grids under the hood has been removed.

## FieldSet

Expand Down
Loading
Loading