I want to use a Python script for interpolation in 1st image to have the result of 2nd one, this result is from simple kriging saga. I tried the code below but I have two main problems.
- Need a lot of time, over 3-4 hours to have
- The result is this in 3d image and it's really different.
from osgeo import gdal,ogr
import numpy as np
import gstools as gs
import time
#read a 2-D array from image
filefolder="C:/Users/fermi/Desktop/SUEWS/SUEWSapp/kriging/"
rasterinput=filefolder+"T2UFfill_T2_01_14.tif"
src_ds = gdal.Open( rasterinput )
b3 = np.array(src_ds.GetRasterBand(1).ReadAsArray())
# the coordinates of x and y
x= np.arange(322580., 335780., 100.)
y = np.arange(3906895., 3913695. , 100.)
start = time.time()
print("hello")
# your input data as x, y, data points
data = b3
cond_pos = gs.generate_grid((x,y))
cond_val = [1.0]
mean = np.mean(data)
# simple gaussian model since your plot looks rather smooth
model = gs.Gaussian(dim=2,len_scale=1)
krige = gs.krige.Simple(model, cond_pos, data)
# output on a structured grid (image) with specified axis
xaxis, yaxis = y,x
field = krige.structured([xaxis, yaxis], return_var=False)
