Interpolation and Regriding

Note

Go to the end to download the full example code.

Interpolation and Regriding#

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

import easyclimate as ecl
import xarray as xr
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import numpy as np

Interpolation from points to grid#

Open sample surface pressure data for the European region

data = ecl.open_tutorial_dataset("PressQFF_202007271200_872.csv")
print(data)

         lat      lon     qff
0    43.8569   4.4064  1016.2
1    56.3265  -3.7273   995.1
2    55.7350  24.4170  1015.5
3    60.3000  -4.2000   997.8
4    57.8500  45.7667  1020.5
..       ...      ...     ...
867  38.8766  22.4360  1011.8
868  68.8493  28.2991  1018.2
869  39.0633  26.5983  1010.9
870  49.4464   2.1272  1009.8
871  61.3667  30.8833  1021.2

[872 rows x 3 columns]

easyclimate.interp.interp_point2mesh enables interpolation from site data to grid point data.

See also

meshdata = ecl.interp.interp_point2mesh(
    data,
    var_name="qff",
    grid_x=37.5,
    grid_y=75.0,
    point=[-26.0, 34.5],
    resolution=32,
    sigma=1,
)
meshdata

<xarray.DataArray 'qff' (lat: 1200, lon: 2400)> Size: 12MB
array([[1023.1991 , 1023.1991 , 1023.1991 , ..., 1002.229  , 1002.22516,
        1002.2213 ],
       [1023.1988 , 1023.1988 , 1023.1988 , ..., 1002.21936, 1002.21545,
        1002.2116 ],
       [1023.1984 , 1023.1984 , 1023.1984 , ..., 1002.2096 , 1002.2056 ,
        1002.2017 ],
       ...,
       [       nan,        nan,        nan, ...,        nan,        nan,
               nan],
       [       nan,        nan,        nan, ...,        nan,        nan,
               nan],
       [       nan,        nan,        nan, ...,        nan,        nan,
               nan]], dtype=float32)
Coordinates:
  * lon      (lon) float64 19kB -25.97 -25.94 -25.91 -25.88 ... 48.94 48.97 49.0
  * lat      (lat) float64 10kB 34.5 34.53 34.56 34.59 ... 71.91 71.94 71.97


Plotting interpolated grid point data and corresponding station locations

fig, ax = plt.subplots(subplot_kw={"projection": ccrs.PlateCarree(central_longitude=0)})

ax.gridlines(draw_labels=["bottom", "left"], color="grey", alpha=0.5, linestyle="--")
ax.coastlines(edgecolor="black", linewidths=0.5)

# Draw interpolation results
meshdata.plot.contourf(
    ax=ax,
    transform=ccrs.PlateCarree(),
    cbar_kwargs={"location": "bottom"},
    cmap="RdBu_r",
    levels=21,
)

# Draw observation stations
ax.scatter(data["lon"], data["lat"], s=1, c="r", transform=ccrs.PlateCarree())
plot interp
<matplotlib.collections.PathCollection object at 0x7f804db31cf0>

Regriding#

Reading example raw grid data

u_data = ecl.tutorial.open_tutorial_dataset("uwnd_202201_mon_mean").sortby("lat").uwnd
u_data

<xarray.DataArray 'uwnd' (time: 2, level: 17, lat: 73, lon: 144)> Size: 1MB
[357408 values with dtype=float32]
Coordinates:
  * time     (time) datetime64[ns] 16B 2022-01-01 2022-02-01
  * lon      (lon) float32 576B 0.0 2.5 5.0 7.5 10.0 ... 350.0 352.5 355.0 357.5
  * lat      (lat) float32 292B -90.0 -87.5 -85.0 -82.5 ... 82.5 85.0 87.5 90.0
  * level    (level) float32 68B 1e+03 925.0 850.0 700.0 ... 50.0 30.0 20.0 10.0
Attributes:
    long_name:     Monthly mean u wind
    units:         m/s
    precision:     2
    var_desc:      u-wind
    level_desc:    Pressure Levels
    statistic:     Mean
    parent_stat:   Other
    dataset:       NCEP Reanalysis Derived Products
    actual_range:  [-68.194824 124.399994]


Define the target grid (only for latitude/longitude and regular grids)

target_grid = xr.DataArray(
    dims=("lat", "lon"),
    coords={
        "lat": np.arange(-89, 89, 6) + 1 / 1.0,
        "lon": np.arange(-180, 180, 6) + 1 / 1.0,
    },
)

easyclimate.interp.interp_point2mesh performs a regridding operation.

See also

EXCITED-CO2/xarray-regrid

regriding_data = ecl.interp.interp_mesh2mesh(u_data, target_grid)
regriding_data

/home/runner/work/easyclimate/easyclimate/src/easyclimate/interp/mesh2mesh.py:93: UserWarning: It seems that the input data longitude range is from -180° to 180°. Currently automatically converted to it from 0° to 360°.
  warnings.warn(

<xarray.DataArray 'uwnd' (time: 2, level: 17, lat: 30, lon: 60)> Size: 490kB
array([[[[ 5.60226893e-01,  1.02667853e+00,  1.49387230e+00, ...,
          -7.23611766e-01, -3.29127944e-01,  1.04678571e-01],
         [-7.13537864e-02,  4.89775395e-01,  1.05180766e+00, ...,
          -1.44522468e+00, -1.06899892e+00, -6.04192466e-01],
         [-3.99954724e+00, -3.65577347e+00, -3.40841788e+00, ...,
          -3.56099902e+00, -4.03412802e+00, -4.18725687e+00],
         ...,
         [ 4.90549397e-01,  6.82065892e-01,  8.77517456e-01, ...,
           3.63861387e+00,  1.77867858e+00,  7.16775246e-01],
         [-1.08903117e+00, -2.07193441e+00, -2.58209565e+00, ...,
           4.92451716e+00,  2.62419448e+00,  4.87097889e-01],
         [-3.48438617e+00, -4.12667664e+00, -4.65919235e+00, ...,
          -1.30906337e+00, -2.02451518e+00, -2.77096647e+00]],

        [[ 5.60516608e-01,  1.02732306e+00,  1.49006493e+00, ...,
          -7.03160739e-01, -3.10322070e-01,  1.14032865e-01],
         [-2.41967231e-01,  3.33032781e-01,  9.23710179e-01, ...,
          -1.57606402e+00, -1.21519297e+00, -7.65322164e-01],
         [-3.85738642e+00, -3.43561254e+00, -3.15858048e+00, ...,
          -3.72483817e+00, -4.09529017e+00, -4.14435428e+00],
...
         [ 3.96105365e+01,  4.09033581e+01,  4.17687487e+01, ...,
           3.27103223e+01,  3.58181055e+01,  3.79991065e+01],
         [ 1.64507136e+01,  1.70246421e+01,  1.75105355e+01, ...,
           1.25291067e+01,  1.42987498e+01,  1.56053565e+01],
         [-1.19253571e+00, -6.11321497e-01, -3.92851830e-03, ...,
          -2.87689324e+00, -2.31967870e+00, -1.76107121e+00]],

        [[-2.16849996e+00, -2.24924996e+00, -2.28517834e+00, ...,
          -1.70846427e+00, -1.90010710e+00, -2.05642876e+00],
         [ 9.25357080e-01,  8.11392877e-01,  6.53821403e-01, ...,
           1.17457153e+00,  1.10617847e+00,  1.01874988e+00],
         [-2.32710742e+00, -2.31753610e+00, -2.39192835e+00, ...,
          -2.25328593e+00, -2.32342877e+00, -2.35374999e+00],
         ...,
         [ 4.42924986e+01,  4.69607846e+01,  4.87700727e+01, ...,
           3.43513566e+01,  3.74099993e+01,  4.09591069e+01],
         [ 1.97691081e+01,  2.20551811e+01,  2.40682136e+01, ...,
           1.48973211e+01,  1.59275005e+01,  1.75919629e+01],
         [-1.64439259e+00, -5.79463816e-01,  4.93000412e-01, ...,
          -4.43860681e+00, -3.60410719e+00, -2.66571409e+00]]]])
Coordinates:
  * time     (time) datetime64[ns] 16B 2022-01-01 2022-02-01
  * level    (level) float32 68B 1e+03 925.0 850.0 700.0 ... 50.0 30.0 20.0 10.0
  * lat      (lat) float64 240B -88.0 -82.0 -76.0 -70.0 ... 68.0 74.0 80.0 86.0
  * lon      (lon) float64 480B 1.0 7.0 13.0 19.0 ... 337.0 343.0 349.0 355.0
Attributes:
    long_name:     Monthly mean u wind
    units:         m/s
    precision:     2
    var_desc:      u-wind
    level_desc:    Pressure Levels
    statistic:     Mean
    parent_stat:   Other
    dataset:       NCEP Reanalysis Derived Products
    actual_range:  [-68.194824 124.399994]


Plotting differences before and after interpolation

fig, ax = plt.subplots(1, 2, figsize=(12, 5))

u_data.sel(level=500).isel(time=0).plot(ax=ax[0])
ax[0].set_title("Before", size=20)

regriding_data.sel(level=500).isel(time=0).plot(ax=ax[1])
ax[1].set_title("After", size=20)
Before, After
Text(0.5, 1.0, 'After')

Total running time of the script: (0 minutes 6.362 seconds)