本文旨在提供一个使用Python从栅格图像中提取多边形区域内NDVI值的实用指南。我们将使用Rasterio和Fiona库，演示如何加载栅格数据和矢量数据，并利用掩膜操作提取特定区域的NDVI平均值。此外，还将介绍如何提取多边形外部的NDVI值，为遥感影像分析提供有效的方法。
从栅格图像中提取多边形区域的NDVI值
本教程将指导您如何使用Python从栅格图像中提取特定多边形区域内的NDVI（归一化植被指数）值。我们将使用 rasterio 和 fiona 这两个强大的库来处理栅格数据和矢量数据。
准备工作
在开始之前，请确保您已安装以下库：
- rasterio: 用于读取和写入栅格数据。
- fiona: 用于读取和写入矢量数据（如Shapefile）。
- numpy: 用于数值计算。
您可以使用 pip 安装这些库：
立即学习“Python免费学习笔记（深入）”；
pip install rasterio fiona numpy
步骤详解
-
导入必要的库
首先，导入所需的库：
import rasterio import fiona import rasterio.mask import numpy as np
-
加载矢量数据（Shapefile）
使用 fiona 库加载包含多边形边界的Shapefile。假设您的Shapefile文件路径存储在 path_to_shapefile 变量中。
path_to_shapefile = "path/to/your/shapefile.shp" # 替换为你的shapefile路径 with fiona.open(path_to_shapefile, "r") as sf: shapes = [feature["geometry"] for feature in sf]这段代码打开Shapefile，并提取所有要素的几何信息，存储在 shapes 列表中。每个要素的几何信息都表示一个多边形。
-
加载栅格数据（TIFF）
使用 rasterio 库加载包含NDVI值的栅格图像（通常是TIFF格式）。假设您的TIFF文件路径存储在 source_tif 变量中。
source_tif = "path/to/your/ndvi.tif" # 替换为你的NDVI栅格图像路径 with rasterio.open(source_tif) as src: out_image, out_transform = rasterio.mask.mask(src, shapes, crop=True)这段代码打开TIFF文件，并使用 rasterio.mask.mask 函数提取多边形区域内的栅格数据。crop=True 参数表示裁剪输出图像到多边形边界。out_image 变量包含裁剪后的栅格数据，out_transform 变量包含裁剪后的图像的地理变换信息。
-
计算NDVI平均值
使用 numpy 库计算提取的NDVI值的平均值。
NDVI_mean = np.mean(out_image) print(f"NDVI 平均值: {NDVI_mean}")这段代码计算 out_image 中所有像素值的平均值，并将结果打印出来。
-
提取多边形外部的NDVI值
rasterio.mask.mask 函数还允许您提取多边形外部的像素值。只需将 invert 参数设置为 True 即可。
with rasterio.open(source_tif) as src: out_image, out_transform = rasterio.mask.mask(src, shapes, crop=True, invert=True) NDVI_mean_outside = np.mean(out_image) print(f"多边形外部 NDVI 平均值: {NDVI_mean_outside}")这段代码提取多边形外部的栅格数据，并计算其平均值。
完整代码示例
以下是一个完整的代码示例，展示了如何从栅格图像中提取多边形区域内和区域外的NDVI平均值：
import rasterio
import fiona
import rasterio.mask
import numpy as np
# 定义文件路径
path_to_shapefile = "path/to/your/shapefile.shp" # 替换为你的shapefile路径
source_tif = "path/to/your/ndvi.tif" # 替换为你的NDVI栅格图像路径
# 加载矢量数据
with fiona.open(path_to_shapefile, "r") as sf:
shapes = [feature["geometry"] for feature in sf]
# 提取多边形内部的NDVI值
with rasterio.open(source_tif) as src:
out_image, out_transform = rasterio.mask.mask(src, shapes, crop=True)
NDVI_mean_inside = np.mean(out_image)
print(f"多边形内部 NDVI 平均值: {NDVI_mean_inside}")
# 提取多边形外部的NDVI值
with rasterio.open(source_tif) as src:
out_image, out_transform = rasterio.mask.mask(src, shapes, crop=True, invert=True)
NDVI_mean_outside = np.mean(out_image)
print(f"多边形外部 NDVI 平均值: {NDVI_mean_outside}")注意事项
- 确保Shapefile和TIFF文件位于正确的路径。
- rasterio.mask.mask 函数返回的 out_image 是一个三维数组，其中第一个维度通常是波段数。如果您的栅格图像只有一个波段（例如，NDVI图像），则可以使用 out_image[0] 来访问实际的像素值。
- 如果您的Shapefile包含多个多边形，rasterio.mask.mask 函数将提取所有多边形内的像素值。如果需要提取特定多边形内的像素值，您需要修改 shapes 列表，只包含您感兴趣的多边形的几何信息。
- 在计算平均值之前，建议检查 out_image 中是否存在无效值（例如，NaN）。可以使用 numpy.nanmean 函数来忽略无效值。
总结
本教程介绍了如何使用Python和 rasterio、fiona 库从栅格图像中提取多边形区域内的NDVI值。通过使用 rasterio.mask.mask 函数，您可以轻松地提取特定区域的像素值，并进行后续分析。这种方法在遥感影像分析、土地覆盖分类等领域具有广泛的应用。
