用 Python 和 GDAL.Warp 完成矢量裁剪栅格（遥感影像掩膜、SHP 裁剪 TIF、图像裁剪）

矢量裁剪栅格：使用 Python 和 GDAL.Warp 的分步指南

什么是矢量裁剪栅格？

在遥感、GIS 和图像处理中，裁剪栅格是指从原始栅格数据中提取特定感兴趣区域（ROI）的过程。而矢量裁剪栅格是一种利用矢量数据（例如 Shapefile）定义 ROI，然后从中裁剪栅格数据的方法。

为什么使用矢量裁剪栅格？

矢量裁剪栅格有几个优点：

• 精确性：矢量数据提供了精确的边界定义，确保裁剪结果的准确性。
• 效率：与像素级裁剪相比，矢量裁剪可以显著提高处理效率，尤其是对于大型栅格数据集。
• 灵活性：矢量数据易于编辑和更新，使您可以轻松调整 ROI 并重新裁剪栅格。

使用 Python 和 GDAL.Warp 进行矢量裁剪栅格

本指南将指导您使用 Python 和 GDAL.Warp 库分步执行矢量裁剪栅格操作。

步骤 1：安装 GDAL

在开始之前，请确保您的计算机上已安装 GDAL。可以在其官方网站上找到适用于不同平台的安装说明。

步骤 2：导入必需的库

在 Python 脚本中，导入 GDAL 和 NumPy 库：

import gdal
import numpy as np

步骤 3：读取栅格和矢量数据

使用 GDAL 打开原始栅格数据和用于裁剪的矢量数据：

raster_path = 'path/to/raster.tif'
raster_ds = gdal.Open(raster_path)

vector_path = 'path/to/vector.shp'
vector_ds = gdal.OpenEx(vector_path, gdal.OF_VECTOR)

步骤 4：定义裁剪范围

从矢量数据中提取裁剪范围：

vector_layer = vector_ds.GetLayer()
extent = vector_layer.GetExtent()

步骤 5：创建输出栅格

创建新栅格以存储裁剪后的数据：

output_path = 'path/to/output.tif'
output_driver = gdal.GetDriverByName('GTiff')
output_ds = output_driver.Create(output_path, width, height, raster_ds.RasterCount, raster_ds.GetRasterBand(1).DataType)

步骤 6：设置输出栅格属性

将裁剪范围、投影和其他属性应用于输出栅格：

output_ds.SetProjection(raster_ds.GetProjection())
output_ds.SetGeoTransform(extent)

步骤 7：执行裁剪操作

使用 GDAL.Warp 工具执行裁剪：

gdal.Warp(output_ds, raster_ds, cutlineDSName=vector_path, cropToCutline=True)

步骤 8：清理

关闭所有打开的数据集：

raster_ds = None
vector_ds = None
output_ds = None

示例代码

以下代码演示了如何使用上述步骤执行矢量裁剪栅格操作：

import gdal
import numpy as np

# 读取栅格数据
raster_path = 'path/to/raster.tif'
raster_ds = gdal.Open(raster_path)

# 读取矢量数据
vector_path = 'path/to/vector.shp'
vector_ds = gdal.OpenEx(vector_path, gdal.OF_VECTOR)

# 定义裁剪范围
vector_layer = vector_ds.GetLayer()
extent = vector_layer.GetExtent()

# 创建输出栅格
output_path = 'path/to/output.tif'
output_driver = gdal.GetDriverByName('GTiff')
output_ds = output_driver.Create(output_path, width, height, raster_ds.RasterCount, raster_ds.GetRasterBand(1).DataType)

# 设置输出栅格属性
output_ds.SetProjection(raster_ds.GetProjection())
output_ds.SetGeoTransform(extent)

# 执行裁剪操作
gdal.Warp(output_ds, raster_ds, cutlineDSName=vector_path, cropToCutline=True)

# 清理
raster_ds = None
vector_ds = None
output_ds = None

常见问题解答

• 问：我可以使用其他库来进行矢量裁剪吗？
  答：是的，除了 GDAL，还有其他库可以执行矢量裁剪，例如 rasterio 和 Shapely。

• 问：如何裁剪包含多个部分的矢量数据？
  答：使用 gdal.Warp 时，您可以通过将 cutlineDSNameOrLayer 设置为代表要裁剪的多部分矢量的矢量层来裁剪包含多个部分的矢量。

• 问：裁剪后的栅格会保留源栅格的属性吗？
  答：是的，裁剪后的栅格将继承源栅格的投影、地理变换和统计信息。

• 问：我可以使用矢量裁剪栅格执行镶嵌吗？
  答：虽然矢量裁剪不直接支持镶嵌，但可以通过分步裁剪并合并结果来实现。

• 问：如何优化矢量裁剪栅格的性能？
  答：对于大型数据集，建议使用分块方法并优化矢量数据的拓扑。