Using PyKrige to interpolate xyz data

[aok289]

I have tried putting my XYZ data into numpy arrays as per the PyKrige documentation

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:/mydata.xyz')
data = np.array(my_data)

However, from this point on I'm not sure how to apply Ordinary Kriging to the data. I'd ideally like the output to be a Tiff image, but don't require the plotting of graph, etc - so do I still need gridx and gridy??

I'm just abit lost beyond this point

[MuellerSeb]

Copying the answer from here: https://gis.stackexchange.com/questions/440450/using-pykrige-for-ordinary-kriging-on-an-xyz-file

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!

[aok289]

Using the above information I tried writing my code

import numpy as np
import pandas as pd
import glob
from pykrige.ok import OrdinaryKriging
from pykrige.kriging_tools import write_asc_grid
import pykrige.kriging_tools as kt
import matplotlib.pyplot as plt
from matplotlib.colors import LinearSegmentedColormap
from matplotlib.patches import Path, PathPatch

datafile='C:/myfile.xyz'
df=pd.read_csv(datafile,delimiter='\t')

print(df)

lons=np.array(df['x'])
lats=np.array(df['y'])
data=np.array(df['z'])

grid_space = 0.5
grid_lon = np.arange(np.amin(lons), np.amax(lons), grid_space) #grid_space is the desired delta/step of the output array
grid_lat = np.arange(np.amin(lats), np.amax(lats), grid_space)

OK = OrdinaryKriging(lons, lats, data, variogram_model='gaussian', verbose=True, enable_plotting=False,nlags=20)
z1, ss1 = OK.execute('grid', grid_lon, grid_lat)

kt.write_asc_grid(grid_lon, grid_lat, z1, filename="C:/users/public/output.asc")

print('job complete')

But now I'm getting memory errors

Traceback (most recent call last):
  File "c:\Users\anonymous\Documents\pykrige_kriginginterpolation.py", line 25, in <module>
    z1, ss1 = OK.execute('grid', grid_lon, grid_lat)
  File "C:\Python10\lib\site-packages\pykrige\ok.py", line 847, in execute
    a = self._get_kriging_matrix(n)
  File "C:\Python10\lib\site-packages\pykrige\ok.py", line 642, in _get_kriging_matrix
    a[:n, :n] = -self.variogram_function(self.variogram_model_parameters, d)
  File "C:\Python10\lib\site-packages\pykrige\variogram_models.py", line 45, in gaussian_variogram_model
    return psill * (1.0 - np.exp(-(d**2.0) / (range_ * 4.0 / 7.0) ** 2.0)) + nugget
numpy.core._exceptions._ArrayMemoryError: Unable to allocate 10.2 GiB for an array with shape (36935, 36935) and data type float64

I'm not even sure I am going about this the right way. I am wanting a python alternative to ArcGIS Ordinary Kriging, but I seem to keep struggling with this module

[LSchueler]

Have you tried reducing the grid_space? The error pretty clearly tells you, that you do not have enough free RAM. 36935² are quite a lot of grid cells. But 10GiB of RAM is not that much nowadays. If you execute the code on a different machine with more RAM, the code should run.