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 python3.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_okTrue) for url in download_urls: filename url.split(/)[-1] with requests.get(url, streamTrue) as r: r.raise_for_status() with open(fGlobal-AI_PET_v3/{filename}, wb) as f: for chunk in r.iter_content(chunk_size8192): f.write(chunk) print(fDownloaded {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__], cropTrue, all_touchedTrue ) 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_colNAME): 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, ascendingFalse).head(10))典型输出示例省份平均AI台湾1.42海南1.38广东1.35......宁夏0.315. 可视化与结果验证5.1 空间分布可视化使用matplotlib绘制区域AI分布import matplotlib.pyplot as plt fig, ax plt.subplots(figsize(12, 8)) china_provinces.merge( province_ai_df, onNAME ).plot( columnmean_AI, legendTrue, axax, cmapRdYlGn, schemequantiles ) ax.set_title(中国省级行政区年平均干旱指数分布) plt.tight_layout() plt.savefig(china_ai_distribution.png, dpi300)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, workers4): with ThreadPoolExecutor(max_workersworkers) 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, maskedTrue) # 此时数据仍在磁盘上 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, chunksoptimal_chunk )