Arrays of interpolated values

[t-ohno0723]

When analyzing the structure of a tropical cyclone, there are cases where we want to calculate the structure averaged in azimuth relative to the center of the tropical cyclone by taking points in concentric circles and averaging the values at those points. How can we obtain an array of interpolated values for an array of arbitrary coordinate points?

[tinaok]

Here is the drawing of what @t-ohno0723 wish to do using xarray.

How can we define point A (black cross in the drawing, not necessarily grid center of healpix cell) and compute the mean values of points (in this drawing, little 10-11 circles in 10km from the center points).

[jmdelouis]

This is a convolution, so it is perfectly well-defined mathematically, with weights corresponding to the intersection area between the circular kernel and each pixel. However, it introduces quantization noise due to the fact that pixel sizes are relatively large compared to the radius of the averaging kernel. Using a sharp threshold to compute a simple mean may not be the best approach.

A better method would be to first convolve the sphere with a Gaussian beam of 10 km and then extract values using bilinear interpolation. This requires just two calls to the healpy library — healpy.smoothing for the convolution, and healpy.get_interp_val for interpolation — if you're working over the full sphere. Since the convolution is performed in spherical harmonic space, it is expected to reduce pixel-induced quantization effects (notably the oscillations due to the sinc function).

If you're working on a local domain in HEALPix, you can use the healpix_convolve library instead.

[tinaok]

Could be nice to create an example and add that in updated tutorial? #138

[tinaok]

@jmdelouis
What if we have use case, requires just one single point? (like comparison with in-situ observation)

How can we retrieve a single value at a given (lat, lon) location — not the center of a HEALPix pixel — by interpolating from neighboring pixels with weights using healpy?

I think it would be very useful to have this as a utility function for geo-users: when the data is stored in Zarr/HEALPix format, the function would extract the appropriate chunk(s), perform the interpolation, and return the value.

Maybe this is something to be implemented in healpy-geo? @benbovy @keewis

[benbovy]

A comment about the API:

we have discussed a bit with @keewis and @scottyhq during our small sprint on Xarray indexes about how to reuse Xarray's interp() with custom indexes. Right now interp() is relying on numpy or scipy.interpolate functions.

There are many ways to interpolate data and the best way highly depends on the application and the kind of data structure. The solutions suggested here above by @jmdelouis are a nice illustration, they work very well with Healpix but couldn't always work with other DGGSs that aren't designed to perform operations in the harmonic space.

A possible way to deal with this may be to provide in Xarray a NamedInterpolator base class. It would be very similar to scipy.interpolate's Interpolator classes but would work with Xarray Dataset and DataArray objects instead of unlabelled array-like objects. Then in xdggs we could provide more-or-less specialized classes like DGGSInterpolator (DGGS-aware but DGGS-agnostic), HealpixInterpolator, etc. (EDIT: a NamedInterpolator base class may not be needed, actually)