python使用开源库进行克里金插值

克里金插值法

import rasterio
import pandas as pd
from rasterio.mask import mask
from rasterio.transform import Affine
from pykrige.ok import OrdinaryKriging
import numpy as np

# ======================================= 输入参数 ==================================================== #

rasterOutPath = "output_tif01.tif"
rasterCropPath = "output_tif02.tif"
raindata = pd.read_csv("RAIN_input.csv")  # 降雨量数据
maskPath = 'mask.shp'

res = 100 #你自己设置的图像分辨率
# ======================================= 输入参数 END ================================================ #

longitude = raindata['coox'].tolist()
latitude = raindata['cooy'].tolist()
elevation = raindata['elevation'].tolist()
# -------------------------------------------------------------------------------------------------------------------- #
# --------------------------------------------       krige插值      --------------------------------------------------- #
# -------------------------------------------------------------------------------------------------------------------- #

longitude = list(map(float,longitude))
latitude = list(map(float,latitude))
# 生成栅格的x、y
right_add = 300
ud_add = 100
# 由于雨量站的坐标不是在广州地图的边缘因此需要更大的面积  这里在栅格上下增加ud_add 在栅格右侧增加right_add
grid_lon = np.linspace(min(longitude)-10, max(longitude)+(10+right_add*100),int(abs(max(longitude) - min(longitude) )/res)+int(right_add*100/res))
grid_lat = np.linspace(min(latitude)-(10+ud_add*100), max(latitude)+(10+ud_add*100),int(abs(max(latitude) - min(latitude))/res)+int(2*ud_add*100/res))
# (min(longitude) - max(longitude))/100

OK = OrdinaryKriging(longitude, latitude, elevation, variogram_model='exponential',nlags=20)
z1, ss1 = OK.execute('grid', grid_lon, grid_lat, mask=maskPath)

transform = Affine.translation(grid_lon[0], grid_lat[0]) * Affine.scale(abs(res), abs(res)) #tif的经纬度生成
# lon 经度，lat 维度 剩下的套用我这句话就可以了

with rasterio.open(
        rasterOutPath,
        'w',
        driver='GTiff',  # 设置数据类型
        height=len(grid_lat),  # 维度的格子数量
        width=len(grid_lon),  # 经度的格子数量
        count=1,  # 波段数量
        dtype=z1.dtype,  # 数据类型
        crs=None,  # 投影
        transform=transform, ) as f:
    meta = f.meta
    f.write(z1.data, 1)  # 想文件中写入数据
# -------------------------------------------------------------------------------------------------------------------- #
# --------------------------------------------       krige插值  end --------------------------------------------------- #
# -------------------------------------------------------------------------------------------------------------------- #


# ******************************************************************************************************************* #
# ***********************************************     加掩膜裁剪      *********************************************** #
# ******************************************************************************************************************* #

from osgeo import gdal

input_raster = rasterOutPath

output_raster = rasterCropPath


input_raster = gdal.Open(input_raster)
ds = gdal.Warp(output_raster, #输出栅格路径
               input_raster, # 输入栅格路径
               format='GTiff', # 影像保存格式
               cutlineDSName=maskPath, # 输入矢量路径
               cropToCutline=False, #(为True时，结果会与输入矢量大小一致。为False时，结果会与带裁剪的输入栅格大小一致)
               dstSRS='EPSG:4326', # 参考系：WGS84
               outputType=gdal.GDT_Float64, #数据类型
               dstNodata=0 )
ds=None

# ******************************************************************************************************************* #
# ***********************************************     加掩膜裁剪  end *********************************************** #
# ******************************************************************************************************************* #

参考自
《地理信息系统算法基础》科学出版社