GDAL（一）：插值和裁剪（python）

在GIS中，经常遇到原始数据是点，但是展示的时候点展示并不好，能区域内连续展示最好了。

这个时候就需要用到插值，把点转成连续的面展示。

大部分的GIS软件中都有插值的工具可以直接使用，不过当遇到批量插值的时候，工具用起来就比较费时间了。

所以就想到写代码，进行批量操作，这样可以运行代码，就不用管了。

既然用代码批量操作，GDAL就是最佳选择。本想是用 .net 的，因为是引用C++的dll，用起来很不方便，而且文档、资料比较少，就选择了 python来写。

一、网格化、算法简介

网上找对应的实现代码是比较好找，但是对这个实现的基本介绍、参数怎么设置等基本没有。

如果只是使用是没问题，但是想用好还是有点问题的。

这里说下自己历程：

• 查找代码实现，这个好找

• 理解参数，这个相对较少

• 算法配置，网上基本没有

最后发现是在官网找到的介绍，比较简单，配置使用是没问题的。

重点：对于成熟的开源项目，看官网！看官网！看官网！这样可以避免很多坑。

这里对于具体的算法也就不拿过来了，大家直接去官网看：

Grid 简介 （对应中文版）

Grid 参数解释、配置（对应中文版）

 

有一个对插值算法解释挺好的，在这里贴出来：（原文地址）

反距离权重插值

'invdist:power=3.6:smoothing=0.2:radius1=0.0:radius2=0.0:angle=0.0:max_points=0:min_points=0:nodata=0.0'

参数解释：

• invdist：表示算法名称 - 反距离权重，即距离越近权重越大，对待计算的网格影像越大；

• power：距离对权重的影响程度）（数字表示指数），默认值为2, 值越大表示距离越近的点对插值的影响程指数倍增；对于该参数的设置，若点较密集且均匀分布，则值设置在2以下，若点相对较少且分布不均，则为了保障插值效果，可将参数设置在3以上；

• smoothing：平滑系数，默认为0，数值越大表示越平滑，同时精度也会越低；

• radius1：表示椭圆x轴方向上的半径；

• radius2：表示椭圆y轴方向上的半径；

• angle：表示椭圆旋转的弧度；

• max_points：表示参与插值的最大点数，默认值0表示搜索到多少就是多少；

• min_points：表示参与插值的最小点数，默认值0表示搜索到多少就是多少；

• nodata：对nodata的填充值

二、代码实现

代码实现起来还是挺简单的：

# field 这里是要插值的字段
def insertRaster(field):
    startTime = ps.to_datetime(datetime.datetime.now())
    print('开始插值：%s' % startTime)
    point_file = 'xxx.shp'
    output_file = 'xxx' + field + '.tif'

    opts = gdal.GridOptions(format="GTiff", outputType=gdal.GDT_Float32, width=2472, height=1000,
                            algorithm="invdist:power=3.5:smothing=0.0:radius=1.0:max_points=12:min_points=0:nodata=0.0", zfield=field)

    gdal.Grid(destName=output_file, srcDS=point_file, options=opts)
    endTime = ps.to_datetime(datetime.datetime.now())
    useTime = endTime - startTime
    print('插值完成！用时：%s s' % useTime.total_seconds())

这里面的 GridOptions 里面的参数都可以看文档，算法是上面给的链接，可以根据里面进行配置。

三、裁剪

为什么要在这里带上裁剪。

插值的输出结果，都是一个矩形。大部分场景下，这和我们需要使用的范围不一致。而且源数据的范围也不一定是矩形，在源数据范围以外，插值的结果都是无效的，所以要裁减掉。

裁剪对应的简介、参数配置也都可以在官网找到。对应直接上代码：

def cutRaster():
    input_shp = 'xxx.shp'
    input_raster = 'xxx-r '.tif'
    result_raster = 'xxx-r'_c.tiff'

    if os.path.exists(input_raster):
        startTime = ps.to_datetime(datetime.datetime.now())
        print('开始裁剪：%s' % startTime)
        ds = gdal.Warp(result_raster, input_raster, format='GTiff',
                       cutlineDSName=input_shp, dstNodata=0)
        endTime = ps.to_datetime(datetime.datetime.now())
        useTime = endTime - startTime
        print('裁剪完成，用时：%s s' % useTime.total_seconds())
    else:
        print('文件不存在：%s' % input_raster)