Regridding: IFS-HRES and ICON-CH2 2 m temperature on a 200 m grid over Valais#

This example illustrates interp2grid() with four methods:

nearest: value of the nearest source grid point;
linear: linear interpolation (bilinear on regular grids, on a triangulation for unstructured grids);
idw: inverse distance weighting of the k nearest grid points, with weights d^(-weight_power);
temp_fixgrad: among the source grid points within r_search grid spacings, the one with the most similar altitude, corrected for the remaining altitude difference with a fixed lapse rate (6 K/km by default).

The 2 m temperature of IFS-HRES (regular 0.1° grid, about 9 km) and of the ICON-CH2 control member (unstructured grid, about 2 km) is regridded to the 200 m digital height model of swisstopo over Valais, a region of deep valleys and high mountains, for the same valid time.

Data: ECMWF open data, © ECMWF, CC BY 4.0. MeteoSwiss open data, CC BY 4.0, source: MeteoSwiss. DHM25/200, © swisstopo. Country borders: Natural Earth.

import contextlib
import datetime as dt
import io
import logging
import os
import re
import time
import timeit
import zipfile
from collections.abc import Callable
from functools import partial
from typing import Any

# import pyproj before earthkit.data to avoid a crash at exit (earthkit-data#1202)
import pyproj  # noqa: F401  # pylint: disable=unused-import
import eccodes_cosmo_resources

# ICON GRIB files need the COSMO definitions of ecCodes
os.environ["ECCODES_DEFINITION_PATH"] = str(eccodes_cosmo_resources.get_definitions_path())
os.environ["ECCODES_VERSION_CHECK_OFF"] = "1"

import cartopy.crs as ccrs
import earthkit.data as ekd
import earthkit.plots as ekp
import geopandas as gpd
import numpy as np
import pandas as pd
import requests
import xarray as xr
from earthkit.data.utils import download_and_cache
from requests.adapters import HTTPAdapter
from urllib3.util import Retry

import grid_process.grid_interp as gri
import grid_process.timeseries as tms
import grid_process.units_conv as uconv
from grid_process import data_model as dm

METHODS: dict[str, dict[str, Any]] = {
    "nearest": {"method": "nearest"},
    "linear": {"method": "linear"},
    "idw": {"method": "idw", "k": 4, "weight_power": 2},
    "temp_fixgrad": {"method": "temp_fixgrad"},
}
LEAD_H_IFS = 24  # lead time [h] of the latest IFS-HRES run defining the common valid time

Visual comparison#

One figure per model, with the same colour scale for both.

lv03 = ccrs.epsg(21781)
t_min = min(float(res[m]["T_2M"].min()) for res in results.values() for m in METHODS)
t_max = max(float(res[m]["T_2M"].max()) for res in results.values() for m in METHODS)
style = ekp.styles.Style(
    levels=np.arange(np.floor(t_min), np.ceil(t_max) + 1, 1), colors="viridis", extend="both", units="°C"
)

for name, res in results.items():
    fig = ekp.Figure(rows=2, columns=2, figsize=(12, 7.5), crs=lv03)
    for i, label in enumerate(METHODS):
        out = res[label]
        subplot = fig.add_map(i // 2, i % 2)
        subplot.pcolormesh(out["T_2M"].squeeze().values, x=out.x.values, y=out.y.values, style=style, transform=lv03)
        subplot.title(f"method={label}")
    fig.borders()
    fig.legend(label="2 m temperature [°C]")
    fig.title(f"{name} 2 m temperature on the DHM25/200 grid, valid {valid_time:%Y-%m-%d %H:%M} UTC")
    fig.show()

../_images/62618716ef14780bf1d52c229ebbfad18cdfa11ed0de2e4f3598ee90bfa02d4e.png

../_images/6dd262a9cf4ab98c74250a5b06b93d4f858b8ab8e196e9807aeb61aad4daf0c8.png

Performance benchmark#

The first call includes building the neighbour search structures and reprojecting the coordinates, which grid_process caches. Subsequent calls on the same grids (e.g. for other variables or runs) reuse them: their runtime is the mean ± standard deviation over 5 repetitions.

def benchmark(
    func: Callable[..., xr.Dataset],
    data: xr.Dataset,
    target: gpd.GeoDataFrame | xr.Dataset,
    methods: dict[str, dict[str, Any]],
    first_call_times: dict[str, float],
    repeat: int = 5,
) -> pd.DataFrame:
    """Runtime of the first call and mean ± standard deviation of the subsequent calls of each method."""
    rows: dict[str, dict[str, Any]] = {}
    for method, kw in methods.items():
        times = np.array(timeit.repeat(partial(func, data, target, **kw), number=1, repeat=repeat))
        rows[method] = {
            "first call [s]": f"{first_call_times[method]:.3f}",
            "subsequent calls [s]": f"{times.mean():.3f} ± {times.std():.3f}",
        }
    return pd.DataFrame.from_dict(rows, orient="index").rename_axis("method")


for name, ds in models.items():
    src_shape = " × ".join(str(ds.sizes[dim]) for dim in dm.get_spatial_dims(ds))
    print(f"{name}: {src_shape} source grid points -> {dst_grid.y.size} × {dst_grid.x.size} cells")
# one column group per model; the try_reg=False rows only exist for the regular IFS grid
tables = {name: benchmark(gri.interp2grid, ds, dst_grid, runs[name], first_call[name]) for name, ds in models.items()}
pd.concat(tables, axis=1).fillna("–")

IFS-HRES: 26 × 45 source grid points -> 496 × 833 cells
ICON-CH2: 283876 source grid points -> 496 × 833 cells

	IFS-HRES		ICON-CH2
	first call [s]	subsequent calls [s]	first call [s]	subsequent calls [s]
method
nearest	0.337	0.065 ± 0.003	0.126	0.058 ± 0.002
nearest (try_reg=False)	0.105	0.060 ± 0.004	–	–
linear	0.090	0.076 ± 0.005	0.108	0.103 ± 0.003
linear (try_reg=False)	0.072	0.065 ± 0.001	–	–
idw	0.168	0.098 ± 0.005	0.212	0.098 ± 0.006
temp_fixgrad	0.288	0.174 ± 0.011	0.313	0.172 ± 0.002

Takeaways#

With IFS-HRES (9 km), nearest shows the source grid boxes. linear and idw are smooth, but like the model they ignore the terrain at 200 m: the Rhône valley gets the temperature of a much higher model altitude, and the summits of a much lower one.
temp_fixgrad restores the terrain structure by correcting with the altitude difference between the source grid point and the 200 m grid: warmer valleys, colder summits. ICON-CH2 (2 km) already resolves the main valleys, so its correction is smaller.
The fixed lapse rate (grad, 6 K/km by default) is a simplification and cannot reproduce temperature inversions.
Every method regrids both models to the 413 000 cells of the destination grid in well under a second. temp_fixgrad is the most expensive.
On the regular IFS grid, the regular path (try_reg=True, the default) and the generic path (try_reg=False) of nearest select the same source grid point, except for a handful of destination cells lying almost on the boundary between two source grid boxes, where the reprojection tips the choice. linear differs slightly everywhere: the regular path is bilinear, while the generic path interpolates linearly on a triangulation of the source grid points. Here the regular path is the faster one.

Regridding: IFS-HRES and ICON-CH2 2 m temperature on a 200 m grid over Valais#

Destination grid#

Source data#

Interpolation#

Visual comparison#

Performance benchmark#

Takeaways#