earthkit-geo - ECMWF internal earthkit training, March 2026

In this notebook you will see how to:

inspect the grid properties of GRIB data
plot the gridpoints of a field
interpolate GRIB data from one grid to another (regridding)
extract the nearest gridpoint from a field
obtain polygons suitable for a polytope request

1. Getting the data¶

The input data is a GRIB file containing 1000 hPa temperature fields on 3 different global grids.

First, fetch the file and list its contents.

import earthkit.data as ekd

ds = ekd.from_source("sample", "grids_3.grib").to_fieldlist()
ds.ls()

2. Select a field¶

Next, select the field you will inspect in the rest of the notebook.

To try another grid type, change the field index below and rerun the subsequent cells of the notebook.

# take one field
print("Possible indexes:", 0, 'to', len(ds)-1)
field_index = 0
f = ds[field_index]
f

2. Inspecting the grid¶

With the geography metadata component, you can inspect all the relevant geographical metadata.

f.geography.grid_type()

Inspect the ‘grid_spec’ - this is a canonical description of the grid used to identify it. See https://earthkit-geo.readthedocs.io/en/release-1.0.0rc0/guide/mir/gridspec.html for more details of this.

f.geography.grid_spec()

The field’s shape is derived from the geography. For grids with regular 2D structure, like regular latitude-longitude grids, the shape is always 2D. Otherwise, like for reduced Gaussian grids, it is 1D.

f.geography.shape()

3. Plotting the gridpoints¶

This example shows you how to plot the original gridpoint positions of the field.

import earthkit.plots as ekp

chart = ekp.Map(size=(7,7))
chart.contourf(f, units="celsius", auto_style=True)

# plot the original grid points
chart.grid_points(f, c="black")

chart.title(f"gridType={f.geography.grid_type()}")
chart.coastlines()
chart.gridlines()
chart.legend()

chart.show()

4. Regridding¶

You can regrid data with earthkit.geo.regrid(). This uses ECMWF’s MIR library to perform the regridding computations. The first time that a new combination of input and output grid is encountered, you may see some extra output logs, and the processing may take a few seconds. Subsequent regriddings with those grids will be faster, as they will use cached interpolation matrices.

import earthkit.geo as ekg

# the target grid is a global 10x10 degree regular latitude-longitude grid
out_grid = {"grid": [10,10]}

# exercise: if you are currently on the 3rd field, i.e. HEALPix, look at the grid_spec printed above, and generate an 'H8' grid
# by creating and using a new grid_spec which is a modification of it

# perform interpolation to produce a new field
f_res = ekg.regrid(f, grid=out_grid, interpolation="linear")

print(f_res.geography.grid_spec())

The next cell plots the original and interpolated fields over a subarea.

figure = ekp.Figure(domain="North Atlantic", rows=2, columns=1)

# the original field
subplot = figure.add_map(0, 0)
subplot.contourf(f, units="celsius", auto_style=True)
subplot.grid_points(f, c="black", s=25)
subplot.title("Original field")

# the interpolated field
subplot = figure.add_map(1, 0)
subplot.contourf(f_res, units="celsius", auto_style=True)
subplot.grid_points(f_res, c="black", s=25)
subplot.title("Interpolated field")

figure.coastlines()
figure.gridlines()

figure.show()

5. Getting the nearest gridpoint¶

In this example you can see how to extract the nearest gridpoint value of a field using earthkit.geo.nearest_point_haversine().

# ref location (lat, lon)
p_ref = (51.45, -0.97)

# get latlon
latlon = f.geography.latlons()
lat = latlon[0]
lon = latlon[1]

# get nearest point index
idx, dist = ekg.nearest_point_haversine(p_ref, (lat, lon))
idx

# extract field value at given index
f.values[idx]

Goto 2¶

Now go back to part 2 of this notebook, increase the field index and re-run the code up to this point

6. Working with NetCDF data¶

It is more natural to work with NetCDF data in xarray representation, which Earthkit also supports.

# download some example NetCDF data and read into an xarray
fs = ekd.from_source("sample", "era5-hourly-2t-20230724T1200Z.nc")
f = fs.to_xarray()
t2 = f["t2m"]
t2

# note that although Eathkit in the future will be able to derive the grid_spec from an xarray, it is currently
# required to explicitly supply this information. From looking at the metadata above, we can see that the xarray
# is a regular 0.25-degree grid, so we add a grid_spec attribute to annotate the DataArray with this information
t2.attrs["ek_grid_spec"] = '{"grid": [0.25, 0.25]}'
t2

# the target grid is a global 10x10 degree regular latitude-longitude grid
out_grid = {"grid": [10,10]}

# perform interpolation
f_res = ekg.regrid(t2, grid=out_grid)
f_res

figure = ekp.Figure(domain="North Atlantic", rows=2, columns=1)

# the original field
subplot = figure.add_map(0, 0)
subplot.contourf(f,  auto_style=True)
subplot.grid_points(f, c="black", alpha=0.3, s=1)
subplot.title("Original field")

# the interpolated field
subplot = figure.add_map(1, 0)
subplot.contourf(f_res,  auto_style=True)
subplot.grid_points(f_res, c="black", s=25)
subplot.title("Interpolated field")

figure.coastlines()
figure.gridlines()

figure.show()

Extra work¶

Extract and plot data for specific regions¶

This example uses earthkit-geo’s country_polygons() function to return the points describing country boundaries. We will then use that in a polytope request that will contact the ECMWF MARS server and request data for just those regions. The result will be in the covjson format, which we will convert to xarray and plot.

from earthkit.geo.cartography import country_polygons

domain = ["France", "Ireland"]

domain_polys = country_polygons(domain)

Quick polytope primer: if the below request fails because of authemtication, you may need to create your auth file by coping your ~/.ecmwfapirc file into ~/.polytopeapirc and renaming ‘key’ to ‘user_key’, and ‘e-mail’ to ‘user_email’.

request = {
    "stream": "oper",
    "levtype": "sfc",
    "param": "2t",
    "step": "12",
    "time": "00",
    "date": "20260314",
    "type": "fc",
    "domain": "g",
    "feature" : {
        "type" : "polygon",
        "shape" : domain_polys,
    },
}

data_from_domain = ekd.from_source("polytope", "ecmwf-mars", request=request, stream=False)

data_from_domain_xr = data_from_domain.to_xarray()
data_from_domain_xr

from shapely.geometry import Polygon

fig = ekp.Figure(rows=1, columns=1)
m = fig.add_map(domain=domain)
m.quickplot(data_from_domain_xr)
m.coastlines()
m.gridlines()
m.legend(location="bottom")
m.title()   # auto‑title if metadata present

fig.show()

Plot grid values¶

Inspect the following code to see how to plot grid point values on the map

The next plot shows a smaller area and displays the grid values at each grid point. We will use the data() method, which returns arrays of latitudes, longitudes and values; this information can then be used to create value labels. And with a bit of cartopy magic, we can shift the labels to sit nicely above the points!

Since the raw data values are actually stored in Kelvin, we will also revert the shading to use Kelvin so we are seeing the ‘real’ data.

import cartopy.crs as ccrs

f = ds[field_index]

chart = ekp.Map(domain="Europe")
chart.contourf(f, units="kelvin", auto_style=True)
# plot the original grid points
chart.grid_points(f, c="black") #marker="+"

# generate grid values
lat, lon, vals = f.data(flatten=True)
labels = [f"{x:.1f}" for x in vals]
for i, lbs in enumerate(labels):
    chart.ax.annotate(lbs, (lon[i], lat[i]), transform=ccrs.Geodetic(), 
                      xytext=(0,5),
                      textcoords="offset pixels", 
                      annotation_clip=True, 
                      fontsize=6, horizontalalignment="center")

chart.title(f"gridType={f.geography.grid_type()}")
chart.coastlines()
chart.gridlines()
chart.legend()

chart.show()