简介本资源是一套面向水文与气候领域初学者及科研人员的GLDAS数据处理入门工具集聚焦水储量估算这一典型应用场景解决遥感水文数据读取难、格式转换繁、多层积分计算复杂等实操痛点。压缩包共3个文件均为MATLAB脚本.m格式包括GLDAS数据读取readgldas.m、标准化水储量转换gldas2TWSt.m及时间序列插值与后处理TWSt2slept.m总大小仅8KB轻量易集成适合作为科研代码模块嵌入现有分析流程。已有1496人学习下载反映出该类基础处理脚本在高校水文建模、干旱监测课程实践及GRACE水储量验证研究中的高频需求。用户可直接调用脚本完成NetCDF格式GLDAS土壤湿度数据的加载、多层湿度加权积分、单位统一转换mm→cm或kg/m²、时间对齐及输出为标准水储量时间序列显著降低从原始数据到可分析结果的技术门槛。1. GLDAS 数据到底是什么能用来干什么如果你做水文、气候、农业遥感或者陆面过程模拟GLDAS 这个名字应该不陌生。中文全称是全球陆地数据同化系统Global Land Data Assimilation System由 NASA 戈达德太空飞行中心联合 NOAA 等多个机构开发。它最大的价值是把卫星遥感、雷达降水、气象站观测这些零散的输入统一喂给陆面过程模型通过数据同化技术输出一套全球范围、时间连续、空间全覆盖的陆面状态变量和通量变量。我最初接触 GLDAS是想算一个流域的水储量变化。当时手里的数据源很杂实测井数据稀疏、遥感土壤湿度有空间缺口想找一个能提供完整水循环分量的产品搜来搜去还是绕回 GLDAS。它一个文件里就能同时给你多层土壤湿度、雪水当量、蒸散发、降水、径流这对做水量平衡的人简直是全家桶。后来用得多了发现很多研究里它也是默认数据源比如干旱监测、蒸散发估算、洪水预报、地下水储量变化反演GLDAS 都是被反复引用的基础数据集。但说实话这个数据用起来远没有下载就能出图那么简单。最常见的问题就是单位。打开 NetCDF 文件变量单位写的是 kg/m²/s怎么换算成我们习惯的毫米四个土壤湿度层都是 kg/m²把它们加起来到底代表什么物理量这些基础问题不搞清楚后面做水储量计算就很容易差出两三个数量级。我见过有人拿 GLDAS 的月平均蒸散发直接画图数值小得离谱就是因为没有把每秒平均通量换算成每月总量。这篇文章我打算按实际工作中的思路来写从数据格式、单位体系、常规处理流程一直讲到拿 GLDAS 算水储量变化的完整代码和注意事项。适合刚开始接触 GLDAS、或者已经下载了 .nc4 文件但不知道怎么下手的读者也适合想用 GLDAS 和 GRACE 卫星数据做水储量变化对比的研究生和工程师。2. 数据格式满屏的 .nc4 到底是什么2.1 NetCDF 的自描述格式怎么理解GLDAS 数据绝大多数是 NetCDF 格式官方产品里常见 .nc 和 .nc4 后缀少部分老版本用 HDF 或者 GRIB。NetCDF 是一种自描述的二进制造型所谓自描述就是文件里不光存了数字矩阵还顺带把每个变量叫什么名字、单位是什么、多少维度、经纬度范围、时间分辨率这些元信息都写进去了。你只要用一个能解析 NetCDF 的工具打开文件就能直接看到这些信息不需要另外找配套文档。这个特点在实际处理里太重要了。比如你拿到一个 GLDAS 文件第一件事就是打印变量列表和单位确认里面到底有哪些字段而不是想当然地按文件名猜。xarray 或者 ncdump 都能干这个事后面会演示。2.2 文件命名规则与网格设计GLDAS 文件命名有一套固定规律以最常见的月平均产品为例GLDAS_NOAH025_M.A201501.021.nc4拆开看就是GLDAS数据集名NOAH025使用的是 Noah 陆面过程模型025 代表 0.25° 分辨率MMonthly表示月平均A201501A 表示数据起始时间为 2015 年 1 月021模型版本.nc4NetCDF4 格式GLDAS-2.1 的 0.25° 产品全球网格是 1440经度×600纬度经度从 -180 到 179.75纬度从 90 到 -60。注意纬度只到南纬 60 度南极洲那一块基本不在覆盖范围内。这个细节做全球制图、叠加海岸线的时候会暴露出来别到时候对着图找奇怪。另外很多老产品或者 GLDAS-1 的 1° 数据经度范围是 0 到 360而 GLDAS-2 的 0.25° 数据用 -180 到 180。把不同来源的数据混在一起用之前第一件事就是统一经度坐标系否则你裁剪出来的区域永远是空白或者图上一半数据跑到太平洋去了。2.3 时间分辨率与变量后缀GLDAS 有几种常见的时间分辨率3 小时一次、月平均部分产品还有 1 小时。3 小时数据适合做连续模拟或者拼日值文件数量很大一年就有 2920 个月平均数据比较轻量适合做长期趋势和季节分析。变量命名里藏着一个关键信息后缀。你会在 GLDAS 文件里看到Tair_f_inst、Rainf_f_tavg、Qs_acc这样的名字。其中_inst表示瞬时值比如某一时刻的气温_tavg表示时间平均比如一段时间内的平均降雨速率_acc表示累积量比如一段时间内的总径流这个区别是单位换算的钥匙。比如 3 小时产品的Rainf_f_tavg虽然是平均降雨速率但如果你只想要这一段时间的总降水直接乘上 3 小时对应的秒数就行。而Qs_acc本身已经是这个过程内的累积径流深单位就是 kg/m²等效毫米不需要再乘时间。很多新手在这里栽跟头后面我详细讲。3. 数据单位最容易翻车的环节3.1 先记住一个基准换算1 kg/m² 1 mm 等效水深GLDAS 的原始单位严格采用国际单位制这对科学计算是严谨的但对日常水文应用来说不太友好。整个单位换算的核心其实只需要记住一句话1 kg/m² 等于 1 mm 等效水深。为什么水的密度按 1000 kg/m³ 算1 mm 水深铺在 1 m² 面积上体积是 0.001 m³质量就是 1 kg。所以 GLDAS 里任何以 kg/m² 为单位的变量都可以直接理解成等效水深毫米。这个关系贯穿所有换算。有了这个基准其他换算统统推得出来蒸散发速率kg/m²/s想换成 mm/d乘以 86400一天秒数月平均蒸散发速率kg/m²/s想换成月总量 mm乘以该月天数再乘 86400降水速率kg/m²/s同样处理注意 2 月要按 28 或 29 天算3.2 土壤湿度分层的单位怎么理解GLDAS-2.1 Noah 模型的土壤湿度变量是分层的名字很直白SoilMoi0_10cm0-10 cm 土壤含水量SoilMoi10_40cm10-40 cm 土壤含水量SoilMoi40_100cm40-100 cm 土壤含水量SoilMoi100_200cm100-200 cm 土壤含水量单位都是 kg/m²也就是单位面积上这一层土壤里水的质量。这里有个容易混淆的点它给的不是体积含水量m³/m³而是等效水深。如果你想算整个 2 米土柱的水量直接把四层加起来得到的就是以毫米为单位的整层土壤有效水储量。这个特性在后续水储量计算里非常方便。如果想把 kg/m² 换算成体积含水量就按层厚除一下。比如 SoilMoi0_10cm 对应厚度 0.1 m那么体积含水量m³/m³ 数值kg/m²/ 0.1m/ 1000kg/m³。做站点对比或者和遥感土壤湿度产品标定时经常需要这个操作。3.3 蒸散发、降水与潜热通量的换算GLDAS 里蒸散发和降水通常以速率形式出现单位是 kg/m²/s也就是每秒蒸发多少毫米水深。月平均数据里这两个变量也是速率千万别当成月总量。换算代码很直接import xarray as xr # 打开一个月平均文件 ds xr.open_dataset(GLDAS_NOAH025_M.A201501.021.nc4) evap_rate ds[Evap_tavg] # 单位 kg/m²/s # 换算成 mm/d evap_mm_day evap_rate * 86400 # 换算成 2015年1月总蒸发量mm evap_mm_month evap_rate * (31 * 24 * 3600)另外GLDAS 能量通量变量里有一个潜热通量单位是 W/m²。如果你手里只有潜热通量想反推蒸散发公式是LE Lv × E其中 LE 是潜热通量W/m²Lv 是汽化潜热约 2.5×10⁶ J/kgE 是蒸散发速率kg/m²/s。所以蒸散发速率 LE / 2.5e6。这个公式在能量平衡闭合检验和蒸发互补法研究里经常用到。3.4 降雨速率与累积量的实践判断GLDAS 文件里Rainf_f_tavg这种变量表面上写着 average实际是该时段内的平均速率不是该时段累积量。3 小时产品的降雨如果算日总降水需要把当天 8 个时次分别乘上 10800 秒3 小时秒数再累加月平均产品则直接乘以时间段的秒数。我踩过的坑第一次拿 GLDAS 月平均降雨数据直接用 xarray 的 sum 沿时间维累加结果是 kg/m²/s 单位下的速率累加根本没有物理意义画出来的年降水总量高达数万毫米。后来养成习惯每次处理带通量性质的变量都先看单位、判断速率还是累积、再乘对应秒数一步都不省。提示判断变量是速率还是累积量最可靠的依据是单位。kg/m²/s 是速率kg/m² 是累积或状态量别只看变量名带没带 avg。4. GLDAS 数据获取与常规处理流程4.1 从 GES DISC 下载注意子集选择和账号GLDAS 官方数据主要在 NASA GES DISC 网站发布。搜索 GLDAS 可以找到 GLDAS-2.1 和 GLDAS-2.2 等集合。下载前建议先用网页端的子集功能选择区域、变量和时间范围这样生成的文件体积小很多。3 小时数据全球全年非常庞大一次拉全国数据动辄几十 GB不做子集选择机器根本带不动。GES DISC 需要注册账号才能下载用 Earthdata 账号登录。如果是在服务器上下载可以配置 .netrc 文件配合 wget 脚本这个网上都有教程不再展开。建议优先下载月平均产品做趋势分析用 3 小时产品做具体过程模拟。4.2 本地环境准备处理 GLDAS 数据我的固定组合是 Python xarray rioxarray。xarray 天然支持多维数组、时间切片、空间裁剪rioxarray 负责把 NetCDF 转成 GeoTIFF方便在 GIS 里展示。如果机器内存不够可以加 dask 实现分块读取。conda create -n gldas python3.10 conda activate gldas conda install -c conda-forge xarray dask netcdf4 h5netcdf rioxarray matplotlib cartopy4.3 读取文件和快速审查元数据拿到 .nc4 文件第一步是打印整个数据集结构看变量、维度、单位、坐标范围import xarray as xr ds xr.open_dataset(GLDAS_NOAH025_M.A201501.021.nc4) print(ds) # 单独查看某个变量的单位 print(ds[SoilMoi40_100cm].attrs[units])输出里会列出所有变量。GLDAS 月平均文件里常见的变量包括变量名含义单位SoilMoi0_10cm0-10cm 土壤含水量kg/m²SoilMoi40_100cm40-100cm 土壤含水量kg/m²SWE雪水当量kg/m²CanopInt冠层截留水kg/m²Evap_tavg蒸散发速率kg/m²/sRainf_f_tavg降雨速率kg/m²/sQs_acc地表径流累积量kg/m²Qsb_acc地下径流累积量kg/m²Tair_f_inst气温瞬时值K查看ds[lat].values能看到纬度是递减排列的从 90 到 -60步长 0.25 度。经度则从 -180 递增到 179.75。这个方向问题在画图时会遇到但 xarray 会自动识别坐标维一般不会出错。4.4 区域裁剪、时间重采样和聚合拿到数据以后最常见的操作就是把一个大范围数据集裁剪到研究区再把时间分辨率聚合到需要的尺度。比如我在处理华北地区时直接用sel方法按经纬度切片# 裁剪华北地区 region ds.sel(latslice(42, 32), lonslice(110, 122))注意latslice(42, 32)表示从北纬 42 到北纬 32因为纬度是递减的所以大的在前。这种写法和常规的 slice 不太一样很多人第一次会搞反。如果是 3 小时数据想合成日值用 xarray 的resample很顺手。对于降水速率类变量合成日总量需要先乘每个时次的秒数再求和对于瞬时温度直接取日均值就行。# 把3小时降水速率转成日累积量 daily_precip (precip * 10800).resample(time1D).sum() # 把3小时气温转成日均温 daily_temp temp.resample(time1D).mean()空间重采样也是常用操作。如果想从 0.25° 降到 1°最直接的办法是区域均值或者双线性插值。xarray 的coarsen方法适合整数倍聚合# 把0.25°聚合到1° coarse ds.coarsen(lat4, lon4, boundarytrim).mean()但要注意简单的块平均没有考虑球面面积权重严格做法是乘上纬度余弦权重后再平均。区域尺度不大的时候误差可以接受全球或大洲尺度建议认真算一下纬度权重。4.5 导出 GeoTIFF 供 GIS 使用水文和遥感团队经常需要把 GLDAS 结果放进 ArcGIS 或 QGIS 里叠加行政区划。用 rioxarray 转 GeoTIFF 非常简单import rioxarray # 给xarray对象设置空间维度 da region[SoilMoi40_100cm].rio.set_spatial_dims(x_dimlon, y_dimlat) da da.rio.write_crs(EPSG:4326) da.rio.to_raster(soil_moi_201501.tif)导出的 GeoTIFF 自带地理坐标GIS 里直接拖进来就能对齐。这里要注意的是如果数据原本纬度是递减的转 GeoTIFF 时 rioxarray 会正确处理不用手动翻转。5. 水储量变化计算GLDAS 的核心应用5.1 从 GLDAS 提取总水储量水文上说的水储量指一个区域所有形态储存在地表和地下的水分总量包括土壤水、雪水、冠层截留水、地表水、地下水等。GLDAS 能直接输出的有土壤水、雪水当量、冠层截留水所以一般用这个公式估算总水储量 TWS 四层土壤水之和 SWE CanopInt所有变量单位都是 kg/m²加起来直接就是毫米等效水深。这个 TWS 序列是后续计算水储量异常的基础。如果研究区在高寒地区SWE 的贡献要特别重视冬春季可能占很大比例。5.2 和 GRACE 卫星数据结合的思路GLDAS 水储量最有价值的应用之一是跟 GRACE 卫星反演的总水储量异常TWSA做对比。GRACE 测的是地球重力场的时变能反演区域总水储量变化但它测到的是所有水分的总和分不清到底是土壤水、地下水还是地表水。而 GLDAS 有过程模型约束能给出浅层土壤水和雪水的分量。很多研究里把 GRACE TWSA 减去 GLDAS 模拟的浅层水储量异常剩下的残差解释为地下水储量变化。这就是GRACE GLDAS 估算地下水变化的经典思路。实际应用时需要注意两者的空间分辨率和时间尺度要先统一通常把 GLDAS 重采样到 GRACE 的网格上或者做区域平均再比较。5.3 完整计算流程从文件列表到时间序列这里我给出一个可以直接改用的完整流程假设你已经下载了 2002 年以来的月平均 GLDAS 数据文件放在当前目录import glob import xarray as xr import numpy as np # 读取所有月平均文件 files sorted(glob.glob(GLDAS_NOAH025_M.A*.021.nc4)) ds xr.open_mfdataset(files, combineby_coords) # 计算总水储量单位kg/m²等效为mm tws (ds[SoilMoi0_10cm] ds[SoilMoi10_40cm] ds[SoilMoi40_100cm] ds[SoilMoi100_200cm] ds[SWE] ds[CanopInt]) # 裁剪到研究区域以华北平原为例 region tws.sel(latslice(40, 32), lonslice(110, 122)) # 区域平均 region_mean region.mean(dim[lat, lon], skipnaTrue) # 计算距平基线取2002-2010年 baseline region_mean.sel(timeslice(2002-01-01, 2010-12-31)).mean(time) anomaly region_mean - baseline # 画出时间序列 anomaly.plot()这段代码的核心思想就是先构造一个逐网格的水储量序列再做区域平均和基线距平。如果要进一步分析趋势可以对 anomaly 做线性回归或者用季节性分解模型把季节项和趋势项分开。5.4 计算水储量变化时的注意事项第一个要注意的是土壤层深度。GLDAS Noah 模型的土壤深度固定是 2 米不包含深层地下水和地表水体。所以在干旱区GLDAS 水储量变化的年际振幅通常比 GRACE 小很多因为它根本没有把深层地下水变化算进去。解读残差时这个系统性偏差要想清楚。第二个是区域平均的纬度权重。如果研究区从南到北跨了好几千公里直接用简单算术平均会高估高纬地区贡献。可以这样加权lat_weight np.cos(np.deg2rad(region[lat])) region_weighted region.weighted(lat_weight).mean(dim[lat, lon])第三个是月平均数据的时间标签。GLDAS 月平均文件的时间坐标通常是每个月 1 号表示这个值代表整月平均而不是 1 号当天瞬时值。画时间序列时把时间戳对齐到每月 1 号没问题但如果你想和站点月值对比要确认站点的月值是自然月平均还是月末瞬时值否则相位会差出半个月。6. 常见问题与排查技巧实录6.1 .nc4 文件打不开怎么办常见原因有两个一是环境里缺少 netCDF4 驱动二是文件本身可能是 HDF 容器。解决办法是在 open_dataset 时指定引擎ds xr.open_dataset(xxx.nc4, engineh5netcdf)或者提前安装好 netcdf4 库。如果两种引擎都试过了还是打不开用命令行 ncdump -h 看一眼文件头确认文件是否完整下载。6.2 画出来的图世界位置不对如果你发现数据在 0-360 经度范围而其他数据是 -180-180画图时美洲会被切到两段。统一经度的写法ds ds.assign_coords(lon(((ds.lon 180) % 360) - 180)).sortby(lon)这一步在合并多套数据源之前做比事后补救省事得多。6.3 蒸散发数值异常偏小或偏大蒸散发偏小十有八九是只用了速率没有乘时间秒数。月平均文件里Evap_tavg是每秒平均的速率量级大概在 0 - 0.0002 kg/m²/s 之间直接画图看起来特别小。换算成 mm/月之后华北地区夏季月蒸发差不多 100 mm 量级才算正常。反过来如果你发现数值偏大可能是把累积量和速率的处理搞反了。6.4 土壤湿度与站点实测对不上土壤湿度实测通常给的是体积含水量 m³/m³GLDAS 给的是 kg/m²。对比前先换算体积含水量 土壤水质量 / 层厚 / 1000。还要注意 GLDAS 每层厚度不一致0-10cm 是 0.1 m10-40cm 是 0.3 m40-100cm 是 0.6 m100-200cm 是 1 m别都按 0.1 m 除。6.5 3 小时数据太多内存爆掉3 小时数据全球范围年份下来就是几十 GB直接全部读进内存必然爆。建议分两种方式一是先裁剪研究区域再载入二是用 dask 分块计算。xarray 配合 dask 的 chunk 参数可以延迟计算只在最后触发计算时占用内存。如果只是要月平均可以先 cdo 或者 xarray 聚合完再存成轻量文件。6.6 时间坐标出现奇怪的数值部分从 FTP 下载的 GLDAS 数据时间坐标可能是一个相对数值而不是可读日期。可以用 pandas 转换import pandas as pd ds[time] pd.to_datetime(ds[time].values, unitD, origin2000-01-01)但 GES DISC 官方子集的 NetCDF 文件一般已经带标准时间轴遇到这种情况先确认数据来源不要盲目替换。6.7 处理批量文件的实用建议批量处理 GLDAS 文件我建议把流程分成三步第一步逐文件检查变量名和单位写成清单第二步统一空间坐标和时间编码第三步做裁剪、换算、聚合。每一步的结果都保存成中间 NetCDF方便排查。不要试图一个脚本从头跑到尾数据量大时中间任何一步出错回头排查的成本都会翻倍。我个人实操中的习惯是每套 GLDAS 文件拿到手先打印 data_vars 和 units再用变量名判断是速率还是累积量最后才开始做计算。这个流程帮我少走了很多弯路。如果后续要往深了做建议把 GLDAS 和 ERA5、CN05.1 等再分析产品在你的研究区做一次交叉对比看看各产品之间差异有多大不用纠结哪个绝对正确而是把这些差异当作不确定性来源写进结论里。这才是水文遥感研究里最有价值的部分。本文还有配套的精品资源点击获取