Least worse way to interpolate a few data points to a gridded region?

[jypeter]

Question criteria

• I added a descriptive title here.
• I searched the xCDAT GitHub Discussions to find a similar question and didn't find it.
• I searched the xCDAT documentation and/or Xarray documentation.

Describe your question

This may not be the right place to ask, but my question is not covered by the Horizontal Regridding example.

I talked with one of our PhD students this morning, and he told me that he was using scipy (scipy.interpolate.griddata I think) to interpolate data from a few scattered observation points to a gridded domain. This did not seem like the correct way to go. Although the target gridded area was small (some sub-region of the French Alps), so the fact that scipy does not know it is interpolating data on a sphere probably does not matter too much.

I feel this may be linked to the upscaling/downscaling/kriging bag of worms that we possibly don't want to open (and that I don't really understand), but maybe having a few hints for students (and others) could be useful

Thanks!

Are there are any possible answers you came across?

No response

Minimal Complete Verifiable Example (MVCE)

No response

Relevant log output

No response

Environment

No response

Anything else we need to know?

No response

[mfwehner]

This paper may be of use. Code is all in R however. Mark Risser and Michael Wehner (2020) The effect of geographic sampling on evaluation of extreme precipitation in high resolution climate models. Advances in Statistical Climatology, Meteorology and Oceanography 6, 115–139, https://doi.org/10.5194/ascmo-6-115-2020

[jypeter]

Thanks @mfwehner for the reference! I will forward it to people who do model-data comparison

But here, I'm more looking for some Python recipe that would allow interns or beginning PhD students to interpolate their sample/observation data points to some grid, in a way that would be slightly less quick and dirty than scipy.interpolate.griddata, without even quantifying the error.

I know it would still be quick and dirty, but it would be an improvement

[tomvothecoder]

Hi @jypeter, I'm tagging @jasonb5, @pochedls, and @chengzhuzhang to see if they might have a solution since I'm not skilled enough in this area.

[pochedls]

I have used cKDTree, which is fast (but I have also found its syntax a little confusing).

[jypeter]

@pochedls Thanks! Finding out the closest model point to a data point (preferably on a sphere) is also a use case we have from time to time. I did that a long time ago in a quick and dirty way (use array syntax to quickly compute the straight-line distance between the data point and all the model points, assuming they are not too far way, and get the index of the point where the distance is smallest)

But my question here is more: interpolating to go from a few scattered data points to some (regular or irregular) grid

[pochedls]

I think ckdtree helps to do this.

This is some pseudo-code I've extracted when I've done this in the past. I believe I looped over each time point and mapped scattered data to a regular lat / lon grid:

# import
from scipy import spatial

# user-data
lat = array-like set of latitude points
lon = array-like set of longitude points
data = array-like set of values

# target grid
LON, LAT = np.meshgrid(np.arange(1.25, 360, 2.5), np.arange(-88.75, 90, 2.5))
shp = LAT.shape
# Create a nearest-neighbor tree for the grid
tree = spatial.cKDTree(list(zip(LAT.flat, LON.flat)))

dis, ind = tree.query(np.array([lat, lon]).T)
n = tree.n
X = np.bincount(ind, weights=data,  minlength=n) # Sum of data in each grid box
cnt = np.bincount(ind, weights=np.ones_like(data), minlength=n) # Total number of matches
with np.errstate(divide='ignore', invalid='ignore'):
    grid = X / cnt # going to get divide by zero here for grid boxes with no data
grid = grid.reshape(shp) # reshape to regular grid

I should note that this discussion is beyond the scope of xcdat, though this seems like it could be a reasonable feature to include in the future (if it doesn't exist already).

[pochedls]

@jypeter : I added an example of this here for E3SM model data. This has limitations – it is nearest neighbor, only works for 2D fields, and does some simple interpolation / extrapolation (you can turn this off). But this could be useful for quick-data looks.