GDAL（一）：插值和裁剪（python）
在GIS中，经常遇到原始数据是点，但是展示的时候点展示并不好，能区域内连续展示最好了。
这个时候就需要用到插值，把点转成连续的面展示。
大部分的GIS软件中都有插值的工具可以直接使用，不过当遇到批量插值的时候，工具用起来就比较费时间了。
所以就想到写代码，进行批量操作，这样可以运行代码，就不用管了。
既然用代码批量操作，GDAL就是最佳选择。本想是用 .net 的，因为是引用C++的dll，用起来很不方便，而且文档、资料比较少，就选择了 python来写。
一、网格化、算法简介
网上找对应的实现代码是比较好找，但是对这个实现的基本介绍、参数怎么设置等基本没有。
如果只是使用是没问题，但是想用好还是有点问题的。
这里说下自己历程：
-
查找代码实现，这个好找
-
理解参数，这个相对较少
-
算法配置，网上基本没有
最后发现是在官网找到的介绍，比较简单，配置使用是没问题的。
重点：对于成熟的开源项目，看官网！看官网！看官网！这样可以避免很多坑。
这里对于具体的算法也就不拿过来了，大家直接去官网看：
有一个对插值算法解释挺好的，在这里贴出来：（原文地址）
反距离权重插值
'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)
