Skip to content

Commit 9f39bea

Browse files
reint-fischerreint-fischer
authored
Merge pull request #2398 from Parcels-code/curvilinear_guide_v4
Add curvilinear how-to guide v4
2 parents e4da651 + 5b6a731 commit 9f39bea

File tree

3 files changed

+317
-218
lines changed

3 files changed

+317
-218
lines changed

‎docs/user_guide/examples/tutorial_nemo_curvilinear.ipynb‎

Lines changed: 316 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,316 @@
1+
{
2+
"cells": [
3+
{
4+
"attachments": {},
5+
"cell_type": "markdown",
6+
"metadata": {},
7+
"source": [
8+
"# Curvilinear grids\n"
9+
]
10+
},
11+
{
12+
"attachments": {},
13+
"cell_type": "markdown",
14+
"metadata": {},
15+
"source": [
16+
"Parcels supports [curvilinear grids](https://www.nemo-ocean.eu/doc/node108.html) such as those used in the [NEMO models](https://www.nemo-ocean.eu/).\n",
17+
"\n",
18+
"We will be using the example dataset `NemoCurvilinear_data`. These fields are a purely zonal flow on an aqua-planet (so zonal-velocity is 1 m/s and meridional-velocity is 0 m/s everywhere, and no land). However, because of the curvilinear grid, the `U` and `V` fields vary for the rotated gridcells north of 20N.\n"
19+
]
20+
},
21+
{
22+
"cell_type": "code",
23+
"execution_count": null,
24+
"metadata": {},
25+
"outputs": [],
26+
"source": [
27+
"import matplotlib.pyplot as plt\n",
28+
"import numpy as np\n",
29+
"import xarray as xr\n",
30+
"\n",
31+
"import parcels"
32+
]
33+
},
34+
{
35+
"attachments": {},
36+
"cell_type": "markdown",
37+
"metadata": {},
38+
"source": [
39+
"We can create a `FieldSet` just like we do for normal grids.\n",
40+
"Note that NEMO is discretised on a C-grid. U and V velocities are not located on the same nodes (see https://www.nemo-ocean.eu/doc/node19.html ).\n",
41+
"\n",
42+
"```\n",
43+
" __V1__\n",
44+
"| |\n",
45+
"U0 U1\n",
46+
"|__V0__|\n",
47+
"```\n",
48+
"\n",
49+
"To interpolate U, V velocities on the C-grid, Parcels needs to read the f-nodes, which are located on the corners of the cells (for indexing details: https://www.nemo-ocean.eu/doc/img360.png).\n",
50+
"\n",
51+
"```{note}\n",
52+
"TODO: add link to grid indexing explanation once implemented in v4\n",
53+
"```\n"
54+
]
55+
},
56+
{
57+
"cell_type": "code",
58+
"execution_count": null,
59+
"metadata": {},
60+
"outputs": [],
61+
"source": [
62+
"data_folder = parcels.download_example_dataset(\"NemoCurvilinear_data\")\n",
63+
"files = data_folder.glob(\"*.nc4\")\n",
64+
"ds_fields = xr.open_mfdataset(\n",
65+
" files, combine=\"nested\", data_vars=\"minimal\", coords=\"minimal\", compat=\"override\"\n",
66+
")\n",
67+
"\n",
68+
"\n",
69+
"# TODO: replace manual fieldset creation with FieldSet.from_nemo() once available\n",
70+
"ds_fields = (\n",
71+
" ds_fields.isel(time_counter=0, drop=True)\n",
72+
" .isel(time=0, drop=True)\n",
73+
" .isel(z_a=0, drop=True)\n",
74+
" .rename({\"glamf\": \"lon\", \"gphif\": \"lat\", \"z\": \"depth\"})\n",
75+
")\n",
76+
"\n",
77+
"import xgcm\n",
78+
"\n",
79+
"xgcm_grid = xgcm.Grid(\n",
80+
" ds_fields,\n",
81+
" coords={\n",
82+
" \"X\": {\"left\": \"x\"},\n",
83+
" \"Y\": {\"left\": \"y\"},\n",
84+
" },\n",
85+
" periodic=False,\n",
86+
" autoparse_metadata=False,\n",
87+
")\n",
88+
"grid = parcels.XGrid(xgcm_grid, mesh=\"spherical\")\n",
89+
"\n",
90+
"U = parcels.Field(\n",
91+
" \"U\", ds_fields[\"U\"], grid, interp_method=parcels.interpolators.XLinear\n",
92+
")\n",
93+
"V = parcels.Field(\n",
94+
" \"V\", ds_fields[\"V\"], grid, interp_method=parcels.interpolators.XLinear\n",
95+
")\n",
96+
"U.units = parcels.GeographicPolar()\n",
97+
"V.units = (\n",
98+
" parcels.GeographicPolar()\n",
99+
") # U and V need GeographicPolar for C-Grid interpolation to work correctly\n",
100+
"UV = parcels.VectorField(\n",
101+
" \"UV\", U, V, vector_interp_method=parcels.interpolators.CGrid_Velocity\n",
102+
")\n",
103+
"fieldset = parcels.FieldSet([U, V, UV])"
104+
]
105+
},
106+
{
107+
"attachments": {},
108+
"cell_type": "markdown",
109+
"metadata": {},
110+
"source": [
111+
"And we can plot the `U` field.\n"
112+
]
113+
},
114+
{
115+
"cell_type": "code",
116+
"execution_count": null,
117+
"metadata": {},
118+
"outputs": [],
119+
"source": [
120+
"fig, ax = plt.subplots(1, 2, figsize=(10, 4))\n",
121+
"pc1 = ds_fields.U.plot(cmap=\"viridis\", ax=ax[0], vmin=0)\n",
122+
"pc2 = ax[1].pcolormesh(\n",
123+
" fieldset.U.grid.lon,\n",
124+
" fieldset.U.grid.lat,\n",
125+
" fieldset.U.data[0, 0, :, :],\n",
126+
" vmin=0,\n",
127+
" vmax=1,\n",
128+
")\n",
129+
"ax[1].set_ylabel(\"Latitude [deg N]\")\n",
130+
"ax[1].set_xlabel(\"Longitude [deg E]\")\n",
131+
"plt.colorbar(pc2, label=\"U\", ax=ax[1])\n",
132+
"plt.tight_layout()\n",
133+
"plt.show()"
134+
]
135+
},
136+
{
137+
"attachments": {},
138+
"cell_type": "markdown",
139+
"metadata": {},
140+
"source": [
141+
"As you see above, the `U` field indeed is 1 m/s south of 20N, but varies with longitude and latitude north of that. We can confirm that Parcels will take care to rotate the `U` and `V` fields by doing a field evaluation at (60N, 50E):"
142+
]
143+
},
144+
{
145+
"cell_type": "code",
146+
"execution_count": null,
147+
"metadata": {},
148+
"outputs": [],
149+
"source": [
150+
"u, v = fieldset.UV.eval(\n",
151+
" np.array([0]), np.array([0]), np.array([20]), np.array([50]), applyConversion=False\n",
152+
") # do not convert m/s to deg/s\n",
153+
"print(f\"(u, v) = ({u[0]:.3f}, {v[0]:.3f})\")\n",
154+
"assert np.isclose(u, 1.0, atol=1e-3)"
155+
]
156+
},
157+
{
158+
"cell_type": "markdown",
159+
"metadata": {},
160+
"source": [
161+
"Now we can run particles as normal."
162+
]
163+
},
164+
{
165+
"cell_type": "code",
166+
"execution_count": null,
167+
"metadata": {
168+
"tags": [
169+
"hide-output"
170+
]
171+
},
172+
"outputs": [],
173+
"source": [
174+
"npart = 20\n",
175+
"lonp = 30 * np.ones(npart)\n",
176+
"latp = np.linspace(-70, 88, npart)\n",
177+
"runtime = np.timedelta64(40, \"D\")\n",
178+
"\n",
179+
"pset = parcels.ParticleSet(fieldset, lon=lonp, lat=latp)\n",
180+
"pfile = parcels.ParticleFile(\n",
181+
" store=\"output_curvilinear.zarr\", outputdt=np.timedelta64(1, \"D\")\n",
182+
")\n",
183+
"\n",
184+
"pset.execute(\n",
185+
" [parcels.kernels.AdvectionEE],\n",
186+
" runtime=runtime,\n",
187+
" dt=np.timedelta64(1, \"D\"),\n",
188+
" output_file=pfile,\n",
189+
")\n",
190+
"np.testing.assert_allclose(pset.lat, latp, atol=1e-1)"
191+
]
192+
},
193+
{
194+
"attachments": {},
195+
"cell_type": "markdown",
196+
"metadata": {},
197+
"source": [
198+
"And then we can plot these trajectories. As expected, all trajectories go exactly zonal and due to the curvature of the earth, ones at higher latitude move more degrees eastward (even though the distance in km is equal for all particles). These particles at high latitudes cross the antimeridian (180 deg E) and keep going east.\n"
199+
]
200+
},
201+
{
202+
"cell_type": "code",
203+
"execution_count": null,
204+
"metadata": {},
205+
"outputs": [],
206+
"source": [
207+
"ds = xr.open_zarr(\"output_curvilinear.zarr\")\n",
208+
"\n",
209+
"plt.plot(ds.lon.T, ds.lat.T, \".-\")\n",
210+
"plt.vlines(np.arange(-180, 901, 360), -90, 90, color=\"r\", label=\"antimeridian\")\n",
211+
"plt.ylabel(\"Latitude [deg N]\")\n",
212+
"plt.xlabel(\"Longitude [deg E]\")\n",
213+
"plt.xticks(np.arange(-180, 901, 90))\n",
214+
"plt.legend(loc=\"lower right\")\n",
215+
"plt.show()"
216+
]
217+
},
218+
{
219+
"cell_type": "markdown",
220+
"metadata": {},
221+
"source": [
222+
"If we want the `particles.lon` to stay within `[-180,180]` (or `[0,360]`), we can either do this in post-processing, or add a periodic boundary Kernel:"
223+
]
224+
},
225+
{
226+
"cell_type": "code",
227+
"execution_count": null,
228+
"metadata": {},
229+
"outputs": [],
230+
"source": [
231+
"# post processing\n",
232+
"ds[\"lon\"] = ds[\"lon\"] % 360\n",
233+
"ds[\"lon\"] = ds[\"lon\"].where(ds[\"lon\"] <= 180, ds[\"lon\"] - 360)"
234+
]
235+
},
236+
{
237+
"cell_type": "code",
238+
"execution_count": null,
239+
"metadata": {
240+
"tags": [
241+
"hide-output"
242+
]
243+
},
244+
"outputs": [],
245+
"source": [
246+
"# with a Kernel\n",
247+
"def periodicBC(particles, fieldset): # pragma: no cover\n",
248+
" particles.dlon = np.where(\n",
249+
" particles.lon + particles.dlon > 180, particles.dlon - 360, particles.dlon\n",
250+
" )\n",
251+
"\n",
252+
"\n",
253+
"pset = parcels.ParticleSet(fieldset, lon=lonp, lat=latp)\n",
254+
"pfile = parcels.ParticleFile(\n",
255+
" store=\"output_curvilinear_periodic.zarr\", outputdt=np.timedelta64(1, \"D\")\n",
256+
")\n",
257+
"\n",
258+
"pset.execute(\n",
259+
" [parcels.kernels.AdvectionEE, periodicBC],\n",
260+
" runtime=runtime,\n",
261+
" dt=np.timedelta64(1, \"D\"),\n",
262+
" output_file=pfile,\n",
263+
")"
264+
]
265+
},
266+
{
267+
"cell_type": "code",
268+
"execution_count": null,
269+
"metadata": {},
270+
"outputs": [],
271+
"source": [
272+
"ds_periodic = xr.open_zarr(\"output_curvilinear_periodic.zarr\")\n",
273+
"\n",
274+
"fig, ax = plt.subplots(1, 2, figsize=(10, 5))\n",
275+
"ax[0].plot(ds.lon.T, ds.lat.T, \".-\")\n",
276+
"ax[0].vlines(np.arange(-180, 360, 360), -90, 90, color=\"r\", label=\"antimeridian\")\n",
277+
"ax[0].set_ylabel(\"Latitude [deg N]\")\n",
278+
"ax[0].set_xlabel(\"Longitude [deg E]\")\n",
279+
"ax[0].set_xticks(np.arange(-180, 181, 45))\n",
280+
"ax[0].set_title(\"in post processing\")\n",
281+
"ax[0].legend(loc=\"lower center\")\n",
282+
"\n",
283+
"ax[1].plot(ds_periodic.lon.T, ds_periodic.lat.T, \".-\")\n",
284+
"ax[1].vlines(np.arange(-180, 360, 360), -90, 90, color=\"r\", label=\"antimeridian\")\n",
285+
"ax[1].set_ylabel(\"Latitude [deg N]\")\n",
286+
"ax[1].set_xlabel(\"Longitude [deg E]\")\n",
287+
"ax[1].set_xticks(np.arange(-180, 181, 45))\n",
288+
"ax[1].set_title(\"with periodic Kernel\")\n",
289+
"ax[1].legend(loc=\"lower center\")\n",
290+
"\n",
291+
"plt.show()"
292+
]
293+
}
294+
],
295+
"metadata": {
296+
"kernelspec": {
297+
"display_name": "test-notebooks",
298+
"language": "python",
299+
"name": "python3"
300+
},
301+
"language_info": {
302+
"codemirror_mode": {
303+
"name": "ipython",
304+
"version": 3
305+
},
306+
"file_extension": ".py",
307+
"mimetype": "text/x-python",
308+
"name": "python",
309+
"nbconvert_exporter": "python",
310+
"pygments_lexer": "ipython3",
311+
"version": "3.11.0"
312+
}
313+
},
314+
"nbformat": 4,
315+
"nbformat_minor": 4
316+
}

0 commit comments

Comments
 (0)