python中利用gdal.warp完成矢量裁剪栅格（遥感影像掩膜、shp裁剪tif、图像裁剪）_python 矢量裁剪栅格

gdal库能够很方便的完成栅格数据的裁剪，代码如下：

import numpy as np
from osgeo import gdal

def extract_by_shp(in_shp_path, inputpath, outputpath):
    input_raster = gdal.Open(inputpath)
        # 利用gdal.Warp进行裁剪
    result = gdal.Warp(
        outputpath,
        input_raster,
        format = 'GTiff',
        cutlineDSName = in_shp_path, # 用于裁剪的矢量
        cropToCutline = False, # 是否使用cutlineDSName的extent作为输出的界线
        dstNodata = 0 # 输出数据的nodata值
        )
    del result
    
if __name__ == "__main__":
    in_shp_path=r"D:\边界.shp"
    inputpath=r"D:\S2_20220812.tif"
    outputpath=r"D:\result.tif"
    extract_by_shp(in_shp_path, inputpath, outputpath)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

但在使用过程中，本人发现其仍然有不足之处，导致裁剪的结果总是不如人意，原因在于gdal.warp方法的倒数第二个参数的设置：

cropToCutline：指输出的裁剪结果按照什么样的边界；
当参数设置为True时，将按照输入的矢量，即“边界.shp”的外接矩阵输出；
当参数设置为False时，将不按照矢量，而按照原来栅格的边界输出。
这其中，无效值的区域将被输出为背景值（下图中黑色部分）。

一般，面矢量与栅格的二维空间关系总共有三种，包含、被包含、相交。下图表示cropToCutline设置为不同值的时候的几种裁剪结果：

可以看出，在很多时候，矢量裁剪的栅格存在很多没必要的背景区域，虽然不影响后续的计算，但是增加行列数和图像的大小。
网上也看到了一些其他的解决方法，没有彻底看懂，加上时间紧迫，只想出了一个比较赘余的笨方法：

1.先将栅格数据转矢量；
2.将栅格转的矢量与输入的边界矢量就求交集；
3.使用上述交集裁剪原始的栅格数据，设置cropToCutline=“True”,得到结果。

其中：
第一步，栅格转矢量：python栅格遥感影像转矢量
第二部：矢量之间求交集：python两个shp面矢量之间求交集
第三步：裁剪按照本文开始所描述方法。

使用这三步，再次裁剪上述的三种情况，得到的结果：

本方法是最笨的一个方法，过程还需要额外生成两个shp文件，希望有看到的同仁，有好的建议不吝赐教！！