简介陆面过程研究离不开再分析数据集而NetCDF格式作为科学数据存储的通用标准承载着大量全球尺度的地表参数。对于从事水文分析、干旱监测或气候变化研究的工程师而言理解数据格式与物理单位是开展工作的第一步。GLDAS作为全球陆面数据同化系统的代表产品提供了土壤湿度、雪水当量、蒸散发等关键变量但其文件名编码、存储单位及缺省值处理常让初学者却步。从数据解压到最终计算水储量异常TWSA每一步都隐藏着容易忽视的细节kg/m²与mm的等价转换、通量变量的时间缩放、经纬度切片的方向以及xarray库的高效应用。掌握这些基础技能才能将原始二进制数据转化为可信赖的区域水文时间序列。本文以实际案例梳理GLDAS数据的完整处理链路帮助读者规避常见陷阱高效完成水储量变化分析。 拿到一个叫GLDAS.zip的压缩包第一反应通常是「解压然后呢」。我最初接触 GLDAS 数据时也是这样压缩包解开后一堆.nc4文件文件名里全是NOAH025_M.A202401.020.nc4这种编码没有接触过 NetCDF 格式的人很容易愣在原地。这篇文章我就以这个 zip 包为起点完整讲清楚 GLDAS 数据的格式结构、单位陷阱、处理流程以及如何用它算水储量变化。适合刚接触陆面数据同化产品的研究生、工程师也适合想用 GLDAS 做区域水文分析但被格式和单位卡住的朋友。GLDASGlobal Land Data Assimilation System全球陆面数据同化系统是目前陆面过程研究里使用最广泛的再分析数据集之一它把观测数据降水、辐射等驱动进多个陆面模型输出土壤湿度、雪水当量、蒸散发、径流等变量。对水储量相关研究来说GLDAS 几乎是绕不开的起点数据。下面我按从拿到 zip 包到最终做分析的完整链路来写。1. 打开GLDAS.zip前必须搞明白的事版本、产品和文件名编码1.1 GLDAS 版本怎么选2.0、2.1 还是 2.2GLDAS 目前常见的有 GLDAS-2.0、GLDAS-2.1 和 GLDAS-2.2 三套产品很多人下载的时候直接在数据门户里搜 GLDAS结果一次下载几十个 GB其实根本用不到这么多。我自己最常用的是 GLDAS-2.1它使用 Noah 陆面模型驱动时间覆盖 2000 年至今空间分辨率 0.25 度大约 25-28 公里时间分辨率有 3 小时、日、月三种。做区域水储量变化、干旱监测、蒸散发估算2.1 是默认选择。GLDAS-2.0 覆盖 1948-2014 年适合做长序列气候分析但它在时间上不延续。GLDAS-2.2 加入了更多卫星观测同化精度理论上更高但数据文件更大处理起来也更费劲。如果你是初学者直接选 GLDAS-2.1 的月平均数据文件名中带M.A或者日平均数据带D.A够用且处理压力小。1.2 文件名编码拆解看见.nc4就知道里面是什么一个典型的 GLDAS-2.1 月数据文件名长这样GLDAS_NOAH025_M.A202401.021.nc4拆开看GLDAS_NOAH025产品名Noah 模型 0.25 度分辨率。M时间平均类型M月平均D日平均3H3小时平均。A202401数据时间2024 年 1 月。021GLDAS 版本号即 2.1。.nc4NetCDF4 格式。日数据文件名类似GLDAS_NOAH025_D.A20240101.021.nc43 小时数据是GLDAS_NOAH025_3H.A2024010100.021.nc4。看到文件名里的模式和版本号不用打开文件你就能判断这个包是不是你要的时间段和分辨率。我每次批量下载前都会先列文件名清单检查有没有缺月份、文件名是否统一避免后续批量处理时因为个别文件名出错而中断。1.3 从官网下载时容易忽略的细节GLDAS 数据托管在 NASA GES DISC需要注册 Earthdata 账号然后通过网站上的数据检索界面勾选变量、时间范围、区域范围最后生成下载链接。这里有几个容易踩的坑不要一次性勾选全部变量文件体积会非常大按需求勾选变量通常选SoilMoi00_10cm、SoilMoi10_40cm、SWE、CanopInt、Evap_tavg等几个核心变量即可。如果做区域分析在检索界面设置好经纬度范围服务端直接裁剪能省大量下载流量和处理时间。下载时注意文件粒度月数据一个月一个文件日数据一天一个文件如果你要 10 年日数据就是 3650 个文件不要用浏览器逐个点用提供的批量下载脚本wget 方式挂机下载。下载完的 zip 包先检查大小GLDAS 月数据单文件大约几 MB日数据单文件约 10-30 MB。如果某个 zip 文件明显偏小大概率下载不完整解压时就会报错。2. 解压和数据格式初探NetCDF4 的结构一句话讲清楚2.1 Linux 下解压 GLDAS.zip 的正确姿势拿到GLDAS.zip后Linux 下直接unzip GLDAS.zip -d GLDAS/如果文件比较多zip 包很大建议先测试完整性unzip -t GLDAS.zipunzip -t会逐个文件测试 CRC 校验如果输出No errors detected in compressed data of GLDAS.zip说明压缩包完好。如果某个文件报错通常是下载过程中网络问题导致文件损坏重新下载该文件即可不必全部重下。Windows 用户可以用 7-Zip 解压但要注意GLDAS 数据文件都是英文名路径不要包含中文否则后面用 Python 读取时会出现一些奇怪的路径编码问题。我碰到过在 Windows 下解压到带中文的目录后xarray 读取时报告文件不存在的情况本质上是系统编码不一致。2.2 NetCDF4 格式维度、变量、属性三要素解开后你会得到一个或多个.nc4文件这是 NetCDF4 格式。NetCDF 本质是一种自描述的二进制科学数据格式意思是文件里不仅存了数据还存了数据的维度、变量名、单位、坐标系等信息。用 Python 打开一个 GLDAS 月文件看结构import xarray as xr ds xr.open_dataset(GLDAS_NOAH025_M.A202401.021.nc4) print(ds)输出会列出类似这样的内容Dimensions: lat: 600 lon: 1440 time: 1 Coordinates: time: 2024-01-16 lat: -59.875 ... 89.875 lon: -179.875 ... 179.875 Data variables: SoilMoi00_10cm(time, lat, lon) float32 ... SoilMoi10_40cm(time, lat, lon) float32 ... SWE(time, lat, lon) float32 ... CanopInt(time, lat, lon) float32 ... ...三个核心概念维度dimensions描述数据轴GLDAS 的常规三维是time、lat、lon。变量variables实际的物理量每个变量都有名字、维度、数据类型和属性。属性attributes变量的元信息包括单位、缺省值、缩放因子等。单位能不能对上全靠这里。新手最容易犯的错误是不看变量属性直接拿原始数值去计算。后面参数、公式全对就是忘了缩放因子或者单位没转换结果算出来的水储量异常大几个数量级。2.3 ncdump 命令检查没有 Python 环境也能快速看结构如果你还没配置 Python 环境可以用 NetCDF 自带的ncdump命令查看文件头ncdump -h GLDAS_NOAH025_M.A202401.021.nc4-h表示只显示头信息包括所有维度和变量属性不打印全部数据。这个命令在服务器上排查问题时特别快不用启动 Python 解释器。比如你想确认这个文件里SWE变量的单位到底是kg/m^2还是mm一条命令就能看到。3. 数据单位里的陷阱kg/m²、mm 和实际物理量的换算3.1 为什么 GLDAS 变量单位不能直接用这是全文最重要的部分。GLDAS 数据里很多变量存储单位是kg/m²但在水文分析中你通常需要的是水柱高度mm或者体积含水量m³/m³。单位不换算后面所有结果都是错的。核心换算规则对水来说1 kg/m² 1 mm的水柱高度。因为水的密度约 1000 kg/m³面积 1 m²、深度 1 mm 的水体积是 0.001 m³质量是 1 kg。所以只要变量表示的是水的等效质量除不除 1000 都可以直接当 mm 用。但土壤湿度不同。SoilMoi00_10cm的含义是 0-10cm 土层中每平方米含有的水质量单位 kg/m²你要得到体积含水量m³/m³需要除以土层厚度米再除以水的密度1000 kg/m³体积含水量 SoilMoi00_10cm / 0.1 / 1000例如某点SoilMoi00_10cm 25 kg/m²那么 0-10cm 土层的平均体积含水量约为25 / 0.1 / 1000 0.25 m³/m³也就是 25% 的体积含水率。这个数值是合理的田间持水量范围。如果你不除以 0.1 和 1000直接用 25 当体积含水量那结果就会比实际夸张 100 倍。3.2 常用变量的单位对照表我把 GLDAS-2.1 Noah 产品里与水储量分析最相关的几个变量列成表方便查阅变量名含义存储单位换算目标换算方式SoilMoi00_10cm0-10cm土壤湿度kg/m²m³/m³数值 / 0.1 / 1000SoilMoi10_40cm10-40cm土壤湿度kg/m²m³/m³数值 / 0.3 / 1000SoilMoi40_100cm40-100cm土壤湿度kg/m²m³/m³数值 / 0.6 / 1000SoilMoi100_200cm100-200cm土壤湿度kg/m²m³/m³数值 / 1.0 / 1000SWE雪水当量kg/m²mm数值不变直接用CanopInt冠层截留水量kg/m²mm数值不变直接用Rainfall降水量通量kg/m²/smm/步长数值 × 步长秒数Evap_tavg蒸散发通量kg/m²/smm/步长数值 × 步长秒数Qs_acc地表径流通量kg/m²/smm/步长数值 × 步长秒数特别注意Rainfall、Evap_tavg这类变量它们在文件里是通量单位单位时间里的质量即 kg/m²/s你要算一天的总量必须乘以时间步长秒数。日数据的时间步长是 86400 秒月数据则要乘以该月的秒数。很多教程会直接忽略这一步结果是降水量的数量级完全不对。3.3 月平均数据做累积量时的隐藏问题如果你用的是月平均文件M.A里面的Rainfall_facc或Evap_tavg这类变量实际上是「该月单位时间内的平均通量」。也就是说文件里某个点的Evap_tavg值代表的是这个月每秒平均蒸散发多少 kg 水。要得到这个月的总蒸散发量mm/月需要乘以当月的总秒数。举例1 月文件里某点Evap_tavg 3.0e-05 kg/m²/s1 月有 31 天总秒数 2678400 秒那么该月总蒸散发为3.0e-05 × 2678400 80.35 kg/m² ≈ 80.35 mm80 mm 的月蒸散发量在湿润区是合理的。如果你忽略时间缩放直接用 3.0e-05 去画图出来的图几乎全是 0完全没法看。所以在写处理脚本之前一定要先print(ds[Evap_tavg].attrs)看一眼原始单位。GLDAS 不同版本、不同变量的单位偶尔会有差异不要迷信任何一个固定公式以文件属性里的units为准。4. 用 Python 把 GLDAS 数据变成能用的时间序列4.1 推荐工具链xarray rioxarray pandas处理 GLDAS 数据我目前最顺手的组合是xarray核心数据处理库天生适配 NetCDF 多维数组。rioxarray在 xarray 基础上增加地理坐标能力方便裁剪和导出 GeoTIFF。pandas时间序列处理和分析。numpy数值计算。安装pip install xarray rioxarray netcdf4 pandas numpynetcdf4是底层引擎xarray 需要它才能读取.nc4文件。不建议用scipy引擎读 NetCDF4对某些压缩变量的支持不好。4.2 读取并查看变量属性读取月文件并打印土壤湿度属性import xarray as xr ds xr.open_dataset(GLDAS_NOAH025_M.A202401.021.nc4) print(ds[SoilMoi00_10cm].attrs)输出通常包含{units: kg/m^2, long_name: Soil moisture content, 0-10 cm layer, _FillValue: 9.96921e36}看到_FillValue了吗这是缺省值的标记。GLDAS 数据在海洋和某些无数据区域会填一个极大值通常是 9.96921e36计算前必须掩膜掉否则均值、总和都会变成天文数字。4.3 区域裁剪和加权平均算一个流域的土壤湿度假设你要算某个流域的平均土壤湿度需要把 GLDAS 网格裁剪到流域边界内然后做面积加权平均。GLDAS 是规则经纬度网格直接用经纬度范围裁剪最简单import xarray as xr import numpy as np ds xr.open_dataset(GLDAS_NOAH025_M.A202401.021.nc4) # 假设流域经纬度范围 lat_min, lat_max 20, 30 lon_min, lon_max 100, 110 ds_sub ds.sel(latslice(lat_max, lat_min), lonslice(lon_min, lon_max)) # 土壤湿度换算为体积含水量 sm_0_10 (ds_sub[SoilMoi00_10cm] / 0.1 / 1000).compute() # 计算区域均值简单算术平均 regional_mean sm_0_10.mean(dim[lat, lon]) print(regional_mean.values)注意sel(latslice(lat_max, lat_min))中纬度必须从大到小因为 GLDAS 的纬度是从南到北递增的如果写成slice(lat_min, lat_max)会得到一个空的数组这个细节坑了我两次。如果要做精确的流域边界裁剪可以用rioxarray配合 shapefileimport rioxarray import geopandas as gpd ds xr.open_dataset(GLDAS_NOAH025_M.A202401.021.nc4) # 选择变量并设置空间参考 da ds[SWE] da da.rio.write_crs(EPSG:4326) # 读取流域边界 gdf gpd.read_file(watershed.shp) # 裁剪 da_clipped da.rio.clip(gdf.geometry, gdf.crs)这种方式比简单的经纬度裁剪更贴合自然流域边界也不会把边界外的网格数据混进来。4.4 提取网格点的长时间序列如果做站点对比需要提取某个经纬度网格的全部时间序列import xarray as xr import pandas as pd ds xr.open_dataset(GLDAS_NOAH025_D.A2024.tar) # 假设整合后的日数据 # 最近邻提取 point ds.sel(lat28.125, lon105.125, methodnearest) ts point[SoilMoi00_10cm].to_dataframe() print(ts.head())这里methodnearest会自动匹配最近的网格中心点避免因为经纬度精度不匹配导致报错。提取前建议先看一下目标位置的网格分辨率0.25 度网格在纬度 28 度附近经度方向的网格面积和纬度方向不一样如果做区域平均务必用面积加权cos(lat)权重否则高纬度网格会被过度代表。4.5 批量处理多个 zip 解压后的文件如果你的 zip 包里有多个文件比如一年的月数据共 12 个.nc4可以用xr.open_mfdataset直接拼接import xarray as xr ds xr.open_mfdataset(GLDAS_NOAH025_M.A2024*.nc4, combineby_coords)open_mfdataset会自动按时间维度拼接。需要注意文件路径里的通配符*要保证能匹配到所有目标文件。如果每个文件都有多余的time维度虽然只有一个时间点combineby_coords通常能自动处理。拼接后立即检查ds.time是否连续用pd.date_range对比一下缺月份是很常见的问题。5. 算水储量变化从各分量到总水储量异常5.1 GLDAS 中各水储量分量怎么加总GLDAS 模拟的总陆地水储量Terrestrial Water Storage, TWS可以近似表示为TWS 土壤含水量(各层之和) 雪水当量 冠层截留水量如果用 GLDAS 变量表达存储单位全部统一为等效水高mm土壤含水量0-200cm 四层总水高四层SoilMoi变量直接相加这里不用转体积含水量因为kg/m²等价于 mm 水柱。雪水当量SWE。冠层截留CanopInt。于是ds[TWS] ( ds[SoilMoi00_10cm] ds[SoilMoi10_40cm] ds[SoilMoi40_100cm] ds[SoilMoi100_200cm] ds[SWE] ds[CanopInt] )注意这里不需要除以 1000因为它本身就是等效水高。这也印证了前面说的kg/m²对水来说就是 mm四层土壤湿度加在一起直接是 0-200cm 的总水柱高度。5.2 水储量异常TWSA计算流程水储量异常Terrestrial Water Storage Anomaly是相对于某个基准期通常取多年平均的差值。GRACE 卫星反演的水储量异常常被用来和 GLDAS 结果对比验证。计算流程读取多年月平均 TWS 数据。计算每个格点的多年月平均气候态比如 2000-2024 年每个月的平均值。用每个月的实际 TWS 减去对应月份的气候态得到 TWSA。# 假设 ds 已经拼接了多年月数据并且包含了 TWS 变量 clim ds[TWS].groupby(time.month).mean(dimtime) anomaly ds[TWS] - climgroupby(time.month)是 xarray 计算气候态最简洁的方式。它会自动按月份分组求出 1-12 月各自的多年平均值然后与原始序列相减。得到的就是每个月相对气候态的偏差。需要特别说明一点GLDAS 模拟的 TWS 不包含地下水层位的完整变化它只模拟到 2 米深度Noah 模型。所以把 GLDAS 的 TWSA 与 GRACE 的 TWSA 比较时两者在绝对量级上会有差异但在季节性波动和干旱事件的相对变化上通常趋势一致。如果你想更精细地分析地下水变化建议结合 GRACE 数据减去 GLDAS 浅层水储量信号来反推。5.3 一个可以复现的完整案例某流域月 TWSA 时间序列我写一段完整可跑的脚本以 2020-2024 年月数据为例import xarray as xr import pandas as pd import matplotlib.pyplot as plt # 1. 读取多年月数据 ds xr.open_mfdataset(GLDAS_NOAH025_M.A202*.nc4, combineby_coords) # 2. 计算总水储量等效水高 mm ds[TWS] ( ds[SoilMoi00_10cm] ds[SoilMoi10_40cm] ds[SoilMoi40_100cm] ds[SoilMoi100_200cm] ds[SWE] ds[CanopInt] ) # 3. 裁剪目标区域比如长江上游某范围内 lat_min, lat_max 25, 35 lon_min, lon_max 90, 110 ds_region ds.sel(latslice(lat_max, lat_min), lonslice(lon_min, lon_max)) # 4. 面积加权区域平均 weights np.cos(np.deg2rad(ds_region[lat])) weights weights / weights.sum() region_tws (ds_region[TWS] * weights).sum(dim[lat, lon]) # 5. 计算气候态和异常 clim region_tws.groupby(time.month).mean(dimtime) anomaly region_tws - clim # 6. 画图 anomaly.plot() plt.show()这段脚本在 8 年以内的月数据量级上运行时间通常不超过 1 分钟。如果你要处理日数据建议先做月平均再算否则计算量会大 30 倍而且日数据的噪声很大趋势反而不容易看出来。6. 实际处理中绕不开的五个坑及排查思路6.1 文件损坏与 zip 解压失败的处理GLDAS 数据包下载频繁出现 zip 包损坏具体表现是解压到一半报错或者解压后发现ncdump -h无法读取文件。排查思路先用unzip -t xxx.zip测试压缩包完整性。如果是某个文件 CRC 校验失败重新下载该文件不要试图用unzip -o强制覆盖损坏的 NetCDF 文件即使强制解压出来也无法使用。如果压缩包是分卷 zip比如.z01.zip需要把所有分卷放在同一目录用 7-Zip 或zip -s合并后再解压。GLDAS 官方一般不提供分卷但有些高校镜像会分卷压缩。解压后批量检查文件是否有效我建议在 Linux 下用一条循环命令for f in GLDAS/*.nc4; do ncdump -h $f /dev/null 21 || echo $f broken; done这条命令会逐个测试 NetCDF 文件的完整性把损坏文件的文件名打印出来。这比写 Python 脚本快得多。6.2file is not a zip file问题的真实原因很多人在拿到 GLDAS 数据时会遇到file is not a zip file的报错。大多数情况下这个文件根本不是 zip 格式而是网页下载失败生成的一个 HTML 错误页面或者是一个 JSON 格式的下载链接列表。你用file命令查看file GLDAS.zip如果输出显示HTML document或ASCII text那就说明下载链接失效或者服务器返回了错误页面。解决办法是回到数据门户重新生成下载链接不要试图用任何 zip 修复工具去修复一个根本不是 zip 的文件。这个排查方法也适用于从各种数据共享平台下载的压缩包。6.3 经纬度范围写反切片后数据全空的排查GLDAS 的纬度坐标是递增的从 -90 到 90所以在用sel(latslice(高纬度, 低纬度))时顺序必须是从大到小。写反了不会报错只会得到空数组很多人在这里浪费几个小时。排查方法是直接查看打印结果ds_region ds.sel(latslice(lat_min, lat_max), lonslice(lon_min, lon_max)) print(ds_region[lat].values.shape) print(ds_region[lon].values.shape)如果打印出来的 shape 是(0,)说明切片方向写反了。解决方法是把lat的 slice 参数从大到小ds_region ds.sel(latslice(lat_max, lat_min), lonslice(lon_min, lon_max))6.4 时间编码问题日期显示成数字或时区偏移GLDAS NetCDF 文件的时间坐标通常是days since 1970-01-01或hours since 2000-01-01这种形式。xarray 一般能自动识别但如果你手动读取或者直接转换 DataFrame时间列可能显示为一串很大的数字。遇到time坐标显示为数字时用 pandas 转换time_index pd.to_datetime(ds[time].values)如果时区不对比如 UTC 时间与北京时间相差 8 小时还需要自行加 8 小时time_beijing time_index pd.Timedelta(hours8)对于日数据和月数据这个差异不影响月/日均值但如果你用 3 小时数据做逐时刻分析必须处理时区问题否则每个时刻都偏移 8 小时分析结果会完全错位。6.5 缺省值_FillValue没有掩膜统计结果全变大GLDAS 的_FillValue是 9.96921e36如果你在计算时忘记掩膜这个极端值会把均值、和、标准差全部拉到一个荒谬的量级。xarray 在大多数情况下会自动识别_FillValue并转为 NaN但以下情况容易被忽视你使用groupby计算气候态时如果某些缺失值一直存在计算 mean 可能仍然返回 NaN。你使用scipy引擎读取时缺省值可能没有被自动转换。你自己构造数组时没有设置mask。最稳妥的方法是读取后立即显式处理ds ds.where(ds ! 9.96921e36)这条语句把所有等于缺省值的元素置为 NaN。执行之后再做任何统计计算都更安全。7. 我把这套流程优化后的个人体会GLDAS 数据处理的门槛不在编程而在对数据本身的理解。很多初学者卡在单位换算和文件结构上其实只要养成先看变量属性、再动笔写处理代码的习惯就能避开大部分坑。我自己的固定动线是解压前检查 zip 完整性解压后用ncdump -h扫一遍文件头进入 Python 后先打印attrs确认单位和_FillValue最后才做裁剪、转换、拼接、统计。关于存储和计算的意见如果你想长期做区域水储量变化分析预先按行政边界或流域把全国 GLDAS 数据裁剪成小文件单独存成 netCDF 或 GeoTIFF会很方便。裁剪后单文件体积能缩小到原来的十分之一后续反复读取和分析的速度提升明显。我的习惯是保留原始经纬度网格不插值只在最终出图时按需重采样因为插值会引入额外误差而水储量分析最忌讳在数据源头引入人为误差。另一个实用技巧是GLDAS 的月数据和日数据可以混用。做趋势分析用月数据即可做极端降水事件对土壤湿度的快速响应分析时再用日数据甚至 3 小时数据。不要一开始就把所有时间分辨率的数据都下载下来存储和处理成本都会成倍增加。先明确科学问题需要的时空分辨率再有针对性地下载处理这才是效率最高的工作方式。最后再分享一个小操作处理完一批 GLDAS 数据后我会把每个 zip 包的来源链接、变量清单、处理脚本版本号记在一个 markdown 文件里和数据文件放在同一个目录。这个习惯在半年后回溯数据时会省下大量时间——尤其是当你发现单位换算想复核、或者处理脚本需要调整参数时不至于对着几 GB 的数据发呆想不起来当时是怎么做出来的。本文还有配套的精品资源点击获取