<jupyter_start><jupyter_text>Interpolation using PyKrige This week we are going to interpolate our soil geochemistry dataset using `pykrige`. So to do that we first need to install it with `pip`<jupyter_code>! pip install pykrige<jupyter_output><empty_output><jupyter_text>This week we are using a lot of the code from last weeks Gaussian process. Mostly we are reusing the grid construction portion and plotting sections. So we will start with the same imports as last week<jupyter_code>import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import geopandas as gpd
from itertools import product
import contextily as ctx

%matplotlib inline<jupyter_output><empty_output><jupyter_text>Then let's read in the data with `geopandas`<jupyter_code>data = gpd.read_file("geochemistry_subset.shp")
data.drop(
    index=[3, 16], inplace=True
)<jupyter_output><empty_output><jupyter_text>And here is the code we used last week to build our grid from the bounds of our soil samples<jupyter_code>x_values = np.linspace(min(data.geometry.x), max(data.geometry.x), num=50)
y_values = np.linspace(min(data.geometry.y), max(data.geometry.y), num=50)
xy = np.array(list(product(x_values, y_values)))
X0, X1 = xy[:, 0].reshape(50, 50), xy[:, 1].reshape(50, 50)<jupyter_output><empty_output><jupyter_text>To use `pykrige` we just need to import `OrdinaryKriging`<jupyter_code>from pykrige.ok import OrdinaryKriging<jupyter_output><empty_output><jupyter_text>We feed in our x and y sample locations, the element concentration (`qvalue`), and pick a `variogram_model`, here I chose `'spherical'` because it fits the data better than other models. I also enable `verbose=True` to see what it's doing, and then `enable_plotting=True` to see what the variogram model looks like<jupyter_code>OK = OrdinaryKriging(
    data.geometry.x,
    data.geometry.y,
    data.qvalue,
    variogram_model="spherical",
    variogram_parameters= {'sill': 350, 'range': 100000, 'nugget': 0.01},
    verbose=True,
    enable_plotting=True,
)<jupyter_output>Plotting Enabled

Adjusting data for anisotropy...
Initializing variogram model...
Coordinates type: 'euclidean' 

Using 'spherical' Variogram Model
Partial Sill: 349.99
Full Sill: 350.0
Range: 100000
Nugget: 0.01<jupyter_text>Now that we have a variogram model we are set to make interpolation predictions on the grid. To do that we use `OK.execute` with `grid` and `x_values` and `y_values`<jupyter_code>predictions, variance = OK.execute("grid", x_values, y_values)<jupyter_output>Executing Ordinary Kriging...<jupyter_text>Use our contextily helper function for plotting on a map<jupyter_code># this cell creates a function that grabs the stamen terrain tiles
def add_basemap(
    ax, zoom, url="http://tile.stamen.com/terrain/tileZ/tileX/tileY.png"
):
    xmin, xmax, ymin, ymax = ax.axis()
    basemap, extent = ctx.bounds2img(
        xmin, ymin, xmax, ymax, zoom=zoom, url=url
    )
    ax.imshow(basemap, extent=extent, interpolation="bilinear")
    # restore original x/y limits
    ax.axis((xmin, xmax, ymin, ymax))<jupyter_output><empty_output><jupyter_text>And then we plot it up on top of the map with the sample locations to see how our model stacked up<jupyter_code>ax = data.plot(
    column="qvalue",
    vmin=0,
    vmax=50,
    figsize=(20, 10),
    legend=True,
    cmap="RdBu_r",
    alpha=0.7,
    edgecolor="k",
    markersize=90,
    zorder=2,
)
add_basemap(ax, zoom=9)  # add our basemap to the plot
im = ax.pcolormesh(
    X0,
    X1,
    predictions.T,
    cmap="RdBu_r",
    alpha=0.5,
    zorder=1,
    edgecolor=(1.0, 1.0, 1.0, 1.0),
    linewidth=0.01,
)
plt.colorbar(im, ax=ax)<jupyter_output><empty_output><jupyter_text>And just like last week we can also plot up the error. In this case its the squared standard deviation.<jupyter_code>ax = data.plot(
    figsize=(20, 10),
    legend=True,
    color="black",
    alpha=0.7,
    markersize=90,
    zorder=2,
)
add_basemap(ax, zoom=9)  # add our basemap to the plot
im = ax.pcolormesh(
    X0,
    X1,
    variance.T,
    cmap="RdBu_r",
    alpha=0.5,
    zorder=1,
    edgecolor=(1.0, 1.0, 1.0, 1.0),
    linewidth=0.01,
)
plt.colorbar(im, ax=ax)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11063423938960419, 'fraction_numerical': 0.02503576537911302, 'mean_word_length': 4.691994572591588, 'pattern_counts': {'":': 0, '<': 39, '<?xml version=': 0, '>': 39, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '231486', 'license': 'Unlicense', 'n_tokens_mistral': 1595, 'n_tokens_neox': 1507, 'n_words': 432, 'path': 'Interpolations/PyKrige Interpolation.ipynb', 'repo_name': 'jessepisel/5minutesofpython', 'size': '4194'}