多边形内提取NDVI值Python教程

本文旨在提供一个使用Python从栅格图像中提取多边形区域内NDVI值的实用指南。我们将使用Rasterio和Fiona库，演示如何加载栅格数据和矢量数据，并利用掩膜操作提取特定区域的NDVI平均值。此外，还将介绍如何提取多边形外部的NDVI值，为遥感影像分析提供有效的方法。

从栅格图像中提取多边形区域的NDVI值

本教程将指导您如何使用Python从栅格图像中提取特定多边形区域内的NDVI（归一化植被指数）值。我们将使用 rasterio 和 fiona 这两个强大的库来处理栅格数据和矢量数据。

准备工作

在开始之前，请确保您已安装以下库：

• rasterio: 用于读取和写入栅格数据。
• fiona: 用于读取和写入矢量数据（如Shapefile）。
• numpy: 用于数值计算。

您可以使用 pip 安装这些库：

pip install rasterio fiona numpy

步骤详解

1. 导入必要的库

   首先，导入所需的库：

   import rasterio
   import fiona
   import rasterio.mask
   import numpy as np

2. 加载矢量数据（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 列表中。每个要素的几何信息都表示一个多边形。

3. 加载栅格数据（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 变量包含裁剪后的图像的地理变换信息。

4. 计算NDVI平均值

   使用 numpy 库计算提取的NDVI值的平均值。

   NDVI_mean = np.mean(out_image)
   print(f"NDVI 平均值: {NDVI_mean}")

   这段代码计算 out_image 中所有像素值的平均值，并将结果打印出来。

5. 提取多边形外部的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 函数，您可以轻松地提取特定区域的像素值，并进行后续分析。这种方法在遥感影像分析、土地覆盖分类等领域具有广泛的应用。