Interpolation on a multi-dimensional curvilinear grid

Maybe it's already doable with Interpolations.jl but I could not figure it out. I have a curvilinear grid of the globe with some data and I want to interpolate it (a fairly common task in Earth sciences). So, for the knots of the grid, I have two m×n arrays, `LAT` and `LON`. For the data I have another m×n array `V`.

I would like to be able to do something like

```julia
itp = interpolate((LAT, LON), V)
```

Here is a non-working minimal example to test (with a plot to illustrate the type of grid):

```julia
lat = -10:10
lon = -20:20
# The curvilinear grid:
LAT = [x + 10cos(y/40) for x in lat, y in lon]
LON = [y + 0.01x^3 for x in lat, y in lon]
# Just a plot to illustrate:
using Plots
lonrows = vcat(([x; NaN] for x in eachrow(LON))...)
loncols = vcat(([x; NaN] for x in eachcol(LON))...)
latrows = vcat(([x; NaN] for x in eachrow(LAT))...)
latcols = vcat(([x; NaN] for x in eachcol(LAT))...)
plot([lonrows; loncols], [latrows; latcols], lab="curvilinear grid")
# The data on that grid:
V = rand(size(LON)...)
# Now I would like to interpolate:
using Interpolations
itp = interpolate(LON, LAT, V) # throws no method matching interpolate(::Matrix, ::Matrix, ::Matrix)
itp = interpolate((LON, LAT), V) # throws no method matching interpolate(::Tuple{Matrix, Matrix}, ::Matrix)
```

The curvilinear grid of this example:

But maybe this is for another package than Interpolations.jl?