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

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)

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

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

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

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

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

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

使用这三步,再次裁剪上述的三种情况,得到的结果:
在这里插入图片描述
本方法是最笨的一个方法,过程还需要额外生成两个shp文件,希望有看到的同仁,有好的建议不吝赐教!!

以下是使用python的ogr和osr库,根据很多坐标点使用gdal.warp函数裁剪栅格数据的示例代码: ```python import gdal import ogr import osr # 设置输入栅格数据路径 input_raster_path = '/path/to/input/raster.tif' # 设置输出栅格数据路径 output_raster_path = '/path/to/output/raster.tif' # 设置裁剪范围的坐标点 # 这里假设有一个shapefile文件,其中包含了多个矩形范围的坐标点 clip_shapefile_path = '/path/to/clip/shapefile.shp' # 打开输入栅格数据 input_raster = gdal.Open(input_raster_path) # 获取输入栅格数据的投影信息 input_projection = input_raster.GetProjection() # 打开裁剪范围的shapefile文件 clip_shapefile = ogr.Open(clip_shapefile_path) # 获取shapefile文件中的第一个图层 clip_layer = clip_shapefile.GetLayer(0) # 获取shapefile文件中每个图形的几何体 clip_geometries = [feature.GetGeometryRef() for feature in clip_layer] # 将输入栅格数据的投影信息转换为裁剪范围的投影信息 clip_spatial_reference = clip_layer.GetSpatialRef() input_spatial_reference = osr.SpatialReference(input_projection) transform = osr.CoordinateTransformation(input_spatial_reference, clip_spatial_reference) # 将裁剪范围的坐标点转换为输入栅格数据的投影坐标系中的坐标点 clip_bounds = [] for clip_geometry in clip_geometries: clip_geometry.Transform(transform) clip_bounds.append(clip_geometry.GetEnvelope()) # 使用gdal.warp函数裁剪栅格数据 gdal.Warp(output_raster_path, input_raster, cutlineDSName=clip_shapefile_path, cropToCutline=True) ``` 在以上代码中,使用gdal.Open函数打开输入栅格数据和裁剪范围的shapefile文件,然后使用gdal.Warp函数裁剪栅格数据。在裁剪栅格数据之前,需要将输入栅格数据的投影信息转换为裁剪范围的投影信息,并将裁剪范围的坐标点转换为输入栅格数据的投影坐标系中的坐标点。最后,使用gdal.Warp函数裁剪栅格数据。
评论 5
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值