Python使用GDAL矢量裁剪栅格，设置背景值为空白（已解决）

一、使用gdal.Warp

gdalwarp 实用程序是一种图像拼接、重投影和扭曲实用程序。该程序可以重新投影到任何支持的投影。如果图像是带有控制信息的“原始”图像，也可以存储原始图像的投影信息。
gdal.warp官网链接：
https://gdal.org/api/python/osgeo.gdal.html#osgeo.gdal.Warp

def extractTiff(input_shape,input_raster,output_raster):#裁剪input_raster = gdal.Open(input_raster)# 开始裁剪result = gdal.Warp(output_raster,input_raster,format='GTiff',cutlineDSName=input_shape,cropToCutline=True)result.FlushCache()del result

gdal.Warp主要参数如下：

gdal.Warp(options = [], format = 'GTiff', outputBounds = None, outputBoundsSRS = one, xRes = None, yRes = None, targetAlignedPixels = False, width = 0, height = 0, srcSRS = None,dstSRS = None, srcAlpha = False, dstAlpha = False, warpOptions = None,errorThreshold = None, warpMemoryLimit = None, creationOptions = None,outputType = GDT_Unknown, workingType = 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, ropToCutline = False, copyMetadata = True, metadataConflictValue = None,setColorInterpretation = False, callback = None, callback_data = None):
其中：
options — 可以是一个字符串数组，一个字符串或者令其为空值，但是使用后面其他的参数来定义。
format — 输出的格式 (例如"GTiff"等)。
outputBounds — 在目标空间参考系统的输出数据集的范围，形式为 (minX,minY, maxX, maxY) 。
outputBoundsSRS — 如果在dstSRS中没有定义的话，使用这个关键字定义输出数据集的边界的空间参考系统。
xRes, yRes — 在目标参考系统中的像元大小。
targetAlignedPixels —是否强制输出边界为输出分辨率的倍数。
width — 输出栅格的像素列数。
height — 输出栅格的像素行数。
srcSRS —源空间参考系统。
dstSRS — 输出空间参考系统。
srcAlpha — 是否强制将输入数据集的最后一个波段作为alpha波段。
dstAlpha — 是否强制创建一个输出数据集的alpha波段。
outputType — 输出类型 (例如gdal.GDT_Byte等)
workingType — working type (gdal.GDT_Byte, etc…)
warpOptions —变形选项列表。
errorThreshold --近似转换的误差阈值(用像素表示) 。
warpMemoryLimit — 工作缓存大小，单位是bytes。
resampleAlg — 重采样模式。
creationOptions — 创建选项列表。
srcNodata — 源数据的nodata值。
dstNodata — 输出数据的nodata值。
multithread — 是否多线程计算和输入输出操作。
tps— 是否使用Thin Plate Spline GCP 转换器。
rpc— 是否使用RPC转换器。
geoloc — 是否使用GeoLocation数组转换器。
polynomialOrder — 多项式GCP插值的阶数。
transformerOptions — 转换参数
cutlineDSName — 剪切线数据集名称。这里的剪切线是指对影像进行剪切的时候所使用的矢量图层。
cutlineLayer — 剪切线图层名称。
cutlineWhere — 剪切线的WHERE语句。
cutlineSQL — 剪切线的SQL 语句。
cutlineBlend — 以像素为单位的剪切线混合距离。
cropToCutline — 是否使用剪切线的extent作为输出的界线。
copyMetadata — 是否拷贝源数据的元数据。
metadataConflictValue — 元数据冲突值。
setColorInterpretation — 是否强制将输入波段的颜色解释赋予输出波段。
callback — 回调函数。
callback_data — 回调函数数据

cropToCutline为True时，使用剪切线的extent作为输出的界线；dstNodata='NULL’时，可将背景值设置为空。此方法非常好用，强烈推荐。
注意cropToCutline默认为False，默认界限为输入矢量图的界限，使用默认设置会导致裁剪后的图片非常大（见下图设置不同的dstNodata后的文件大小）。另外还可以参考这篇文章的说明：python批量遥感影像裁剪的三种方法

二、使用PIL结合GDAL

基本原理可以参考：GDAL和 Pillow 的互操作
网上公开的代码是通过PIL与GDAL栅格数据互操作，将栅格数据转换为PIL中的对象，使用PyShp读取矢量边界，裁剪PIL，最后通过数组，转回栅格数据。以下是完整代码

#!/usr/bin/env python
# -*- encoding: utf-8 -*-
'''
@File    :   k_factor.py
@Time    :   2022/09/26 08:45:09
@Author  :   Junlai
@Version :   1.0
@Desc    :   None
'''
import os,shapefile
from osgeo import gdal, gdal_array, ogr
from PIL import Image, ImageDrawdef image2Array(i):"""将一个Python图像库的数组转换为一个gdal_array图片"""a = gdal_array.numpy.frombuffer(i.tobytes(), 'b')a.shape = i.im.size[1], i.im.size[0]return adef world2Pixel(geoMatrix, x, y):"""使用GDAL库的geomatrix对象((gdal.GetGeoTransform()))计算地理坐标的像素位置"""ulx = geoMatrix[0]uly = geoMatrix[3]xDist = geoMatrix[1]yDist = geoMatrix[5]rtnX = geoMatrix[2]rtnY = geoMatrix[4]pixel = int((x - ulx) / xDist)line = int((uly - y) / abs(yDist))return (pixel, line)def extractFromShp(input_raster,input_shape,output_raster):# 用于裁剪的栅格数据raster = input_raster# 裁剪后的栅格数据output = output_raster# 将数据源作为gdal_array载入srcArray = gdal_array.LoadFile(raster)# 同时载入gdal库的图片从而获取geotransformsrcImage = gdal.Open(raster)geoTrans = srcImage.GetGeoTransform()geoProj = srcImage.GetProjection()# # 使用PyShp库打开shp文件shp = input_shaper = shapefile.Reader("{}.shp".format(shp))# # 将图层扩展转换为图片像素坐标minX, minY, maxX, maxY = r.bboxulX, ulY = world2Pixel(geoTrans, minX, maxY)lrX, lrY = world2Pixel(geoTrans, maxX, minY)# 计算新图片的尺寸pxWidth = int(lrX - ulX)pxHeight = int(lrY - ulY)clip = srcArray[ulY:lrY, ulX:lrX]# 为图片创建一个新的geomatrix对象以便附加地理参照数据geoTrans = list(geoTrans)geoTrans[0] = minXgeoTrans[3] = maxY# 在一个空白的8字节黑白掩膜图片上把点映射为像元# 边界线pixels = []for p in r.shape(0).points:pixels.append(world2Pixel(geoTrans, p[0], p[1]))rasterPoly = Image.new("L", (pxWidth, pxHeight), 1)# 使用PIL创建一个空白图片用于绘制多边形rasterize = ImageDraw.Draw(rasterPoly)rasterize.polygon(pixels, 0)# 使用PIL图片转换为Numpy掩膜数组mask = image2Array(rasterPoly)# 根据掩膜图层对图像进行裁剪clip = gdal_array.numpy.choose(mask, (clip, 0)).astype(gdal_array.numpy.uint8)print(clip.max())# 保存为tiff文件gdal_array.SaveArray(clip, "{}.tif".format(output),format="GTiff")if __name__ == '__main__':pj_path = os.path.abspath(os.path.dirname(os.path.abspath(os.path.dirname(__file__))))fpath = os.path.join(pj_path, 'data')input_shape = f'{fpath}\\boundary\\input'output_raster = f'{fpath}\\out'# tif输入路径，打开文件input_raster = f'{fpath}\\lu\\input.tif'# 矢量裁剪栅格extractFromShp(input_raster,input_shape,output_raster)

该方法存在一个小bug，当矢量文件属性表存在中文，会出现以下报错，

解决方法可以简单粗暴，直接修改PyShp源码：
self.encoding = kwargs.pop('encoding', 'utf-8')修改为self.encoding = kwargs.pop('encoding', 'gbk')
笔者尝试直接用ogr替代PyShp，完成上述操作，完全可行，结果一致，在此不表，感兴趣的可以自行尝试。

三、使用ArcTools裁剪

过于简单粗暴，在此不表

四、BUG

有个小的BUG，暂时还没摸清楚原因，MARK一下：
使用GDAL结合PIL裁剪后的栅格图，会比ArcTool少一行，GDAL.WRAP裁剪时会少两行一列。

使用ArcTool裁剪
使用PIL结合GDAL裁剪时
使用GDAL.WRAP裁剪时