Skip to content

Commit 5f41847

Browse files
reint-fischerreint-fischer
authored
Merge branch 'v4-dev' into diffusion_guide_v4
2 parents 11fb5c3 + 09bb041 commit 5f41847

File tree

3 files changed

+129
-61
lines changed

3 files changed

+129
-61
lines changed

‎docs/conf.py‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -521,7 +521,7 @@ def linkcode_resolve(domain, info):
521521
# -- Options for MyST parser ----------------------------------------------
522522
myst_heading_anchors = 3
523523

524-
myst_enable_extensions = ["substitution"]
524+
myst_enable_extensions = ["substitution", "amsmath", "dollarmath"]
525525

526526
# -- Options for MyST-nb --------------------------------------------------
527527
nb_execution_mode = "cache"

docs/user_guide/examples_v3/tutorial_unitconverters.ipynb renamed to docs/user_guide/examples/tutorial_unitconverters.ipynb

Lines changed: 127 additions & 59 deletions
Original file line numberDiff line numberDiff line change
@@ -5,7 +5,7 @@
55
"cell_type": "markdown",
66
"metadata": {},
77
"source": [
8-
"# Spherical grids and UnitConversion\n"
8+
"# Spherical grids and unit converters\n"
99
]
1010
},
1111
{
@@ -15,15 +15,15 @@
1515
"source": [
1616
"In most applications, Parcels works with `spherical` meshes, where longitude and latitude are given in degrees, while depth is given in meters. But it is also possible to use `flat` meshes, where longitude and latitude are given in meters (note that the dimensions are then still called `longitude` and `latitude` for consistency reasons).\n",
1717
"\n",
18-
"In all cases, velocities are given in m/s. So Parcels seamlessly converts between meters and degrees, under the hood. For transparency, this tutorial explain how this works.\n"
18+
"In all cases, velocities are given in m/s. Parcels seamlessly converts between meters and degrees, under the hood. For transparency, this guide explains how this works.\n"
1919
]
2020
},
2121
{
2222
"attachments": {},
2323
"cell_type": "markdown",
2424
"metadata": {},
2525
"source": [
26-
"Let's first import the relevant modules, and create dictionaries for the `U`, `V` and `temp` data arrays, with the velocities 1 m/s and the temperature 20C.\n"
26+
"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"
2727
]
2828
},
2929
{
@@ -34,6 +34,7 @@
3434
"source": [
3535
"import matplotlib.pyplot as plt\n",
3636
"import numpy as np\n",
37+
"import xarray as xr\n",
3738
"\n",
3839
"import parcels"
3940
]
@@ -44,26 +45,29 @@
4445
"metadata": {},
4546
"outputs": [],
4647
"source": [
47-
"xdim, ydim = (10, 20)\n",
48-
"data = {\n",
49-
" \"U\": np.ones((ydim, xdim), dtype=np.float32),\n",
50-
" \"V\": np.ones((ydim, xdim), dtype=np.float32),\n",
51-
" \"temp\": 20 * np.ones((ydim, xdim), dtype=np.float32),\n",
52-
"}\n",
53-
"dims = {\n",
54-
" \"lon\": np.linspace(-15, 5, xdim, dtype=np.float32),\n",
55-
" \"lat\": np.linspace(35, 60, ydim, dtype=np.float32),\n",
56-
"}"
48+
"from parcels._datasets.structured.generated import simple_UV_dataset\n",
49+
"\n",
50+
"nlat = 10\n",
51+
"nlon = 18\n",
52+
"ds = simple_UV_dataset(dims=(1, 1, nlat, nlon), mesh=\"spherical\").isel(time=0, depth=0)\n",
53+
"ds[\"temperature\"] = ds[\"U\"] + 20 # add temperature field of 20 deg\n",
54+
"ds[\"U\"].data[:] = 1.0 # set U to 1 m/s\n",
55+
"ds[\"V\"].data[:] = 1.0 # set V to 1 m/s\n",
56+
"ds"
5757
]
5858
},
5959
{
6060
"attachments": {},
6161
"cell_type": "markdown",
6262
"metadata": {},
6363
"source": [
64-
"We can convert these data and dims to a FieldSet object using `FieldSet.from_data`. We add the argument `mesh='spherical'` (this is the default option) to signal that all longitudes and latitudes are in degrees.\n",
64+
"To create a `parcels.FieldSet` object, we define the `parcels.Field`s and the structured grid (`parcels.XGrid`) the fields are defined on. We add the argument `mesh='spherical'` to the `parcels.XGrid` to signal that all longitudes and latitudes are in degrees.\n",
65+
"\n",
66+
"```{note}\n",
67+
"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",
68+
"```\n",
6569
"\n",
66-
"Plotting the `U` field indeed shows a uniform 1m/s eastward flow.\n"
70+
"Plotting the `U` field indeed shows a uniform 1 m/s eastward flow.\n"
6771
]
6872
},
6973
{
@@ -72,8 +76,25 @@
7276
"metadata": {},
7377
"outputs": [],
7478
"source": [
75-
"fieldset = parcels.FieldSet.from_data(data, dims, mesh=\"spherical\")\n",
76-
"plt.pcolormesh(fieldset.U.lon, fieldset.U.lat, fieldset.U.data[0, :, :], vmin=0, vmax=1)\n",
79+
"grid = parcels.XGrid.from_dataset(ds, mesh=\"spherical\")\n",
80+
"U = parcels.Field(\"U\", ds[\"U\"], grid, interp_method=parcels.interpolators.XLinear)\n",
81+
"V = parcels.Field(\"V\", ds[\"V\"], grid, interp_method=parcels.interpolators.XLinear)\n",
82+
"UV = parcels.VectorField(\"UV\", U, V)\n",
83+
"temperature = parcels.Field(\n",
84+
" \"temperature\", ds[\"temperature\"], grid, interp_method=parcels.interpolators.XLinear\n",
85+
")\n",
86+
"fieldset = parcels.FieldSet([U, V, UV, temperature])\n",
87+
"\n",
88+
"plt.pcolormesh(\n",
89+
" fieldset.U.grid.lon,\n",
90+
" fieldset.U.grid.lat,\n",
91+
" fieldset.U.data[0, 0, :, :],\n",
92+
" vmin=0,\n",
93+
" vmax=1,\n",
94+
" shading=\"gouraud\",\n",
95+
")\n",
96+
"plt.ylabel(\"Latitude\")\n",
97+
"plt.xlabel(\"Longitude\")\n",
7798
"plt.colorbar()\n",
7899
"plt.show()"
79100
]
@@ -92,8 +113,12 @@
92113
"metadata": {},
93114
"outputs": [],
94115
"source": [
95-
"print(fieldset.UV[0, 0, 40, -5])\n",
96-
"print(fieldset.temp[0, 0, 40, -5])"
116+
"time = np.array([0])\n",
117+
"z = np.array([0])\n",
118+
"lat = np.array([40])\n",
119+
"lon = np.array([-5])\n",
120+
"print(fieldset.UV[time, z, lat, lon])\n",
121+
"print(fieldset.temperature[time, z, lat, lon])"
97122
]
98123
},
99124
{
@@ -103,7 +128,7 @@
103128
"source": [
104129
"While the temperature field indeed is 20C, as we defined, these printed velocities are much smaller.\n",
105130
"\n",
106-
"This is because Parcels converts under the hood from m/s to degrees/s. This conversion is done with a `UnitConverter` object, which is stored in the `.units` attribute of each Field. Below, we print these\n"
131+
"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"
107132
]
108133
},
109134
{
@@ -112,7 +137,7 @@
112137
"metadata": {},
113138
"outputs": [],
114139
"source": [
115-
"for fld in [fieldset.U, fieldset.V, fieldset.temp]:\n",
140+
"for fld in [fieldset.U, fieldset.V, fieldset.temperature]:\n",
116141
" print(f\"{fld.name}: {fld.units}\")"
117142
]
118143
},
@@ -121,7 +146,7 @@
121146
"cell_type": "markdown",
122147
"metadata": {},
123148
"source": [
124-
"So the U field has a `GeographicPolar` UnitConverter object, the V field has a `Geographic` Unitconverter and the `temp` field has a `UnitConverter` object.\n",
149+
"So the U field has a `GeographicPolar` UnitConverter object, the V field has a `Geographic` UnitConverter and the `temp` field has a `UnitConverter` object.\n",
125150
"\n",
126151
"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"
127152
]
@@ -132,7 +157,7 @@
132157
"metadata": {},
133158
"outputs": [],
134159
"source": [
135-
"u, v = fieldset.UV[0, 0, 40, -5]\n",
160+
"u, v = fieldset.UV[time, z, lat, lon]\n",
136161
"print(v * 1852 * 60)"
137162
]
138163
},
@@ -150,7 +175,15 @@
150175
"metadata": {},
151176
"outputs": [],
152177
"source": [
153-
"print(fieldset.UV.eval(0, 0, 40, -5, applyConversion=False))"
178+
"print(\n",
179+
" fieldset.UV.eval(\n",
180+
" time,\n",
181+
" z,\n",
182+
" lat,\n",
183+
" lon,\n",
184+
" applyConversion=False,\n",
185+
" )\n",
186+
")"
154187
]
155188
},
156189
{
@@ -166,7 +199,7 @@
166199
"cell_type": "markdown",
167200
"metadata": {},
168201
"source": [
169-
"If longitudes and latitudes are given in meters, rather than degrees, simply add `mesh='flat'` when creating the FieldSet object.\n"
202+
"If longitudes and latitudes are given in meters, rather than degrees, simply add `mesh='flat'` when creating the XGrid object.\n"
170203
]
171204
},
172205
{
@@ -175,20 +208,38 @@
175208
"metadata": {},
176209
"outputs": [],
177210
"source": [
178-
"fieldset_flat = parcels.FieldSet.from_data(data, dims, mesh=\"flat\")\n",
211+
"ds_flat = simple_UV_dataset(dims=(1, 1, nlat, nlon), mesh=\"flat\").isel(time=0, depth=0)\n",
212+
"ds_flat[\"temperature\"] = ds_flat[\"U\"] + 20 # add temperature field of 20 deg\n",
213+
"ds_flat[\"U\"].data[:] = 1.0 # set U to 1 m/s\n",
214+
"ds_flat[\"V\"].data[:] = 1.0 # set V to 1 m/s\n",
215+
"grid = parcels.XGrid.from_dataset(ds_flat, mesh=\"flat\")\n",
216+
"U = parcels.Field(\"U\", ds_flat[\"U\"], grid, interp_method=parcels.interpolators.XLinear)\n",
217+
"V = parcels.Field(\"V\", ds_flat[\"V\"], grid, interp_method=parcels.interpolators.XLinear)\n",
218+
"UV = parcels.VectorField(\"UV\", U, V)\n",
219+
"temperature = parcels.Field(\n",
220+
" \"temperature\",\n",
221+
" ds_flat[\"temperature\"],\n",
222+
" grid,\n",
223+
" interp_method=parcels.interpolators.XLinear,\n",
224+
")\n",
225+
"fieldset_flat = parcels.FieldSet([U, V, UV, temperature])\n",
179226
"\n",
180227
"plt.pcolormesh(\n",
181-
" fieldset_flat.U.lon,\n",
182-
" fieldset_flat.U.lat,\n",
183-
" fieldset_flat.U.data[0, :, :],\n",
228+
" fieldset_flat.U.grid.lon,\n",
229+
" fieldset_flat.U.grid.lat,\n",
230+
" fieldset_flat.U.data[0, 0, :, :],\n",
184231
" vmin=0,\n",
185232
" vmax=1,\n",
233+
" shading=\"gouraud\",\n",
186234
")\n",
187235
"plt.colorbar()\n",
188236
"plt.show()\n",
189237
"\n",
190-
"print(\"Velocities:\", fieldset_flat.UV[0, 0, 40, -5])\n",
191-
"for fld in [fieldset_flat.U, fieldset_flat.V, fieldset_flat.temp]:\n",
238+
"print(\n",
239+
" \"Velocities:\",\n",
240+
" fieldset_flat.UV[time, z, lat, lon],\n",
241+
")\n",
242+
"for fld in [fieldset_flat.U, fieldset_flat.V, fieldset_flat.temperature]:\n",
192243
" print(f\"{fld.name}: {fld.units}\")"
193244
]
194245
},
@@ -225,23 +276,35 @@
225276
"kh_zonal = 100 # in m^2/s\n",
226277
"kh_meridional = 100 # in m^2/s\n",
227278
"\n",
228-
"fieldset.add_field(\n",
229-
" parcels.Field(\n",
230-
" \"Kh_zonal\",\n",
231-
" kh_zonal * np.ones((ydim, xdim), dtype=np.float32),\n",
232-
" grid=fieldset.U.grid,\n",
233-
" )\n",
279+
"ds[\"Kh_zonal\"] = xr.DataArray(\n",
280+
" data=kh_zonal * np.ones((nlat, nlon), dtype=np.float32), dims=[\"YG\", \"XG\"]\n",
234281
")\n",
235-
"fieldset.add_field(\n",
236-
" parcels.Field(\n",
237-
" \"Kh_meridional\",\n",
238-
" kh_meridional * np.ones((ydim, xdim), dtype=np.float32),\n",
239-
" grid=fieldset.U.grid,\n",
240-
" )\n",
282+
"\n",
283+
"kh_zonal_field = parcels.Field(\n",
284+
" \"Kh_zonal\",\n",
285+
" ds[\"Kh_zonal\"],\n",
286+
" grid=fieldset.U.grid,\n",
287+
" interp_method=parcels.interpolators.XLinear,\n",
288+
")\n",
289+
"\n",
290+
"fieldset.add_field(kh_zonal_field)\n",
291+
"\n",
292+
"ds[\"Kh_meridional\"] = xr.DataArray(\n",
293+
" data=kh_meridional * np.ones((nlat, nlon), dtype=np.float32), dims=[\"YG\", \"XG\"]\n",
241294
")\n",
242295
"\n",
296+
"kh_meridional_field = parcels.Field(\n",
297+
" \"Kh_meridional\",\n",
298+
" ds[\"Kh_meridional\"],\n",
299+
" grid=grid,\n",
300+
" interp_method=parcels.interpolators.XLinear,\n",
301+
")\n",
302+
"\n",
303+
"fieldset.add_field(kh_meridional_field)\n",
304+
"\n",
243305
"for fld in [fieldset.Kh_zonal, fieldset.Kh_meridional]:\n",
244-
" print(f\"{fld.name}: {fld[0, 0, 40, -5]:e} {fld.units}\")"
306+
" val = fld[time, z, lat, lon]\n",
307+
" print(f\"{fld.name}: {val} {fld.units}\")"
245308
]
246309
},
247310
{
@@ -261,7 +324,7 @@
261324
"outputs": [],
262325
"source": [
263326
"deg_to_m = 1852 * 60\n",
264-
"print(fieldset.Kh_meridional[0, 0, 40, -5] * deg_to_m**2)"
327+
"print(fieldset.Kh_meridional[time, z, lat, lon] * deg_to_m**2)"
265328
]
266329
},
267330
{
@@ -281,11 +344,11 @@
281344
"\n",
282345
"| Field name | Converter object | Conversion for `mesh='spherical'` | Conversion for `mesh='flat'` |\n",
283346
"| ---------------- | ----------------------- | --------------------------------------------------------- | ---------------------------- |\n",
284-
"| 'U' | `GeographicPolar` | $1852 \\cdot 60 \\cdot \\cos(lat \\cdot \\frac{\\pi}{180})$ | 1 |\n",
285-
"| 'V' | `Geographic` | $1852 \\cdot 60$ | 1 |\n",
286-
"| 'Kh_zonal' | `GeographicPolarSquare` | $(1852 \\cdot 60 \\cdot \\cos(lat \\cdot \\frac{\\pi}{180}))^2$ | 1 |\n",
287-
"| 'Kh_meridional' | `GeographicSquare` | $(1852 \\cdot 60)^2$ | 1 |\n",
288-
"| All other fields | `UnitConverter` | 1 | 1 |\n",
347+
"| `\"U\"` | `GeographicPolar` | $1852 \\cdot 60 \\cdot \\cos(lat \\cdot \\frac{\\pi}{180})$ | 1 |\n",
348+
"| `\"V\"` | `Geographic` | $1852 \\cdot 60$ | 1 |\n",
349+
"| `\"Kh_zonal\"` | `GeographicPolarSquare` | $(1852 \\cdot 60 \\cdot \\cos(lat \\cdot \\frac{\\pi}{180}))^2$ | 1 |\n",
350+
"| `\"Kh_meridional\"` | `GeographicSquare` | $(1852 \\cdot 60)^2$ | 1 |\n",
351+
"| All other fields | `UnitConverter` | 1 | 1 |\n",
289352
"\n",
290353
"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",
291354
"\n",
@@ -298,12 +361,19 @@
298361
"metadata": {},
299362
"outputs": [],
300363
"source": [
364+
"ds[\"Ustokes\"] = xr.DataArray(\n",
365+
" data=np.ones((nlat, nlon), dtype=np.float32), dims=[\"YG\", \"XG\"]\n",
366+
")\n",
367+
"\n",
301368
"fieldset.add_field(\n",
302369
" parcels.Field(\n",
303-
" \"Ustokes\", np.ones((ydim, xdim), dtype=np.float32), grid=fieldset.U.grid\n",
370+
" \"Ustokes\",\n",
371+
" ds[\"Ustokes\"],\n",
372+
" grid=fieldset.U.grid,\n",
373+
" interp_method=parcels.interpolators.XLinear,\n",
304374
" )\n",
305375
")\n",
306-
"print(fieldset.Ustokes[0, 0, 40, -5])"
376+
"print(fieldset.Ustokes[time, z, lat, lon])"
307377
]
308378
},
309379
{
@@ -320,17 +390,15 @@
320390
"metadata": {},
321391
"outputs": [],
322392
"source": [
323-
"from parcels.tools.converters import GeographicPolar\n",
324-
"\n",
325-
"fieldset.Ustokes.units = GeographicPolar()\n",
326-
"print(fieldset.Ustokes[0, 0, 40, -5])\n",
327-
"print(fieldset.Ustokes[0, 0, 40, -5] * 1852 * 60 * np.cos(40 * np.pi / 180))"
393+
"fieldset.Ustokes.units = parcels.GeographicPolar()\n",
394+
"print(fieldset.Ustokes[time, z, lat, lon])\n",
395+
"print(fieldset.Ustokes[time, z, lat, lon] * 1852 * 60 * np.cos(40 * np.pi / 180))"
328396
]
329397
}
330398
],
331399
"metadata": {
332400
"kernelspec": {
333-
"display_name": "Python 3",
401+
"display_name": "test-notebooks",
334402
"language": "python",
335403
"name": "python3"
336404
},
@@ -344,7 +412,7 @@
344412
"name": "python",
345413
"nbconvert_exporter": "python",
346414
"pygments_lexer": "ipython3",
347-
"version": "3.11.6"
415+
"version": "3.11.0"
348416
}
349417
},
350418
"nbformat": 4,

‎docs/user_guide/index.md‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -24,6 +24,7 @@ The tutorials written for Parcels v3 are currently being updated for Parcels v4.
2424
:titlesonly:
2525
2626
examples/tutorial_nemo_curvilinear.ipynb
27+
examples/tutorial_unitconverters.ipynb
2728
```
2829

2930
<!-- examples/documentation_indexing.ipynb -->
@@ -32,7 +33,6 @@ examples/tutorial_nemo_curvilinear.ipynb
3233
<!-- examples/tutorial_timevaryingdepthdimensions.ipynb -->
3334
<!-- examples/tutorial_periodic_boundaries.ipynb -->
3435
<!-- examples/tutorial_interpolation.ipynb -->
35-
<!-- examples/tutorial_unitconverters.ipynb -->
3636

3737
```{toctree}
3838
:caption: Create ParticleSets

0 commit comments

Comments
 (0)