1. ERA5-Land数据简介与获取方法ERA5-Land是欧洲中期天气预报中心ECMWF推出的高分辨率陆地表面再分析数据集空间分辨率达到9公里0.1°时间跨度从1950年延续至今。这个数据集特别适合需要精细分析地表气候特征的研究场景比如农业气象监测、流域水文模拟或区域气候变化研究。我处理过不少气象数据集ERA5-Land最让我印象深刻的是它的数据完整性。记得有次做青藏高原的气候分析其他数据集在高原边缘总有缺失但ERA5-Land的数据覆盖非常完整。数据集包含三个核心变量总降水量total precipitation月累计降水毫米数2米气温2m temperature近地表月平均气温开尔文单位地表向下太阳辐射Surface solar radiation downwards单位面积接收的短波辐射量获取数据最便捷的方式是通过Copernicus Climate Data StoreCDS的API。第一次使用时需要注册账号并在用户目录下创建.cdsapirc配置文件。这里分享一个我常用的Python下载脚本模板import cdsapi c cdsapi.Client() def download_era5land(year, month, variables, target_file): c.retrieve( reanalysis-era5-land-monthly-means, { format: netcdf, variable: variables, year: year, month: month, time: 00:00, product_type: monthly_averaged_reanalysis }, target_file)实际使用时要注意CDS的配额限制。我建议按年度分批下载比如先获取2010-2020年的目标区域数据。如果研究区域不大可以先用area参数限定下载范围北/西/南/东的经纬度这能显著减小文件体积。2. 数据预处理实战技巧2.1 格式转换与空间参考处理原始下载的NetCDF格式虽然包含丰富的元数据但很多GIS软件对它的支持有限。我习惯先用GDAL将其转为GeoTIFF这个过程中会遇到三个典型问题维度顺序混乱有时变量数据的维度排列是(time, lat, lon)有时却是(time, lon, lat)无效值处理海洋区域的填充值可能干扰后续计算坐标精度损失直接转换可能导致小数点后精度降低这是我优化过的转换代码加入了维度自动检测和异常值过滤from osgeo import gdal import numpy as np def safe_nc_to_tif(nc_path, var_name, output_dir): ds gdal.Open(fNETCDF:{nc_path}:{var_name}) band ds.GetRasterBand(1) arr band.ReadAsArray() # 处理无效值 arr[arr 1e10] np.nan # 自动检测维度顺序 if arr.shape[0] ds.RasterYSize: arr arr.T # 创建输出文件 driver gdal.GetDriverByName(GTiff) out_ds driver.Create( f{output_dir}/{var_name}.tif, ds.RasterXSize, ds.RasterYSize, 1, gdal.GDT_Float32) # 复制空间参考 out_ds.SetGeoTransform(ds.GetGeoTransform()) out_ds.SetProjection(ds.GetProjection()) out_ds.GetRasterBand(1).WriteArray(arr) out_ds.FlushCache()2.2 研究区裁剪与重投影做区域分析时经常需要将全球数据裁剪到研究区范围。这里有个效率技巧先投影再裁剪比先裁剪再投影要快3-5倍。以中国区域为例我推荐使用Albers等面积投影EPSG:102025这个投影能保持面积准确性特别适合降水总量这类需要面积计算的分析。用GDAL做投影转换时重采样方法的选择很关键气温数据用双线性插值bilinear降水数据用最邻近法nearest避免产生负值辐射数据用三次卷积cubic保持梯度平滑# 使用gdalwarp进行投影转换的示例 gdalwarp -s_srs EPSG:4326 -t_srs EPSG:102025 \ -r bilinear -tr 10000 10000 \ -of GTiff input.tif output_proj.tif3. 关键变量的单位换算ERA5-Land的原始数据单位可能不符合常规分析需求必须进行标准化处理3.1 降水数据转换原始降水单位是米/月要转为更常用的毫米/月。这里有个坑需要根据当月天数调整转换系数。我整理过各月份的换算系数表月份天数换算系数1,3,5,7,8,10,1231天×1000×314,6,9,1130天×1000×302月平年28天×1000×282月闰年29天×1000×29Python实现时可以用calendar模块自动判断闰年import calendar def convert_precip(month, year, raster): days_in_month calendar.monthrange(year, month)[1] return raster * 1000 * days_in_month3.2 温度数据转换2米气温原始单位是开尔文K转为摄氏度℃的公式很简单℃ K - 273.15但在处理多年数据时要注意ERA5-Land的温度缺省值是32767转换前需要先过滤temp_c np.where(temp_k 1000, np.nan, temp_k - 273.15)3.3 辐射数据转换地表向下太阳辐射的原始单位是焦耳/平方米J/m²转为更直观的千瓦时/平方米/天kWh/m²/day需要两步转换将月累计值转为日平均值÷当月天数单位转换÷3600000因为1kWh3.6×10⁶Jdef convert_radiation(month, year, raster): days calendar.monthrange(year, month)[1] return raster / (days * 3600000)4. 月度数据整合与分析4.1 时间序列构建处理多年数据时我推荐使用xarray库构建时间序列数据集比单独处理每个文件高效得多import xarray as xr # 合并多个年份的nc文件 ds xr.open_mfdataset(era5_*.nc, combineby_coords) # 提取中国区域 china ds.sel( latitudeslice(55, 15), longitudeslice(70, 140)) # 计算区域平均值 temp_mean china[t2m].mean(dim(latitude, longitude))4.2 典型分析场景农业干旱监测结合降水、温度和辐射数据可以计算标准化降水蒸散发指数SPEI。我常用的计算流程是用Thornthwaite方法计算潜在蒸散发PET计算降水与PET的差值P-PET进行标准化处理def calculate_spei(precip, temp, radiation): # 计算PET简化版 pet 16 * (10 * temp / radiation)**0.5 # 计算水分盈亏 balance precip - pet # 标准化 spei (balance - balance.mean()) / balance.std() return spei城市热岛分析比较城市与周边农村的温度差异时建议提取城市边界缓冲区如半径50km选取周边3-5个农村对照点计算月平均温差urban_mask urban_area.buffer(50000) rural_points [(lat1, lon1), (lat2, lon2)...] urban_temp temp.sel(urban_mask).mean() rural_temp temp.sel_points(rural_points).mean() heat_island urban_temp - rural_temp5. 常见问题解决方案在实际项目中遇到过几个典型问题这里分享我的解决方法问题1数据出现异常条纹原因通常是投影转换时的插值问题解决改用最邻近法重采样或使用GDAL的-et 0.1参数提高精度阈值问题2边缘区域出现NaN值原因WGS84与投影坐标系的范围不完全匹配解决转换时适当扩大处理范围完成后用研究区掩膜精确裁剪问题3跨时区的时间标注错误原因ERA5-Land使用UTC时间直接转换时可能出错解决用pandas的时区转换功能统一处理import pandas as pd # 将UTC时间转为本地时区 df[local_time] df[time].dt.tz_localize(UTC).dt.tz_convert(Asia/Shanghai)处理大型数据集时内存管理很关键。我的经验是对于超过5GB的数据改用分块处理chunking使用Dask延迟加载大数据文件及时关闭不再使用的文件句柄# 使用dask分块处理 import dask.array as da chunked da.from_array(big_array, chunks(1000, 1000)) result chunked.mean().compute()
从数据获取到分析应用：基于ERA5-Land月尺度降水、气温与辐射数据的全流程实践
1. ERA5-Land数据简介与获取方法ERA5-Land是欧洲中期天气预报中心ECMWF推出的高分辨率陆地表面再分析数据集空间分辨率达到9公里0.1°时间跨度从1950年延续至今。这个数据集特别适合需要精细分析地表气候特征的研究场景比如农业气象监测、流域水文模拟或区域气候变化研究。我处理过不少气象数据集ERA5-Land最让我印象深刻的是它的数据完整性。记得有次做青藏高原的气候分析其他数据集在高原边缘总有缺失但ERA5-Land的数据覆盖非常完整。数据集包含三个核心变量总降水量total precipitation月累计降水毫米数2米气温2m temperature近地表月平均气温开尔文单位地表向下太阳辐射Surface solar radiation downwards单位面积接收的短波辐射量获取数据最便捷的方式是通过Copernicus Climate Data StoreCDS的API。第一次使用时需要注册账号并在用户目录下创建.cdsapirc配置文件。这里分享一个我常用的Python下载脚本模板import cdsapi c cdsapi.Client() def download_era5land(year, month, variables, target_file): c.retrieve( reanalysis-era5-land-monthly-means, { format: netcdf, variable: variables, year: year, month: month, time: 00:00, product_type: monthly_averaged_reanalysis }, target_file)实际使用时要注意CDS的配额限制。我建议按年度分批下载比如先获取2010-2020年的目标区域数据。如果研究区域不大可以先用area参数限定下载范围北/西/南/东的经纬度这能显著减小文件体积。2. 数据预处理实战技巧2.1 格式转换与空间参考处理原始下载的NetCDF格式虽然包含丰富的元数据但很多GIS软件对它的支持有限。我习惯先用GDAL将其转为GeoTIFF这个过程中会遇到三个典型问题维度顺序混乱有时变量数据的维度排列是(time, lat, lon)有时却是(time, lon, lat)无效值处理海洋区域的填充值可能干扰后续计算坐标精度损失直接转换可能导致小数点后精度降低这是我优化过的转换代码加入了维度自动检测和异常值过滤from osgeo import gdal import numpy as np def safe_nc_to_tif(nc_path, var_name, output_dir): ds gdal.Open(fNETCDF:{nc_path}:{var_name}) band ds.GetRasterBand(1) arr band.ReadAsArray() # 处理无效值 arr[arr 1e10] np.nan # 自动检测维度顺序 if arr.shape[0] ds.RasterYSize: arr arr.T # 创建输出文件 driver gdal.GetDriverByName(GTiff) out_ds driver.Create( f{output_dir}/{var_name}.tif, ds.RasterXSize, ds.RasterYSize, 1, gdal.GDT_Float32) # 复制空间参考 out_ds.SetGeoTransform(ds.GetGeoTransform()) out_ds.SetProjection(ds.GetProjection()) out_ds.GetRasterBand(1).WriteArray(arr) out_ds.FlushCache()2.2 研究区裁剪与重投影做区域分析时经常需要将全球数据裁剪到研究区范围。这里有个效率技巧先投影再裁剪比先裁剪再投影要快3-5倍。以中国区域为例我推荐使用Albers等面积投影EPSG:102025这个投影能保持面积准确性特别适合降水总量这类需要面积计算的分析。用GDAL做投影转换时重采样方法的选择很关键气温数据用双线性插值bilinear降水数据用最邻近法nearest避免产生负值辐射数据用三次卷积cubic保持梯度平滑# 使用gdalwarp进行投影转换的示例 gdalwarp -s_srs EPSG:4326 -t_srs EPSG:102025 \ -r bilinear -tr 10000 10000 \ -of GTiff input.tif output_proj.tif3. 关键变量的单位换算ERA5-Land的原始数据单位可能不符合常规分析需求必须进行标准化处理3.1 降水数据转换原始降水单位是米/月要转为更常用的毫米/月。这里有个坑需要根据当月天数调整转换系数。我整理过各月份的换算系数表月份天数换算系数1,3,5,7,8,10,1231天×1000×314,6,9,1130天×1000×302月平年28天×1000×282月闰年29天×1000×29Python实现时可以用calendar模块自动判断闰年import calendar def convert_precip(month, year, raster): days_in_month calendar.monthrange(year, month)[1] return raster * 1000 * days_in_month3.2 温度数据转换2米气温原始单位是开尔文K转为摄氏度℃的公式很简单℃ K - 273.15但在处理多年数据时要注意ERA5-Land的温度缺省值是32767转换前需要先过滤temp_c np.where(temp_k 1000, np.nan, temp_k - 273.15)3.3 辐射数据转换地表向下太阳辐射的原始单位是焦耳/平方米J/m²转为更直观的千瓦时/平方米/天kWh/m²/day需要两步转换将月累计值转为日平均值÷当月天数单位转换÷3600000因为1kWh3.6×10⁶Jdef convert_radiation(month, year, raster): days calendar.monthrange(year, month)[1] return raster / (days * 3600000)4. 月度数据整合与分析4.1 时间序列构建处理多年数据时我推荐使用xarray库构建时间序列数据集比单独处理每个文件高效得多import xarray as xr # 合并多个年份的nc文件 ds xr.open_mfdataset(era5_*.nc, combineby_coords) # 提取中国区域 china ds.sel( latitudeslice(55, 15), longitudeslice(70, 140)) # 计算区域平均值 temp_mean china[t2m].mean(dim(latitude, longitude))4.2 典型分析场景农业干旱监测结合降水、温度和辐射数据可以计算标准化降水蒸散发指数SPEI。我常用的计算流程是用Thornthwaite方法计算潜在蒸散发PET计算降水与PET的差值P-PET进行标准化处理def calculate_spei(precip, temp, radiation): # 计算PET简化版 pet 16 * (10 * temp / radiation)**0.5 # 计算水分盈亏 balance precip - pet # 标准化 spei (balance - balance.mean()) / balance.std() return spei城市热岛分析比较城市与周边农村的温度差异时建议提取城市边界缓冲区如半径50km选取周边3-5个农村对照点计算月平均温差urban_mask urban_area.buffer(50000) rural_points [(lat1, lon1), (lat2, lon2)...] urban_temp temp.sel(urban_mask).mean() rural_temp temp.sel_points(rural_points).mean() heat_island urban_temp - rural_temp5. 常见问题解决方案在实际项目中遇到过几个典型问题这里分享我的解决方法问题1数据出现异常条纹原因通常是投影转换时的插值问题解决改用最邻近法重采样或使用GDAL的-et 0.1参数提高精度阈值问题2边缘区域出现NaN值原因WGS84与投影坐标系的范围不完全匹配解决转换时适当扩大处理范围完成后用研究区掩膜精确裁剪问题3跨时区的时间标注错误原因ERA5-Land使用UTC时间直接转换时可能出错解决用pandas的时区转换功能统一处理import pandas as pd # 将UTC时间转为本地时区 df[local_time] df[time].dt.tz_localize(UTC).dt.tz_convert(Asia/Shanghai)处理大型数据集时内存管理很关键。我的经验是对于超过5GB的数据改用分块处理chunking使用Dask延迟加载大数据文件及时关闭不再使用的文件句柄# 使用dask分块处理 import dask.array as da chunked da.from_array(big_array, chunks(1000, 1000)) result chunked.mean().compute()
