Skip to content
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

Address additional comments for regridding section #12

Merged
merged 3 commits into from
Jan 30, 2023
Merged

Address additional comments for regridding section #12

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
3b3e46e
address additional comments
stephenworsley Jan 25, 2023
418d258
jupytext sync
stephenworsley Jan 25, 2023
aa93e71
remove unnecessary imports
stephenworsley Jan 26, 2023
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
89 changes: 43 additions & 46 deletions notebooks/.notebook_shadow_copies/Sec_05_Regridding.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -17,12 +17,7 @@ jupyter:
## Set up

```python
from pathlib import Path
import numpy as np

from esmf_regrid.experimental.unstructured_scheme import MeshToGridESMFRegridder, GridToMeshESMFRegridder
import iris
from iris import load, load_cube
from iris.coords import DimCoord
from iris.cube import Cube
```
Expand Down Expand Up @@ -104,6 +99,7 @@ import matplotlib.pyplot as plt
from testdata_fetching import um_temp
grid_temp = um_temp()

# We slice the cube to make sure it is 2D for plotting.
iqplt.pcolormesh(grid_temp[0, 0])
plt.gca().coastlines()
plt.show()
Expand All @@ -113,11 +109,12 @@ plt.gca().coastlines()
plt.show()
```

We can then plot the difference between the UM data and the data regridded from LFRic. Since our data is now on a latlon grid we can do this with matplotlib as normal.
We can then plot the difference between the UM data and the data regridded from LFRic. Since all our data is now on a latlon grid we can subtract to find the difference between the regridded LFRic data and equivalent UM data and plot this with matplotlib as normal.

```python
temp_diff = result_2 - grid_temp

# We choose a colormap that makes it clear where the differences are.
iqplt.pcolormesh(temp_diff[0, 0], vmin=-4,vmax=4, cmap="seismic")
plt.gca().coastlines()
plt.show()
Expand Down Expand Up @@ -206,10 +203,10 @@ lat_bands = DimCoord(
lon_full = DimCoord(0, bounds=[[-180, 180]], standard_name="longitude", units="degrees")
```

**Step 3:** Create a single celled cube (i.e. `Cube([[0]])`) and attach the latitude and longitude coordinates to it.
**Step 3:** Create a six celled cube (i.e. `Cube([[0, 0, 0, 0, 0, 0]])`) and attach the latitude and longitude coordinates to it.

```python
lat_band_cube = Cube(np.zeros((1,) + lat_bands.shape))
lat_band_cube = Cube([[0, 0, 0, 0, 0, 0]])
lat_band_cube.add_dim_coord(lat_bands, 1)
lat_band_cube.add_dim_coord(lon_full, 0)
lat_band_cube
Expand All @@ -227,7 +224,7 @@ Initialise a `MeshToGridESMFRegridder` with `mesh_cube` and your single celled c
lat_band_mean_calculator_10 = MeshToGridESMFRegridder(mesh_cube, lat_band_cube, resolution=10)
```

**Step 5:** Apply this regridder to `mesh_cube` and print the data from this result (i.e. `print(result_cube.data)`).
**Step 5:** Apply this regridder to `mesh_cube` and print the data from this result (i.e. `print(result_cube.data)`) and plot with `iqplt.pcolormesh`.

```python
lat_band_mean_10 = lat_band_mean_calculator_10(mesh_cube)
Expand All @@ -239,7 +236,7 @@ plt.show()

**Step 6:** Repeat step 4 and 5 for `resolution=100`.

Note the difference in value. Also note that it takes more time to initialise a regridder with higher resolution.
Note the difference in value. Also note that it takes more time to initialise a regridder with higher resolution. Higher resolutions ought to be more accurate but there is a tradeoff between performance and accuracy.

```python
lat_band_mean_calculator_100 = MeshToGridESMFRegridder(mesh_cube, lat_band_cube, resolution=100)
Expand Down Expand Up @@ -274,64 +271,68 @@ print(lon_band_mean_10.data)

## Exercise 3: Hovmoller plots

If we have data on aditional dimensions, we can use the same approach as exercise 2 to produce a Hovmoller diagram. That is, if we have data that varies along time we can take the area weighted mean over latitude bands and plot the data aginst longitude and time (or similarly, we can plot against latitude and time).

**Step 1:** Load a temperature cube using the `testdata_fetching` function `lfric_temp`. Extract a single pressure slice using the `cube.extract` method with a constraint `iris.Constraint(pressure=850)` as the argument (we choose this level because it has noticable details).
If we have data on aditional dimensions, we can use the same approach as exercise 2 to produce a Hovmoller diagram. That is, if we have data that varies along time we can take the area weighted mean over latitude bands and plot the data aginst latitude and time (or similarly, we can plot against longitude and time).

**Step 1:** Load a cube with humidity data using the `testdata_fetching` function `lfric_rh_alltimes_3d`.

from testdata_fetching import fric_temp
mesh_temp = lfric_temp()

temp_slice = mesh_temp.extract(iris.Constraint(pressure=850))
temp_slice
```python
from testdata_fetching import lfric_rh_alltimes_3d

humidity_cube = lfric_rh_alltimes_3d()
humidity_cube
```

**Step 2:** Create a target cube whose longitude coordinate is derived from the UM cube loaded from `um_orography` and whose latitude coordinate has bounds `[[-180, 180]]`. This can be done by slicing a cube derived from `um_orography` (using the slice `[:1]` so that this dimension isnt collapsed), removing the latitude coordinate and adding a latitude coordinate with bounds `[[-180, 180]]` (you can reuse the coordinate from exercise 2).
**Step 2:** Create a target cube whose latitude coordinate is derived from the UM cube loaded from `um_orography` and whose longitude coordinate has bounds `[[-180, 180]]`. This can be done by slicing a cube derived from `um_orography` (using the slice `[:, :1]` so that this dimension isnt collapsed), removing the longitude coordinate and adding a longitude coordinate with bounds `[[-180, 180]]` (you can reuse the coordinate from exercise 2).

```python
target_cube_lons = grid_cube[:1]
target_cube_lons.remove_coord("latitude")
target_cube_lons.add_dim_coord(lat_full, 0)
target_cube_lons
target_cube_lats = grid_cube[:,:1]
target_cube_lats.remove_coord("longitude")
target_cube_lats.add_dim_coord(lon_full, 1)
target_cube_lats
```

```python
# We also can do the same thing for bands of constant latitude.
# We also can do the same thing for bands of constant longitude.

# target_cube_lats = grid_cube[:,:1]
# target_cube_lats.remove_coord("longitude")
# target_cube_lats.add_dim_coord(lon_full, 1)
# target_cube_lats
# target_cube_lons = grid_cube[:1]
# target_cube_lons.remove_coord("latitude")
# target_cube_lons.add_dim_coord(lat_full, 0)
# target_cube_lons
```

**Step 3:** Create a `MeshToGridESMFRegridder` regridder from the slice of the temperature cube onto the target cube. Set the resolution keyword to 2 (this should be sufficient since these are bands of constant longitude). Use this regridder to create a resulting cube.
**Step 3:** Create a `MeshToGridESMFRegridder` regridder from the slice of the humidity cube onto the target cube. Set the resolution keyword to 500 (this should be good balance of accuracy and performance). Use this regridder to create a resulting cube.

```python
um_lon_band_mean_calculator = MeshToGridESMFRegridder(temp_slice, target_cube_lons, resolution=2)
um_lon_bound_means = um_lon_band_mean_calculator(temp_slice)
um_lon_bound_means
um_lat_band_mean_calculator = MeshToGridESMFRegridder(humidity_cube, target_cube_lats, resolution=500)
um_lat_band_means = um_lat_band_mean_calculator(humidity_cube)
um_lat_band_means
```

```python
# um_lat_band_mean_calculator = MeshToGridESMFRegridder(temp_slice, target_cube, resolution=500)
# um_lat_bound_means = um_lat_band_mean_calculator(temp_slice)
# um_lat_bound_means
# Continuing for bands of constant longitude.
# Note: this code takes about 2 minutes to run. I think this is due to with the way ESMF handles cells
# with unusual shapes. See https://github.com/SciTools-incubator/iris-esmf-regrid/issues/234.

# um_lon_band_mean_calculator = MeshToGridESMFRegridder(humidity_cube, target_cube_lons, resolution=2)
# um_lon_band_means = um_lon_band_mean_calculator(humidity_cube)
# um_lon_band_means
```

**Step 4:** Plot the data in the resulting cube. This can be done with `iqplt.pcolormesh`. Note that the resulting cube will have an unnecessary dimension which will have to be sliced (using `[:, 0]`). Note that extra steps can be taken to format the dates for this plot.
**Step 4:** Plot the data in the resulting cube. This can be done with `iqplt.pcolormesh`. Note that the resulting cube will have an unnecessary dimension which will have to be sliced (using `[:, :, 0]`). Note that extra steps can be taken to format the dates for this plot.

```python
import matplotlib.dates as mdates

iqplt.pcolormesh(um_lon_bound_means[:, 0, :])
plt.gca().yaxis.set_major_formatter(mdates.DateFormatter("%D"))
# We use a colormap which highlights fine details.
iqplt.pcolormesh(um_lat_band_means[:, :, 0])
plt.gca().xaxis.set_major_formatter(mdates.DateFormatter("%D"))
plt.show()
```

```python
# import matplotlib.dates as mdates
# iqplt.pcolormesh(um_lat_bound_means[:, :, 0])
# plt.gca().xaxis.set_major_formatter(mdates.DateFormatter("%D"))
# Continuing for bands of constant longitude.

# iqplt.pcolormesh(um_lon_band_means[:, 0])
# plt.gca().yaxis.set_major_formatter(mdates.DateFormatter("%D"))
# plt.show()
```

Expand All @@ -342,7 +343,3 @@ plt.show()
# Use this regridder to compare how well bilinear regridding and area weighted
# regridding preserve area weighted mean after round tripping.
```

```python

```
Loading