Interpolation: AbstractGrid -> arbitrary grid

[milankl]

@hottad I felt today like writing an Interpolator to get from any Grid<:AbstractGrid to an arbitrary grid that's just defined by two vectors of latitude/longitude points. This should allow us to have a semi-Lagrangian advection scheme and also allows us to get from one grid to another without going through spectral space. E.g.

julia> A = zeros(HEALPix4Grid,24)     # start with a P24 HEALPix4Grid 
julia> A.data .= rand(length(A.data)) # and some random data on it

julia> npoints = 100                  # lets interpolate onto 100 points
julia> θs = 20*randn(npoints)         # some latitudes in [-90˚,90˚N]
julia> λs = 360*rand(npoints)         # some longitudes in [0˚,360˚E]
julia> SpeedyWeather.interpolate(θs,λs,A)
100-element Vector{Float64}:
 0.5777433235930121
 0.9188505312368747
 0.702817220557382
 0.15673643528695175
 0.7337283563851287
 ...

Done. 🥳

What happens under the hood is that we create I::Interpolator which knows about the geometry of A::AbstractGrid and contains the arrays that define which grid points of A have to be used to interpolate onto all the θs,λs and with which weights

julia> I = Interpolator(A,npoints)    # initialize an Interpolator struct
julia> SpeedyWeather.update_locator!(I,θs,λs)    # translate θ,λ to grid indices and weights to interpolate from

Now the Interpolator can be used but also reused in case several different fields should be interpolated onto the same grid

julia> B = zeros(npoints)                  # array to write into
julia> SpeedyWeather.interpolate!(B,A,I)   # interpolate A onto the arbitrary grid in I and store in B
100-element Vector{Float64}:
 0.5777433235930121
 0.9188505312368747
 0.702817220557382
 0.15673643528695175
 0.7337283563851287
...

I still need to test that everything works and it should already be type-stable and flexible. We could even put a interpolate!(B::AbstractGrid,A::AbstractGrid) on top of these functions and can then conveniently jump between all our grids.

[hottad]

Great work, Milan!
You implemented bilinear interpolation based on anvil-shaped averaging, which will be very useful, but depending on how we use it, we may need to implement more elaborate interpolation methods.
For example, models like IFS use quasi-cubic interpolation (cubic interpolation with 4 stencil points in the longitudinal direction for the nearest two latitudes, linear interpolation with 2 stencil points in the longitudinal direction for the farther two latitudes, and then cubic interpolation in the latitude direction, which uses 12 stencil points in total).
JMA-GSM even uses quintic interpolation for the winds for the sake of accuracy.

Also, for computational performance, instead of anvil averaging, you could consider using the two-step interpolation (i.e., doing longitude interpolations for the adjacent rings, then doing latitudinal interpolation using the interpolated value that align on the same longitude).

Perhaps we can introduce AbstractInterpolator type and make the current bilinear interpolator a default one, but we can do this later when it becomes necessary.

[milankl]

We may need to implement more elaborate interpolation methods.

Yes, absolutely. I thought I start easy though and I like your idea to introduce AbstractInterpolator and then we can generalise this in the future. I'll keep the name Interpolator for now, although we should probably give it names like LinearInterpolator, CubicInterpolator etc. in the future.

Also, for computational performance, instead of anvil averaging, you could consider using the two-step interpolation (i.e., doing longitude interpolations for the adjacent rings, then doing latitudinal interpolation using the interpolated value that align on the same longitude).

The latter is what I've implemented here, I've tried to sketch this in the docstring

        0..............1    # fraction of distance Δab between a,b
        |<  Δab   >|
0^      a -------- o - b    # anvil-shaped average of a,b,c,d at location x
.Δy                |
.                  |
.v                 x 
.                  |
1         c ------ o ---- d
          |<  Δcd >|
          0...............1 # fraction of distance Δcd between c,d

So first we interpolate onto the os which share the same longitude determined by x, then those are interpolated onto x. I've just called it "anvil interpolation" as it looks like an anvil to me. Is that already an otherwise used term?

[milankl]

Re: Generalisation. We now have

abstract type AbstractLocator{NF} end
struct AnvilLocator{NF<:AbstractFloat} <: AbstractLocator{NF}
abstract type AbstractInterpolator{NF,G} end
struct AnvilInterpolator{NF<:AbstractFloat,G<:AbstractGrid} <: AbstractInterpolator{NF,G}
    geometry::GridGeometry{G}
    locator::AnvilLocator{NF}
end

So only a NewLocator <: AbstractLocator would need to be defined for higher order interplation scheme. As well as a new method for the find_grid_indices function, which is essentially the longitude points. And the actual averaging function of course. Easy 😆

[hottad]

Thank you, Milan! and sorry that I commented without carefully reading the code.
I am not aware of any use of "anvil" interpolator etc. Nor a standard terminology for the two-step interpolation method that you implemented here.

I am always impressed with the elegance of the code that you write in Julia. Thank you!

[milankl]

I am always impressed with the elegance of the code that you write in Julia. Thank you!

Ha! Thanks but I feel a lot of elegance credit also goes to Julia.