如何用开源工具在重叠多边形内计算栅格景观类别占比（免全栅格加载）

针对你提到的全球大栅格+带自相交问题的多边形分类占比统计需求，我推荐用Python的rasterio + geopandas组合来实现——这套方案不需要一次性加载整个超大栅格到内存，还能处理多边形的自相交问题，完全开源免费，我之前处理过类似的大尺度空间统计需求，亲测高效可靠。

核心思路 #

1. 先修复多边形自相交：自相交是很多空间分析工具的“拦路虎”，用geopandas的小技巧就能快速修复；
2. 分块读取栅格：利用rasterio的分块读取能力，只处理当前块内与多边形重叠的区域，避免内存溢出；
3. 逐年度栅格统计：遍历1992-2015年的每一幅年度栅格，对每个ABC_ID对应的多边形，统计其覆盖范围内各栅格类别的占比；
4. 整合结果生成表格：把所有年度的统计结果汇总，输出为CSV或Excel格式的标准化统计表格。

具体步骤与代码示例 #

1. 安装依赖库

先确保你安装了必要的工具库：

pip install rasterio geopandas pandas numpy shapely

2. 修复多边形的自相交问题

import geopandas as gpd

# 读取矢量文件
gdf = gpd.read_file("你的多边形文件.shp")

# 修复自相交：用buffer(0)是修复自相交多边形的经典技巧，能快速修正大部分拓扑问题
gdf["geometry"] = gdf["geometry"].apply(lambda geom: geom.buffer(0) if not geom.is_valid else geom)

# 过滤掉修复后仍无效的多边形（如果有的话）
gdf = gdf[gdf["geometry"].is_valid]

3. 分块处理栅格并统计占比

import rasterio
from rasterio.features import geometry_mask
import pandas as pd
import numpy as np

# 定义23幅年度栅格的路径列表（按年份排序）
raster_paths = ["1992.tif", "1993.tif", ..., "2015.tif"]

# 初始化结果容器
all_results = []

for raster_path in raster_paths:
    # 从文件名提取年份
    year = int(raster_path.split(".")[0])
    
    with rasterio.open(raster_path) as src:
        # 按栅格内部块大小分块读取，避免加载整个全球栅格
        for _, window in src.block_windows(1):
            # 读取当前块的栅格数据
            raster_data = src.read(1, window=window)
            # 获取当前块的地理范围
            window_bounds = src.window_bounds(window)
            
            # 筛选出与当前块重叠的多边形，提升处理效率
            overlapping_gdf = gdf[gdf.intersects(gpd.GeoDataFrame(
                {'geometry': [window_bounds]}, 
                crs=gdf.crs
            ).geometry.iloc[0])]
            
            if overlapping_gdf.empty:
                continue
            
            # 遍历每个重叠的多边形，统计类别占比
            for idx, row in overlapping_gdf.iterrows():
                abc_id = row["ABC_ID"]
                geom = row["geometry"]
                
                # 生成当前多边形在该栅格块内的掩码
                mask = geometry_mask(
                    [geom],
                    out_shape=raster_data.shape,
                    transform=src.window_transform(window),
                    invert=True
                )
                
                # 提取掩码覆盖的栅格值，排除NoData（如果有）
                masked_values = raster_data[mask & (raster_data != src.nodata)]
                if len(masked_values) == 0:
                    continue
                
                # 统计各类别的像素数量
                value_counts = np.bincount(masked_values, minlength=221)  # 类别范围10-220，确保覆盖所有可能值
                total_pixels = len(masked_values)
                
                # 计算占比并整理结果
                for category in range(10, 221):
                    count = value_counts[category]
                    if count == 0:
                        continue
                    proportion = count / total_pixels
                    all_results.append({
                        "ABC_ID": abc_id,
                        "Year": year,
                        "Category": category,
                        "Pixel_Count": count,
                        "Proportion": proportion
                    })

# 将结果转换为DataFrame并保存为CSV
result_df = pd.DataFrame(all_results)
result_df.to_csv("栅格类别占比统计.csv", index=False)

替代方案：GDAL命令行工具 #

如果你不想写代码，也可以用GDAL的命令行工具组合实现：

1. 修复多边形自相交：用ogr2ogr命令修正拓扑：

   ogr2ogr -f "ESRI Shapefile" 修复后的多边形.shp 原多边形.shp -dialect sqlite -sql "SELECT *, ST_Buffer(geometry, 0) AS geometry FROM 原图层名"
   

2. 对每个栅格，用gdal_rasterize将多边形栅格化为掩码，再结合gdalinfo提取统计信息，最后用脚本（bash/Python）汇总占比。不过这个方法灵活性不如Python代码，适合简单场景。

注意事项 #

• 全球栅格分块处理是核心，rasterio的block_windows会自动适配栅格的内部块大小，最大化内存效率；
• 若多边形数量极多，可以手动创建空间索引进一步提速：gdf.sindex；
• 如果栅格存在NoData值，一定要在提取值时过滤，避免影响占比计算。

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