Global-AI_PET_v3 数据集实战：Python 加载 5.67GB 栅格数据并计算区域平均 AI

Global-AI_PET_v3 数据集实战：Python 加载 5.67GB 栅格数据并计算区域平均 AI

当面对全球尺度的气候数据时，如何高效处理海量栅格文件成为许多研究者的首要挑战。Global-AI_PET_v3 作为当前最权威的干旱指数数据集之一，其5.67GB的体积和30弧秒的高分辨率既带来了丰富的信息量，也对数据处理技术提出了更高要求。本文将手把手带你用Python实现从数据下载到空间分析的全流程，特别针对计算行政区划内平均干旱指数这一典型需求，提供可复用的代码方案。

1. 环境准备与数据获取

1.1 安装必要的Python库

处理地理空间数据需要特定的工具链。推荐使用conda创建独立环境，避免依赖冲突：

conda create -n geo_env python=3.9 conda activate geo_env conda install -c conda-forge rasterio xarray geopandas dask netCDF4 matplotlib

关键库的作用：

• rasterio：高效读取GeoTIFF格式的栅格数据
• xarray：处理NetCDF格式的多维数据集
• geopandas：行政区划矢量数据操作
• dask：大数据分块处理

1.2 下载数据集

Global-AI_PET_v3 数据集存储在Figshare平台，可通过Python自动下载：

import requests import os download_urls = [ "https://figshare.com/ndownloader/files/.../Global-AI_monthly_v3.zip", "https://figshare.com/ndownloader/files/.../Global-AI_ET0_annual_v3.zip" ] os.makedirs("Global-AI_PET_v3", exist_ok=True) for url in download_urls: filename = url.split('/')[-1] with requests.get(url, stream=True) as r: r.raise_for_status() with open(f"Global-AI_PET_v3/{filename}", 'wb') as f: for chunk in r.iter_content(chunk_size=8192): f.write(chunk) print(f"Downloaded {filename}")

注意：实际URL需替换为Figshare提供的真实下载链接。完整数据集包含月度/年度AI数据、ET0数据及说明文档。

2. 数据加载与内存优化

2.1 读取NetCDF格式数据

Global-AI_PET_v3 的年平均数据采用NetCDF格式，xarray是处理此类数据的理想选择：

import xarray as xr # 加载年度干旱指数数据 ds = xr.open_dataset("Global-AI_PET_v3/Global-AI_annual_v3.nc") print(ds)

典型输出显示数据维度信息：

<xarray.Dataset> Dimensions: (lat: 21600, lon: 43200) Coordinates: * lat (lat) float64 89.99 89.98 89.97 ... -89.98 -89.99 * lon (lon) float64 -180.0 -179.99 -179.98 ... 179.98 179.99 Data variables: AI (lat, lon) float32 ...

2.2 分块处理大文件

对于5.67GB的大型数据集，直接加载可能导致内存溢出。Dask提供了优雅的解决方案：

import dask.array as da # 创建分块读取的xarray数据集 ds_chunked = xr.open_dataset( "Global-AI_PET_v3/Global-AI_annual_v3.nc", chunks={'lat': 1000, 'lon': 1000} ) # 查看分块情况 print(ds_chunked.AI.data)

此时数据并未实际加载，只有在计算时才会按需读取相应块，显著降低内存压力。

3. 空间分析与区域统计

3.1 准备行政区划边界

以中国省级行政区为例，需要获取对应的矢量边界数据：

import geopandas as gpd # 加载省级行政区划（示例使用Natural Earth数据） china_provinces = gpd.read_file("path_to_shapefile/china_provinces.shp") # 选择特定省份（如云南省） yunnan = china_provinces[china_provinces.NAME == "Yunnan"]

3.2 区域平均AI计算

结合rasterio和geopandas实现精确的空间统计：

from rasterio.mask import mask import numpy as np def calculate_mean_ai(raster_path, geometry): with rasterio.open(raster_path) as src: # 裁剪到目标区域 clipped, transform = mask( src, [geometry.__geo_interface__], crop=True, all_touched=True ) clipped = clipped.squeeze() # 计算有效像素平均值 valid_pixels = clipped[clipped != src.nodata] return np.mean(valid_pixels) # 计算云南省平均AI yunnan_ai = calculate_mean_ai( "Global-AI_PET_v3/Global-AI_annual_v3.tif", yunnan.geometry.iloc[0] ) print(f"云南省年平均干旱指数: {yunnan_ai:.2f}")

4. 完整工作流示例

将上述步骤整合为可复用的工作流：

import pandas as pd def batch_calculate_ai(raster_path, gdf, id_col="NAME"): results = [] for _, row in gdf.iterrows(): mean_ai = calculate_mean_ai(raster_path, row.geometry) results.append({ id_col: row[id_col], "mean_AI": mean_ai }) return pd.DataFrame(results) # 计算全国各省平均AI province_ai_df = batch_calculate_ai( "Global-AI_PET_v3/Global-AI_annual_v3.tif", china_provinces ) # 结果排序输出 print(province_ai_df.sort_values("mean_AI", ascending=False).head(10))

典型输出示例：

5. 可视化与结果验证

5.1 空间分布可视化

使用matplotlib绘制区域AI分布：

import matplotlib.pyplot as plt fig, ax = plt.subplots(figsize=(12, 8)) china_provinces.merge( province_ai_df, on="NAME" ).plot( column="mean_AI", legend=True, ax=ax, cmap="RdYlGn", scheme="quantiles" ) ax.set_title("中国省级行政区年平均干旱指数分布") plt.tight_layout() plt.savefig("china_ai_distribution.png", dpi=300)

5.2 数据质量检查

验证计算结果合理性：

# 检查全球AI值范围 global_ai = xr.open_dataarray("Global-AI_PET_v3/Global-AI_annual_v3.nc") print(f"全球AI最小值: {global_ai.min().values:.2f}") print(f"全球AI最大值: {global_ai.max().values:.2f}") print(f"全球AI平均值: {global_ai.mean().values:.2f}") # 对比文献值 literature_value = 0.65 # 参考值 deviation = abs(province_ai_df["mean_AI"].mean() - literature_value) print(f"与文献值偏差: {deviation:.2%}")

6. 性能优化技巧

6.1 并行计算加速

利用多核CPU加速区域统计：

from concurrent.futures import ThreadPoolExecutor def parallel_calculate_ai(raster_path, gdf, workers=4): with ThreadPoolExecutor(max_workers=workers) as executor: futures = [] for _, row in gdf.iterrows(): futures.append( executor.submit( calculate_mean_ai, raster_path, row.geometry ) ) return [f.result() for f in futures] # 测试并行版本 %timeit -n 3 parallel_calculate_ai(raster_path, china_provinces)

6.2 内存映射技术

对于超大规模数据处理，使用内存映射避免内存爆炸：

# 使用rasterio的内存映射模式 with rasterio.open("Global-AI_PET_v3/Global-AI_annual_v3.tif") as src: # 创建内存映射数组 ai_array = src.read(1, masked=True) # 此时数据仍在磁盘上 mean_global = ai_array.mean() print(f"内存映射计算的平均值: {mean_global:.2f}")

7. 常见问题解决方案

7.1 坐标系不一致问题

当栅格数据与矢量数据的坐标系不匹配时：

# 统一坐标系 china_provinces = china_provinces.to_crs("EPSG:4326") # WGS84 # 或者动态转换 with rasterio.open(raster_path) as src: target_crs = src.crs gdf_reprojected = gdf.to_crs(target_crs)

7.2 缺失值处理

针对数据中的nodata值：

def safe_mean_ai(clipped_array, nodata): masked = np.ma.masked_equal(clipped_array, nodata) if masked.count() == 0: return np.nan return masked.mean() # 在calculate_mean_ai函数中使用 valid_mean = safe_mean_ai(clipped, src.nodata)

7.3 大文件分块处理策略

优化分块大小以平衡I/O和计算效率：

# 根据内存大小调整分块 optimal_chunk = { 'lat': 2000, 'lon': 2000, 'band': 1 } ds_optimized = xr.open_dataset( "Global-AI_PET_v3/Global-AI_monthly_v3.nc", chunks=optimal_chunk )