Using PyKrige for ordinary kriging on an XYZ file

Question

I've previously used ArcGIS for Ordinary Kriging, but I now need a Python alternative.

I found the PyKrige module which i decided to give a try, but I am bit confused by it.

So far, I have read my XYZ file into a numpy array like so

import matplotlib.pyplot as plt
import numpy as np
import pykrige.kriging_tools as kt
from pykrige.ok import OrdinaryKriging
from numpy import genfromtxt

np.set_printoptions(suppress=True)

my_data = genfromtxt('C:/myfile.xyz')
data = np.array(my_data)

print(data)

which prints my data like so

[[626769.27 233220.76     14.66]
 [626768.91 233221.02     14.65]
 [626768.58 233221.26     14.64]
 ...

However, the next part of the PyKrige example found at https://geostat-framework.readthedocs.io/projects/pykrige/en/stable/examples/00_ordinary.html#sphx-glr-examples-00-ordinary-py uses the following code after reading the data into an array

gridx = np.arange(0.0, 5.5, 0.5)
gridy = np.arange(0.0, 5.5, 0.5)

OK = OrdinaryKriging(
data[:, 0],
data[:, 1],
data[:, 2],
variogram_model="linear",
verbose=False,
enable_plotting=False,
)

z, ss = OK.execute("grid", gridx, gridy)

kt.write_asc_grid(gridx, gridy, z, filename="pykrigoutput.asc")
plt.imshow(z)
plt.show()

This is what's confusing me. I don't know how the gridx and gridy are determined? Similarly, is z, ss the variable name or what is this?

My output will be a GeoTiff (hopefully that's supported) and I hope to set an output cell size of 0.5

But can anyone advise me how to successfully use this module as the documentation gives no explanation of the example.

Answer 1

1. gridx and gridy are the axes of the output grid (the grid where you want to generate the kriged values on)
2. z are the output values at the specified grid points
3. ss stands for "sigma-square" and is the kriging variance at the specified grid points

In order to have the output as a geotiff file, you would need to write one and specify your desired grid beforehand:

gridx = np.arange(0.0, 5.5, 0.5)
gridy = np.arange(0.0, 5.5, 0.5)

Already defines a grid with a step-size of 0.5. You still need to define ranges for x and y (in this case both range from 0 to 5 [since 5.5 is excluded]).

For writing a GeoTIFF file from a numpy array, have a look here: https://here.isnew.info/how-to-save-a-numpy-array-as-a-geotiff-file-using-gdal.html

Hope that helps!