Python提取多边形内外NDVI均值方法

本文详细介绍了如何利用Python高效提取栅格数据（如NDVI图像）中多边形区域内外NDVI均值，为遥感数据分析提供了一种可重复且精准的方法。文章重点讲解了如何运用`rasterio`和`fiona`库读取栅格和矢量数据，并利用`rasterio.mask.mask`函数实现基于多边形的栅格裁剪与掩膜操作。此外，还阐述了如何处理多个多边形，以及计算裁剪区域内外NDVI均值，同时考虑了数据类型一致性和坐标系一致性等关键问题。通过本文的学习，读者能够掌握使用Python进行遥感数据分析的核心技术，为更复杂的遥感应用奠定基础。

本文介绍了如何使用Python编程提取栅格图像（如Landsat影像生成的NDVI图像）中，特定多边形区域内部和外部的NDVI均值。通过rasterio和fiona库，可以轻松实现对栅格数据的裁剪、掩膜以及统计计算，为遥感数据分析提供了一种高效且可重复的方法。

准备工作

在开始之前，需要确保安装了以下Python库：

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

可以使用pip命令安装这些库：

pip install rasterio fiona numpy

核心步骤

1. 读取矢量数据 (Shapefile)

   首先，需要使用fiona库读取包含多边形信息的shapefile。以下代码展示了如何读取shapefile并提取几何信息：

   import fiona
   
   shapefile_path = "path/to/your/shapefile.shp"  # 替换为你的shapefile路径
   
   with fiona.open(shapefile_path, "r") as sf:
       shapes = [feature["geometry"] for feature in sf]

   这段代码打开指定的shapefile，并提取每个要素的几何信息，存储在shapes列表中。shapes列表将包含代表多边形的几何对象。

2. 读取栅格数据 (TIFF)

   接下来，使用rasterio库读取栅格数据（例如，NDVI图像）。以下代码展示了如何读取TIFF图像：

   import rasterio
   
   raster_path = "path/to/your/ndvi.tif"  # 替换为你的NDVI图像路径
   
   with rasterio.open(raster_path) as src:
       raster_image = src.read(1) # 读取第一个波段，假设NDVI数据存储在第一个波段
       raster_meta = src.meta.copy() # 复制元数据，后续写入使用

   这段代码打开指定的TIFF图像，并将图像数据读取到raster_image变量中。同时，复制图像的元数据，以便后续写入处理后的数据。

3. 使用多边形裁剪栅格数据

   使用rasterio.mask.mask函数，可以根据多边形裁剪栅格数据。以下代码展示了如何裁剪栅格数据，并计算裁剪区域内的NDVI均值：

   import rasterio.mask
   import numpy as np
   
   with rasterio.open(raster_path) as src:
       out_image, out_transform = rasterio.mask.mask(src, shapes, crop=True)
       out_image = out_image.astype(src.dtypes[0]) # 确保数据类型一致
   
   # 计算裁剪区域内的NDVI均值，忽略无效值
   NDVI_mean = np.nanmean(out_image)
   print(f"多边形内部NDVI均值: {NDVI_mean}")

   这段代码使用shapes列表中的多边形裁剪栅格图像。crop=True参数表示裁剪到最小范围。np.nanmean用于计算均值，并忽略NaN值。

4. 提取多边形外部的NDVI值

   可以使用rasterio.mask.mask函数，通过设置invert=True参数来提取多边形外部的像素值。以下代码展示了如何提取多边形外部的NDVI值：

   with rasterio.open(raster_path) as src:
       out_image, out_transform = rasterio.mask.mask(src, shapes, crop=True, invert=True)
       out_image = out_image.astype(src.dtypes[0]) # 确保数据类型一致
   
   # 计算裁剪区域外的NDVI均值，忽略无效值
   NDVI_mean_outside = np.nanmean(out_image)
   print(f"多边形外部NDVI均值: {NDVI_mean_outside}")

   通过设置invert=True，rasterio.mask.mask函数将返回多边形外部的像素值。

5. 分别处理多个多边形

   如果需要分别处理多个多边形，可以循环遍历shapes列表，并对每个多边形执行裁剪和均值计算操作。

   for i, shape in enumerate(shapes):
       with rasterio.open(raster_path) as src:
           out_image, out_transform = rasterio.mask.mask(src, [shape], crop=True)
           out_image = out_image.astype(src.dtypes[0]) # 确保数据类型一致
   
       NDVI_mean = np.nanmean(out_image)
       print(f"多边形 {i+1} 内部NDVI均值: {NDVI_mean}")
   
       with rasterio.open(raster_path) as src:
           out_image, out_transform = rasterio.mask.mask(src, [shape], crop=True, invert=True)
           out_image = out_image.astype(src.dtypes[0]) # 确保数据类型一致
   
       NDVI_mean_outside = np.nanmean(out_image)
       print(f"多边形 {i+1} 外部NDVI均值: {NDVI_mean_outside}")

   这段代码循环遍历shapes列表，并对每个多边形分别计算内部和外部的NDVI均值。

完整示例代码

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图像路径

# 读取矢量数据
with fiona.open(shapefile_path, "r") as sf:
    shapes = [feature["geometry"] for feature in sf]

# 循环处理每个多边形
for i, shape in enumerate(shapes):
    # 多边形内部
    with rasterio.open(raster_path) as src:
        out_image, out_transform = rasterio.mask.mask(src, [shape], crop=True)
        out_image = out_image.astype(src.dtypes[0]) # 确保数据类型一致

    NDVI_mean = np.nanmean(out_image)
    print(f"多边形 {i+1} 内部NDVI均值: {NDVI_mean}")

    # 多边形外部
    with rasterio.open(raster_path) as src:
        out_image, out_transform = rasterio.mask.mask(src, [shape], crop=True, invert=True)
        out_image = out_image.astype(src.dtypes[0]) # 确保数据类型一致

    NDVI_mean_outside = np.nanmean(out_image)
    print(f"多边形 {i+1} 外部NDVI均值: {NDVI_mean_outside}")

注意事项

• 确保shapefile和栅格图像的坐标系一致。如果坐标系不一致，需要进行投影转换。
• rasterio.mask.mask函数返回的out_image是一个三维数组，第一个维度的大小为1，需要使用out_image[0]来访问二维图像数据。
• 在计算均值时，使用np.nanmean可以忽略无效值（例如，NaN值）。
• 根据实际需求，可以修改代码以提取其他统计信息，例如最大值、最小值、标准差等。
• 需要确保栅格数据类型与计算时的数据类型一致，如果数据类型不一致，可能会导致计算结果错误，使用 out_image = out_image.astype(src.dtypes[0]) 可以保证数据类型一致。

总结

本文介绍了如何使用Python的rasterio和fiona库提取栅格数据中多边形内外NDVI均值。通过这些方法，可以高效地进行遥感数据分析，并提取有用的信息。 掌握这些技术，可以为更复杂的遥感应用打下坚实的基础。

今天关于《Python提取多边形内外NDVI均值方法》的内容介绍就到此结束，如果有什么疑问或者建议，可以在golang学习网公众号下多多回复交流；文中若有不正之处，也希望回复留言以告知！