Interpolate grid data to station data

This example demonstrates how to interpolate 2m temperature data from a regular grid (ERA5 reanalysis data) to specific point locations (cities in Eastern China) using the easyclimate.interp.interp_mesh2point function from the easyclimate package.

The workflow includes:

1. Loading and visualizing the grid data

2. Creating a DataFrame with point locations

3. Performing the interpolation

4. Visualizing the results with interpolated values

Before proceeding with all the steps, first import some necessary libraries and packages

import easyclimate as ecl
import xarray as xr
import cartopy.crs as ccrs
import numpy as np
import pandas as pd

Load 2m temperature data from a NetCDF file (ERA5 reanalysis data) The data contains temperature values on a regular latitude-longitude grid

t2m_data = xr.open_dataset("js_t2m_ERA5_2025052000.nc").t2m
t2m_data

<xarray.DataArray 't2m' (lat: 81, lon: 121)> Size: 39kB
[9801 values with dtype=float32]
Coordinates:
    number      int64 8B ...
    valid_time  datetime64[ns] 8B ...
  * lat         (lat) float64 648B 45.0 44.75 44.5 44.25 ... 25.5 25.25 25.0
  * lon         (lon) float64 968B 100.0 100.2 100.5 100.8 ... 129.5 129.8 130.0
    expver      <U4 16B ...
Attributes: (12/32)
    GRIB_paramId:                             167
    GRIB_dataType:                            an
    GRIB_numberOfPoints:                      1038240
    GRIB_typeOfLevel:                         surface
    GRIB_stepUnits:                           1
    GRIB_stepType:                            instant
    ...                                       ...
    GRIB_totalNumber:                         0
    GRIB_units:                               K
    long_name:                                2 metre temperature
    units:                                    degC
    standard_name:                            unknown
    GRIB_surface:                             0.0

xarray.DataArray

't2m'

• lat: 81
• lon: 121

• ...

  [9801 values with dtype=float32]

• Coordinates: (5)
  • number
    ()
    int64
    ...

    long_name :

    ensemble member numerical id

    units :

    1

    standard_name :

    realization

    [1 values with dtype=int64]

  • valid_time
    ()
    datetime64[ns]
    ...

    long_name :

    time

    standard_name :

    time

    [1 values with dtype=datetime64[ns]]

  • lat
    (lat)
    float64
    45.0 44.75 44.5 ... 25.5 25.25 25.0

    units :

    degrees_north

    standard_name :

    latitude

    long_name :

    latitude

    stored_direction :

    decreasing

    array([45.  , 44.75, 44.5 , 44.25, 44.  , 43.75, 43.5 , 43.25, 43.  , 42.75,
           42.5 , 42.25, 42.  , 41.75, 41.5 , 41.25, 41.  , 40.75, 40.5 , 40.25,
           40.  , 39.75, 39.5 , 39.25, 39.  , 38.75, 38.5 , 38.25, 38.  , 37.75,
           37.5 , 37.25, 37.  , 36.75, 36.5 , 36.25, 36.  , 35.75, 35.5 , 35.25,
           35.  , 34.75, 34.5 , 34.25, 34.  , 33.75, 33.5 , 33.25, 33.  , 32.75,
           32.5 , 32.25, 32.  , 31.75, 31.5 , 31.25, 31.  , 30.75, 30.5 , 30.25,
           30.  , 29.75, 29.5 , 29.25, 29.  , 28.75, 28.5 , 28.25, 28.  , 27.75,
           27.5 , 27.25, 27.  , 26.75, 26.5 , 26.25, 26.  , 25.75, 25.5 , 25.25,
           25.  ])

  • lon
    (lon)
    float64
    100.0 100.2 100.5 ... 129.8 130.0

    units :

    degrees_east

    standard_name :

    longitude

    long_name :

    longitude

    array([100.  , 100.25, 100.5 , 100.75, 101.  , 101.25, 101.5 , 101.75, 102.  ,
           102.25, 102.5 , 102.75, 103.  , 103.25, 103.5 , 103.75, 104.  , 104.25,
           104.5 , 104.75, 105.  , 105.25, 105.5 , 105.75, 106.  , 106.25, 106.5 ,
           106.75, 107.  , 107.25, 107.5 , 107.75, 108.  , 108.25, 108.5 , 108.75,
           109.  , 109.25, 109.5 , 109.75, 110.  , 110.25, 110.5 , 110.75, 111.  ,
           111.25, 111.5 , 111.75, 112.  , 112.25, 112.5 , 112.75, 113.  , 113.25,
           113.5 , 113.75, 114.  , 114.25, 114.5 , 114.75, 115.  , 115.25, 115.5 ,
           115.75, 116.  , 116.25, 116.5 , 116.75, 117.  , 117.25, 117.5 , 117.75,
           118.  , 118.25, 118.5 , 118.75, 119.  , 119.25, 119.5 , 119.75, 120.  ,
           120.25, 120.5 , 120.75, 121.  , 121.25, 121.5 , 121.75, 122.  , 122.25,
           122.5 , 122.75, 123.  , 123.25, 123.5 , 123.75, 124.  , 124.25, 124.5 ,
           124.75, 125.  , 125.25, 125.5 , 125.75, 126.  , 126.25, 126.5 , 126.75,
           127.  , 127.25, 127.5 , 127.75, 128.  , 128.25, 128.5 , 128.75, 129.  ,
           129.25, 129.5 , 129.75, 130.  ])

  • expver
    ()
    <U4
    ...

    [1 values with dtype=<U4]

• Indexes: (2)
  • lat
    PandasIndex

    PandasIndex(Index([ 45.0, 44.75,  44.5, 44.25,  44.0, 43.75,  43.5, 43.25,  43.0, 42.75,
            42.5, 42.25,  42.0, 41.75,  41.5, 41.25,  41.0, 40.75,  40.5, 40.25,
            40.0, 39.75,  39.5, 39.25,  39.0, 38.75,  38.5, 38.25,  38.0, 37.75,
            37.5, 37.25,  37.0, 36.75,  36.5, 36.25,  36.0, 35.75,  35.5, 35.25,
            35.0, 34.75,  34.5, 34.25,  34.0, 33.75,  33.5, 33.25,  33.0, 32.75,
            32.5, 32.25,  32.0, 31.75,  31.5, 31.25,  31.0, 30.75,  30.5, 30.25,
            30.0, 29.75,  29.5, 29.25,  29.0, 28.75,  28.5, 28.25,  28.0, 27.75,
            27.5, 27.25,  27.0, 26.75,  26.5, 26.25,  26.0, 25.75,  25.5, 25.25,
            25.0],
          dtype='float64', name='lat'))

  • lon
    PandasIndex

    PandasIndex(Index([ 100.0, 100.25,  100.5, 100.75,  101.0, 101.25,  101.5, 101.75,  102.0,
           102.25,
           ...
           127.75,  128.0, 128.25,  128.5, 128.75,  129.0, 129.25,  129.5, 129.75,
            130.0],
          dtype='float64', name='lon', length=121))

• Attributes: (32)

  GRIB_paramId :

  167

  GRIB_dataType :

  an

  GRIB_numberOfPoints :

  1038240

  GRIB_typeOfLevel :

  surface

  GRIB_stepUnits :

  1

  GRIB_stepType :

  instant

  GRIB_gridType :

  regular_ll

  GRIB_uvRelativeToGrid :

  0

  GRIB_NV :

  0

  GRIB_Nx :

  1440

  GRIB_Ny :

  721

  GRIB_cfName :

  unknown

  GRIB_cfVarName :

  t2m

  GRIB_gridDefinitionDescription :

  Latitude/Longitude Grid

  GRIB_iDirectionIncrementInDegrees :

  0.25

  GRIB_iScansNegatively :

  0

  GRIB_jDirectionIncrementInDegrees :

  0.25

  GRIB_jPointsAreConsecutive :

  0

  GRIB_jScansPositively :

  0

  GRIB_latitudeOfFirstGridPointInDegrees :

  90.0

  GRIB_latitudeOfLastGridPointInDegrees :

  -90.0

  GRIB_longitudeOfFirstGridPointInDegrees :

  0.0

  GRIB_longitudeOfLastGridPointInDegrees :

  359.75

  GRIB_missingValue :

  3.4028234663852886e+38

  GRIB_name :

  2 metre temperature

  GRIB_shortName :

  2t

  GRIB_totalNumber :

  0

  GRIB_units :

  K

  long_name :

  2 metre temperature

  units :

  degC

  standard_name :

  unknown

  GRIB_surface :

  0.0

Create a basemap plot focused on Eastern China region

fig, ax = ecl.plot.quick_draw_spatial_basemap(figsize=(5, 5))
ax.set_extent([116, 123, 30, 36], crs = ccrs.PlateCarree())

# Select and plot temperature data for the region of interest
draw_data = t2m_data.sel(lon = slice(100, 140), lat = slice(45, 25))
draw_data.plot.contourf(
    ax = ax,
    transform=ccrs.PlateCarree(),
    cbar_kwargs = {'location': 'bottom'},
    levels = np.linspace(17, 29, 21),
    cmap = "plasma"
)

ax.set_title("2m temperature")

Text(0.5, 1.0, '2m temperature')

Create a DataFrame containing city locations (longitude and latitude) These are the points where we want to interpolate temperature values

data = {
    "Site": ["Nanjing (南京)", "Suzhou (苏州)", "Shanghai (上海)", "Chuzhou (滁州)", "Changzhou (常州)", "Xuzhou (徐州)", "Yancheng (盐城)"],
    "lon": [118.7788631, 120.6212881, 121.4700152, 118.3139455, 119.9691539, 117.1810431, 120.1577019],
    "lat": [32.0438284, 31.311123, 31.2312707, 32.3027377, 31.8122623, 34.2665258, 33.349559]
}
df = pd.DataFrame(data)
df

	Site	lon	lat
0	Nanjing (南京)	118.778863	32.043828
1	Suzhou (苏州)	120.621288	31.311123
2	Shanghai (上海)	121.470015	31.231271
3	Chuzhou (滁州)	118.313946	32.302738
4	Changzhou (常州)	119.969154	31.812262
5	Xuzhou (徐州)	117.181043	34.266526
6	Yancheng (盐城)	120.157702	33.349559

Use interp_mesh2point to interpolate grid values to point locations Parameters:

• t2m_data: Input grid data (xarray DataArray)

• df: DataFrame with point locations

• lon/lat_dim_mesh: Name of lon/lat dimensions in grid data

• lon/lat_dim_df: Name of lon/lat columns in DataFrame

df_interp = ecl.interp.interp_mesh2point(
    t2m_data, df,
    lon_dim_mesh = "lon",
    lat_dim_mesh = "lat",
    lon_dim_df = "lon",
    lat_dim_df = "lat"
)
df_interp

	Site	lon	lat	interpolated_value
0	Nanjing (南京)	118.778863	32.043828	27.434105
1	Suzhou (苏州)	120.621288	31.311123	27.056907
2	Shanghai (上海)	121.470015	31.231271	25.860559
3	Chuzhou (滁州)	118.313946	32.302738	27.139358
4	Changzhou (常州)	119.969154	31.812262	27.489566
5	Xuzhou (徐州)	117.181043	34.266526	28.625858
6	Yancheng (盐城)	120.157702	33.349559	26.331501

Create a combined plot showing both the grid data and interpolated points

proj_trans = ccrs.PlateCarree() # Coordinate reference system for transformations

fig, ax = ecl.plot.quick_draw_spatial_basemap(figsize=(5, 5))
ax.set_extent([116, 123, 30, 36], crs = proj_trans)

draw_data = t2m_data.sel(lon = slice(100, 140), lat = slice(45, 25))

# Plot the grid data again for reference
draw_data.plot.contourf(
    ax = ax,
    transform=ccrs.PlateCarree(),
    cbar_kwargs = {'location': 'bottom'},
    levels = np.linspace(17, 29, 21),
    cmap = "plasma"
)

# Plot the point locations as red dots
ax.scatter(
    df_interp["lon"],
    df_interp["lat"],
    transform = proj_trans,
    color = 'r',
    s = 5
)

# Add temperature values as text labels near each point
for i, row in df_interp.iterrows():
    ax.text(
        row["lon"],
        row["lat"],
        str(np.round(row["interpolated_value"], decimals=2)), # Rounded to 2 decimal places
        transform=proj_trans,
        fontsize=10,
        ha='center',  # Horizontal alignment
        va='bottom',  # Vertical alignment
        color='blue'
)

ax.set_title("2m temperature (Points)")

Text(0.5, 1.0, '2m temperature (Points)')