[FR] LonLat to LonLat interpolation

[Datseris]

Would be nice to have interpolation from one LonLatGrid to another type of LonLatGrid. I believe this will be relatively easy using https://github.com/JuliaMath/Interpolations.jl .

We first create a new copy array, whose other dimensions are the same as in ClimArray while the Lon Lat dims are the new ones according to user input.
We can iterate over other dimensions of the ClimArray using oidxs = otheridxs(A, (Lon(), Lat())) same way as it is done in the source code of spacemean. For each Lon x Lat slice we do a bilinear/bicubic interpolation using Interpolations.jl.

Done!

[Datseris]

@Balinus what kind of Julia package did you use for implementing interpolation in the current ClimateTools.jl?

[Balinus]

There's code already in ClimateTools for that. Initially, I used Interpolations.jl but I had some problems with irregular grids (interpolations from a global model onto a regional climate models for example). So, I added a Python dependency scipy.griddata (that was straightforward and rather quick) so that I was able to continue my work while looking for another exhaustive approach.

Hence, my next goal was to properly implement GeoStats.jl. There's also code with that in ClimateTools. The idea was to define a custom solver (e.g. Kriging) wikth the ability to let the user use a custom solver from the large collection of methods from GeoStats.jl. The implementation of GeoStats.jl is not fully implemented but functional if I recall correctly.

See lines between 170-427.
https://github.com/JuliaClimate/ClimateTools.jl/blob/master/src/functions.jl

@juliohm

edit- note that in functions.jl, there's useful code for polygon clipping based on shapefiles.

[Datseris]

OKay, I see, thanks. This is too advanced for me so I think I need to leave it to you to adjust the codebase you already have to the new ClimArray, which shouldn't be too much work.

[juliohm]

Thank you @Balinus for considering the integration with the GeoStats.jl stack. We are constantly evolving the API to handle more complex domains. In the latest release we support any domain as defined by Meshes.jl, which is already a nice step forward. If you feel that a domain is missing for lat-lon coordinates, Meshes.jl is the appropriate place to submit a PR.

Regarding the GeoStats.jl solvers, we have various methods implemented as you pointed out. It would be nice to know what you felt was missing in the solvers. Contributions are very welcome and we also value issues. We cannot fix what is not in our radar :)

Please let me know if can I still be of any help here. I try to document the GeoStats.jl stack as much as possible, but it is really challenging to implement the features and add the docs at the same time as you know.

[Balinus]

Thanks @juliohm for the answers. I do not feel that there is a missing solver. It's just that I changed jobs last year and I never finished the integration with GeoStats :)

Interesting to know the integrations with Meshes.jl. Perhaps there is some fusion possible with https://github.com/JuliaClimate/MeshArrays.jl

MeshArrays.jl seems focused on "Gridded earth variables, domain decomposition, and climate model C-grid support".

[juliohm]

It would be nice to figure out if MeshArrays.jl can follow the Meshes.jl API or if we could migrate their mesh types to Meshes.jl for greater visibility, maintainability, etc. GeoStats.jl already supports the Meshes.jl API and we are working actively there to make sure that things integrate more broadly.

[Balinus]

Pinging @gaelforget, the developer of MeshArrays.jl

[gaelforget]

Revisiting the internals of MeshArrays.jl and the inner array types especially has been on my to do list for a while, and I am hoping to get around to it soon... Have not looked at Meshes.jl yet but it's been on my radar in that context.

By the way, the interpolation scheme in MeshArrays.jl avoids using lon and lat since that does not work for distorted grids, or even for a cube sphere. Instead I use nearest neighbors in 3D, grid cells as polygons, and local stereo projections taking advantage of the C-grid setup. Not sure how relevant this is to this discussion, but seems like it could be.

Side question for @juliohm if I may : I might create a wrapper for https://github.com/sahrk/DGGRID , or just artifact it, for ocean stats that I tend to compute on cube sphere grids and such (see MeshArrays.jl). Is there already a Julia wrapper / build for that? A pure Julia version? A good place / package to put something like this within JuliaGeo or JuliaEarth packages?

ps. thanks for ping @Balinus

[juliohm]

Thank you for pinging @gaelforget, and sorry for the delay. I was in a workshop last week and couldn't find the time to reply.

It would be really nice if we could use/evolve a common mesh interface here. GeoStats.jl v0.24 already uses Meshes.jl's meshes, which means that you can do interpolation, simulation, learning directly on these domain discretizations. If you could port MeshArrays.jl to Meshes.jl or at least implement the traits of that package, we would be a in a much better situation regarding data manipulation and statistics on spherical/geoid domains.

In this long-term vision we want to avoid different APIs for different types of domains. We would like users to define geospatial problems on any kind of mesh, no matter if it is a global grid of hexagons, a Cartesian grid of quadrangles, a Bezier curve embedded in 3D, etc. We are working actively on these projects and more contributors are more than welcome.

Check our meshes.jl and geostats.jl streams in Zulip if you want to discuss further details.

[gaelforget]

Thank you for pinging @gaelforget, and sorry for the delay. I was in a workshop last week and couldn't find the time to reply.

It would be really nice if we could use/evolve a common mesh interface here. GeoStats.jl v0.24 already uses Meshes.jl's meshes, which means that you can do interpolation, simulation, learning directly on these domain discretizations. If you could port MeshArrays.jl to Meshes.jl or at least implement the traits of that package, we would be a in a much better situation regarding data manipulation and statistics on spherical/geoid domains.

In this long-term vision we want to avoid different APIs for different types of domains. We would like users to define geospatial problems on any kind of mesh, no matter if it is a global grid of hexagons, a Cartesian grid of quadrangles, a Bezier curve embedded in 3D, etc. We are working actively on these projects and more contributors are more than welcome.

Check our meshes.jl and geostats.jl streams in Zulip if you want to discuss further details.

Sounds like a lot of good ideas but also possibly a fair amount of work. My todo list for this spring isn't looking like I can add this to it. Maybe I could suggest this as an undergraduate project here this fall. Need to think about this more and spend more time on zulip too. Anyway, thanks a lot for suggesting this and sorry for the delay in response

[Balinus]

Here's a working example for interpolations using GeoStats.jl :

JuliaEarth/GeoStats.jl#22