Kriging Interpolation

[hazardgoat]

Description of the desired feature
I would like to request that an ordinary kriging interpolation method be added to PyGMT.

While several python libraries exist that can perform kriging interpolation, they are confusing to implement for geospatial data and/or are limited in the file formats they can export to. A kriging interpolation method would allow for easy visualization of an interpolation of readings at scattered stations, and if the ability to export to something like netcdf, geotiff, or shapefile is also implemented, then that would be ideal.

Are you willing to help implement and maintain this feature? Yes/No
I am not technically skilled enough to implement or maintain this feature, but I would do my best to help however I can. I aspire to be knowledgeable enough to assist meaningfully.

[weiji14]

Hi @hazardgoat, so the thing is, kriging isn't implemented in the GMT C library, so it will be tricky for PyGMT to include this feature 😅 Not that we can't do it, but it would require either:

1. Have GMT C implement kriging at https://github.com/GenericMappingTools/gmt, and we can then wrap it in PyGMT, or
2. Implement this purely in Python/PyGMT, which I'm not against, but it's worth raising with @GenericMappingTools/pygmt-maintainers if this is something to bring into the scope of the project.

While several python libraries exist that can perform kriging interpolation, they are confusing to implement for geospatial data and/or are limited in the file formats they can export to.

Could you link to some of those projects? I found https://stackoverflow.com/questions/45175201/how-can-i-interpolate-station-data-with-kriging-in-python which mentions pykriging and scikit-learn's GaussianProcessRegressor (apparently another term for kriging).

However, if you're interested in just 'gridding', GMT does offer a few options at https://docs.generic-mapping-tools.org/6.2/modules.html#gridding. The status of these functions in PyGMT are:

• greenspline (no PR) | Interpolate using Green’s functions for splines in 1-3 dimensions
• nearneighbor (Wrap nearneighbor #1379) | Grid table data using a "Nearest neighbor" algorithm
• sphinterpolate (Wrap sphinterpolate #1418) | Spherical gridding in tension of data on a sphere
• surface (Wrap surface #245) | Grid table data using adjustable tension continuous curvature splines
• triangulate (Wrap triangulate #731) | Delaunay triangulation or Voronoi partitioning and gridding of Cartesian data
• grdinterpolate (no PR) | Interpolate a 3-D cube, 2-D grids or 1-D series from a 3-D data cube or stack of 2-D grids

So far, only 1 is available (surface) in PyGMT v0.4.1, but 3 are near completion (nearneighbor, sphinterpolate, triangulate) and should be available in PyGMT v0.5.0.

I guess the main question is, what do you want to achieve? Each of these functions do slightly different things, and depending on whether your input data is Cartesian/Geographic, and what type of output you want to get, the choice of which algorithm to use can be different.

P.S. You can read @leouieda's paper to get a brief overview of gridding:

• Uieda, (2018). Verde: Processing and gridding spatial data using Green’s functions. Journal of Open Source Software, 3(30), 957, https://doi.org/10.21105/joss.00957

[hazardgoat]

@weiji14

Thanks for such a thorough reply! I've come across those examples you linked, as well as these ones for gstools and xarray. While I've been able to get the gstools one to work, it only exports to .vtr files, so I'm not sure how to get a raster or shapefile out of it which I need for plotting in PyGMT and investigating with other tools.

Honestly, I'm certain my own ignorance is largely to blame, but it's also what motivates me to request the feature in PyGMT. The way methods are implemented in PyGMT, as well as how helpful and friendly the PyGMT community is, is non-trivial with respect to ease of learning. My rational is that if kriging interpolation were a method in PyGMT then I would be able to use it in the way I need.

I admit, I'm not certain what griding is. The paper you linked describes creating a regular grid, but does that include all the in-between points, or just repositioning the existing points?

I've been using the grdtrack method for data analysis, so my view of PyGMT is one of a mapping tool that is also capable of investigating datasets. The motivating task of my request is that I have a bunch of scattered stations at geographic coordinates with z-direction observations. I want to make, visualize, and export an interpolation of those observations that interpolates all the empty space between them. The reason I want to use kriging interpolation is that I've read kriging is ideal for statistical analysis as it can produce an error map.

I hope I was able to clarify and answer your question :)

[leouieda]

My 2 cents here: kriging is not a simple thing to implement and it would be a considerable effort to include it in PyGMT. It would also mean maintaining a large chunk of numerical code, which is quite different from the rest of the package. So my view is that it’s not worth doing it here.

There are Python packages exclusively for kriging out there as well. PyKrige seems to be the most mature and well developed. I don’t think we could easily match their capabilities and frankly it would be best to improve their software through PRs than to try to reimplement everything. Kriging is also complex to use and requires a good understating of the underlying method to do so correctly. Much more so than the interpolation we have in GMT/PyGMT.

[seisman]

Based on the above discussions, I think PyGMT won't implement Kriging interpolation unless GMT does. So the options are:

• 🎉 Close the feature request
• 👀 Transfer the feature request to GMT

Please place your vote on this comment.

[seisman]

Closing the feature request as it won't be implemented.