Regridding HEALPix data to 1x1 degree

For most analyses, we recommend to stay on the HEALPix grid as long as possible. However, the comparison to other data products or the application of some algorithms may require an interpolation to a regular grid. In this example, we show how the HEALPix hierarchy can also help in this regard.

Usually, regridding of storm-resolving simulation output is computationally heavy. By making use of the pre-computed grid hierarchy, we can reduce this significantly. For our example of the 1x1 degree grid, the zoom-level 7 in HEALPix is rather close to our target resolution. Regridding from this zoom level is already 64 times faster than the original model output. In the following, we briefly outline how to use Python and CDO to perform efficient regriddings.

Using Python

import intake
import matplotlib.pyplot as plt
import healpy as hp
import cmocean
import xarray as xr
import numpy as np

cat = intake.open_catalog("https://data.nextgems-h2020.eu/catalog.yaml")
data = cat["ICON"]["ngc3028"](zoom=7, time="P1D").to_dask()

# Find the HEALPix pixels that are closest to the 1x1 degree grid points.
lon = np.arange(0, 360, 1)
lat = np.arange(90, -91, -1)
pix = xr.DataArray(
    hp.ang2pix(data.crs.healpix_nside, *np.meshgrid(lon, lat), nest=True, lonlat=True),
    coords=(("lat", lat), ("lon", lon)),
)

# Plot the 1x1 world map for a given variable and time.
data.tas.isel(cell=pix).sel(time="2023-05-30").plot(
    vmin=253, vmax=303, cmap=cmocean.cm.thermal
)

<matplotlib.collections.QuadMesh at 0x7fff7d472cb0>

Using CDO

We can also do the same nearest-neighbor interpolation using CDO. By default, the CDO global grid is defined from -180 to 180 degrees. To be consistent with FESOM, we can store our own grid definition in a text file:

%%file global_1x1.txt
#
# Global 1x1 deg grid
#
gridtype = lonlat
gridsize = 65160
xsize = 360
ysize = 181
xname = lon
xlongname = "longitude"
xunits = "degrees_east"
yname = lat
ylongname = "latitude"
yunits = "degrees_north"
xfirst = 0
xinc = 1
yfirst = 90
yinc = -1

Overwriting global_1x1.txt

Now we can use CDO to perform a nearest-neighbor interpolation (rempann) onto this grid and plot the result for comparison.

For a more in-depth documentation on how to handle HEALPix with CDO have a look at the CDO documentation.

%%bash
module use /work/k20200/k202134/hsm-tools/outtake/module
module load hsm-tools/unstable
module load cdo

DATAPATH=$(query_yaml.py ICON ngc3028 -s time=P1D zoom=7 --cdo)

cdo \
  -f nc4 \
  -remapnn,global_1x1.txt \
  [ -select,name=tas,date="2023-05-30" ${DATAPATH} ] \
  ngc3028_tas.nc

cdo    remapnn: Nearest neighbor weights from unstructured (196608) to lonlat (360x181) grid

xr.open_dataset("ngc3028_tas.nc").tas.plot(vmin=253, vmax=303, cmap=cmocean.cm.thermal)

<matplotlib.collections.QuadMesh at 0x7fff791f2fb0>