如何在Python中对大尺寸栅格进行分区统计？

问题背景 #

需要用同一个shapefile对数百个大尺寸栅格（单栅格尺寸13665×30641）做分区统计，结果导出为CSV。测试小数据时代码正常，但处理实际数据时，无论是用zonal_stats还是设置block_size的gen_zonal_stats，都出现内存溢出错误。

原代码 #

##### This code now works for small but not large size problem

import rasterio
import geopandas as gpd
from rasterstats import zonal_stats
from rasterstats import gen_zonal_stats
import pandas as pd
import glob

# Open the shapefile
polygons = gpd.read_file('/test/LargeDomainBoundary.shp')

# List of rasters to process
# Get a list of all raster files in a directory
raster_list = glob.glob('/test1/*.tif')
print(raster_list)

# Creating an empty dataframe to store the results
results_df = pd.DataFrame()

# Loop through the rasters
for raster_path in raster_list:

    # Open the raster file
    with rasterio.open(raster_path) as src:
      raster = src.read(1)
      print(raster.shape)
      trans = src.transform  # --> 这里用src.transform代替src.affine
      print(trans)

    # Perform zonal statistics
    stats = zonal_stats(polygons, raster, stats=['mean', 'sum', 'min', 'count', 'max'], affine=trans, tolerance=0.1, nodata=-32768) #nodata=-9999
#     stats = list(
#         gen_zonal_stats(polygons, raster, affine=trans, stats=['mean', 'sum', 'min', 'count', 'max'], nodata=-32768,
#                         block_size=(1000, 1000)))


    # Creating a dataframe from the stats
    df = pd.DataFrame(stats)
    # Adding the name of the raster processed as a column
    df["Raster"] = raster_path
    # Concatenating the result to the results dataframe
    results_df = pd.concat([results_df, df])

# Saving the results to a CSV file
results_df.to_csv('/globalhome/test/zonal_stats_results.csv', index=False)

错误信息 #

Traceback (most recent call last):
File "C:\GDAL\zonalusingsameshapefile.py", line 32, in 
stats = list(
File "C:\Users\Anaconda3\envs\venvgdal\lib\site-packages\rasterstats\main.py", line 155, in gen_zonal_stats
fsrc = rast.read(bounds=geom_bounds)
File "C:\Users\Anaconda3\envs\venvgdal\lib\site-packages\rasterstats\io.py", line 298, in read
new_array = boundless_array(
File "C:\Users\Anaconda3\envs\venvgdal\lib\site-packages\rasterstats\io.py", line 179, in boundless_array
out = np.ones(shape=window_shape) * nodata
File "C:\Users\Anaconda3\envs\venvgdal\lib\site-packages\numpy\core\numeric.py", line 204, in ones
a = empty(shape, dtype, order)
ValueError: array is too big; `arr.size * arr.dtype.itemsize` is larger than the maximum possible size.

解决思路与修改后代码 #

核心问题是原代码一次性把整个大栅格读入内存，导致内存不足。修改方向是避免全量加载栅格，同时优化矢量和结果存储逻辑：

修改后代码 #

import rasterio
import geopandas as gpd
from rasterstats import gen_zonal_stats
import pandas as pd
import glob

# 加载矢量边界
polygons = gpd.read_file('/test/LargeDomainBoundary.shp')

# 获取所有待处理的栅格文件
raster_list = glob.glob('/test1/*.tif')

# 用列表收集结果，避免多次拼接DataFrame浪费内存
results_list = []

for raster_path in raster_list:
    print(f"正在处理: {raster_path}")
    # 打开栅格但不读取全部数据，直接传入数据源对象给统计函数
    with rasterio.open(raster_path) as src:
        # 用栅格原生的块大小作为分块参数，效率更高
        block_size = src.block_shapes[0]
        # 使用gen_zonal_stats，直接传入src而非已读取的raster数组
        stats = list(gen_zonal_stats(
            polygons,
            src,
            stats=['mean', 'sum', 'min', 'count', 'max'],
            nodata=-32768,
            block_size=block_size
        ))
    
    # 将当前栅格的统计结果转为DataFrame，添加栅格路径标识
    df = pd.DataFrame(stats)
    df["Raster"] = raster_path
    results_list.append(df)

# 一次性合并所有结果，减少内存碎片
results_df = pd.concat(results_list, ignore_index=True)

# 保存最终结果到CSV
results_df.to_csv('/globalhome/test/zonal_stats_results.csv', index=False)

关键修改点 #

1. 取消全量加载栅格：移除src.read(1)，直接将rasterio的src对象传给gen_zonal_stats，rasterstats会自动按需读取栅格块，不会一次性占用大量内存。
2. 匹配原生栅格块大小：用src.block_shapes[0]获取栅格存储时的原生块大小，比手动设置100×100更高效，减少IO次数。
3. 优化结果存储：用列表收集每个栅格的结果DataFrame，最后一次性拼接，避免多次concat导致的内存浪费和碎片。

额外优化（针对超大单个多边形） #

如果你的shapefile是单个超大多边形，即使分块读取栅格仍可能内存溢出，建议先把大面拆分成多个小网格面，分批次统计后再合并结果：

from shapely.geometry import box

# 获取原多边形的边界范围
xmin, ymin, xmax, ymax = polygons.total_bounds
# 拆分网格数量（可根据实际情况调整）
x_split = 10
y_split = 10
# 计算网格步长
x_step = (xmax - xmin) / x_split
y_step = (ymax - ymin) / y_split

# 生成网格多边形
grid_polygons = []
for i in range(x_split):
    for j in range(y_split):
        grid_polygons.append(box(
            xmin + i*x_step, 
            ymin + j*y_step, 
            xmin + (i+1)*x_step, 
            ymin + (j+1)*y_step
        ))

# 转换为GeoDataFrame并裁剪到原多边形范围
grid_gdf = gpd.GeoDataFrame(geometry=grid_polygons, crs=polygons.crs)
clipped_grid = gpd.overlay(grid_gdf, polygons, how='intersection')

之后用clipped_grid代替原polygons进行统计，最后对每个栅格的小网格统计结果做聚合（比如sum字段累加，mean字段按面积加权平均）。

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