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?
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.))