<jupyter_start><jupyter_text>X lines of Python Gridding map dataThis notebook goes with the [Agile Scientific blog post from 8 March 2019](https://agilescientific.com/blog/2019/3/8/x-lines-of-python-gridding-map-data).I'm using a small dataset originally from [**<NAME>**](http://people.ku.edu/~gbohling/) at the Kansas Geological Survey. I can no longer find the data online.We will look at four ways to do this:- Using SciPy with `scipy.interpolate.Rbf`- Using SciPy with `scipy.griddata()`- Using the Scikit-Learn machine learning library with `sklearn.gaussian_process`- Using the [Verde](https://github.com/fatiando/verde) spatial gridding library with `verde.Spline` (thank you to Leo Uieda for contributing this!) Load and inspect the data<jupyter_code>import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt<jupyter_output><empty_output><jupyter_text>We'll load the data from the web; if you are offline but have the repo, the data file is also in `../data/ZoneA.dat`.<jupyter_code># Essential line number 1.
df = pd.read_csv('https://www.dropbox.com/s/6dyfc4fl5slhgry/ZoneA.dat?raw=1',
                 sep=' ',
                 header=9,
                 usecols=[0, 1, 2, 3],
                 names=['x', 'y', 'thick', 'por'],
                 dtype="float64"
                )
df.head()
df.describe()
sns.distplot(df.por)<jupyter_output><empty_output><jupyter_text>This looks a bit unpleasant, but we're just getting out min and max values for the x and y columns.<jupyter_code># Line 2.
extent = x_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]<jupyter_output><empty_output><jupyter_text>Later on, we'll see a nicer way to do this using the Verde library.Now we can plot the data:<jupyter_code>fig, ax = plt.subplots(figsize=(10,6))

ax.scatter(df.x, df.y, c=df.por)
ax.set_aspect(1)
ax.set_xlim(*extent[:2])
ax.set_ylim(*extent[2:])
ax.set_xlabel('Easting [m]')
ax.set_ylabel('Northing [m]')
ax.set_title('Porosity %')
ax.grid(c='k', alpha=0.2)

plt.show()<jupyter_output><empty_output><jupyter_text>Make a regular gridWe must make a grid, which represents the points we'd like to predict.<jupyter_code># Line 3.
grid_x, grid_y = np.mgrid[x_min:x_max:500, y_min:y_max:500]

# Use *shape* argument to specify the *number* (not size) of bins:
# grid_x, grid_y = np.mgrid[x_min:x_max:100j, y_min:y_max:100j]
plt.figure(figsize=(10,6))
plt.scatter(grid_x, grid_y, s=10)<jupyter_output><empty_output><jupyter_text>That was easy! Interpolation with radial basis functionNow we make an interpolator and use it to predict into the grid 'cells'.<jupyter_code>from scipy.interpolate import Rbf

# Make an n-dimensional interpolator. This is essential line number 4.
rbfi = Rbf(df.x, df.y, df.por)

# Predict on the regular grid. Line 5.
di = rbfi(grid_x, grid_y)<jupyter_output><empty_output><jupyter_text>Let's plot the result. First, we'll need the min and max of the combined sparse and gridded data, so we can plot them with the same colourmap ranges:<jupyter_code>mi = np.min(np.hstack([di.ravel(), df.por.values]))
ma = np.max(np.hstack([di.ravel(), df.por.values]))<jupyter_output><empty_output><jupyter_text>Notice the transpose and the `origin='lower'`, to keep everything matched up with the original dataset.<jupyter_code>plt.figure(figsize=(15,15))
c1 = plt.imshow(di.T, origin="lower", extent=extent, vmin=mi, vmax=ma)
c2 = plt.scatter(df.x, df.y, s=60, c=df.por, edgecolor='#ffffff66', vmin=mi, vmax=ma)

plt.colorbar(c1, shrink=0.67)
plt.show()<jupyter_output><empty_output><jupyter_text>The circles (the data) are the same colour as the grid (the model), so we can see that the error on this prediction is almost zero. In fact, the default parameters force the model to pass through all the data points (interpolation, as opposed to estimation or approximation). The `Rbf()` interpolator has [a few options](https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.Rbf.html). The most important one is probably `smooth`, which is the thing to increase if you end up with a singular matrix (because it can't converge on a solution). Anything above 0 relaxes the constraint that the surface must pass through every point. If you get an error, you probably need to change the smoothing.You can also change the `function` (default is `multiquadric`, which also has an `epsilon` parameter to vary the range of influence of each point).<jupyter_code>rbfi = Rbf(df.x, df.y, df.por, smooth=0.2)
di = rbfi(grid_x, grid_y)

plt.imshow(di.T, origin="lower", extent=extent)
plt.scatter(df.x, df.y, s=2, c='w')
plt.show()<jupyter_output><empty_output><jupyter_text>We can also make a histogram and kernel density estimation of the errors, by making predictions at the original input locations:<jupyter_code>por_hat = rbfi(df.x, df.y)

sns.distplot(por_hat - df.por)<jupyter_output><empty_output><jupyter_text>With the smoothing set to 0.2, we end up with a smoother surface, but pay for it with larger errors.  Interpolation with `scipy.griddata()` The `Rbf()` interpolator is the one to know about, because it has lots of useful parameters. It's probably the only one you need to know. But there is also `scipy.griddata()`. For example see [this SciPy recipe](https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.interpolate.griddata.html).The interface is slightly different &mdash; we have to pass a single array of the coordinates (the (x, y) locations of the points we know). We also pass the values to interpolate, and the grids.The function will not accept Pandas `Series` objects, so we'll use the `Series.values` attribute to get at the NumPy array representation.First, let's make the 2D array of coordinates:<jupyter_code>points = df[['x', 'y']].values<jupyter_output><empty_output><jupyter_text>The grdding step is easy. We'll try three different algorithms:<jupyter_code>from scipy.interpolate import griddata

grid_z0 = griddata(points, df.por.values, (grid_x, grid_y), method='nearest')
grid_z1 = griddata(points, df.por.values, (grid_x, grid_y), method='linear')
grid_z2 = griddata(points, df.por.values, (grid_x, grid_y), method='cubic')<jupyter_output><empty_output><jupyter_text>Inspect the results.<jupyter_code>fig, axs = plt.subplots(ncols=3, figsize=(15, 5))

ax = axs[0]
ax.imshow(grid_z0, origin='lower', extent=extent)
ax.scatter(df.x, df.y, s=2, c='w')
ax.set_title('Nearest')

ax = axs[1]
ax.imshow(grid_z1, origin='lower', extent=extent)
ax.scatter(df.x, df.y, s=2, c='w')
ax.set_title('Linear')

ax = axs[2]
ax.imshow(grid_z2, origin='lower', extent=extent)
ax.scatter(df.x, df.y, s=2, c='w')
ax.set_title('Cubic')

plt.show()<jupyter_output><empty_output><jupyter_text>I don't particularly like any of these results.  Using `verde.Spline` One of the options in scipy's `Rbf` interpolator is the "thin-plate" kernel. This is what the `verde.Spline` interpolator is based on but with a few modifications, like damping regularization to smooth the solution. It's similar to the `RBF` and `GaussianProcessRegressor` approach but Verde provides a more conenient API for gridding tasks.For example, we now have a nicer way to define `extent`, using `vd.pad_region()`:<jupyter_code>import verde as vd

extent = x_min, x_max, y_min, y_max = vd.pad_region(vd.get_region((df.x, df.y)), pad=1000)
spline = vd.Spline(mindist=2000, damping=1e-4)
spline.fit((df.x, df.y), df.por)<jupyter_output><empty_output><jupyter_text>To make a grid, use the `.grid` method of the spline:<jupyter_code>grid = spline.grid(region=extent, spacing=500)
grid<jupyter_output><empty_output><jupyter_text>This returns an [`xarray.Dataset`](http://xarray.pydata.org/en/stable/generated/xarray.Dataset.htmlxarray.Dataset) which can be easily plotted, saved to disk as netCDF, or used for computations. The coordinates for the grid are automatically generated and populated based on the desired region and spacing. The spacing is adjusted to fit the desired region exactly. Optionally, you can set `adjust="region"` to adjust the size of the region so that the spacing is exact. And plot the predicted grid with the input data using the same colourmap:<jupyter_code># Compute min and max of all the data:
mi = np.min(np.hstack([grid.scalars.values.ravel(), df.por.values]))
ma = np.max(np.hstack([grid.scalars.values.ravel(), df.por.values]))

# Plot it all.
plt.figure(figsize=(15,15))
im = plt.imshow(grid.scalars, origin='lower', extent=extent, vmin=mi, vmax=ma)
pts = plt.scatter(df.x, df.y, c=df.por, s=80, edgecolor='#ffffff66', vmin=mi, vmax=ma)
plt.colorbar(im, shrink=0.67)
plt.show()<jupyter_output><empty_output><jupyter_text>As before, we can compute the error by making a prediction on the original (x, y) values and comparing to the actual measured porosities at those locations:<jupyter_code>por_hat = spline.predict((df.x, df.y))

sns.distplot(por_hat - df.por)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12369182810064573, 'fraction_numerical': 0.014807392562903586, 'mean_word_length': 5.212309820193638, 'pattern_counts': {'":': 0, '<': 70, '<?xml version=': 0, '>': 70, 'https://': 5, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '17', 'dir': 'jupyter-structured-clean-dedup', 'id': '39122', 'license': 'Apache-2.0', 'n_tokens_mistral': 3260, 'n_tokens_neox': 3108, 'n_words': 955, 'path': 'Gridding_map_data.ipynb', 'repo_name': 'EvanBianco/geocomputing_demos', 'size': '8989'}