简介GLDAS.zip 面向气候研究、水文分析与遥感应用的学习者聚焦全球陆地数据同化系统GLDAS的读取、处理与水储量估算。资源包内含 3 个 .m 脚本文件压缩包约 8KB体量轻便主要提供 MATLAB 环境下的数据处理与转换代码便于在本地快速搭建 GLDAS 数据读取与计算流程。已有 1502 人学习下载说明该方向具备一定关注度。通过脚本可衔接 NetCDF 格式数据的读取、土壤湿度分层积分与总水储量TWS估算等关键环节帮助读者理解 GLDAS 各变量单位差异如毫米、摄氏度、毫米/日及时间空间重采样、裁剪等预处理思路并为进一步分析干旱、洪水与水文循环动态提供可复用的代码基础适合具备一定 MATLAB 与遥感基础的中高级用户参考。1. GLDAS 水储量数据落地从 zip 包到可用的月度异常图拿到一个叫GLDAS.zip的压缩包解压出来一堆.nc4文件变量名是TWS_tavg、单位写着kg m-2时间维度还是「距 2000-01-01 的小时数」——这是很多人第一次接触 GLDAS 水储量数据的真实场景。GLDASGlobal Land Data Assimilation System把 Noah、CLSM、VIC 等陆面模型和卫星观测同化到一起输出全球陆地表面的土壤湿度、雪水当量、冠层水、地下水等分量把它们按列相加就得到陆地水储量Terrestrial Water StorageTWS的近似值。这套数据常被用来做流域干旱监测、地下水亏损评估、重力卫星 GRACE 的对比验证。本文不聊概念史只讲一件事怎么把 GLDAS 的 zip 包变成一张能进论文、能进业务系统的水储量异常时间序列中间的单位、格式、缺测、尺度四个坑怎么绕过去。2. GLDAS 数据格式与单位先搞清楚 nc4 里到底存了什么2.1 为什么 GLDAS 是 nc4 而不是 GeoTIFFGLDAS 官方分发的是 NetCDF4.nc4格式不是 GeoTIFF也不是 CSV。原因在于它本质是四维数据时间 × 纬度 × 经度 × 变量。一个GLDAS_NOAH025_M.A202001.021.nc4文件里同时装着SoilMoi0_10cm_inst、SoilMoi10_40cm_inst、SWE_inst、CanopInt_inst、TWS_inst等十几个变量每个变量共享同一套经纬度网格。GeoTIFF 只能表达二维栅格硬转会把变量维度压掉后面做水储量分量求和时就得反复对齐文件反而更麻烦。NetCDF4 的另一个好处是自带元数据。用ncdump -h看一眼头信息单位、缺测值、时间基准全在里面ncdump -h GLDAS_NOAH025_M.A202001.021.nc4 | head -40输出里会看到类似TWS_inst:units kg m-2、time:units hours since 2000-01-01 00:00:00、_FillValue -9999.f这样的行。这三行决定了后面所有处理逻辑单位是面密度不是体积时间是小时偏移不是日期字符串缺测是 -9999 不是 NaN。2.2 水储量单位换算kg m-2 到 mm 到 km³GLDAS 水储量相关变量的单位是kg m-2。水的密度是 1000 kg m-3所以 1 kg m-2 的水层厚度就是 1 mm。这一步换算系数是 1不用乘任何东西但很多人会在这里翻车——把 kg m-2 当成 kg 直接乘面积结果量纲全错。真正需要换算的是从「等效水高mm」到「体积km³」import numpy as np # 假设 tws 是 (time, lat, lon) 的 numpy 数组单位 kg m-2 即 mm # 网格分辨率 0.25 度计算每个格点的实际面积 lat np.arange(-59.875, 90, 0.25) # GLDAS 纬度范围 lon np.arange(-179.875, 180, 0.25) R 6371.0 # 地球半径 km dlat np.deg2rad(0.25) dlon np.deg2rad(0.25) # 每个纬度圈上的格点面积 km2 area np.zeros((len(lat), len(lon))) for i, la in enumerate(lat): area[i, :] R**2 * dlat * dlon * np.cos(np.deg2rad(la)) # mm - km31 mm 水层 × 1 km2 1e-3 km3 tws_km3 tws * area[np.newaxis, :, :] * 1e-3这段代码的关键在np.cos(np.deg2rad(la))。高纬度格点的实际面积比赤道小如果直接用等经纬度面积乘青藏高原和亚马逊的体量会被系统性高估。参数上R6371.0是地球平均半径dlat、dlon是弧度制的网格步长1e-3是 mm·km² 到 km³ 的换算系数。注意GLDAS 的纬度是从北到南还是从南到北不同版本可能不一样。用ncdump看lat:units和实际数值顺序别默认。2.3 时间维解析小时偏移怎么变成 pandas 日期time:units hours since 2000-01-01 00:00:00意味着时间变量里存的是从 2000 年 1 月 1 日零点起算的小时数。Noah 月度产品里这个值通常是月初的小时数比如 2020 年 1 月对应 175320 小时左右。用 xarray 读的时候它会自动解码但如果你用netCDF4库手动读就得自己转import netCDF4 as nc import pandas as pd ds nc.Dataset(GLDAS_NOAH025_M.A202001.021.nc4) time_var ds.variables[time] times nc.num2date(time_var[:], time_var.units, only_use_cftime_datetimesFalse) df_time pd.to_datetime([str(t) for t in times]) print(df_time)nc.num2date的only_use_cftime_datetimesFalse参数很关键设成 True 会返回 cftime 对象pandas 处理起来多一层转换。输出应该是DatetimeIndex([2020-01-01], dtypedatetime64[ns])。如果输出是 1999 年或 2000 年 1 月 2 日说明时间基准理解错了回去检查time:units字符串。3. 用 xarray 批量处理 GLDAS.zip从解压到水储量分量求和3.1 解压与文件命名规律GLDAS 的 zip 包解压后文件名通常遵循GLDAS_模型_分辨率_时间频率.年月.版本.nc4的格式。比如GLDAS_NOAH025_M.A202001.021.nc4表示 Noah 模型、0.25 度、月度、2020 年 1 月、2.1 版本。批量处理时不要手动一个个读用glob按模式匹配import glob import xarray as xr files sorted(glob.glob(GLDAS_NOAH025_M.A*.021.nc4)) print(f共找到 {len(files)} 个月度文件) # 输出示例共找到 240 个月度文件2000-01 到 2019-12sorted不能省。glob 返回的顺序依赖文件系统不排序的话时间轴会乱后面算异常值时月度气候态就全错了。3.2 用 xarray 的 open_mfdataset 合并时间维xarray 的open_mfdataset能把多个 nc4 文件按时间维拼成一个 Dataset比循环netCDF4快一个数量级ds xr.open_mfdataset( files, combineby_coords, parallelTrue, chunks{time: 12, lat: 200, lon: 200} ) print(ds)combineby_coords让 xarray 根据时间坐标自动对齐不需要手动指定 concat_dim。chunks参数启用 dask 分块time: 12表示一次处理 12 个月内存不够时可以调小。parallelTrue在多文件场景下能明显加速但要注意 dask 的线程数别超过 CPU 核数。读进来之后检查一下变量列表print([v for v in ds.data_vars if TWS in v or SoilMoi in v or SWE in v])Noah 模型里水储量分量通常包括SoilMoi0_10cm_inst、SoilMoi10_40cm_inst、SoilMoi40_100cm_inst、SoilMoi100_200cm_inst、SWE_inst、CanopInt_inst。有些版本直接提供TWS_inst有些需要自己加。如果TWS_inst存在优先用它省去分量求和的麻烦。3.3 分量求和得到 TWS注意缺测传播如果没有现成的TWS_inst就手动求和soil_vars [SoilMoi0_10cm_inst, SoilMoi10_40cm_inst, SoilMoi40_100cm_inst, SoilMoi100_200cm_inst] tws ds[soil_vars[0]] for v in soil_vars[1:]: tws tws ds[v] tws tws ds[SWE_inst] ds[CanopInt_inst] tws.name TWS这里有个隐蔽的坑如果某个分量在某个格点上是_FillValue-9999直接相加会把 -9999 传播到结果里导致该格点 TWS 变成 -30000 多。正确做法是先做掩膜tws tws.where(tws -9000) # 把 -9999 变成 NaNwhere的条件阈值设 -9000 而不是 0因为有些干旱区土壤湿度确实接近 0但不会到 -9999 这个量级。这一步做完再算区域平均或画图时就不会被异常值拉偏。4. GLDAS 水储量避坑五个让结果全错的细节4.1 现象月度气候态算出来是 12 条平行线原因时间轴没有按月份分组或者分组时用了错误的time.month。xarray 的groupby(time.month)依赖时间坐标是 datetime 类型如果时间坐标还是数值型的小时偏移分组会失败或产生无意义结果。解决在读入后立刻检查ds.time.dtype如果不是datetime64[ns]用ds xr.decode_cf(ds)解码。然后clim tws.groupby(time.month).mean(dimtime) anom tws.groupby(time.month) - climanom就是水储量异常单位还是 mm。这一步是后续所有分析的基础气候态错了后面全错。4.2 现象青藏高原区域平均水储量比预期低一个量级原因纬度顺序反了。GLDAS 的纬度有些版本是从 90 到 -60 递减有些是从 -60 到 90 递增。如果代码里假设了递增顺序去算面积权重高纬度格点会被分配到错误的cos(lat)。解决读入后打印ds.lat.values[:5]和ds.lat.values[-5:]确认顺序。如果是递减用ds ds.sortby(lat)翻转。翻转后面积权重的计算才正确。4.3 现象某个月份全球平均 TWS 突然跳变原因GLDAS 不同版本之间可能有基准差异或者某个月份的文件缺失导致open_mfdataset跳过了那个月时间轴不连续。跳变往往出现在版本切换点或数据补丁月。解决合并后检查时间轴连续性time_diff ds.time.to_series().diff().dt.days print(time_diff[time_diff 40]) # 月度数据间隔应约 30 天如果有大于 40 天的间隔说明中间缺月。缺月不能简单插值要么补下对应月份的文件要么在分析时把该时段标记为缺测。4.4 现象用TWS_inst和手动求和结果不一致原因TWS_inst可能包含了手动求和没算进去的分量比如GroundWater_inst或SnowDepth_inst。不同 GLDAS 版本对 TWS 的定义略有差异。解决先看ncdump -h里TWS_inst的long_name和comment属性确认它包含哪些分量。如果定义不明确以手动求和为准并在论文里写清楚求和公式。两者差异超过 5% 时不要混用。4.5 现象算出来的水储量趋势在赤道附近是正的中纬度是负的原因没有去掉季节循环就直接做线性回归。水储量的季节振幅在中纬度很大直接回归会把季节信号混进趋势里。解决先算异常anom再对anom做趋势from scipy import stats import numpy as np # anom 是 (time, lat, lon) t np.arange(len(anom.time)) slope np.apply_along_axis( lambda y: stats.linregress(t, y)[0] if not np.isnan(y).all() else np.nan, axis0, arranom.values )stats.linregress返回的斜率单位是 mm/月乘以 12 就是 mm/年。这一步做完赤道附近的趋势应该接近 0中纬度干旱区应该出现负值这才符合物理预期。5. 从异常到趋势GLDAS 水储量的进阶验证与出图技巧5.1 用 GRACE 做交叉验证的快速方法GLDAS 的 TWS 异常可以和 GRACE 的等效水高做对比。GRACE 的空间分辨率粗约 300 km需要先降尺度到 GLDAS 的 0.25 度网格或者把 GLDAS 升尺度到 GRACE 的网格。我一般用后者因为升尺度比降尺度稳定# 把 GLDAS 0.25 度升尺度到 1 度便于和 GRACE 对比 tws_1deg tws.coarsen(lat4, lon4, boundarytrim).mean()coarsen的boundarytrim会丢掉边缘不完整的格点避免引入偏差。升尺度后两者的相关系数在大部分流域应该超过 0.6如果低于 0.4检查时间对齐和单位是否一致。5.2 出图时的色标与投影选择水储量异常图常用RdBu色标0 在中间红色负异常、蓝色正异常。但RdBu的红色端在打印成灰度时会偏暗如果论文要黑白印刷换成BrBG或PuOr。投影用PlateCarree就行GLDAS 是等经纬度网格不需要复杂投影。import matplotlib.pyplot as plt import cartopy.crs as ccrs fig plt.figure(figsize(10, 5)) ax plt.axes(projectionccrs.PlateCarree()) anom_mean anom.mean(dimtime) anom_mean.plot(axax, transformccrs.PlateCarree(), cmapRdBu, vmin-100, vmax100, cbar_kwargs{label: TWS anomaly (mm)}) ax.coastlines() plt.savefig(tws_anomaly.png, dpi300, bbox_inchestight)vmin、vmax设 ±100 mm 是经验值大部分区域的水储量异常在这个范围内。如果某个流域异常超过 200 mm单独调整色标范围不要用全局统一色标把其他区域压成一片白。5.3 一个我踩过的坑dask 惰性计算导致的内存爆炸用open_mfdataset读 240 个月的文件时如果直接.values转 numpydask 会一次性把所有数据拉进内存16 GB 内存的机器直接卡死。正确做法是分块计算或者先做时间平均再取 values# 错误tws.values 会触发全量加载 # 正确先降维再取 tws_clim tws.groupby(time.month).mean(dimtime).compute().compute()放在降维之后dask 只需要同时持有 12 个月的数据块内存占用降到原来的二十分之一。这个习惯我用了三年每次处理 GLDAS 都提醒自己一遍。5.4 参数速查表参数常用值说明网格分辨率0.25° / 1.0°Noah 是 0.25°CLSM 有 1.0°时间频率M月度/ 3H3 小时水储量趋势用月度就够缺测值-9999所有变量统一水储量单位kg m-2 mm换算系数 1体积换算× 面积 × 1e-3mm·km² 到 km³气候态分组time.month需要 datetime 时间坐标这张表我贴在显示器边上每次新建脚本先抄一遍省得反复查文档。做 GLDAS 水储量处理最深的教训是单位、时间基准、纬度顺序这三样每换一个数据版本就要重新确认一遍别信「上次就是这么做的」。我现在的习惯是拿到新 nc4 文件先跑三行检查——ncdump -h看单位print(ds.time)看时间print(ds.lat.values[:3])看纬度顺序。这三行花不了一分钟但能省掉后面三天的返工。希望帮到你。本文还有配套的精品资源点击获取