Global-AI_PET_v3 数据集实战：1km水文气候数据的Python全流程解析
当我们需要在全球尺度上研究水资源分布、农业灌溉需求或生态系统脆弱性时，高精度的潜在蒸散量(PET)和干旱指数(AI)数据变得至关重要。Global-AI_PET_v3作为目前最新版本的数据集，以1km分辨率提供了1970-2000年间的月度和年度平均值，为相关领域的研究者带来了前所未有的分析可能性。本文将带你从数据获取到初步分析，完整走通这个强大数据集的Python处理流程。
1. 环境准备与数据获取
1.1 安装必要的Python库
处理地理空间数据需要特定的工具链。推荐使用conda创建独立环境以避免依赖冲突：
conda create -n hydro python=3.9 conda activate hydro conda install -c conda-forge rasterio xarray numpy pandas matplotlib cartopy核心库功能说明：
| 库名称 | 用途 | 关键特性 |
|---|---|---|
| rasterio | GeoTIFF读取 | 高效处理大型栅格文件 |
| xarray | 多维数据分析 | 标签化维度，方便时间序列处理 |
| cartopy | 地理可视化 | 专业级地图投影支持 |
1.2 从Figshare下载数据集
数据集包含四个主要文件，总计约5.67GB。我们可以编写自动化下载脚本：
import requests from tqdm import tqdm files = [ "Global-AI_monthly_v3.zip", "Global-ET0_monthly_v3.zip", "Global-AI_ET0_annual_v3.zip", "Global-AI_PET_v3_Readme.pdf" ] base_url = "https://figshare.com/ndownloader/files/" file_ids = ["34271745", "34271748", "34271751", "34271754"] # 实际ID需确认 for name, fid in zip(files, file_ids): url = f"{base_url}{fid}" print(f"下载 {name}...") response = requests.get(url, stream=True) with open(name, "wb") as f: for chunk in tqdm(response.iter_content(chunk_size=8192)): if chunk: f.write(chunk)注意：大规模下载前请确认Figshare的使用条款，商业用途可能需要额外授权。
2. 数据读取与初步探索
2.1 解压与文件结构
下载完成后，使用Python标准库解压文件：
import zipfile import os for zip_file in [f for f in os.listdir() if f.endswith(".zip")]: with zipfile.ZipFile(zip_file, 'r') as zip_ref: zip_ref.extractall("data")解压后的典型目录结构：
data/ ├── Global-AI_monthly_v3/ │ ├── ai_01.tif │ ├── ai_02.tif │ └── ... (12个月度文件) ├── Global-ET0_monthly_v3/ │ ├── et0_01.tif │ ├── et0_02.tif │ └── ... (12个月度文件) └── Global-AI_ET0_annual_v3/ ├── ai_year.tif └── et0_year.tif2.2 使用rasterio读取单个月份数据
让我们先加载一月份的干旱指数数据：
import rasterio import matplotlib.pyplot as plt with rasterio.open("data/Global-AI_monthly_v3/ai_01.tif") as src: ai_data = src.read(1) profile = src.profile print(f"数据尺寸: {ai_data.shape}") print(f"无效值: {src.nodata}") print(f"坐标参考系统(CRS): {src.crs}") print(f"地理变换矩阵:\n{src.transform}")典型输出示例：
数据尺寸: (21600, 43200) 无效值: -3.4028235e+38 坐标参考系统(CRS): EPSG:4326 地理变换矩阵: | 0.00, 0.00,-180.00| | 0.00,-0.00, 90.00| | 0.00, 0.00, 0.00|2.3 使用xarray处理时间序列
对于多时间步长的分析，xarray提供了更便捷的接口：
import xarray as xr import glob # 创建月度AI数据集 ai_files = sorted(glob.glob("data/Global-AI_monthly_v3/ai_*.tif")) ds_ai = xr.concat([xr.open_rasterio(f) for f in ai_files], dim=xr.DataArray(range(1,13), dims=["month"], name="month")) ds_ai.name = "aridity_index" # 创建月度ET0数据集 et0_files = sorted(glob.glob("data/Global-ET0_monthly_v3/et0_*.tif")) ds_et0 = xr.concat([xr.open_rasterio(f) for f in et0_files], dim=xr.DataArray(range(1,13), dims=["month"], name="month")) ds_et0.name = "potential_evapotranspiration" # 合并数据集 ds = xr.merge([ds_ai, ds_et0])3. 数据分析与可视化
3.1 基本统计计算
计算全球年平均干旱指数和潜在蒸散量：
# 年度平均计算 annual_ai = ds_ai.mean(dim="month") annual_et0 = ds_et0.mean(dim="month") print(f"全球年平均干旱指数统计:") print(f" 均值: {annual_ai.mean().values:.2f}") print(f" 标准差: {annual_ai.std().values:.2f}") print(f" 最小值: {annual_ai.min().values:.2f}") print(f" 最大值: {annual_ai.max().values:.2f}") # 按纬度带分析 lat_bins = range(-90, 91, 10) lat_labels = [f"{i}°-{i+10}°" for i in lat_bins[:-1]] ds['lat_bin'] = xr.DataArray( pd.cut(ds.y.values, bins=lat_bins, labels=lat_labels), dims=['y'], coords={'y': ds.y} ) zonal_stats = ds.groupby('lat_bin').mean().compute()3.2 地理空间可视化
使用cartopy创建专业地图：
import cartopy.crs as ccrs import cartopy.feature as cfeature def plot_global_map(data, title, cmap='viridis'): fig = plt.figure(figsize=(12, 6)) ax = fig.add_subplot(1, 1, 1, projection=ccrs.Robinson()) # 添加地理要素 ax.add_feature(cfeature.LAND, facecolor='lightgray') ax.add_feature(cfeature.COASTLINE, linewidth=0.5) ax.add_feature(cfeature.BORDERS, linestyle=':', linewidth=0.5) # 绘制数据 im = ax.imshow(data, transform=ccrs.PlateCarree(), cmap=cmap, origin='upper') # 添加色标 plt.colorbar(im, ax=ax, orientation='horizontal', pad=0.05, aspect=50, label=title) plt.title(title) plt.tight_layout() return fig # 绘制七月干旱指数 plot_global_map(ds_ai.sel(month=7)[0], "七月全球干旱指数", cmap='RdYlBu_r') plt.savefig('global_ai_july.png', dpi=300, bbox_inches='tight')3.3 时间序列分析
提取特定区域的时间变化模式：
# 定义感兴趣区域（以经纬度边界表示） regions = { "Sahara": {"lon": (-20, 40), "lat": (15, 30)}, "Amazon": {"lon": (-80, -50), "lat": (-20, 10)}, "East_Asia": {"lon": (100, 140), "lat": (20, 50)} } # 区域提取函数 def extract_region(ds, region): mask_lon = (ds.x >= region["lon"][0]) & (ds.x <= region["lon"][1]) mask_lat = (ds.y >= region["lat"][0]) & (ds.y <= region["lat"][1]) return ds.where(mask_lon & mask_lat, drop=True) # 计算各区域月均值 regional_data = {} for name, bounds in regions.items(): regional_data[name] = extract_region(ds, bounds).mean(dim=['x', 'y']).compute() # 绘制季节变化曲线 fig, axes = plt.subplots(2, 1, figsize=(10, 8)) for name, data in regional_data.items(): axes[0].plot(data.month, data.aridity_index, label=name) axes[1].plot(data.month, data.potential_evapotranspiration, label=name) axes[0].set_title("干旱指数季节变化") axes[0].set_ylabel("AI") axes[1].set_title("潜在蒸散量季节变化") axes[1].set_ylabel("ET0 (mm)") for ax in axes: ax.legend() ax.grid(True) ax.set_xlabel("月份") plt.tight_layout()4. 高级分析与应用示例
4.1 干旱热点区域识别
基于干旱指数的分类标准：
# 干旱程度分类 def classify_aridity(ai): conditions = [ (ai < 0.03), # 极端干旱 (ai < 0.2), # 干旱 (ai < 0.5), # 半干旱 (ai < 0.65), # 半湿润 (ai >= 0.65) # 湿润 ] classes = [0, 1, 2, 3, 4] return np.select(conditions, classes) # 应用分类 annual_ai_class = xr.apply_ufunc( classify_aridity, annual_ai, dask='parallelized', output_dtypes=[np.int8] ) # 统计各类面积占比 class_area = np.histogram(annual_ai_class, bins=5, range=(0,4))[0] total_pixels = class_area.sum() class_percent = class_area / total_pixels * 100 print("全球干旱类型面积占比:") for i, p in enumerate(class_percent): print(f" 类别{i}: {p:.1f}%")4.2 与植被指数的相关性分析
假设我们有MODIS NDVI数据，可以探索干旱指数与植被生长的关系：
# 示例代码 - 需要实际NDVI数据 def calculate_correlation(ai_data, ndvi_data): # 重采样到相同分辨率 ai_resampled = ai_data.coarsen(x=10, y=10).mean() # 计算相关系数 corr = xr.corr(ai_resampled, ndvi_data, dim=['x', 'y']) # 可视化 fig, ax = plt.subplots(figsize=(10, 6)) corr.plot(ax=ax, cmap='coolwarm', vmin=-1, vmax=1) ax.set_title("干旱指数与NDVI的空间相关性") return fig4.3 数据导出与共享
将处理结果导出为NetCDF格式以便后续使用：
# 设置编码以压缩数据 encoding = { "aridity_index": {"zlib": True, "complevel": 5}, "potential_evapotranspiration": {"zlib": True, "complevel": 5} } # 保存年度数据集 ds_annual = xr.merge([ annual_ai.rename("annual_ai"), annual_et0.rename("annual_et0") ]) ds_annual.to_netcdf("output/global_ai_pet_annual.nc", encoding=encoding) # 保存处理后的月度数据 ds_monthly = ds.stack(sample=('month', 'y', 'x')) ds_monthly.to_netcdf("output/global_ai_pet_monthly.nc", encoding=encoding)5. 性能优化与大规模处理技巧
处理全球1km分辨率数据时，内存管理至关重要。以下是几个实用技巧：
5.1 分块处理策略
使用dask进行惰性计算和分块处理：
import dask.array as da # 使用dask打开数据集 ai_dask = xr.open_mfdataset("data/Global-AI_monthly_v3/ai_*.tif", chunks={'x': 1000, 'y': 1000}, engine='rasterio', combine='nested', concat_dim='month') # 定义分块计算任务 def annual_mean(chunk): return chunk.mean(dim='month') # 并行计算 result = ai_dask.map_blocks(annual_mean) annual_ai = result.compute() # 只在需要时执行计算5.2 内存映射技术
对于单机处理，可以使用内存映射减少内存占用：
# 创建内存映射文件 def create_memmap(filename, shape, dtype): if not os.path.exists(filename): fp = np.memmap(filename, dtype=dtype, mode='w+', shape=shape) fp.flush() return np.memmap(filename, dtype=dtype, mode='r+', shape=shape) # 示例：逐月处理并累积结果 annual_sum = create_memmap("temp/annual_sum.dat", (21600, 43200), np.float32) for month in range(1, 13): with rasterio.open(f"data/Global-AI_monthly_v3/ai_{month:02d}.tif") as src: data = src.read(1) valid_mask = (data != src.nodata) annual_sum[valid_mask] += data[valid_mask]5.3 区域提取最佳实践
当只需要研究特定区域时，先裁剪再处理能显著提升效率：
from shapely.geometry import box import geopandas as gpd def clip_to_region(input_file, output_file, bbox): # 创建地理边界 geometry = box(*bbox) gdf = gpd.GeoDataFrame({'geometry': [geometry]}, crs="EPSG:4326") # 使用GDAL进行裁剪 cmd = f"gdalwarp -te {bbox[0]} {bbox[1]} {bbox[2]} {bbox[3]} " \ f"-of GTiff {input_file} {output_file}" os.system(cmd) # 示例：裁剪东南亚区域 clip_to_region( "data/Global-AI_monthly_v3/ai_01.tif", "output/se_asia_ai.tif", (90, -15, 150, 30) # minx, miny, maxx, maxy )在实际项目中，根据数据使用场景的不同，这些技术可以组合使用。例如，先使用GDAL进行初步的区域裁剪，然后在处理阶段应用dask进行并行计算，最后将结果保存为压缩的NetCDF格式以便团队成员共享使用。
