x lines of Python: Gridding map data
Difficulty rating: moderate.
Welcome to the latest in the X lines of Python series. You probably thought it had died, gawn to ‘eaven, was an x-series. Well, it’s back!
Today we’re going to fit a regularly sampled surface — a grid — to an irregular set of points in (x, y) space. The points represent porosity, measured in volume percent.
Here’s what we’re going to do; it all comes to only 9 lines of code!
Load the data from a text file (needs 1 line of code).
Compute the extents and then the coordinates of the new grid (2 lines).
Make a radial basis function interpolator using SciPy (1 line).
Perform the interpolation (1 line).
Make a plot (4 lines).
As usual, there’s a Jupyter Notebook accompanying this blog post, and you can run it right now without installing anything.
Run the accompanying notebook in MyBinder
Run the notebook in Google Colaboratory
Just the juicy bits
The notebook goes over the workflow in a bit more detail — with more plots and a few different ways of doing the interpolation. For example, we try out triangulation and demonstrate using scikit-learn’s Gaussian process model to show how we might use kriging (turns out kriging was machine learning all along!).
If you don’t have time for all that, and just want the meat of the notebook, here it is:
This results in the following plot, in which the points are the original data, plotted with the same colourmap as the surface itself (so they should be the same colour, more or less, as their background).