Interpolation gridded data to geographical point location

Question

I am a big fan of MetPy and had a look at their interpolation functions (https://unidata.github.io/MetPy/latest/api/generated/metpy.interpolate.html) but could not find what I was looking for.

I am looking for a function to interpolate a gridded 2D (lon and lat) or 3D (lon, lat and vertical levels) climate data field to a specific geographic location (lat/lon).

The function would take 5 arguments: a 2D/3D data variable and associated latitude and longitude variables, as well as the two desired latitude and longitude coordinate values. Returned is either a single value (for 2D field) or a vertical profile (for 3D field).

I am basically looking for an equivalent to the old Basemap function bm.interp(). Cartopy does not have an equivalent. The CDO (Climate Data Operators) operator 'remapbil,lon=/lat=' does the same thing but works directly on netCDF files from the command line, I'm looking for a Python solution.

I think such a function would be a useful addition to the MetPy library as it allows for comparing gridded data (e.g., model or satellite data) with point observations such as from weather stations or radiosonde profiles (treated as just a vertical profile here).

Can you point me in the right direction?

Answer 1

I think what you're looking for already exists in scipy.interpolate (scipy is one of MetPy's dependencies). Here we can use interpn to interpolate linearly in n dimensions:

import numpy as np
from scipy.interpolate import interpn

# Array of synthetic grid to interpolate--ordered z,y,x
a = np.arange(24).reshape(2, 3, 4)

# Locations of grid points along each dimension
z = np.array([1.5, 2.5])
y = np.array([-1., 0., 1.])
x = np.array([-3.5, -1, 1, 3.5])

interpn((z, y, x), a, (2., 0.5, 2.))

Answer 2

This can be done easily with my nctoolkit package (https://nctoolkit.readthedocs.io/en/latest/). It uses CDO as a backend, and defaults to bilinear interpolation. The following would regrid a .nc file to a single grid point and then convert it to an xarray dataset.

import nctoolkit as nc
import pandas as pd
data = nc.open_data("example.nc")
grid = pd.DataFrame({"lon":[0], "lat":[50]})
data.regrid(grid)
ds = data.to_xarray()

Answer 3

To add one more solution, if you're already using multidimensional netCDF files and want a Python solution: check out xarray's interpolation tools. They support multidimensional, label-based interpolation with usage similar to xarray's indexing interface. These are built on top of the same scipy.interpolate otherwise mentioned, and xarray is also a MetPy dependency.