简介GLDAS全球陆地数据同化系统由NASA开发提供土壤湿度、地表温度、降水等网格化陆面数据是气候与水循环研究的重要数据源。这份资源面向需要下载并处理GLDAS数据、进一步估算水储量的科研人员和高年级本科生包含3个MATLAB脚本压缩包整体仅8KB。脚本覆盖了从NetCDF数据读取、变量提取到水储量计算的关键环节用户可直接调用或在此基础上修改无需从零编写全套代码。已有1503人学习下载。借助这批脚本读者可以理解GLDAS数据格式与各变量的单位含义学会通过土壤湿度分层积分估算陆地水储量并掌握结合降水、蒸散等要素分析水文动态的基本思路适合作为水文学、遥感或GIS相关课题的入门工具与代码参考。1. GLDAS 数据看着不难单位换算和格网才是最劝退的地方做干旱监测、流域水量平衡或者 GRACE 卫星水储量验证的人基本都要跟 GLDAS 数据打交道。但很多第一次拿到 GLDAS 数据包的人都会栽在“水储量”这一步明明按教程读出了土壤湿度、雪水当量一求和出来的数值却有十万毫米跟水位实测完全对不上。这不是 GLDAS 数据本身是黑匣子而是它的数据单位、分层深度、经纬度范围以及 zip 包解压后的文件格式每个环节都藏着约定俗成的“惯例”。这篇笔记会从 GLDAS 数据格式和单位讲起再把水储量计算的完整流程和常见坑拆开最后用 GRACE 数据做验证给出一套可以直接复现的处理路径。2. 读懂 GLDAS 数据格式从 zip 里的 NetCDF 到变量单位2.1 文件名拆解3 小时文件、月平均文件和你真正要的那个GLDAS 最常见的分发产品是 Noah 陆面模式驱动的 0.25 度全球数据时间分辨率分为 3 小时和月平均两套。下载站点一般允许你按变量、年份批量下载到了本地经常是打包成 zip 或 tar 的压缩包。文件名长得很像但用途完全不同例如文件命名时间分辨率常见用途GLDAS_NOAH025_3H.A20040101.0000.021.nc43 小时降水过程、土壤湿度日变化、蒸散发日尺度模拟GLDAS_NOAH025_M.A200401.021.nc4月平均月尺度水储量反演、GRACE 对比、长期趋势我做月尺度水储量时通常只用_M那套月平均文件。3 小时文件虽然能还原更多过程细节但单个年份就有 2920 个文件做全时段水储量序列时文件数量大、读取慢而且 3 小时数据里Rainf_tavg这类通量变量必须按时间步长累加才能变成“毫米”一步错就全错。GLDAS-2.1 版本开始Noah 产品默认是 NetCDF4 格式变量名也比 GLDAS-2.0 的 grib 风格规范很多。GLDAS-2.0 里土壤湿度还是SoilMoi00_10cm_inst这样的“层索引 瞬时值”命名到了 2.1 统一成SoilMoist_tavg单位也从kg m-2改成m3/m3这是最容易让人踩坑的一个变化。2.2 水储量相关变量单位对照表如果你要算“水储量”核心变量只有三个土壤湿度、雪水当量、冠层截留水。GLDAS Noah 的 NetCDF 变量属性里写得很清楚但单位并不统一变量名物理含义标准单位算水储量时的换算SoilMoist_tavg土壤体积含水量分 4 层m3/m3乘以每层厚度米再乘 1000得到毫米SWE_tavg雪水当量kg/m21 kg/m2 约等于 1 mm 水深CanopInt_tavg冠层截留水量kg/m21 kg/m2 约等于 1 mm 水深另外还要注意通量类变量Rainf_tavg和Evap_tavg的单位是kg m-2 s-1这是“速率”而不是“量”。有人直接把月平均文件里的Rainf_tavg当月降水量用出来的数值会小得离谱。正确做法是乘以对应时间段的秒数比如 3 小时步长的降水rain_mm Rainf_tavg * 10800如果是月平均文件里的降水速率就乘以该月天数再乘 86400。还有几个非水储量变量也常用到Tair_tavg单位是 KPsurf_f_inst单位是 PaPotEvap_tavg单位是 W/m2。做蒸发和干旱指数时经常要判断“这到底是不是瞬时值”命名里的_f_inst表示瞬时_tavg表示时间平均这个后缀比变量名本身更值得信任。2.3 格网坐标0~360 经度和陆表掩膜GLDAS 0.25 度网格的经度坐标范围经常是 0 到 360而不是我们熟悉的 -180 到 180。也就是说中国区域的经纬度会落在 73°E~135°E 附近这本身体感上没毛病但当你把 GLDAS 和 GRACE、CPC 降水或者站点数据对齐时不少工具默认经度是 -180~180直接画图会出现中国区域被劈到图幅两边的情况再做区域裁剪可能整片是 NaN。处理方式很简单先判断再转换if float(ds.lon.max()) 180: ds ds.assign_coords(lon(((ds.lon 180) % 360) - 180)).sortby(lon)这一点我会在后面避坑章节再展开。另一个和格网相关的属性是掩膜GLDAS 虽然覆盖全球但海洋、大湖这些格点通常用_FillValue填充。NetCDF 读取时如果不处理填充值计算流域平均值时这些格点会变成很大的数污染整条水储量序列。3. GLDAS 数据处理的复现路径zip 解压、xarray 读取与重采样3.1 先给 zip 包做体检testzip 和 EOCD 报错从网上下载的 GLDAS 压缩包最容易遇到的是“下载了一半”或“服务器返回了 HTML 错误页”导致的 zip 损坏。Python 在解压时常见报错是could not find EOCD或者解压到一半抛出BadZipFile。我一般不会直接zfile.extractall()而是先做一遍完整校验。import zipfile from pathlib import Path zip_path Path(GLDAS_NOAH025_M_2004.zip) out_dir Path(gldas_raw) out_dir.mkdir(exist_okTrue) with zipfile.ZipFile(zip_path) as zf: bad_file zf.testzip() if bad_file is not None: raise RuntimeError(fzip 内文件损坏: {bad_file}) zf.extractall(out_dir) nc_files sorted(out_dir.glob(GLDAS_NOAH025_M*.nc4)) print(f解压出 {len(nc_files)} 个月文件)testzip()会逐个读取压缩包内文件的 CRC 校验值返回第一个损坏的文件名如果所有文件完好返回None。这一步在批量下载几十 GB 数据时很值得做能避免后续xr.open_mfdataset在读取到一半时才报错。如果testzip()提示文件损坏常规做法是用下载工具的断点续传重新拉取比如curl -C -接着下载或者回到数据服务端重新生成订单。别用本地“修复 zip”的工具硬解NetCDF 文件被截断后即使能解出来变量维度也是残的后面时间段对齐会非常痛苦。3.2 xarray 读取 NetCDF变量清单和 FillValue 处理GLDAS-2.1 的 nc4 文件可以用 xarray 直接读取。因为是一堆月文件我会用open_mfdataset按时间维度合并。合并前先打开单个文件打印变量列表确认层维度名和填充值属性。import xarray as xr # 先看单文件结构 sample xr.open_dataset(nc_files[0], enginenetcdf4) print(sample) print(sample[SoilMoist_tavg].attrs) print(sample[SoilMoist_tavg].dims)打印结果里最有用的两个信息一个是SoilMoist_tavg的维度名是depth还是layer另一个是_FillValue和scale_factor。xarray 默认会读取_FillValue并转成 NaN但前提是你不要传mask_and_scaleFalse。有些从旧 grib 转出来的文件没有标准_FillValue而是用负值如-9999表示无效这种要手动替换。确认后再批量打开ds xr.open_mfdataset( nc_files, combineby_coords, enginenetcdf4, parallelTrue, ) print(ds.time.values[:3], ds.time.values[-3:])parallelTrue会按文件分块读取适合年份较多的场景。如果内存紧张可以只选择需要的变量读取比如variables[SoilMoist_tavg, SWE_tavg, CanopInt_tavg]这一步能省掉一半以上内存。3.3 把 3 小时数据统一到月尺度如果你手里只有 3 小时文件要把它合成月序列需要注意变量类型。土壤湿度、雪水当量这类“状态变量”应该按月做平均降水、蒸散这类“通量变量”应该按月做累加。ds3h xr.open_mfdataset( sorted(Path(gldas_raw_3h).glob(GLDAS_NOAH025_3H*.nc4)), combineby_coords, parallelTrue, ) # 状态变量月平均 sm_monthly ds3h[SoilMoist_tavg].resample(time1MS).mean(dimtime) # 通量变量按当月秒数累加为 mm secs_per_month ds3h[time].dt.days_in_month * 86400 rain_mm (ds3h[Rainf_tavg] * secs_per_month).resample(time1MS).sum(dimtime)这里的关键是resample的时间频率参数月平均应该用1MS每月 1 日或1ME每月末。如果错用1M部分 pandas 版本会提示频率歧义也有版本直接按自然月包含处理容易丢掉最后一个月。我用1MS是因为 GLDAS 月文件的时间戳大多落在每月 1 日 00:00便于和 GRACE 月数据对齐。4. 从 GLDAS 积出水储量层厚、单位换算与基线4.1 总水储量公式为什么只算到 2 米GLDAS Noah 的土壤湿度分四层层边界分别是 0~10cm、10~40cm、40~100cm、100~200cm。总水储量的“常规算法”是 2 米土柱里的水分加上雪水当量和冠层截留水TWS SMsurface~2m SWE CanopInt写成带单位的公式SMsurface~2m (mm) Σ (θ_i × thickness_i) × 1000其中 θ_i 是第 i 层土壤体积含水量单位m3/m3thickness_i 是第 i 层厚度单位 m。Noah 四层厚度分别取 0.1、0.3、0.6、1.0 米。如果你下载的 NetCDF 里depth坐标写的是层底深度 0.1/0.4/1.0/2.0先差分再使用import numpy as np depth_bottom np.array([0.0, 0.1, 0.4, 1.0, 2.0]) # 单位m thickness_m np.diff(depth_bottom) # [0.1, 0.3, 0.6, 1.0]为什么只积到 2 米一方面 GLDAS Noah 模式本身就只模拟到 2 米深没有深层地下水分量另一方面GRACE 卫星反演的总水储量包含深层地下水所以 GLDAS 与 GRACE 对比必然存在系统性差异。这不是数据错误而是物理含义不同。做验证时要清楚这个差异才不会把两者数值差当成异常。4.2 用 2004—2009 基线计算月距平水储量在不同地区的绝对值差异很大为了和 GRACE 对比通常要转成“距平”也就是去掉月气候态。取哪段时间做气候态是常见选择GRACE 数据处理里常用 2004—2009 年作为基线GLDAS 也取同样窗口这样两条序列的零点定义一致。# 土壤湿度m3/m3 - mm layer_thickness xr.DataArray( [0.1, 0.3, 0.6, 1.0], dimsdepth, coords{depth: ds[SoilMoist_tavg][depth]}, ) sm_mm (ds[SoilMoist_tavg] * layer_thickness * 1000).sum(dimdepth) # 雪水和冠层截留kg/m2 约等于 mm swe_mm ds[SWE_tavg] canop_mm ds[CanopInt_tavg] tws_mm sm_mm swe_mm canop_mm # 气候态均值 clim ( tws_mm.sel(timeslice(2004-01-01, 2009-12-31)) .groupby(time.month) .mean(dimtime) ) twsa tws_mm.groupby(time.month) - clim这里SoilMoist_tavg的维度如果叫layer而不是depth把dimsdepth改成dimslayer即可。代码里我把雪水当量直接当作毫米因为 1 kg/m2 的水深约等于 1 mm这对水储量序列的量级影响可以忽略。4.3 流域平均与保存结果区域尺度分析通常要输出流域平均曲线。简单做法是用经纬度矩形范围做sel但流域边界经常是不规则的。这里我推荐regionmask它支持直接把 shapefile 转成掩膜import geopandas as gpd import regionmask basin gpd.read_file(basin.shp) regions regionmask.from_geopandas(basin, namesbasin) mask regions.mask(twsa.lon, twsa.lat) # 返回区域编号 twsa_basin twsa.where(mask 0).mean(dim(lat, lon)) twsa_basin.to_netcdf(twsa_basin.nc)regions.mask(twsa.lon, twsa.lat)会生成一个二维数组落在流域内的格点编号为 0外部为 NaN。之后用where(mask 0)把其他区域滤掉再做空间平均。这里有一个细节如果流域面积太小比如小于一个 0.25 度格点平均结果会由少数格点主导甚至全是 NaN这种情况我会先做双线性插值到 0.05 度再计算。5. GLDAS 数据处理避坑5 个最容易翻车的现场5.1 把 m3/m3 直接当成 mm 叠加现象土壤湿度明明是 0.2~0.4直接加上雪水当量后总水储量只有几十毫米和 GRACE 的几百毫米量级差很远。原因SoilMoist_tavg的单位是体积含水量不是一个等量水深。只有乘以层厚才能转成毫米。解决先打印变量属性确认单位。如果是m3/m3必须乘以对应层的厚度单位米再乘 1000。这一步少做了后续所有分析和对比都没有意义。5.2 zip 解压报 “could not find EOCD”现象zipfile.ZipFile(GLDAS.zip)直接抛BadZipFile: File is not a zip file或者报could not find EOCD。原因压缩包下载不完整。很多数据订单是先排队再生成压缩包的下载工具断线后生成的 zip 只有几十字节到几百字节拿这种文件去解压只会得到错误。解决先看文件大小是否和网页上标注的一致。如果差太多用支持断点续传的工具重新下载。我常用的是curl -C - -O 下载地址续传完成后一定跑一遍testzip()再删除源文件。5.3 经纬度 0~360 未转换导致区域裁剪全 NaN现象用 shapefile 裁剪 GLDAS 数据后mean()结果全是 NaN画图时中国区域被分成两半。原因GLDAS 数据的经度范围是 0 到 360而 shapefile 和多数 GIS 工具使用 -180 到 180。解决按前面 2.3 节的代码做判断并转换。注意转换后要用sortby(lon)重排否则后面画图还是会出现锯齿状的经度不连续。5.4 月平均文件里的通量变量被当成总量现象月降水量算出来只有零点几毫米或者蒸散发为负值。原因GLDAS 月平均文件里的Rainf_tavg、Evap_tavg等变量属性写的单位是kg m-2 s-1是速率。必须乘以时间秒数而不是直接当作月总量。解决对月平均文件用days_in_month * 86400作乘法对 3 小时文件用固定的 10800 秒。注意days_in_month要取数据本身所在月份的天数不要用 365/12 这种平均天数否则 2 月的偏差会很明显。5.5 _FillValue 把海洋格点变成 1e20 量级现象水储量序列里突然出现一个巨大的尖峰且尖峰总是出现在固定格点或固定月份。原因GLDAS 的陆面变量在海洋、大湖区域填了_FillValue类似 1e20。流域平均时如果只做了mean()而没有mask_and_scale的过滤这些无效格点就会污染平均结果。解决用xr.open_dataset(..., mask_and_scaleTrue)或者手动where过滤掉大于 1e10 的值。无效值的边界一般用陆表掩膜变量也就是Land或mask字段把它们乘进去landmask ds[Land] # 陆地为1非陆地为0 twsa_clean twsa.where(landmask 1)6. 验证 GLDAS 水储量和 GRACE 对线的三个步骤拿到一条 GLDAS 水储量序列后先不要急着往报告里放。我会用 GRACE 或 GRACE-FO 的月水储量异常数据做一次交叉验证重点看“季节振幅对得上、趋势方向一致”。GRACE 反演的是总水储量异常TWSAGLDAS 少了深层地下水分量所以两者的绝对值差距不重要重要的是相位和振幅。第一个步骤是对齐基线和时间窗口。GRACE 处理中常用的基线是 2004—2009 年每个月的平均如果你的 GLDAS 序列也是这个窗口做的距平那么二者直接相减不会出现整体偏移。如果 GLDAS 用的是全时段气候态而 GRACE 用的 2004—2009 基线对比图里就会出现一条接近常数的位移这个位移不算物理信号。第二个步骤是去掉季节循环再看相关。GLDAS 和 GRACE 在月尺度上都以年周期主导直接算 Pearson 相关系数经常能到 0.9 以上但这并不代表模型好因为两者都“自带季节”。更严格的做法是做一个 12 个月滑动平均或者月度去趋势再看剩余信号的相关性。常用的验证代码gldas_twsa_12 twsa_basin.rolling(time12, centerTrue).mean() grace_twsa_12 grace_twsa.rolling(time12, centerTrue).mean() from scipy import stats r, p stats.pearsonr( gldas_twsa_12.dropna(dimtime), grace_twsa_12.dropna(dimtime), )第三个步骤是看散点图的斜率。理想情况下散点应围绕 1:1 线旁边分布但因为有地下水分量GLDAS 的振幅一般小于 GRACE散点回归斜率多在 0.6~0.9。如果斜率远小于 0.5先检查是不是只积了第一层土壤湿度如果斜率大于 1再检查雪水当量单位是否被重复放大了。我这几年的习惯是任何 GLDAS 水储量结果出来先打印变量单位再做testzip()最后用 GRACE 对一条相关。流程虽多但每一步都能挡住一次“结果看起来合理、实际全错”的返工。希望这个习惯能帮到你。本文还有配套的精品资源点击获取