2D interpolation routines
Closed this issue · 3 comments
Currently 2D interpolation uses the radial basis function method of scipy.interpolate.Rbf which offers several out of the box options:
'multiquadric': sqrt((r/self.epsilon)**2 + 1)
'inverse': 1.0/sqrt((r/self.epsilon)**2 + 1)
'gaussian': exp(-(r/self.epsilon)**2)
'linear': r
'cubic': r**3
'quintic': r**5
'thin_plate': r**2 * log(r)And also allows for user defined functions to be passed.
The problem is that most of these methods can result in gradients being overestimated and large local extrema. I have been experimenting with its smoothing parameter which is not well documented. It seems that the "linear" option with smoothing=-1e-3 is the best combination so far testing on the UC and example data. There may be a way to scale the smoothing and epsilon parameters of RBF methods adaptively to deal with different datasets, i.e. sparse or dense.
For now I am going to change the default to Rbf "linear", add the smoothing parameter option with default -1e-3. Results are reasonable.
inverse distance weighting
Apparently the "linear" Rbf behaves somewhat similar to inverse distance weighting but can be quite a bit different. Unfortunately there is no IDW in Python as far as I have found but I have found example codes I am trying to test and modify for our use, see: here and here.
Also I gdal has a command line tool. On the surface (haven't tested it out yet) it has some major advantages: it interpolates and makes the raster in one step with parallelization options, it has several other useful methods like IDW K nearest neighbors. Although it would require a couple preliminary steps like writing a CSV file and a VRT meta file but they look straightforward.
other options and comparisons
I have also tested scipy.interpolate.griddata and it is fairly straightforward however it has a few disadvantages:
- only interpolates within inner bounding from points (convex hull)
- only has 3 methods: nearest, linear, and cubic
see:
https://stackoverflow.com/questions/50816375/scipy-griddata-with-linear-and-cubic-yields-nan
This post made me initially lean towards RBFs:
https://stackoverflow.com/questions/37872171/how-can-i-perform-two-dimensional-interpolation-using-scipy?rq=1
Thanks for the info John, this is great. I like the idea of using the GDAL command line tool. Ya - that's what I don't like much about splines and the like that overestimate the gradient at extremes. Ultimately we'll probably make the decision based on what looks best - nothing like human eyes and mind!
The last commit a7e9cc8 added a Python version of inverse distance weighting, modified from a post by Roger Rovira who also happens to have a lot of great blogs about geospatial stuff. However the gdal_grid option I also got working (windows and Linux) and it has huge advantages. I'm building a class to interface between gridwxcomp and command line usage of gdal_grid. The main two advantages are speed (incredibly faster even without parallelization) and it has several other interpolation methods and options. I think we should keep this Python version around as a backup, it works fine for small regions or with a large scale factor.
Last commit added all of gdal interpolation methods, still needs to be merged into the spatial.py main workflow, at that point the Python version for inverse distance weighting to a power will be removed. However the gdal algorithms can be used now within Python and were tested on multiple datasets. Example usage of the new class gridwxcomp.interpgdal.InterpGdal can be found in the docstrings in interpgdal.py.