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.

import intake


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

/builds/easy/gems/.venv/lib/python3.13/site-packages/intake_xarray/base.py:21: FutureWarning: The return type of `Dataset.dims` will be changed to return a set of dimension names in future, in order to be more consistent with `DataArray.dims`. To access a mapping from dimension names to lengths, please use `Dataset.sizes`.
  'dims': dict(self._ds.dims),

One of the defining properties of HEALPix 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:

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:

from easygems.show import map_show

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

<Figure size 640x480 with 0 Axes>

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:

from easygems.show import map_contour

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

<Figure size 640x480 with 0 Axes>

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).

Nearest-neighbor resampling

Next, we will load an ICON dataset on its native grid. The ICON grid 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.

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 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:

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:

map_show(ds["ts"].isel(time=0), resampler=kd_resampler);

<Figure size 640x480 with 0 Axes>

Delaunay triangulation

Finally, there is a specialised resampler for plotting smoothed maps. It will perform a Delaunay triangulation and return a weighted average of the surrounding grid points.

from easygems.resample import DelaunayResampler

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

This can be especially useful when creating regional maps.

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);

/builds/easy/gems/.venv/lib/python3.13/site-packages/cartopy/io/__init__.py:241: DownloadWarning: Downloading: https://naturalearth.s3.amazonaws.com/10m_physical/ne_10m_coastline.zip
  warnings.warn(f'Downloading: {url}', DownloadWarning)