基于Python提取Landsat5 NDVI影像多边形内外均值并分析差异

以下是用Python完成需求的具体流程，依赖geopandas（矢量处理）、rasterio（栅格处理）、numpy（数值计算）、scipy（统计检验）这几个库。

1. 环境准备 #

先安装所需依赖：

pip install geopandas rasterio numpy scipy

2. 加载矢量与栅格数据 #

读取包含红、蓝多边形的shapefile，同时获取NDVI影像的地理范围：

import geopandas as gpd
import rasterio
from rasterio.mask import mask
import numpy as np
from scipy.stats import f_oneway
from shapely.geometry import box

# 读取shapefile，自行替换文件路径和区分红/蓝的属性列（比如'poly_color'）
gdf = gpd.read_file("你的矢量文件路径.shp")
red_poly = gdf[gdf['poly_color'] == 'red'].geometry.unary_union
blue_poly = gdf[gdf['poly_color'] == 'blue'].geometry.unary_union

# 读取任意一幅NDVI影像，获取影像的地理边界多边形
with rasterio.open("某幅NDVI影像.tif") as src:
    bounds = src.bounds
    extent_poly = box(bounds.left, bounds.bottom, bounds.right, bounds.top)
    # 确保矢量与栅格坐标系一致，不一致则转换矢量
    if gdf.crs != src.crs:
        gdf = gdf.to_crs(src.crs)
        red_poly = gdf[gdf['poly_color'] == 'red'].geometry.unary_union
        blue_poly = gdf[gdf['poly_color'] == 'blue'].geometry.unary_union

3. 构建「外部」多边形 #

用影像范围减去红、蓝多边形的并集，得到外部区域：

# 合并红、蓝多边形，再用范围多边形做差集
combined_poly = red_poly.union(blue_poly)
outer_poly = extent_poly.difference(combined_poly)

4. 提取各区域NDVI均值 #

编写通用函数，对每幅NDVI影像提取三个区域的均值：

def get_ndvi_mean(ndvi_path, polygon):
    with rasterio.open(ndvi_path) as src:
        # 掩膜提取区域内的NDVI数据，自动裁剪并忽略nodata
        out_img, _ = mask(src, [polygon], crop=True, nodata=src.nodata)
        # 过滤无效值后计算均值
        valid_data = out_img[out_img != src.nodata]
        return valid_data.mean() if len(valid_data) > 0 else np.nan

# 遍历3幅NDVI影像，批量提取均值
ndvi_files = ["ndvi_1.tif", "ndvi_2.tif", "ndvi_3.tif"]
region_means = []
for file in ndvi_files:
    means = {
        "影像文件名": file,
        "红色区域均值": get_ndvi_mean(file, red_poly),
        "蓝色区域均值": get_ndvi_mean(file, blue_poly),
        "外部区域均值": get_ndvi_mean(file, outer_poly)
    }
    region_means.append(means)

# 打印提取结果
for res in region_means:
    print(res)

5. 验证均值差异（单因素方差分析） #

对每幅影像的三个区域均值做方差分析，判断是否存在显著差异：

for idx, res in enumerate(region_means, 1):
    data = [res["红色区域均值"], res["蓝色区域均值"], res["外部区域均值"]]
    # 跳过含无效值的情况
    if any(np.isnan(d) for d in data):
        print(f"第{idx}幅影像存在无效数据，无法完成统计检验")
        continue
    # 执行单因素方差分析
    f_stat, p_val = f_oneway(*data)
    print(f"\n第{idx}幅影像统计结果：")
    print(f"F值={f_stat:.3f}，P值={p_val:.3f}")
    print("三个区域NDVI均值存在显著差异" if p_val < 0.05 else "三个区域NDVI均值无显著差异")

关键注意事项 #

• 必须确保矢量文件与NDVI影像的坐标系完全一致，否则掩膜提取会出错。
• 如果shapefile中红/蓝多边形是多个，用unary_union合并为单个多边形，避免掩膜时出现多个独立区域的问题。
• 若NDVI影像无明确nodata值，可手动指定nodata参数（比如设为-9999），避免将背景值计入均值计算。

内容的提问来源于stack exchange，提问作者Fabio Giacomelli