Skip to content

Commit 27b8dfc

Browse files
reint-fischerreint-fischer
reint-fischer
authored and
reint-fischer
committed
move tutorial to user guide
1 parent ba56c6b commit 27b8dfc

3 files changed

Lines changed: 250 additions & 895 deletions

File tree

‎docs/user_guide/examples/tutorial_interpolation.ipynb‎

Lines changed: 250 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,250 @@
1+
{
2+
"cells": [
3+
{
4+
"attachments": {},
5+
"cell_type": "markdown",
6+
"metadata": {},
7+
"source": [
8+
"# Configure Field Interpolation\n"
9+
]
10+
},
11+
{
12+
"attachments": {},
13+
"cell_type": "markdown",
14+
"metadata": {},
15+
"source": [
16+
"Parcels comes with a number of different interpolation methods for tracer fields, such as temperature. Here, we will show how these work, in an idealised example. For more guidance on the sampling of such fields, check out the [sampling guide](./tutorial_sampling).\n",
17+
"\n",
18+
"We first import the relevant modules"
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+
"from matplotlib import cm\n",
30+
"\n",
31+
"import parcels"
32+
]
33+
},
34+
{
35+
"attachments": {},
36+
"cell_type": "markdown",
37+
"metadata": {},
38+
"source": [
39+
"We create a small 2D grid where `P` is a tracer that we want to interpolate. In each grid cell, `P` has a random value between 0.1 and 1.1. We then set `P[1,1]` to `0`, which for Parcels specifies that this is a land cell\n"
40+
]
41+
},
42+
{
43+
"cell_type": "code",
44+
"execution_count": null,
45+
"metadata": {},
46+
"outputs": [],
47+
"source": [
48+
"from parcels._datasets.structured.generated import simple_UV_dataset\n",
49+
"\n",
50+
"ds = simple_UV_dataset(dims=(1, 1, 5, 4), mesh=\"flat\").isel(time=0, depth=0)\n",
51+
"ds[\"lat\"][:] = np.linspace(0.0, 1.0, len(ds.YG))\n",
52+
"ds[\"lon\"][:] = np.linspace(0.0, 1.0, len(ds.XG))\n",
53+
"dx, dy = 1.0 / len(ds.XG), 1.0 / len(ds.YG)\n",
54+
"ds[\"P\"] = ds[\"U\"] + np.random.rand(5, 4) + 0.1\n",
55+
"ds[\"P\"][1, 1] = 0\n",
56+
"ds"
57+
]
58+
},
59+
{
60+
"cell_type": "markdown",
61+
"metadata": {},
62+
"source": [
63+
"From this dataset we create a `parcels.FieldSet`. Parcels requires an interpolation method to be set for each `parcels.Field`, which we will later adapt to see the effects of the different interpolators. A common interpolator for fields on structured grids is (tri)linear, implemented in `parcels.interpolators.XLinear`."
64+
]
65+
},
66+
{
67+
"cell_type": "code",
68+
"execution_count": null,
69+
"metadata": {},
70+
"outputs": [],
71+
"source": [
72+
"grid = parcels.XGrid.from_dataset(ds, mesh=\"flat\")\n",
73+
"U = parcels.Field(\"U\", ds[\"U\"], grid, interp_method=parcels.interpolators.XLinear)\n",
74+
"V = parcels.Field(\"V\", ds[\"V\"], grid, interp_method=parcels.interpolators.XLinear)\n",
75+
"UV = parcels.VectorField(\"UV\", U, V)\n",
76+
"P = parcels.Field(\"P\", ds[\"P\"], grid, interp_method=parcels.interpolators.XLinear)\n",
77+
"fieldset = parcels.FieldSet([U, V, UV, P])"
78+
]
79+
},
80+
{
81+
"attachments": {},
82+
"cell_type": "markdown",
83+
"metadata": {},
84+
"source": [
85+
"We create a Particle class that can sample this field `P`."
86+
]
87+
},
88+
{
89+
"cell_type": "code",
90+
"execution_count": null,
91+
"metadata": {},
92+
"outputs": [],
93+
"source": [
94+
"SampleParticle = parcels.Particle.add_variable(\n",
95+
" parcels.Variable(\"p\", dtype=np.float32, initial=np.nan)\n",
96+
")\n",
97+
"\n",
98+
"\n",
99+
"def SampleP(particles, fieldset):\n",
100+
" particles.p = fieldset.P[particles]"
101+
]
102+
},
103+
{
104+
"attachments": {},
105+
"cell_type": "markdown",
106+
"metadata": {},
107+
"source": [
108+
"Now, we perform four different interpolations on `P`, which we can control by setting `fieldset.P.interp_method`. Note that this can always be done _after_ the `FieldSet` creation. We store the results of each interpolation method in an entry in the dictionary `pset`."
109+
]
110+
},
111+
{
112+
"cell_type": "code",
113+
"execution_count": null,
114+
"metadata": {},
115+
"outputs": [],
116+
"source": [
117+
"pset = {}\n",
118+
"from parcels.interpolators import (\n",
119+
" CGrid_Tracer,\n",
120+
" XLinear,\n",
121+
" XLinearInvdistLandTracer,\n",
122+
" XNearest,\n",
123+
")\n",
124+
"\n",
125+
"interp_methods = [XLinear, XLinearInvdistLandTracer, XNearest, CGrid_Tracer]\n",
126+
"for p_interp in interp_methods:\n",
127+
" fieldset.P.interp_method = (\n",
128+
" p_interp # setting the interpolation method for fieldset.P\n",
129+
" )\n",
130+
"\n",
131+
" print(fieldset.P.interp_method.__name__)\n",
132+
" xv, yv = np.meshgrid(np.linspace(0, 1, 8), np.linspace(0, 1, 8))\n",
133+
" pset[p_interp.__name__] = parcels.ParticleSet(\n",
134+
" fieldset, pclass=SampleParticle, lon=xv.flatten(), lat=yv.flatten()\n",
135+
" )\n",
136+
" pset[p_interp.__name__].execute(\n",
137+
" SampleP, runtime=np.timedelta64(1, \"D\"), dt=np.timedelta64(1, \"D\")\n",
138+
" )"
139+
]
140+
},
141+
{
142+
"attachments": {},
143+
"cell_type": "markdown",
144+
"metadata": {},
145+
"source": [
146+
"And then we can show each of the four interpolation methods, by plotting the interpolated values on the `Particles` locations (circles) on top of the `Field` values (background colors)\n"
147+
]
148+
},
149+
{
150+
"cell_type": "code",
151+
"execution_count": null,
152+
"metadata": {},
153+
"outputs": [],
154+
"source": [
155+
"fig, ax = plt.subplots(1, 4, figsize=(18, 6))\n",
156+
"for i, p in enumerate(pset.keys()):\n",
157+
" data = np.copy(fieldset.P.data[0, 0, :, :])\n",
158+
" data[1, 1] = np.nan\n",
159+
" x = np.linspace(-dx / 2, 1 + dx / 2, len(ds.XG) + 1)\n",
160+
" y = np.linspace(-dy / 2, 1 + dy / 2, len(ds.YG) + 1)\n",
161+
" if p == \"CGrid_Tracer\":\n",
162+
" for lat in fieldset.P.grid.lat:\n",
163+
" ax[i].axhline(lat, color=\"k\", linestyle=\"--\")\n",
164+
" for lon in fieldset.P.grid.lon:\n",
165+
" ax[i].axvline(lon, color=\"k\", linestyle=\"--\")\n",
166+
" pc = ax[i].pcolormesh(x, y, data, vmin=0.1, vmax=1.1)\n",
167+
" ax[i].scatter(\n",
168+
" pset[p].lon, pset[p].lat, c=pset[p].p, edgecolors=\"k\", s=50, vmin=0.1, vmax=1.1\n",
169+
" )\n",
170+
" xp, yp = np.meshgrid(fieldset.P.grid.lon, fieldset.P.grid.lat)\n",
171+
" ax[i].plot(xp, yp, \"kx\")\n",
172+
" ax[i].set_title(f\"Using interp_method='{p}'\")\n",
173+
"plt.colorbar(pc, ax=ax, orientation=\"horizontal\", shrink=0.5)\n",
174+
"plt.show()"
175+
]
176+
},
177+
{
178+
"attachments": {},
179+
"cell_type": "markdown",
180+
"metadata": {},
181+
"source": [
182+
"The white box is here the 'land' point where the tracer is set to zero and the crosses are the locations of the grid points. As you see, the interpolated value is always equal to the field value if the particle is exactly on the grid point (circles on crosses).\n",
183+
"\n",
184+
"For `interp_method='XNearest'`, the particle values are the same for all particles in a grid cell. They are also the same for `interp_method='CGrid_Tracer'`, but the grid cells have then shifted. That is because in a C-grid, the tracer grid cell is on the top-right corner (black dashed lines in right-most panel).\n",
185+
"\n",
186+
"For `interp_method='XLinearInvdistLandTracer'`, we see that values are the same as `interp_method='XLinear'` for grid cells that don't border the land point. For grid cells that do border the land cell, the `XLinearInvdistLandTracer` interpolation method gives higher values, as also shown in the difference plot below.\n"
187+
]
188+
},
189+
{
190+
"cell_type": "code",
191+
"execution_count": null,
192+
"metadata": {},
193+
"outputs": [],
194+
"source": [
195+
"plt.scatter(\n",
196+
" pset[\"XLinear\"].lon,\n",
197+
" pset[\"XLinear\"].lat,\n",
198+
" c=pset[\"XLinearInvdistLandTracer\"].p - pset[\"XLinear\"].p,\n",
199+
" edgecolors=\"k\",\n",
200+
" s=50,\n",
201+
" cmap=cm.bwr,\n",
202+
" vmin=-0.05,\n",
203+
" vmax=0.05,\n",
204+
")\n",
205+
"xp, yp = np.meshgrid(fieldset.P.grid.lon, fieldset.P.grid.lat)\n",
206+
"plt.plot(xp, yp, \"kx\")\n",
207+
"plt.colorbar()\n",
208+
"plt.title(\n",
209+
" \"Difference between 'interp_method=XLinear' \"\n",
210+
" \"and 'interp_method=XLinearInvdistLandTracer'\"\n",
211+
")\n",
212+
"plt.show()"
213+
]
214+
},
215+
{
216+
"attachments": {},
217+
"cell_type": "markdown",
218+
"metadata": {},
219+
"source": [
220+
"So in summary, Parcels has four different interpolation schemes for tracers:\n",
221+
"\n",
222+
"1. `interp_method=parcels.interpolators.XLinear`: compute linear interpolation\n",
223+
"2. `interp_method=parcels.interpolators.XLinearInvdistLandTracer`: compute linear interpolation except near land (where field value is zero). In that case, inverse distance weighting interpolation is computed, weighting by squares of the distance.\n",
224+
"3. `interp_method=parcels.interpolators.XNearest`: return nearest field value\n",
225+
"4. `interp_method=parcels.interpolators.CGridTracer`: return nearest field value supposing C cells\n"
226+
]
227+
}
228+
],
229+
"metadata": {
230+
"kernelspec": {
231+
"display_name": "test-notebooks",
232+
"language": "python",
233+
"name": "python3"
234+
},
235+
"language_info": {
236+
"codemirror_mode": {
237+
"name": "ipython",
238+
"version": 3
239+
},
240+
"file_extension": ".py",
241+
"mimetype": "text/x-python",
242+
"name": "python",
243+
"nbconvert_exporter": "python",
244+
"pygments_lexer": "ipython3",
245+
"version": "3.11.0"
246+
}
247+
},
248+
"nbformat": 4,
249+
"nbformat_minor": 4
250+
}

0 commit comments

Comments
 (0)