本文介绍了如何使用Python编程提取栅格图像（如Landsat影像生成的NDVI图像）中，特定多边形区域内部和外部的NDVI均值。通过rasterio和fiona库，可以轻松实现对栅格数据的裁剪、掩膜以及统计计算，为遥感数据分析提供了一种高效且可重复的方法。
准备工作
在开始之前，需要确保安装了以下Python库：
- rasterio: 用于读取和写入栅格数据。
- fiona: 用于读取矢量数据（例如，shapefile）。
- numpy: 用于数值计算，例如计算均值。
可以使用pip命令安装这些库：
pip install rasterio fiona numpy
核心步骤
-
读取矢量数据 (Shapefile)
立即学习“Python免费学习笔记（深入）”；
首先，需要使用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列表将包含代表多边形的几何对象。
-
读取栅格数据 (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变量中。同时，复制图像的元数据，以便后续写入处理后的数据。
-
使用多边形裁剪栅格数据
使用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值。
-
提取多边形外部的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函数将返回多边形外部的像素值。
-
分别处理多个多边形
如果需要分别处理多个多边形，可以循环遍历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均值。通过这些方法，可以高效地进行遥感数据分析，并提取有用的信息。 掌握这些技术，可以为更复杂的遥感应用打下坚实的基础。
