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 0x7faa16da4130>

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.

regriding_data = ecl.interp.interp_mesh2mesh(u_data, target_grid)
regriding_data
/home/runner/work/easyclimate/easyclimate/src/easyclimate/interp/mesh2mesh.py:94: 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')

Interpolation from model layers to altitude layers#

Suppose the following data are available

uwnd_data = xr.DataArray(
    np.array([  -15.080393 , -10.749071 ,  -4.7920494,  -2.3813725,  -1.4152431,
                -0.6465206,  -7.8181705, -14.718096 , -14.65539  , -14.948015 ,
                -13.705519 , -11.443476 ,  -8.865583 ,  -7.9528713,  -8.329103 ,
                -7.2445316,  -6.7150173,  -5.5189686,  -4.139448 ,  -3.2731838,
                -2.2931194,  -1.0041752,  -1.8983078,  -2.3674374,  -2.8203106,
                -3.2940865,  -3.526329 ,  -3.8654022,  -4.164995 ,  -4.2834396,
                -4.2950516,  -4.3438225,  -4.3716908,  -4.7688255,  -4.6155453,
                -4.5528393,  -4.4831676,  -4.385626 ,  -4.2950516,  -4.0953217]),
    dims=("model_level",),
    coords={
        "model_level": np.array([   36,  44,  51,  56,  60,  63,  67,  70,  73,  75,  78,  81,  83,
                                    85,  88,  90,  92,  94,  96,  98, 100, 102, 104, 105, 107, 109,
                                    111, 113, 115, 117, 119, 122, 124, 126, 129, 131, 133, 135, 136,
                                    137])
    },
)

p_data = xr.DataArray(
    np.array([  2081.4756,   3917.6995,   6162.6455,   8171.3506,  10112.652 ,
                11811.783 ,  14447.391 ,  16734.607 ,  19317.787 ,  21218.21  ,
                24357.875 ,  27871.277 ,  30436.492 ,  33191.027 ,  37698.96  ,
                40969.438 ,  44463.73  ,  48191.92  ,  52151.29  ,  56291.098 ,
                60525.63  ,  64770.367 ,  68943.69  ,  70979.66  ,  74908.17  ,
                78599.67  ,  82012.95  ,  85122.69  ,  87918.29  ,  90401.77  ,
                92584.94  ,  95338.72  ,  96862.08  ,  98165.305 ,  99763.38  ,
                100626.21  , 101352.69  , 101962.28  , 102228.875 , 102483.055 ]),
    dims=("model_level",),
    coords={
        "model_level": np.array([   36,  44,  51,  56,  60,  63,  67,  70,  73,  75,  78,  81,  83,
                                    85,  88,  90,  92,  94,  96,  98, 100, 102, 104, 105, 107, 109,
                                    111, 113, 115, 117, 119, 122, 124, 126, 129, 131, 133, 135, 136,
                                    137])
    },
)

Now we interpolate the data located on the mode plane to the isobaric plane. Note that the units of the given isobaric surface are consistent with pressure_data.

result = ecl.interp.interp1d_vertical_model2pressure(
    pressure_data=p_data,
    variable_data=uwnd_data,
    vertical_input_dim="model_level",
    vertical_output_dim="plev",
    vertical_output_level=np.array(
        [100000, 92500, 85000, 70000, 60000, 50000, 40000, 30000, 20000, 10000]
    ),
)
result
<xarray.DataArray (plev: 10)> Size: 80B
array([ -4.59829518,  -4.29460496,  -3.85226099,  -2.14340702,
        -2.41097966,  -4.87545826,  -7.55671981,  -9.28851613,
       -14.76362259,  -1.46601784])
Coordinates:
  * plev     (plev) int64 80B 100000 92500 85000 70000 ... 30000 20000 10000


Simply calibrate the interpolation effect.

fig, ax = plt.subplots()
ax.plot(p_data, uwnd_data, label="Original")
ax.plot(result.plev, result.data, "o", label="Interpolated")
ax.invert_xaxis()
ax.set_xlabel("Pressure [Pa]")
ax.set_ylabel("Zonal Wind [m/s]")
plt.legend()
plot interp
<matplotlib.legend.Legend object at 0x7faa17168f70>

Interpolation from pressure layers to altitude layers#

Suppose the following data are available

z_data = xr.DataArray(
    np.array([  214.19354,   841.6774 ,  1516.871  ,  3055.7097 ,  4260.5806 ,
                5651.4194 ,  7288.032  ,  9288.193  , 10501.097  , 11946.71   ,
                13762.322  , 16233.451  , 18370.902  , 20415.227  , 23619.033  ,
                26214.322  , 30731.807  ]),
    dims=("level"),
    coords={
        "level": np.array([ 1000.,  925.,  850.,  700.,  600.,  500.,  400.,  300.,  250.,
                            200.,  150.,  100.,   70.,   50.,   30.,   20.,   10.])
    },
)
uwnd_data = xr.DataArray(
    np.array([  -2.3200073, -3.5700073, -2.5800018,  8.080002 , 14.059998 ,
                22.119995 , 33.819992 , 49.339996 , 57.86     , 64.009995 ,
                62.940002 , 49.809998 , 31.160004 , 16.59999  , 10.300003 ,
                10.459991 ,  9.880005 ]),
    dims=("level"),
    coords={
        "level": np.array([ 1000.,  925.,  850.,  700.,  600.,  500.,  400.,  300.,  250.,
                            200.,  150.,  100.,   70.,   50.,   30.,   20.,   10.])
    },
)

Then we need to interpolate the zonal wind data (located on the isobaric surface) to the altitude layers.

target_heights = np.linspace(0, 10000, 101)

result = ecl.interp.interp1d_vertical_pressure2altitude(
    z_data=z_data,
    variable_data=uwnd_data,
    target_heights=target_heights,
    vertical_input_dim="level",
    vertical_output_dim="altitude",
)
result
<xarray.DataArray (altitude: 101)> Size: 808B
array([-1.621444  , -1.9583984 , -2.27662079, -2.57157517, -2.83872555,
       -3.07353594, -3.27147033, -3.42799274, -3.53856716, -3.5986576 ,
       -3.60372807, -3.54924256, -3.43066509, -3.24345966, -2.98309027,
       -2.64502092, -2.22597052, -1.73105426, -1.16883252, -0.54787616,
        0.12324395,  0.83595696,  1.58169199,  2.35187819,  3.13794469,
        3.93132063,  4.72343514,  5.50571737,  6.26959645,  7.00650152,
        7.70786171,  8.36531168,  8.97677015,  9.54745024, 10.08297407,
       10.58896375, 11.07104138, 11.53482909, 11.98594898, 12.43002317,
       12.87267376, 13.31952288, 13.77619264, 14.2482397 , 14.73858718,
       15.2466844 , 15.77174317, 16.31297529, 16.86959254, 17.44080673,
       18.02582965, 18.62387309, 19.23414887, 19.85586876, 20.48824458,
       21.13048811, 21.78181115, 22.44143654, 23.10885874, 23.7838541 ,
       24.46621207, 25.15572206, 25.85217352, 26.55535589, 27.26505859,
       27.98107108, 28.70318277, 29.4311831 , 30.16486152, 30.90400746,
       31.64841035, 32.39785963, 33.15214473, 33.91105503, 34.67432735,
       35.44155111, 36.21229007, 36.98610797, 37.76256858, 38.54123564,
       39.32167292, 40.10344416, 40.88611311, 41.66924355, 42.4523992 ,
       43.23514385, 44.01704122, 44.79765509, 45.5765492 , 46.35328731,
       47.12743317, 47.89855053, 48.66620316, 49.42995455, 50.18915041,
       50.94252257, 51.68869548, 52.42629358, 53.15394131, 53.8702631 ,
       54.57388341])
Coordinates:
  * altitude  (altitude) float64 808B 0.0 100.0 200.0 ... 9.8e+03 9.9e+03 1e+04


Now we can check the interpolation results.

plt.plot(z_data[:9], uwnd_data[:9], "o", label="Original")
plt.plot(result.altitude, result.data, label="Interpolated")
plt.xlabel("Altitude [m]")
plt.ylabel("Zonal Wind [m/s]")
plt.legend()
plot interp
<matplotlib.legend.Legend object at 0x7faa1730f6a0>

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