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.
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:
import numpy as np import pandas as pd import matplotlib.pyplot as plt from scipy.interpolate import Rbf # Load the data. df = pd.read_csv('../data/ZoneA.dat', sep=' ', header=9, usecols=[0, 1, 2, 3], names=['x', 'y', 'thick', 'por'] ) # Build a regular grid with 500-metre cells. extent = x_min, x_max, y_min, y_max = [df.x.min()-1000, df.x.max()+1000, df.y.min()-1000, df.y.max()+1000] grid_x, grid_y = np.mgrid[x_min:x_max:500, y_min:y_max:500] # Make the interpolator and do the interpolation. rbfi = Rbf(df.x, df.y, df.por) di = rbfi(grid_x, grid_y) # Make the plot. plt.figure(figsize=(15, 15)) plt.imshow(di.T, origin="lower", extent=extent) cb = plt.scatter(df.x, df.y, s=60, c=df.por, edgecolor='#ffffff66') plt.colorbar(cb, shrink=0.67) plt.show()
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).