使用 rasterstats 库进行栅格与矢量数据的空间分析

在地理信息系统（GIS）领域，栅格数据和矢量数据是两类常见的数据类型。栅格数据通常代表像素网格，如遥感影像或土地利用图，而矢量数据则通常表示具体的地理实体，如行政区划或土地边界。如何有效地结合这两类数据进行空间分析是许多地理研究中的关键问题。

rasterstats 是一个用于处理栅格和矢量数据的 Python 库，提供了便捷的工具来实现栅格统计、空间叠加分析等。本文将介绍如何使用 rasterstats 库来进行栅格与矢量数据的空间分析，特别是利用 zonal statistics 功能进行不同区域的土地利用统计。

1. rasterstats 库简介

rasterstats 是一个轻量级的 Python 库，主要用于从栅格数据中提取基于矢量区划的统计信息。它可以进行区域统计（zonal statistics）、分类统计、面状统计等功能，简化了栅格与矢量数据的交互。

主要功能包括：

• zonal_stats：针对矢量区域的统计功能，从栅格数据中提取统计信息，如均值、最小值、最大值、方差等。
• point_query：用于从栅格中提取给定点的数值信息。
• gen_zonal_stats：生成基于多边形的统计结果，可用于大规模的并行处理。

2. 环境设置与安装

在开始之前，确保你的环境中已安装了 rasterstats。可以通过以下命令安装：

pip install rasterstats

同时，还需要 geopandas、rasterio 等 GIS 相关库的支持，建议安装以下依赖库：

pip install geopandas rasterio

3. 基础示例：土地利用的区域统计

假设我们有一个代表城市的土地利用栅格数据（如土地利用分类图），以及一个表示行政区划的矢量数据。我们想要统计每个行政区内不同土地利用类型的面积。

3.1 加载数据

首先，我们需要加载矢量数据（行政区划）和栅格数据（土地利用图）：

import geopandas as gpd
from rasterstats import zonal_stats

# 加载行政区划矢量数据
admin_zones = gpd.read_file("path/to/行政区划.shp")

# 土地利用栅格文件路径
landuse_raster = "path/to/Landuse_2020_mainCity1.tif"

3.2 执行 zonal_stats 统计

接下来，使用 zonal_stats 函数对行政区划进行分类统计。我们会统计每个行政区内不同土地利用类型的像素数量，并将其转换为面积：

from rasterstats import zonal_stats

# 使用分类统计 (categorical=True)
stats = zonal_stats(admin_zones, landuse_raster, categorical=True, nodata=0)

3.3 计算面积

土地利用图中的每个像素代表某个具体的面积，比如 10 平方米。我们可以基于像素分辨率计算每种土地利用类型的总面积：

import rasterio

# 打开栅格以获取分辨率
with rasterio.open(landuse_raster) as src:
    pixel_area = src.res[0] * src.res[1]  # 每个像素的面积，例如平方米

# 计算每种土地利用类型的面积
for zone_index, zone_stats in enumerate(stats):
    print(f"行政区划 {admin_zones.iloc[zone_index]['name1']} 的土地利用统计：")
    for landuse_type, pixel_count in zone_stats.items():
        area = pixel_count * pixel_area  # 计算面积
        print(f"土地利用类型 {landuse_type}: 面积 = {area:.2f} 平方米")

3.4 输出结果到 Shapefile

我们还可以将这些统计结果写入到新的 shapefile 中，以便后续的可视化和分析：

# 创建新的 GeoDataFrame 列来存储统计结果
for landuse_type in landuse_types:
    admin_zones[f'landuse_{landuse_type}_area'] = [zone_stats.get(landuse_type, 0) * pixel_area for zone_stats in stats]

# 保存到新的 Shapefile
admin_zones.to_file("output/admin_zones_with_landuse_stats.shp", encoding='utf-8')

4. 高级功能

rasterstats 还提供了其他功能，例如可以进行更多的高级统计分析：

• 统计均值、最大值、最小值等：通过设置 stats=['mean', 'max', 'min'] 等参数，可以提取更为丰富的统计信息。
• 面状统计：可以对多边形进行分析，计算面积、周长等属性。

# 使用更多的统计指标
stats = zonal_stats(admin_zones, landuse_raster, stats="mean max min")
print(stats)

5. 性能与优化

在处理大型数据集时，可能需要优化计算性能。rasterstats 支持并行处理，可以使用 multiprocessing 参数来加速区域统计。

# 使用 multiprocessing 提高性能
stats = zonal_stats(admin_zones, landuse_raster, categorical=True, nodata=0, multiprocessing=True)

6. 总结

rasterstats 是一个功能强大且易于使用的 Python 库，能够简化栅格与矢量数据的结合与分析。无论是进行土地利用统计，还是其他类型的区域统计，rasterstats 都能帮助我们轻松实现这些任务。在实际应用中，我们可以将统计结果与 GIS 数据结合，以进行更深入的空间分析和决策支持。

通过本文的介绍，希望你对如何使用 rasterstats 进行栅格数据的区域统计有了清晰的了解。可以根据自己的需求，结合更多的功能，将这一工具应用于实际项目中。