GDAL与ArcGIS空间插值分析及代码实现（上）

朱雀 2021-12-01 13:50 493阅读 0赞

空间插值，这是通过离散点表征区域特征值的一种方式，空间插值具体的概念及原理请自己百度。实验中对某一时刻的数据进行插值，我们可以使用ArcMap中的工具进行手动处理。但是项目中使用的空间插值，如对每分钟/小时/天的AQI数据进行插值，就不可能手动进行处理了，这时需要考虑开发插值服务解决问题。

空间插值属于GIS的范围，GIS中常见的解决方案包括开源和商业，其中GDAL和ArcGIS分别是各自的优秀代表，下面就对GDAL和ArcGIS空间插值进行说明。

**1、GDAL空间插值(以python为例[https://gdal.org/python/][https_gdal.org_python]）**

（1）GDAL中进行空间插值的函数是 **Grid**，下面是官网对其介绍：

Grid(destName, srcDS, **kwargs)
    
    Create raster from the scattered data.
    Arguments are :
      destName --- Output dataset name
      srcDS --- a Dataset object or a filename
    Keyword arguments are :
      options --- return of gdal.GridOptions(), string or array of strings
      other keywords arguments of gdal.GridOptions()
    If options is provided as a gdal.GridOptions() object, other keywords are ignored.

它的参数包括：

destName：输出数据集（一般为tif文件）的路径（包括文件名）；

srcDS：输入数据集对象或文件名称，可以使用OpenEx打开geojson格式的数据，或使用Open打开shapefile或直接输入shapefile路径；

\*\*kwargs：插值参数，可以使用 key = value 关键字参数的形式进行设置，也可以使用 option = GridOptions()的形式进行一次性设置，

（2）下面对**GridOptions**的内容进行说明：

GridOptions(options=[], format=None, outputType=GDT_Unknown, width=0, height=0, creationOptions=None, outputBounds=None, outputSRS=None, noData=None, algorithm=None, layers=None, SQLStatement=None, where=None, spatFilter=None, zfield=None, z_increase=None, z_multiply=None, callback=None, callback_data=None)
     
    Create a GridOptions() object that can be passed to gdal.Grid()
    Keyword arguments are :
      options --- can be be an array of strings, a string or let empty and filled from other keywords.
      format --- output format ("GTiff", etc...)
      outputType --- output type (gdal.GDT_Byte, etc...)
      width --- width of the output raster in pixel
      height --- height of the output raster in pixel
      creationOptions --- list of creation options
      outputBounds --- assigned output bounds: [ulx, uly, lrx, lry]
      outputSRS --- assigned output SRS
      noData --- nodata value
      algorithm --- e.g "invdist:power=2.0:smoothing=0.0:radius1=0.0:radius2=0.0:angle=0.0:max_points=0:min_points=0:nodata=0.0"
      layers --- list of layers to convert
      SQLStatement --- SQL statement to apply to the source dataset
      where --- WHERE clause to apply to source layer(s)
      spatFilter --- spatial filter as (minX, minY, maxX, maxY) bounding box
      zfield --- Identifies an attribute field on the features to be used to get a Z value from. This value overrides Z value read from feature geometry record.
      z_increase --- Addition to the attribute field on the features to be used to get a Z value from. The addition should be the same unit as Z value. The result value will be Z value + Z increase value. The default value is 0.
      z_multiply - Multiplication ratio for Z field. This can be used for shift from e.g. foot to meters or from  elevation to deep. The result value will be (Z value + Z increase value) * Z multiply value.  The default value is 1.
      callback --- callback method
      callback_data --- user data for callback

**GridOptions**中的参数，需要以 key=value 关键字参数的形式进行设置，基本的参数包括：

format(输出文件格式)、outputType（输出数据类型，如字节GDT\_Byte、32位浮点型gdal.GDT\_Float32等）、width和height（输出栅格宽度和高度）、outputBounds（插值范围）、outputSrs（输出坐标系）、noData（没有值的标记）、algorithm（插值算法，下面会做详细介绍）、SQLStatement（通过sql语句对输入数据进行过滤）、where（通过条件对参与插值的数据进行过滤）、zfield（插值字段。如aqi、co等）、z\_increase和z\_multiply（）对插值字段进行处理的算子），callback（插值完成后的回调函数）、callback\_data（回调函数的用户数据）。

（3）[插值算法][Link 1]

最常用的是 反距离权重算法，下面对该算法的参数进行说明：

'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的填充值。