59 lines
1.9 KiB
Python
59 lines
1.9 KiB
Python
import numpy as np
|
|
import pykrige.kriging_tools as kt
|
|
from pykrige import OrdinaryKriging
|
|
import rasterio as rio
|
|
from ..util.utils import *
|
|
|
|
|
|
class KrigingSpiData(object):
|
|
|
|
def __init__(self, data, save_path, XY):
|
|
self.save_path = save_path
|
|
self.data = data
|
|
self.XY = XY
|
|
self.get_conf()
|
|
|
|
def get_minmax(self):
|
|
delta = 0.1
|
|
miny = self.XY['Y'].min() - delta
|
|
minx = self.XY['X'].min() - delta
|
|
maxy = self.XY['Y'].max() + delta
|
|
maxx = self.XY['X'].max() + delta
|
|
# print(miny, minx, maxy, maxx)
|
|
return (minx, miny, maxx, maxy)
|
|
|
|
def get_conf(self):
|
|
(minx, miny, maxx, maxy) = self.get_minmax()
|
|
res = 0.003
|
|
self.gridx = np.arange(minx, maxx, res)
|
|
self.gridy = np.arange(miny, maxy, res)
|
|
self.height = len(self.gridy)
|
|
self.width = len(self.gridx)
|
|
from rasterio.transform import from_origin
|
|
self.transform = from_origin(minx - res / 2, maxy + res / 2, res, res)
|
|
|
|
def write_tif(self, kriging_data):
|
|
out_meta = {"driver": "GTiff",
|
|
"height": self.height,
|
|
"width": self.width,
|
|
"transform": self.transform}
|
|
from rasterio.crs import CRS
|
|
with rio.open(self.save_path, "w", count=1, dtype='float64', crs=CRS(init='epsg:4326'), **out_meta) as dest:
|
|
dest.write_band(1, kriging_data)
|
|
|
|
def kriging(self):
|
|
OK = OrdinaryKriging(np.asarray(self.XY['X']), np.asarray(self.XY['Y']), self.data, variogram_model='spherical', verbose=False,
|
|
enable_plotting=False)
|
|
z, ss = OK.execute('grid', self.gridx, self.gridy)
|
|
self.write_tif(z)
|
|
cut_raster_by_geo(self.save_path, self.save_path, get_region_geo())
|
|
# copy to geoserver
|
|
from shutil import copyfile
|
|
copyfile(self.save_path, os.path.join(GEOSERVER_SPI_PATH, 'temp.tif'))
|
|
# reproject_by_gdal(self.save_path, os.path.join(TEMP_PATH, 'temp_3857.tif'))
|
|
|
|
|
|
if __name__ == '__main__':
|
|
kd = KrigingSpiData("D:/work/3.scala-gis/geotrellis-through/geotrellis/bj_197901.txt", "D://")
|
|
kd.kriging()
|