把1km栅格数据按县级边界拆成一个个县的独立数据，这叫什么操作？GIS软件和Python各怎么实现？

### **1. 问题解构** 用户的问题核心可以拆解为： 1. **操作名称**：涉及一种将栅格数据按照矢量面（县级行政区）边界进行裁切、分割、属性提取或重采样的空间操作。 2. **GIS桌面软件操作**：使用如ArcGIS、QGIS等图形界面软件完成此任务的工作流。 3. **Python编程操作**：使用Python及相关地理空间库（主要是GDAL和Rasterio）实现自动批处理的方法。结合问题描述（“处理成每个县每个县的数据”），该操作通常被称为 **“按掩膜提取”** 或 **“分区统计”**。两者的目标略有不同： * **按掩膜提取 (Extract by Mask)**：以每个县的矢量边界为掩膜，从原始的1km栅格中裁剪出完全落在该县范围内的栅格数据块。结果是得到一系列独立的、仅包含特定县域内像元的栅格文件。 * **分区统计 (Zonal Statistics)**：以每个县的矢量边界为分区，计算落在该分区内的所有栅格像元的统计值（如平均值、总和、最大值等）。结果是得到一个属性表，其中每个县对应一条记录，包含一个或多个计算出的统计字段，但不会生成新的栅格文件[ref_1]。在实际项目中，这两种操作常常结合使用[ref_3]。 ### **2. 方案推演** 要完成“每个县每个县的数据”处理，一个典型的流程如下： 1. **数据准备**：确保栅格数据和县级SHP文件使用相同的坐标参考系（Coordinate Reference System, CRS）。如果不同，需要先进行投影转换，使其空间范围对齐[ref_4]。 2. **核心处理**：根据最终需求选择上述两种主要方法之一或其组合。 * **生成分县栅格文件**：执行`按掩膜提取`。 * **获取分县统计值**：执行`分区统计`。 * **两者都要**：先提取分县栅格，再对每个分县栅格进行统计计算。 3. **自动化实现**：在Python中，通过遍历县级SHP文件中的每个要素（即每个县），循环执行上述核心操作，实现批处理。 ### **3. 具体操作方法** #### **3.1 GIS桌面软件操作（以QGIS为例）** 在QGIS中，可以使用图形化工具或处理模型实现。 | 操作目标 | 核心工具/插件 | 简要步骤 | | :--- | :--- | :--- | | **按掩膜提取** | **处理工具箱 → GDAL → 栅格提取 → 按掩膜图层裁剪栅格** | 1. 加载栅格和县界SHP。 2. 打开上述工具。 3. 输入栅格选择你的1km栅格。 4. 掩膜图层选择县界SHP。 5. 勾选“为每个要素创建单独的裁剪栅格”。 6. 设置输出目录和文件名模板（如 `result_{县名字段}.tif`）。 7. 运行。 | | **分区统计** | **处理工具箱 → Zonal Statistics** 或 **Zonal Statistics Plugin** | 1. 加载栅格和县界SHP。 2. 打开分区统计工具。 3. 输入矢量图层选择县界SHP。 4. 输入栅格图层选择你的1km栅格。 5. 选择要计算的统计量（如mean, sum）。 6. 输出列的前缀自定义。 7. 运行后，县界SHP的属性表会增加新的统计字段。 | > **注意**：在ArcGIS中，对应的工具分别是 **“按掩膜提取 (Extract by Mask)”**（在Spatial Analyst工具箱）和 **“分区统计以表格显示 (Zonal Statistics as Table)”**（在Spatial Analyst工具箱）[ref_6]。 #### **3.2 Python编程操作** 使用Python的`rasterio`和`geopandas`库可以更灵活地实现。以下代码演示了两种核心操作的自动化流程。 ```python import os import geopandas as gpd import rasterio import rasterio.mask import numpy as np # 定义文件路径 raster_path = "path/to/your/1km_raster.tif" county_shp_path = "path/to/your/county_boundaries.shp" output_dir = "path/to/output/folder" # 1. 读取数据 with rasterio.open(raster_path) as src: raster_data = src.read(1) # 读取第一个波段，根据实际情况调整 raster_profile = src.profile # 获取栅格的元数据（仿射变换、CRS等） raster_crs = src.crs counties_gdf = gpd.read_file(county_shp_path) # 确保矢量数据与栅格数据CRS一致 if counties_gdf.crs != raster_crs: print("正在进行坐标转换...") counties_gdf = counties_gdf.to_crs(raster_crs) # 假设矢量数据中有一个字段叫'NAME'用于标识县名 name_field = 'NAME' # 2. 方法A：循环执行“按掩膜提取”，为每个县生成独立的栅格文件 print("开始按县裁剪栅格...") for index, county in counties_gdf.iterrows(): county_name = county[name_field] # 获取当前县的几何形状 geom = county.geometry # 使用rasterio.mask进行裁剪 try: out_image, out_transform = rasterio.mask.mask(src, [geom], crop=True, filled=False) except ValueError: # 如果该县与栅格区域无交集，则跳过 print(f"县 {county_name} 与栅格数据无交集，已跳过。") continue # 更新输出栅格的元数据 out_profile = raster_profile.copy() out_profile.update({ "driver": "GTiff", "height": out_image.shape[1], "width": out_image.shape[2], "transform": out_transform }) # 构建输出文件路径 output_raster_path = os.path.join(output_dir, f"{county_name}.tif") # 写入新的栅格文件 with rasterio.open(output_raster_path, "w", **out_profile) as dst: dst.write(out_image) print(f"已生成: {output_raster_path}") # 3. 方法B：计算“分区统计”（以计算县域内栅格平均值和总和为例） print("\n开始计算分区统计...") stats_results = [] for index, county in counties_gdf.iterrows(): county_name = county[name_field] geom = county.geometry try: # 裁剪出该县范围内的栅格数据 out_image, out_transform = rasterio.mask.mask(src, [geom], crop=True, filled=False) # 剔除无效值（例如nodata） valid_pixels = out_image[0][~np.isnan(out_image[0])] # 假设第一个波段 if len(valid_pixels) > 0: mean_val = np.mean(valid_pixels) sum_val = np.sum(valid_pixels) else: mean_val = sum_val = None except (ValueError, rasterio.errors.RasterioIOError): mean_val = sum_val = None stats_results.append({ name_field: county_name, 'raster_mean': mean_val, 'raster_sum': sum_val }) # 将统计结果转换为DataFrame并与原始矢量属性表合并 import pandas as pd stats_df = pd.DataFrame(stats_results) counties_with_stats_gdf = counties_gdf.merge(stats_df, on=name_field) # 保存带有统计结果的矢量文件 output_stats_shp_path = os.path.join(output_dir, "counties_with_statistics.shp") counties_with_stats_gdf.to_file(output_stats_shp_path) print(f"分区统计结果已保存至: {output_stats_shp_path}") ``` **代码关键点说明**： 1. **坐标系统一**：代码首先检查并统一栅格与矢量的CRS，这是所有空间操作准确的前提[ref_4]。 2. **按掩膜提取**：核心是`rasterio.mask.mask`函数，它利用矢量几何对栅格进行裁剪。`crop=True`参数确保输出栅格的范围紧贴目标几何边界。 3. **分区统计**：在裁剪的基础上，使用`numpy`对裁剪出的像元数组进行统计计算（如`np.mean`, `np.sum`）。这里处理了无数据（NaN）的情况。 4. **循环与批处理**：通过`for`循环遍历`GeoDataFrame`中的每个县，实现了全自动的批处理。 5. **结果输出**：方法A输出一系列TIFF栅格文件；方法B输出一个增加了统计字段的新SHP文件。用户可以根据最终的数据需求（是需要分县的栅格图，还是只需要汇总的统计表格），选择使用方法A、方法B，或依次运行两者。对于处理全国数千个县的数据，Python脚本在效率和可重复性上远胜于手动图形界面操作[ref_1]。

创作声明：本文部分内容由AI辅助生成（AIGC），仅供参考

下一篇 C语言的void func(void)和Python的def func()看起来都表示无参函数，但为什么实际用法差别这么大？