---
file_format: mystnb
kernelspec:
  name: python3
  display_name: Python 3
execution:
  timeout: 180
---

# Plotting geospatial data on a map

## Plotting in pixel-space

When working with geospatial data, it is common to remap datasets before plotting them on a map.
While this approach is convenient, it not only obscures where the data actually originates on Earth but also requires substantial computational resources, especially for high-resolution fields.
Moreover, this workflow performs an often overlooked extra step: even after remapping to a new grid, the data is *again* resampled into pixel-space during plotting.
In other words, traditional approaches effectively remap the data twice.

In this tutorial, we take a different route: we work directly in pixel-space.
Instead of remapping the data upfront, we first create a Cartopy map axis and then determine the geographic coordinates associated with each pixel on that axis.
By sampling our original dataset at those coordinates, we construct the map directly from the source data.
This method avoids unnecessary interpolation, reduces overhead, and provides a clearer view of how data and map projections interact.

## HEALPix resampling

As a first example, we will plot data on the HEALPix grid.
For this purpose, we will use a HEALPix-ified version of the ERA5 reanalysis dataset.

```{code-cell} python3
import intake


cat = intake.open_catalog("https://data.nextgems-h2020.eu/online.yaml")
ds = cat.ERA5(zoom=6).to_dask()
```

One of the defining properties of [HEALPix](healpix/index) is that the grid points closest to a given latitude and longitude can be queried efficiently based on analytical equations.
All that is required is knowledge of the defining parameters of the HEALPix grid, i.e. its refinement level and nesting scheme.
The `easygems` package provides a `HEALPixResampler` helper class that can be initialised using these parameters:

```{code-cell} python3
from easygems.resample import HEALPixResampler

hp_resampler = HEALPixResampler(nside=ds.crs.healpix_nside)
```

Now we can pass both our datasets and the resampler to the `map_show()` function, which will create a default world map and fill it based on our resampler:

```{code-cell} python3
from easygems.show import map_show

map_show(ds["2t"].isel(time=0), resampler=hp_resampler);
```

The resampler is designed to simply return values for a given coordinate.
Therefore, it can also be used to create different types of map, such as contour lines:

```{code-cell} python3
from easygems.show import map_contour

map_contour(ds["2t"].isel(time=0), resampler=hp_resampler);
```

```{tip}
Both `map_show()` and `map_contour()` can be used with any resampler.
For smooth contour lines, it is recommended to use an interpolating resampler (see, e.g., [Delaunay triangulation](#delaunay-triangulation)).
```

## Nearest-neighbor resampling

Next, we will load an ICON dataset on its native grid.
[The ICON grid](playing_with_triangles/index) is based on an icosahedron that is bisected to create different resolutions, which are then further processed for the numerical stability of simulations.
While visually still close to the original icosahedron, the resulting grid points are technically a rather unstructured set of latitudes and longitudes.

```{code-cell} python3
import xarray as xr

from easygems.show import map_show


ds = xr.open_dataset(
    "https://eerie.cloud.dkrz.de/datasets/icon-esm-er.eerie-control-1950.v20240618.atmos.native.2d_6h_inst/kerchunk",
    engine="zarr",
    chunks={},
).squeeze()
```

In order to resample the grid, we will construct a [_k_-d tree](https://en.wikipedia.org/wiki/K-d_tree) to efficiently query nearest-neighbor points.
This querying is best done in _xyz_ space to find the actual closest points at the poles and date line.
The `easygems` package does provide a convenience class for this:

```{code-cell} python3
from easygems.resample import KDTreeResampler

kd_resampler = KDTreeResampler(lat=ds.lat, lon=ds.lon)
```

Now we can pass both our datasets and the resampler to the `map_show()` function, which will create a default world map and fill it based on our resampler:

```{code-cell} python3
map_show(ds["ts"].isel(time=0), resampler=kd_resampler);
```

## Delaunay triangulation

Finally, there is a specialised resampler for plotting smoothed maps.
It will perform a [Delaunay triangulation](https://en.wikipedia.org/wiki/Delaunay_triangulation) and return a weighted average of the surrounding grid points.

```{code-cell} python3
from easygems.resample import DelaunayResampler

tri_resampler = DelaunayResampler(lat=ds.lat, lon=ds.lon)
```

This can be especially useful when creating regional maps.


```{code-cell} python3
import cartopy.crs as ccrs
import matplotlib.pyplot as plt


fig, ax = plt.subplots(subplot_kw={"projection": ccrs.PlateCarree()})
ax.set_extent([-12, 2, 46, 60])
ax.coastlines()
map_show(ds["uas"].isel(time=0), resampler=tri_resampler);
```