GIS之arcgis系列09：arcpy实现克里金差值

矢量点数据经过克里金差值后可以转换成栅格数据，那么就需要了解一下什么是克里金差值。

什么是克里金法?

　　IDW(反距离加权法)和样条函数法插值工具被称为确定性插值方法，因为这些方法直接基于周围的测量值或确定生成表面的平滑度的指定数学公式。第二类插值方法由地统计方法(如克里金法)组成，该方法基于包含自相关(即，测量点之间的统计关系)的统计模型。因此，地统计方法不仅具有产生预测表面的功能，而且能够对预测的确定性或准确性提供某种度量。

　　克里金法假定采样点之间的距离或方向可以反映可用于说明表面变化的空间相关性。克里金法工具可将数学函数与指定数量的点或指定半径内的所有点进行拟合以确定每个位置的输出值。克里金法是一个多步过程;它包括数据的探索性统计分析、变异函数建模和创建表面，还包括研究方差表面。当您了解数据中存在空间相关距离或方向偏差后，便会认为克里金法是最适合的方法。该方法通常用在土壤科学和地质中。

完整版代码如下：

代码均已经过测试，可直接微调后再工具箱内启动。

import arcpy
import os
import numpy as np

input_folder = arcpy.GetParameterAsText(0)
field_name = arcpy.GetParameterAsText(1)
pixel_size = arcpy.GetParameterAsText(2)
output_folder = arcpy.GetParameterAsText(3)
name = field_name

if not os.path.exists(output_folder):
    os.makedirs(output_folder)

def check_all_zero(input_shp, field_name):
    with arcpy.da.SearchCursor(input_shp, [field_name]) as cursor:
        for row in cursor:
            if row[0] != 0:
                return False
    return True

for root, dirs, files in os.walk(input_folder):
    for file in files:
        if file.endswith(".shp"):
            input_shp = os.path.join(root, file)

            if arcpy.Exists(input_shp):
                projected_shp = os.path.join(output_folder, "projected.shp")
                arcpy.Delete_management(projected_shp)
                arcpy.Project_management(input_shp, projected_shp, arcpy.SpatialReference(32651))

                relative_path = os.path.relpath(root, input_folder)
                output_subfolder = os.path.join(output_folder, relative_path)
                if not os.path.exists(output_subfolder):
                    os.makedirs(output_subfolder)

                output_raster = os.path.join(output_subfolder, os.path.splitext(file)[0] + "_{}.tif".format(name))

                if check_all_zero(input_shp, field_name):
                    desc = arcpy.Describe(projected_shp)
                    extent = desc.extent
                    sr = arcpy.SpatialReference(32651)
                    array = np.zeros((int((extent.YMax - extent.YMin) / int(pixel_size)),
                                      int((extent.XMax - extent.XMin) / int(pixel_size))), dtype=np.float32)
                    raster = arcpy.NumPyArrayToRaster(array, extent.lowerLeft, int(pixel_size), int(pixel_size), -9999)
                    raster.save(output_raster)


                    arcpy.DefineProjection_management(output_raster, sr)
                else:
                    arcpy.CheckOutExtension("Spatial")
                    Output_variance_of_prediction_raster = ""
                    arcpy.gp.Kriging_sa(projected_shp,field_name, output_raster, "Spherical 2548.447994", pixel_size,
                                        "VARIABLE 12",Output_variance_of_prediction_raster)
                    arcpy.CheckInExtension("Spatial")

                arcpy.Delete_management(projected_shp)