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 python3.9 conda activate hydro conda install -c conda-forge rasterio xarray numpy pandas matplotlib cartopy核心库功能说明库名称用途关键特性rasterioGeoTIFF读取高效处理大型栅格文件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, streamTrue) with open(name, wb) as f: for chunk in tqdm(response.iter_content(chunk_size8192)): 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.4028235e38 坐标参考系统(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], dimxr.DataArray(range(1,13), dims[month], namemonth)) 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], dimxr.DataArray(range(1,13), dims[month], namemonth)) ds_et0.name potential_evapotranspiration # 合并数据集 ds xr.merge([ds_ai, ds_et0])3. 数据分析与可视化3.1 基本统计计算计算全球年平均干旱指数和潜在蒸散量# 年度平均计算 annual_ai ds_ai.mean(dimmonth) annual_et0 ds_et0.mean(dimmonth) 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}°-{i10}° for i in lat_bins[:-1]] ds[lat_bin] xr.DataArray( pd.cut(ds.y.values, binslat_bins, labelslat_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, cmapviridis): fig plt.figure(figsize(12, 6)) ax fig.add_subplot(1, 1, 1, projectionccrs.Robinson()) # 添加地理要素 ax.add_feature(cfeature.LAND, facecolorlightgray) ax.add_feature(cfeature.COASTLINE, linewidth0.5) ax.add_feature(cfeature.BORDERS, linestyle:, linewidth0.5) # 绘制数据 im ax.imshow(data, transformccrs.PlateCarree(), cmapcmap, originupper) # 添加色标 plt.colorbar(im, axax, orientationhorizontal, pad0.05, aspect50, labeltitle) plt.title(title) plt.tight_layout() return fig # 绘制七月干旱指数 plot_global_map(ds_ai.sel(month7)[0], 七月全球干旱指数, cmapRdYlBu_r) plt.savefig(global_ai_july.png, dpi300, bbox_inchestight)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, dropTrue) # 计算各区域月均值 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, labelname) axes[1].plot(data.month, data.potential_evapotranspiration, labelname) 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, daskparallelized, output_dtypes[np.int8] ) # 统计各类面积占比 class_area np.histogram(annual_ai_class, bins5, 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(x10, y10).mean() # 计算相关系数 corr xr.corr(ai_resampled, ndvi_data, dim[x, y]) # 可视化 fig, ax plt.subplots(figsize(10, 6)) corr.plot(axax, cmapcoolwarm, vmin-1, vmax1) 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, encodingencoding) # 保存处理后的月度数据 ds_monthly ds.stack(sample(month, y, x)) ds_monthly.to_netcdf(output/global_ai_pet_monthly.nc, encodingencoding)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}, enginerasterio, combinenested, concat_dimmonth) # 定义分块计算任务 def annual_mean(chunk): return chunk.mean(dimmonth) # 并行计算 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, dtypedtype, modew, shapeshape) fp.flush() return np.memmap(filename, dtypedtype, moder, shapeshape) # 示例逐月处理并累积结果 annual_sum create_memmap(temp/annual_sum.dat, (21600, 43200), np.float32) for month in range(1, 13): with rasterio.open(fdata/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]}, crsEPSG:4326) # 使用GDAL进行裁剪 cmd fgdalwarp -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格式以便团队成员共享使用。
