2

I have an xyz file in which I'd like to apply the interpolation method of kriging to via python.

In order to do this, I have come across the module PyKrige and duplicated this example perfectly on my own machine https://colab.research.google.com/drive/1mZDVJ-RioBc4QHHpe2HBNJrx8sJ4-1Lw?pli=1#scrollTo=z0XaQWJq_gJ7

This gave me hope that the module could be easily implemented on my own data - sadly not.

My xyz file looks like so, and there are 36,000 sets of xyz points (973KB)

626814.86   233370.97   14.56
626814.43   233371.28   14.55
626814.01   233371.56   14.54
626813.57   233371.85   14.53
626813.10   233372.15   14.51
626812.66   233372.44   14.49
626812.23   233372.72   14.46
626811.78   233373.01   14.43

I have then altered the script from the example above to suit my data, by reading them into their 3 separate arrays and what not

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

x1 = np.loadtxt("c:/test.xyz", usecols=0, dtype='float')
y1 = np.loadtxt("c:/test.xyz", usecols=1, dtype='float')
z1 = np.loadtxt("c:/test.xyz", usecols=2, dtype='float')

x = list(x1) #may not be a necessary step, but I added it so that after each x coordinate there was a comma like the examples array. Script wasn't working even before doing this
y = list(y1)
z = list(z1)


cax = plt.scatter(x, y, c=z)
cbar = plt.colorbar(cax, fraction=0.03)
plt.title('Test Points')


OK = OrdinaryKriging(
    x, 
    y, 
    z, 
    variogram_model='gaussian',
    verbose=True,
    enable_plotting=False,
    nlags=5,
)

OK.variogram_model_parameters

gridx = np.arange(626000, 629094, 5, dtype='float64')
gridy = np.arange(231500, 234500, 5, dtype='float64')
zstar, ss = OK.execute("grid", gridx, gridy)

print(zstar.shape)
print(ss.shape)

cax = plt.imshow(zstar, extent=(626000, 629094, 231500, 234500), origin='lower')
plt.scatter(x, y, c='z', marker='.')
cbar=plt.colorbar(cax)
plt.title('Test')
savefig( 'C:/users/public/test.png' )
print('image saved')

My issue is that it gets to the execution of the Ordinary Kriging and then complains

numpy.core._exceptions._ArrayMemoryError: Unable to allocate 10.2 GiB for an array with shape (36936, 36936) and data type float64

Understandably, I thought that the issue was the 36,000 points being too many. So, I decided to try a smaller snippet of the data. Now when I say smaller, I mean 20 points (1KB!!!). THE SAME ERROR!!!!

Can someone please tell me what I'm doing wrong? Realistically, I only want the interpolated PNG, not the points and scatterplot

0

1 Answer 1

1

Fix this line gridx = np.arange(62600, 629094, 5, dtype='float64')

You're missing a digit in the start parameter - "62600", should have 6 digits eg "626000". So your current script is trying to interpolate onto a very large array.

9
  • oh yes I see that now. Howvever, even with that fixed I'm still getting the error numpy.core._exceptions._ArrayMemoryError: Unable to allocate 10.2 GiB for an array with shape (36936, 36936) and data type float64
    – hi2938021
    Sep 14, 2022 at 22:32
  • 1
    Works for me with your code and those example points. It only blows out the memory with gridx = np.arange(62600, 629094, 5, dtype='float64')
    – user2856
    Sep 15, 2022 at 0:00
  • Yes I got it working for the 8 points, but the full 36,000-point dataset is proving to be an issue
    – hi2938021
    Sep 15, 2022 at 9:56
  • you can try with dtype="float32" if you are sure you don't need the float64 precision otherwise you might have to do it in batches Sep 15, 2022 at 10:21
  • @LouisCottereau Realistically, I only want the raster image from the background. I don't even need the scattergraph/plots - so I'm wondering would it help if I were to try and remove them
    – hi2938021
    Sep 15, 2022 at 10:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.