Python gdal warp函数 矢量裁剪 遥感影像投影变换

gdal.Warp可以对栅格数据依据输入的矢量进行裁剪，还可以进行栅格重采样，重投影等操作，关键在于options的参数设置。Options的参数设置可以参见网址：https://gdal.org/python/，其中在左侧可以选择Warp和WarpOptions两个模块查看参数。

1. Warp参数解释

warp是扭曲或者变形的意思，该函数是对一个数据或者多个数据集进行变形，用法如下：

Warp(destNameOrDestDS, srcDSOrSrcDSTab, **kwargs)

参数包括：

• destNameOrDestDS --- 输出数据集名称或对象
• srcDSOrSrcDSTab --- 数据集对象或文件名的数组，或数据集对象或文件名的数组

关键字参数是：

• options --- 返回gdal.WarpOptions()，字符串或字符串数组
• gdal.WarpOptions()的其他关键词参数

如果选项是作为gdal.WarpOptions()对象提供的，则会忽略其他关键字。

2. WarpOptions参数解释（可以略过，直接看代码实例）

warp工具实现的关键在于参数options的设置， WarpOptions的功能即是创建一个可以传递给gdal.Warp()的WarpOptions()对象。下面分析WarpOptions参数，其参数较多，如下所示：

WarpOptions(options=None, format=None, outputBounds=None, outputBoundsSRS=None, xRes=None, yRes=None, targetAlignedPixels=False, width=0, height=0, srcSRS=None, dstSRS=None, coordinateOperation=None, srcAlpha=False, dstAlpha=False, warpOptions=None, errorThreshold=None, warpMemoryLimit=None, creationOptions=None, outputType=gdalconst.GDT_Unknown, workingType=gdalconst.GDT_Unknown, resampleAlg=None, srcNodata=None, dstNodata=None, multithread=False, tps=False, rpc=False, geoloc=False, polynomialOrder=None, transformerOptions=None, cutlineDSName=None, cutlineLayer=None, cutlineWhere=None, cutlineSQL=None, cutlineBlend=None, cropToCutline=False, copyMetadata=True, metadataConflictValue=None, setColorInterpretation=False, overviewLevel='AUTO', callback=None, callback_data=None)

挨个对参数进行下解释。

• options --- 可以是一个字符串数组，一个字符串或让空和填充从其他关键字。
• format --- 输出格式(“GTiff”等)
• outputBounds --- 在输出的空间参照系统中的输出边界(minX, minY, maxX, maxY)
• outputBoundsSRS --- 表示输出边界的SRS以防它们没有在dstSRS中表示，其中SRS是spatial reference system的缩写，即空间参照系统
• xRes, yRes --- 目标SRS的空间分辨率
• targetAlignedPixels --- 是否强制输出边界为输出分辨率的倍数
• width --- 输出光栅的像素宽度
• height --- 输出光栅的像素高度
• srcSRS --- 输入的 SRS
• dstSRS --- 输出的 SRS
• coordinateOperation -- 坐标操作是 PROJ string 或者是 WKT string
• srcAlpha --- 是否强制将输入数据集的最后一个频带视为alpha频带
• dstAlpha --- 是否强制创建输出波段
• outputType --- 输出类型 (gdalconst.GDT_Byte, etc...)
• workingType --- 操作类型 (gdalconst.GDT_Byte, etc...)
• warpOptions --- warping options 列表
• errorThreshold --- 近似转换器的误差阈值 (in pixels)
• warpMemoryLimit --- 工作缓冲区的大小，以MB为单位
• resampleAlg --- 重采样模式
• creationOptions --- creation options 列表
• srcNodata --- source nodata value(s)
• dstNodata --- output nodata value(s)
• multithread --- 是否多线程计算和I/O操作
• tps --- 是否使用薄板样条GCP转换器
• rpc --- 是否使用RPC转换器
• geoloc --- 是否使用GeoLocation阵列转换器
• polynomialOrder --- 多项式GCP插值阶
• transformerOptions --- transformer options 列表
• cutlineDSName --- 切割线的数据集的名称
• cutlineLayer --- 切割线的图层名称
• cutlineWhere --- 切割线的 WHERE 从句
• cutlineSQL --- 切割线的 SQL 语句
• cutlineBlend --- 以像素为单位的切割线混合距离
• cropToCutline --- 是否使用切割线范围作为输出边界
• copyMetadata --- 是否拷贝源元数据
• metadataConflictValue --- 元数据数据冲突值
• setColorInterpretation --- 是否强制输入波段的颜色解释到输出波段
• overviewLevel --- 指定必须使用的源文件的概述级别
• callback --- 回调函数
• callback_data --- 用于回调的用户数据

3. 栅格数据矢量裁剪（批量）

在矢量裁剪中，主要用到的是cutline相关的参数，包括切割后的数据集（cutlineDSName），切割图层（cutlineLayer ），是否使用切割线范围作为输出边界（cropToCutline）等，下面展示实例代码。

from osgeo import gdal

ds = gdal.Warp(tif_out, tif_input, format='GTiff', cutlineDSName=shp_boundary, cropToCutline=True, dstNodata=0)

4. 遥感影像重投影

对坐标系进行重投影也很简单，主要是坐标系的设置，如果是非标准的自定义的坐标系，可以参见方法ESRI中WKT格式的空间投影坐标转换为GDAL中的Proj4格式转化为proj4的格式。具体实例代码如下：

from osgeo import gdal

ds = gdal.Warp(tif_out, tif_input, srcSRS='EPSG:4326', dstSRS=srs_proj4) # srs_proj4是自定义的坐标系

5. 遥感影像重采样（修改分辨率）

栅格数据或者说遥感影像的重采样也比较简单，重采样方法的确定在gdal中可以通过ReprojectImage()函数查看，参见https://gdal.org/python/。

采用Warp函数的话，主要修改的resampleAlg参数，默认的最近邻重采样方法，查看源码发现其他可选参数如下。具体重采样方法的差异不再赘述。

      if resampleAlg is not None: 

              mapMethodToString = { 
                  gdalconst.GRA_NearestNeighbour: 'near', 
                   gdalconst.GRA_Bilinear: 'bilinear', 
                  gdalconst.GRA_Cubic: 'cubic', 
                  gdalconst.GRA_CubicSpline: 'cubicspline', 
                  gdalconst.GRA_Lanczos: 'lanczos', 
                  gdalconst.GRA_Average: 'average', 
                  gdalconst.GRA_RMS: 'rms', 
                  gdalconst.GRA_Mode: 'mode', 
                  gdalconst.GRA_Max: 'max', 
                  gdalconst.GRA_Min: 'min', 
                  gdalconst.GRA_Med: 'med', 
                  gdalconst.GRA_Q1: 'q1', 
                  gdalconst.GRA_Q3: 'q3', 
                  gdalconst.GRA_Sum: 'sum', 
              }

代码如下：

from osgeo import gdal

ds = gdal.Warp(tif_out, tif_input, resampleAlg=gdalconst.GRA_NearestNeighbour, xRes=500, yRes=500)

至此，栅格数据（如遥感影像）的裁剪、重投影、重采样都可以用Warp函数实现，而且可以在同一个函数下实现所有的操作，只需要把各个options放在一起即可。

初次使用gdal完毕，撒花~