结合ArcGIS与Arcpy的批量克里金插值与反距离权重插值

尝试了通过Arcpy将shp文件中多字段批量克里金插值的办法，也尝试了直接用Arcpy调用excel再插值的方法，只不过大失败 QAQ 有时间再尝试吧！

准备shp

首先通过ArcGIS界面将excel表（xls格式）导入。【文件——添加数据——添加XY数据——选择sheet表——选择XY字段】，这里我的表格中有两套坐标，分别是用经纬度表示的lat和lon，以及在上一篇博客提到的UTM48N投影坐标系下的投影坐标。如果选择地理坐标系，那么在之后可能需要与上一篇博客类似的操作，将坐标系转为投影坐标。此处我选择投影坐标系下的坐标，与ArcGIS中图层的坐标系一致，然后通过【数据管理工具——要素——要素转点】或者右键点击该事件选择【导出】，将其存为拥有坐标+多个字段的shp文件。

这里需要注意的是，坐标最好与待插值的字段邻近，否则在后续插值时可能照顺序处理不需要插值的解释性字段，出现错误。

生成一个base文件（可选）

由于插值涉及到空间分辨率信息，可以在插值时选择空间分辨率cellsize为固定值，也可以选择基准文件。详细内容见arcpy官方文档： 

克里金法 (3D Analyst)—ArcMap | 文档

在这里为了方便我先用其中一个字段进行了克里金插值（【空间分析工具——插值分析——克里金法】），选择普通克里金法，设置了半变异函数形式（球面函数、高级参数中步长为想要的空间分辨率1000m），网格大小设置为1000m。

同时也尝试用反距离权重插值。【空间分析工具——插值分析——反距离权重法】，参数使用默认参数，网格大小同样也1000m。

两种插值中设置环境，范围均为之前根据气象站位置框定的boundary的utm投影范围。

生成了两个base文件，命名为base.tif（克里金法）和idw.tif（反距离权重法）。

反距离权重法官方文档：

反距离权重法 (地统计分析)—ArcMap | 文档

后续尝试径向基函数插值法，参考官方文档：

径向基函数(RBF)插值法 (地统计分析)—ArcMap | 文档

通过Arcpy批量插值

# -*- encoding: utf-8 -*-
import pandas as pd
import arcpy
from arcpy.sa import *

# 设置环境变量，以便输出栅格数据与原始Excel表格在同一路径下.
arcpy.env.workspace = r"E:\IM_PRO\Meteorological data(1959-2020)\kriging\rainfall\tif\/"
#设置处理范围，这里设置为与生成的base文件一致
arcpy.env.extent = r'E:\IM_PRO\Meteorological data(1959-2020)\kriging\rainfall\tif/base.tif'
#允许覆盖
arcpy.env.overwriteOutput = True

#设置输入shp
inFeatures = r'E:\IM_PRO\Meteorological data(1959-2020)\kriging\rainfall\station_shp/rain47_removezero_utm.shp'

#读取字段
fields = arcpy.ListFields(inFeatures)
field_names = [field.name for field in fields][len(arcpy.Describe(inFeatures).shapeFieldName.split("\\\\")):]
# 设置列表的范围（只包括坐标和待插值字段，去除所有解释性字段）
field_names = field_names[14:-1] 
field_names = [field.encode('ascii') for field in field_names]

#设置输出路径
output_dir = r"E:\IM_PRO\Meteorological data(1959-2020)\kriging\rainfall\tif\idw\\"
def krig_multifields(input,field,output_dir):
    out_raster = output_dir + field +'.tif'
    #设置空间分辨率，此处设置与base文件相同。也可以直接设置像元大小
    cellSize = r'E:\IM_PRO\Meteorological data(1959-2020)\kriging\rainfall\tif/base.tif'
    #设置半变异函数形式，这里都选用默认参数
    kModel = 'Spherical'
    kRadius = "Variable 12"
    #检查是否可以调用扩展模块
    arcpy.CheckOutExtension("3D")
    arcpy.CheckOutExtension("Spatial")
    # 执行克里金插值
    print r'执行克里金插值，插值字段为：%s，输出像元大小为： %s' % (field, cellSize)
    arcpy.Kriging_3d(input, field, out_raster, kModel,
                     cellSize, kRadius)
    print 'save: %s' % out_raster

def idw_multifields(input,field,output_dir):
    out_raster = output_dir + field + '.tif'
    cellSize = r'E:\IM_PRO\Meteorological data(1959-2020)\kriging\rainfall\tif/idw.tif'
    #设置插值所需参数
    power = 2
    searchRadius = "Variable 12"
    arcpy.CheckOutExtension("3D")
    #执行反距离权重插值
    print r'执行反距离权重插值，插值字段为：%s，输出像元大小为： %s' % (field, cellSize)
    arcpy.Idw_3d(input, field, out_raster, cellSize,
                     power, searchRadius)
    print 'save: %s' % out_raster

def rbf_multifields(input,field,output_dir):
    out_raster = output_dir + field + '.tif'
    cellSize = 1000 #直接设置为需要的空间分辨率
    outLayer = "outRBF"
    # 设置插值所需参数
    rbf = "COMPLETELY_REGULARIZED_SPLINE"  #根据官方文档，选择不同的插值函数
    smallscaleParam = "" 
    #使用ArcGIS上的默认值
    majSemiaxis = 235000
    minSemiaxis = 235000 
    angle = 0
    maxNeighbors = 15
    minNeighbors = 10
    sectorType = "ONE_SECTOR"
    searchNeighbourhood = arcpy.SearchNeighborhoodStandard(majSemiaxis,minSemiaxis,
                                                           angle, maxNeighbors,
                                                           minNeighbors, sectorType)

    # Check out the ArcGIS Geostatistical Analyst extension license
    arcpy.CheckOutExtension("GeoStats")
    # 执行插值
    print r'执行径向基函数插值，插值字段为：%s，输出像元大小为： %s' % (field, cellSize)
    arcpy.RadialBasisFunctions_ga(input, field, outLayer, out_raster,
                                  cellSize, searchNeighbourhood, rbf, smallscaleParam)
    print 'save: %s' % out_raster


for field in field_names:
    idw_multifields(inFeatures,field,output_dir)
    #krig_multifields(inFeatures,field)
    #rbf_multifields(inFeatures,field, output_dir)

print('ok')

这里由于shp文件、base文件均为utm投影坐标系啊，没有设置输出坐标系。查阅官方文档，设置输出坐标系的代码为：

arcpy.env.outputCoordinateSystem = arcpy.SpatialReference("WGS 1984 UTM Zone 48N")
arcpy.env.outputCoordinateSystem = arcpy.SpatialReference(4326)

或者是坐标系全称，或者使用EPSG代码都可以。

通过pykrige包插值

import pandas as pd
from pykrige.ok import OrdinaryKriging

def kriging_interpolate(file,time): 
    df = pd.read_excel(file, sheet_name='Sheet')
    Lon = df['Lon']
    Lat = df['Lat']

    data1 = df[time]
    Krig = OrdinaryKriging(Lon,Lat,data1,variogram_model='gaussian')

    lon_range = np.arange(106.3927, 112.959076, 0.0083333)
    lat_range = np.arange(37.633824, 44.440166, 0.0083333)

    data,ssl = Krig.execute('grid', lon_range, lat_range)

    plt.contourf(lon_range, lat_range, data)
    plt.show()

这里只包括对一个字段插值的代码，因为我的格式是坐标+时间（字符串格式），所以使用pykrige包时循环表格中的时间列即可实现批量插值。但这个方法得到的结果与通过arcgis得到的结果相差较大，所以没有再深入探索这个包。后面有时间再说吧！