NDVI（归一化植被指数）是遥感分析中最常用的植被健康度指标，通过近红外波段与红光波段的反射率差异来量化植被覆盖程度。在农业监测、生态评估、城市规划等实际项目中，我们往往需要从特定地理边界（如地块、行政区、保护区）中提取 NDVI 的统计值，而不仅仅是单像素的读取。将矢量多边形与栅格遥感影像结合，实现分区统计，是这一需求的标准解决路径。
NDVI 计算原理与数据准备
NDVI 的计算公式为 (NIR - Red) / (NIR + Red)，其中 NIR 代表近红外波段反射率，Red 代表红光波段反射率。健康植被在近红外波段具有高反射率，在红光波段由于叶绿素吸收而反射率较低，因此 NDVI 值趋近于 1；裸土和水体则接近 0 甚至负值。理解这一原理有助于我们在后续处理中正确设置有效值范围，通常将 NDVI 限制在 -1 到 1 之间，超出此范围的值往往是传感器噪声或云遮挡导致的异常数据。
在数据准备阶段，我们需要两类核心数据：一是包含近红外和红光波段的栅格影像，常见来源包括 Landsat、Sentinel-2 或 MODIS；二是定义感兴趣区域的矢量多边形文件，通常是 GeoJSON 或 Shapefile 格式。以 Sentinel-2 为例，其 Band 4 对应红光波段，Band 8 对应近红外波段，空间分辨率为 10 米，非常适合中尺度植被监测。矢量数据则需要包含每个多边形的唯一标识字段，便于后续统计结果与原始地块一一对应。
数据准备中一个容易被忽视的环节是坐标参考系统（CRS）的统一。栅格影像和矢量多边形可能采用不同的投影坐标系，例如影像常用 UTM 投影，而矢量数据可能是 WGS84 地理坐标系。如果两者 CRS 不一致，直接进行空间操作会导致位置偏移甚至报错。因此，在正式提取之前，必须将矢量数据重投影到与栅格影像一致的 CRS 上，这一步骤使用 geopandas 的 to_crs 方法即可完成。
使用 rasterio 与 geopandas 读取并处理数据
rasterio 是 Python 生态中处理栅格数据的主流库，它基于 GDAL 提供了简洁的 Pythonic 接口。geopandas 则用于矢量数据的读取与操作，两者配合可以高效完成栅格与矢量的联合处理。下面演示如何读取 Sentinel-2 的 Band 4 和 Band 8 影像，并计算 NDVI 矩阵。读取栅格数据时需要注意 profile 属性中保存的变换矩阵和 CRS 信息，这些元数据在后续空间操作中至关重要。
import rasterio
import numpy as np
import geopandas as gpd
# 读取红光波段（Band 4）
with rasterio.open('B04.tif') as red_src:
red = red_src.read(1).astype(float)
profile = red_src.profile
raster_crs = red_src.crs
# 读取近红外波段（Band 8）
with rasterio.open('B08.tif') as nir_src:
nir = nir_src.read(1).astype(float)
# 计算 NDVI，注意分母为零的防护
np.seterr(divide='ignore', invalid='ignore')
ndvi = (nir - red) / (nir + red)
ndvi = np.nan_to_num(ndvi, nan=0.0)
# 将 NDVI 矩阵写入新的栅格文件
profile.update(dtype=rasterio.float32, count=1)
with rasterio.open('ndvi.tif', 'w', **profile) as dst:
dst.write(ndvi.astype(rasterio.float32), 1)
# 读取矢量多边形并重投影到栅格CRS
polygons = gpd.read_file('fields.geojson')
polygons = polygons.to_crs(raster_crs)
print(f'共加载 {len(polygons)} 个多边形')
上述代码中，计算 NDVI 时使用了 np.seterr 来抑制除零警告，并通过 nan_to_num 将无效值替换为 0。实际项目中，更严谨的做法是将无效像素设为 np.nan，这样在后续统计中可以自动排除。此外，将 NDVI 结果写入栅格文件是一个可选但推荐的步骤，方便后续多次复用，避免重复计算。矢量数据的重投影通过 to_crs 方法完成，参数传入栅格影像的 CRS 对象即可。
在数据量较大的场景下，直接将整个影像读入内存可能导致内存不足。rasterio 支持窗口读取（windowed reading），可以按需读取特定区域的像素。如果多边形分布分散且覆盖范围远小于影像范围，可以先计算所有多边形的联合边界框，然后只读取该范围内的栅格数据，显著降低内存占用和计算时间。
基于 rasterstats 执行多边形分区统计
rasterstats 库专门用于在栅格数据上执行矢量多边形的分区统计（Zonal Statistics），它封装了底层的几何运算和像素遍历逻辑，提供了简洁的 API 接口。相比手动实现掩膜提取和统计计算，rasterstats 在正确性和效率上都有明显优势，尤其适合批量处理大量多边形的场景。它支持均值、中位数、标准差、最小值、最大值等多种统计指标，还可以自定义统计函数。
from rasterstats import zonal_stats
import json
# 执行分区统计
stats = zonal_stats(
polygons,
'ndvi.tif',
stats=['mean', 'median', 'min', 'max', 'std'],
geojson_out=True,
nodata=np.nan
)
# 将统计结果合并到原始GeoDataFrame
result_gdf = gpd.GeoDataFrame.from_features(stats)
print(result_gdf[['name', 'mean', 'median', 'min', 'max', 'std']].head(10))
# 导出结果为GeoJSON文件
result_gdf.to_file('ndvi_stats.geojson', driver='GeoJSON')
上述代码中，zonal_stats 函数的第一个参数是矢量多边形，可以是 GeoDataFrame、GeoJSON 字符串或文件路径。第二个参数是栅格文件路径。stats 参数指定需要计算的统计指标。geojson_out=True 表示输出结果包含原始几何信息和属性字段，方便直接导出为完整的矢量文件。nodata 参数用于指定栅格中的无效值，确保统计时自动排除。
统计结果默认以字典列表形式返回，每个多边形对应一个包含统计指标的字典。设置 geojson_out=True 后，结果变为 GeoJSON FeatureCollection 格式，可以直接转换为 GeoDataFrame。最终结果中，每个多边形除了原始属性外，还会新增 mean、median、min、max、std 等字段，分别代表该多边形内 NDVI 的各项统计值。导出为 GeoJSON 后，可以在 QGIS 或 ArcGIS 中可视化查看各区域的植被状况。
结果验证与常见问题排查
完成统计后，务必对结果进行验证。一个有效的验证方法是随机选取几个多边形，在 GIS 软件中手动绘制其 NDVI 直方图，对比 Python 计算的统计值是否一致。常见的问题包括：统计值为 None，通常是因为多边形与栅格范围无交集，或者 CRS 不一致导致几何位置偏移；统计值全部相同，可能是栅格读取时波段选择错误，导致 NDVI 计算公式输入了相同的波段。
另一个需要注意的问题是像元中心对齐。rasterstats 默认采用像元中心法判断一个像素是否属于某个多边形，这意味着位于多边形边界上的像素可能被排除或包含。如果多边形面积较小且仅覆盖少数像素，这种边界效应会显著影响统计结果。在这种情况下，可以通过 all_touched=True 参数改为包含所有被多边形触及的像素，代价是可能引入边界像素的噪声。选择哪种策略取决于项目对精度的要求。
# 使用all_touched策略处理小面积多边形
stats_all_touched = zonal_stats(
polygons,
'ndvi.tif',
stats=['mean', 'count'],
all_touched=True,
nodata=np.nan
)
# 检查是否有统计结果为None的多边形
none_count = sum(1 for s in stats_all_touched if s['mean'] is None)
print(f'统计结果为None的多边形数量: {none_count}')
# 自定义统计函数：计算NDVI大于0.5的像素占比
def healthy_vegetation_ratio(arr):
valid = arr[~np.isnan(arr)]
if len(valid) == 0:
return None
return np.sum(valid > 0.5) / len(valid)
custom_stats = zonal_stats(
polygons,
'ndvi.tif',
add_stats={'healthy_ratio': healthy_vegetation_ratio},
nodata=np.nan
)
自定义统计函数是 rasterstats 的高级功能，通过 add_stats 参数传入一个字典，键为统计指标名称，值为接收一维数组的函数。上面的示例计算了每个多边形内 NDVI 大于 0.5 的像素占比，这一指标在实际农业监测中非常有用，可以量化健康植被的覆盖比例。通过灵活运用自定义统计函数，可以根据具体业务需求提取更有针对性的指标，而不仅限于基础的均值和极值。