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