GeoStat-Framework / PyKrige

Kriging Toolkit for Python
https://pykrige.readthedocs.io
BSD 3-Clause "New" or "Revised" License
740 stars 186 forks source link

Kriging for geotif image #210

Closed Empanag closed 1 year ago

Empanag commented 2 years ago

Hello, i m trying to use kriging to tiff image. the image is that below with number 1 that can be used as an array, it is regular grid. the 2nd image is the one that saga simple kriging module generates. I 'd like to produce a similar product using python. Is it possible somehow PyKrige? I tried with the example "ordinary kriging" that there is in the site had but the product was the same with input image.

image

Empanag commented 2 years ago

i tried the code below but the product is the same with the input image


import gstools as gs
import pykrige.kriging_tools as kt
from pykrige.ok import OrdinaryKriging
import matplotlib.pyplot as plt
from osgeo import gdal,ogr

filefolder="C:/Users/fermi/Desktop/SUEWS/SUEWSapp/kriging/"
rasterpath=filefolder+"T2UFfill_T2_01_14.tif"
print(rasterpath)
raster=gdal.Open(rasterpath)
print(raster)
rasterBand=raster.GetRasterBand(1)
rasterArray= rasterBand.ReadAsArray()
print(rasterArray.shape)

lstarray=[]
gridx = np.arange(3225.80, 3357.80, 1.0)
gridy = np.arange(39068.95, 39136.95, 1.0)
for i in range(len(gridx)):
    for j in range(len(gridy)):
        lstarray.append([gridx[i],gridy[j],rasterArray[j,i]])
rasterarray=np.array(lstarray)

cov_model = gs.Gaussian(dim=2, len_scale=0.1,  var=0.01)

OK2 = OrdinaryKriging(rasterarray[:, 0], rasterarray[:, 1], rasterarray[:, 2], variogram_model=cov_model)
z1, ss1 = OK2.execute("grid", gridx, gridy)

plt.imshow(np.flip(z1,0), origin="lower")
plt.show()
MuellerSeb commented 1 year ago

See: https://github.com/GeoStat-Framework/GSTools/discussions/199