GDAL与ArcGIS空间插值分析及代码实现(上)
空间插值,这是通过离散点表征区域特征值的一种方式,空间插值具体的概念及原理请自己百度。实验中对某一时刻的数据进行插值,我们可以使用ArcMap中的工具进行手动处理。但是项目中使用的空间插值,如对每分钟/小时/天的AQI数据进行插值,就不可能手动进行处理了,这时需要考虑开发插值服务解决问题。空间插值属于GIS的范围,GIS中常见的解决方案包括开源和商业,其中GDA...
空间插值,这是通过离散点表征区域特征值的一种方式,空间插值具体的概念及原理请自己百度。实验中对某一时刻的数据进行插值,我们可以使用ArcMap中的工具进行手动处理。但是项目中使用的空间插值,如对每分钟/小时/天的AQI数据进行插值,就不可能手动进行处理了,这时需要考虑开发插值服务解决问题。
空间插值属于GIS的范围,GIS中常见的解决方案包括开源和商业,其中GDAL和ArcGIS分别是各自的优秀代表,下面就对GDAL和ArcGIS空间插值进行说明。
1、GDAL空间插值(以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)插值算法
最常用的是 反距离权重算法,下面对该算法的参数进行说明:
'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的填充值。
2、ArcGIS空间插值(以ArcPy为例说明)
下次再进行补充...
开放原子开发者工作坊旨在鼓励更多人参与开源活动,与志同道合的开发者们相互交流开发经验、分享开发心得、获取前沿技术趋势。工作坊有多种形式的开发者活动,如meetup、训练营等,主打技术交流,干货满满,真诚地邀请各位开发者共同参与!
更多推荐
所有评论(0)