【PyGIS】GDAL矢量裁剪栅格，设置背景值为空白_python 使用矢量裁剪

一、使用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
1
2
3
4
5
6
7
8
9
10
11

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 — 回调函数数据
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51

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, ImageDraw

def image2Array(i):
    """
    将一个Python图像库的数组转换为一个gdal_array图片
    """
    a = gdal_array.numpy.frombuffer(i.tobytes(), 'b')
    a.shape = i.im.size[1], i.im.size[0]
    return a

def 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库的图片从而获取geotransform
    srcImage = gdal.Open(raster)
    geoTrans = srcImage.GetGeoTransform()
    geoProj = srcImage.GetProjection()
    # # 使用PyShp库打开shp文件
    shp = input_shape
    r = shapefile.Reader("{}.shp".format(shp))
    # # 将图层扩展转换为图片像素坐标
    minX, minY, maxX, maxY = r.bbox
    ulX, 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] = minX
    geoTrans[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)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88

该方法存在一个小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裁剪时