本文介绍如何使用Python提取栅格图像（例如Landsat5影像生成的NDVI图像）中，特定多边形区域内和区域外的NDVI均值。我们将使用Rasterio和Fiona这两个强大的库，通过读取Shapefile文件获取多边形边界，然后使用掩膜操作提取指定区域的NDVI值，最后计算均值。
准备工作
在开始之前，请确保已安装以下Python库：
- Rasterio: 用于读取和写入栅格数据。
- Fiona: 用于读取和写入矢量数据，例如Shapefile。
- Numpy: 用于数值计算，例如计算均值。
可以使用pip安装这些库：
pip install rasterio fiona numpy
提取多边形内部NDVI均值
以下代码演示了如何提取Shapefile文件定义的多边形内部的NDVI均值。
立即学习“Python免费学习笔记（深入）”；
import rasterio
import fiona
import rasterio.mask
import numpy as np
# 定义输入文件路径
shapefile_path = "path/to/your/shapefile.shp" # 替换为你的Shapefile文件路径
raster_path = "path/to/your/ndvi.tif" # 替换为你的NDVI栅格文件路径
try:
# 加载Shapefile文件，读取多边形几何信息
with fiona.open(shapefile_path, "r") as sf:
shapes = [feature["geometry"] for feature in sf]
# 打开栅格文件
with rasterio.open(raster_path) as src:
# 使用掩膜提取多边形内部的像素值
out_image, out_transform = rasterio.mask.mask(src, shapes, crop=True)
# 计算NDVI均值
NDVI_mean = np.mean(out_image)
print(f"多边形内部NDVI均值: {NDVI_mean}")
except fiona.errors.DriverError as e:
print(f"Fiona 错误: 无法打开 Shapefile 文件。请检查文件路径和文件格式。\n{e}")
except rasterio.RasterioIOError as e:
print(f"Rasterio 错误: 无法打开栅格文件。请检查文件路径和文件格式。\n{e}")
except Exception as e:
print(f"发生未知错误: {e}")代码解释：
- 导入必要的库: 导入 rasterio, fiona, rasterio.mask 和 numpy。
- 定义文件路径: shapefile_path 变量存储Shapefile文件的路径，raster_path 变量存储NDVI栅格文件的路径。请根据实际情况修改这些路径。
- 加载Shapefile文件: 使用 fiona.open() 函数打开Shapefile文件，并读取所有要素的几何信息。 shapes 变量是一个包含所有多边形几何对象的列表。
- 打开栅格文件: 使用 rasterio.open() 函数打开NDVI栅格文件。
- 使用掩膜提取像素值: rasterio.mask.mask() 函数使用Shapefile中的多边形作为掩膜，提取栅格图像中多边形内部的像素值。 crop=True 参数表示裁剪输出图像到掩膜的范围。 out_image 变量存储提取的像素值，out_transform 变量存储输出图像的地理变换信息。
- 计算NDVI均值: 使用 numpy.mean() 函数计算提取的像素值的均值。
- 打印结果: 打印计算得到的NDVI均值。
- 异常处理: 使用 try...except 块来捕获可能发生的异常，例如文件无法打开或文件格式错误。
提取多边形外部NDVI均值
要提取多边形外部的NDVI均值，可以使用 rasterio.mask.mask() 函数的 invert=True 参数。
import rasterio
import fiona
import rasterio.mask
import numpy as np
# 定义输入文件路径
shapefile_path = "path/to/your/shapefile.shp" # 替换为你的Shapefile文件路径
raster_path = "path/to/your/ndvi.tif" # 替换为你的NDVI栅格文件路径
try:
# 加载Shapefile文件，读取多边形几何信息
with fiona.open(shapefile_path, "r") as sf:
shapes = [feature["geometry"] for feature in sf]
# 打开栅格文件
with rasterio.open(raster_path) as src:
# 使用掩膜提取多边形外部的像素值，设置 invert=True
out_image, out_transform = rasterio.mask.mask(src, shapes, crop=True, invert=True)
# 计算NDVI均值
NDVI_mean = np.mean(out_image)
print(f"多边形外部NDVI均值: {NDVI_mean}")
except fiona.errors.DriverError as e:
print(f"Fiona 错误: 无法打开 Shapefile 文件。请检查文件路径和文件格式。\n{e}")
except rasterio.RasterioIOError as e:
print(f"Rasterio 错误: 无法打开栅格文件。请检查文件路径和文件格式。\n{e}")
except Exception as e:
print(f"发生未知错误: {e}")注意：
- invert=True 参数指示 rasterio.mask.mask() 函数提取多边形 外部 的像素值。
- 确保Shapefile文件和栅格文件位于正确的路径，并且具有正确的格式。
同时提取多个多边形区域的NDVI均值
如果要提取多个多边形区域的NDVI均值，只需在Shapefile文件中包含多个要素（多边形）。 rasterio.mask.mask() 函数会自动处理Shapefile文件中的所有多边形。
总结
本文介绍了如何使用Python和Rasterio、Fiona库提取栅格图像中多边形区域内和区域外的NDVI均值。 这些技术可以应用于各种遥感分析任务，例如土地覆盖分类、植被监测和环境评估。 通过灵活运用 rasterio.mask.mask() 函数，可以方便地提取指定区域的栅格数据，并进行进一步的分析。
