尧图网页设计 
简介本资源是一套面向地球物理、水文遥感及GRACE重力卫星数据应用研究者的MATLAB工具集聚焦于结合GLDAS陆面模型数据解算区域总水储量变化解决水文质量迁移引起的重力扰动建模与反演难题。压缩包共16个文件包含9个核心MATLAB脚本如main.m主流程、readPotentialCoefficients.m读取球谐系数、gravityDisturbance_fast.m高效重力扰动计算、totalWaterStorage_fast.m水储量快速反演、2个GRACE球谐系数gfc文件ITG-Grace2010系列、2份球谐函数理论PDF讲义、1个海岸线dat数据及辅助脚本整体仅1.52MB轻量实用。已有796人学习下载提供从GLDAS数据读取、勒让德函数计算、Love数加载、高斯滤波到大地水准面修正的完整处理链代码模块清晰、注释充分特别适合初涉GRACE水储量反演的研究者快速上手并理解各环节物理含义与数值实现细节。 看到“GRACE水储量解算_read_gldas_GLDAS_IWant!IWant_use”这个工程目录名我第一反应是亲切——这种命名方式十有八九是同一个脚本改了十几版、被各种read错误折磨到凌晨、最后终于“能用了”的产物。我刚做GRACE水储量解算那会儿也建过类似名字的文件夹。GRACE水储量解算的核心产物是月尺度的陆地水储量异常TWSA而要完成这套解算几乎绕不开GLDAS全球陆面数据同化系统这套陆面再分析数据。GLDAS提供了土壤水、雪水当量、冠层水等分量无论你拿它做地下水储量分离还是做泄漏误差校正第一步都得先把GLDAS数据准确、完整地读进来。这篇就把我从写read_gldas脚本到把它真正跑进GRACE解算流水线的过程捋一遍重点说清楚为什么read这一步最容易拖后腿以及那些年我遇到过的read系列报错到底该怎么治。1. GRACE水储量解算流水线里read_gldas为何卡在最前面1.1 卫星重力怎么变成“水有多深”GRACE和后续的GRACE-FO是两颗相距大约220公里的卫星星间用K波段微波测距持续测量距离变化。当地下水资源增加比如流域遭遇大范围降水入渗或者大型水库蓄水卫星下方的重力就会增强双星距离会跟着发生极其微小的变化。把这种变化按月份累积起来做球谐展开就能得到全球重力场模型的月解通常以球谐系数形式发布比如CSR、JPL、GFZ三家机构的RL06产品。但球谐系数不是水文人员直接能用的“水深”。重力场的变化是大气、海洋、陆水、固体地球内部质量变化的共同结果数据预处理时大气和海洋非潮汐部分已经被一并用去混叠产品扣掉了剩下的信号里陆地水储量变化占主导。要把球谐系数变成空间格网的等效水高需要代入负荷勒夫数公式做一次球面合成输出的单位是毫米水柱这才是TWSA月场。这个过程在数学上不复杂但实际处理里有几个绕不开的坎高阶球谐系数噪声大必须先做去相关滤波常用P4M6再做300公里左右的高斯平滑或者直接用DDK滤波否则TWSA图上全是南北向条带根本没法用。这也是为什么后面做泄漏误差校正时需要一个独立的陆面模型数据来配合——GLDAS就登场了。1.2 GLDAS在解算里的三个硬性用途GLDAS不是随便拿来对比的“参考数据”它在GRACE解算流程里有三个非常具体的用途。第一地下水储量分离。GRACE的TWS是整根土柱的水储量变化包括土壤水、地下水、雪水、地表水、植被冠层截留水。如果研究目标是地下水就必须把土壤水和雪水从GRACE里扣掉GLDAS恰恰提供了这些分量的月尺度估计。经典做法是GWS ≈ TWS_GRACE − TWS_GLDAS其中TWS_GLDAS取到根系深度或200厘米深度的土壤水加上雪水当量加冠层截留。这个式子看着简单但GLDAS分层定义、单位、时间分辨率都会直接影响结果。第二泄漏误差校正。GRACE格网数据经过滤波平滑后信号会从强水储量区“糊”到周边区域比如长江流域汛期的强信号会糊到海里这叫泄漏误差。恢复真实振幅的常用办法是用GLDAS或者水文模型生成一个先验TWS场对它做与GRACE完全相同的滤波和平滑然后比较原始场和平滑场的振幅比得到一个尺度因子再乘回GRACE数据。这个尺度因子的质量直接取决于GLDAS场是不是合理。第三独立交叉验证。GRACE反演出的季节振幅和相位和GLDAS模拟结果能否对得上是判断GRACE解算参数是否合理的第一道检查。如果某流域GRACE夏季TWS异常和GLDAS差得离谱且不是泄漏误差可以解释的那多半是滤波参数或者数据版本出了问题。所以read_gldas在整条流水线里的地位是“地基中的地基”——读错一个变量后面所有分析全白搭。2. GLDAS数据结构和read_gldas脚本的编写思路2.1 数据版本、变量命名和坐标体系我平时用得最多的是GLDAS-2.1陆面模式是NOAH空间分辨率0.25度时间范围从2000年延伸到当前提供3小时步长和月平均两种产品。月尺度水储量解算直接用月平均文件最省事省掉自己聚合几百个3小时文件的麻烦。GLDAS-2.1的变量命名有一个特点看着直观但一个不留神就会写错。土壤水是按深度分层命名的SoilMoi00_10cm_inst代表0到10厘米土壤湿度SoilMoi10_40cm_inst代表10到40厘米后面还有40到100厘米、100到200厘米单位都是kg/m²_inst表示瞬时值。雪水当量是SWE_inst冠层拦截水是CanopInt_inst。坐标体系上0.25度网格的lat数组是从北到南排列的也就是89.875、89.625……一路降到−89.875这个顺序问题很多人第一次画图时才意识到图怎么上下颠倒了。还需要注意一点GLDAS-2.2的CLM模型处理方式不同变量命名变成SoilMoist_SFC_inst表层、SoilMoist_RZ_inst根区、SoilMoist_PROF_inst深层剖面并非NOAH模型逐层命名的风格。你的脚本如果在2.1上跑通了直接拿去读2.2大概率会报变量找不到。所以写脚本一开始就别把变量名焊死在代码里。2.2 用NetCDF4还是xarray两个都要读取NetCDF格式Python生态里最底层的是netCDF4库再上层一点是xarray。我的实际建议是数据读入和重采样用xarray文件健康校验用netCDF4。xarray的好处在于标签索引不用手动管理纬度顺序和掩膜。比如我要把0.25度的GLDAS插值到1度网格直接用.interp(lat..., lon...)就行代码量很少。但xarray在碰到损坏文件时报错信息往往比较绕。有一次我批量处理时看到一个异常xarray报了一长串内部坐标错误后来用netCDF4.Dataset(file)单独一读才发现文件本身就是个只有0KB的坏文件。所以两个库配合用各管一摊效率最高。2.3 先定“整块读”还是“逐层读”写read_gldas脚本有个工程决策一次性把整个NetCDF读进内存还是按需读取变量。我一开始图省事直接用xarray.open_dataset(file_path)让它把整个文件载入结果一个月平均文件就把内存吃掉了不少连续处理几十个文件后Python进程越来越慢最后只能重启。后来改成“按需、逐变量、算完即关”的模式只读目标变量每个变量单独提取算完TWS后立刻close()释放文件句柄。对这种长时间序列批量处理来说这个模式最稳。核心读取函数大概长这样import xarray as xr import numpy as np NOAH_MONTHLY_VARS [ SoilMoi00_10cm_inst, SoilMoi10_40cm_inst, SoilMoi40_100cm_inst, SoilMoi100_200cm_inst, SWE_inst, CanopInt_inst, ] def read_gldas_monthly_twsa(file_path): 读取单个GLDAS-2.1月平均文件返回总水储量项(mm)。 直接使用xarray变量缺失时打印警告避免整个流程中断。 ds xr.open_dataset(file_path, enginenetcdf4) ds ds.sortby(lat) # 统一纬度升序避免和GRACE网格方向不一致 tws None for var in NOAH_MONTHLY_VARS: if var not in ds.variables: print(f[WARN] {var} not found in {file_path}) continue da ds[var].squeeze(dimtime, dropTrue).astype(float64) # 处理填充值GLDAS的无效值不能直接参与累加 fill_value ds[var].attrs.get(_FillValue, -9999.0) da da.where(da ! fill_value) tws da if tws is None else tws da ds.close() return tws.rename(tws_gldas)这个函数里有一个很多人容易忽略的点_FillValue。GLDAS某些格点是无效值值为−9999或者NaN直接累加会把无效污染到整块区域。用.where(da ! fill_value)把无效值置成NaN后面累加时NaN会自然传递你画图或统计时再统一处理缺失区域。3. 下载GLDAS数据read系列报错的根源与排查3.1 retrying连环刷屏连接中断与wget参数调优写批处理脚本下载GLDAS时最烦的就是终端刷屏warning: retrying (retry(total2, connectnone, readnone, redirectnone, statusnone))看着像是无害的警告其实它背后是服务器端连接被重置或者限流。retry(total2)意味着默认只重试两次两次失败后文件就下载失败了。我做2003到2016年月平均文件清单时差不多有五分之一的文件就是因为这个原因缺的。后来我把wget参数调成了这样wget -c -t 10 --waitretry30 --random-wait --timeout60 \ --load-cookies ~/.urs_cookies \ --save-cookies ~/.urs_cookies \ --auth-no-challenge \ -i download_list.txt几个关键参数解释一下-c断点续传。read中断后下次运行能从已下载的部分继续不用重头来。-t 10重试次数从默认2次提到10次。--waitretry30每次重试间隔30秒避免连续请求把限流触发得更厉害。--random-wait让请求间隔随机化降低整体请求频率特征。--timeout60防止连接长期卡在半死状态。另外NASA GES DISC的数据下载需要Earthdata账号认证。如果认证没配好服务器返回的其实是一个HTML登录页面你却把它存成了.nc4这种“文件”一读就是unrecognized format。所以下载之后第一件事不是急着读而是先检查文件大小和文件头。3.2 curl 18 RPC failed大文件传输中断的另一种典型场景如果你是从Git仓库同步GRACE/GLDAS配套脚本或者测试数据可能碰到这样的报错error: rpc failed; curl 18 transfer closed with outstanding read data remain这个错常见于git clone一个携带大二进制文件的仓库。Git传输大文件时默认的HTTP缓冲区相对有限中间如果经过网关或者代理单次请求传输的数据量过大就可能被掐断于是出现curl 18 transfer closed。我试过有效的处理方式git config http.postBuffer 524288000 git config http.version HTTP/1.1第一行把HTTP post缓冲区调大到500MB第二行强制使用HTTP/1.1而不是HTTP/2。HTTP/2的多路复用特性在这种大文件传输下更容易被中间设备判定为异常强制降级到HTTP/1.1反而更稳定。如果这样还是失败那就不要从Git仓库拉二进制数据直接去原始数据服务端下载压缩包绕开Git的传输层问题。3.3 cannot read file. unrecognized file format如何快速定位坏文件读取GLDAS时最直接的报错是cannot read file. unrecognized file format或者拼写略有不同cannot read file. unreconized file format这种报错出现时先别怀疑代码十有八九是文件坏了。GLDAS是标准NetCDF4格式只要文件头完整库就能正常打开。报unrecognized file format通常有三种情况文件是0KB下载彻底失败。文件是HTML登录页认证没通过却保存成了.nc4。文件下载到一半被截断NetCDF头不完整。我写了一个批量健康检查函数每次下载完先跑一遍把坏文件单独拿出来重下from pathlib import Path import netCDF4 as nc def check_gldas_files(folder): bad_files [] nc4_files sorted(Path(folder).glob(*.nc4)) for f in nc4_files: try: with nc.Dataset(f, r) as ds: _ ds.variables[time][:1] except Exception as e: bad_files.append((f.name, str(e))) return bad_files这个函数读每个文件的时间变量第一个值如果连这个都读不出来直接判定为坏文件。实际用下来它能抓出大多数下载残缺的.nc4比我以前一个个手工ncdump效率高太多。4. 从GLDAS到TWSA单位、分层与重采样三大陷阱4.1 单位换算的冷知识kg/m²和mm是同一个数很多教程在读取GLDAS时只说一句“单位换算一下”具体怎么算含糊带过。我第一次读GLDAS时也纠结过SoilMoi00_10cm_inst的单位是kg/m²是不是要乘以一个系数才能变成毫米水柱其实不用。推导一遍1 kg/m²表示1平方米面积上站着1千克水。水的密度是1000 kg/m³1千克水体积是0.001 m³铺在1平方米上厚度就是0.001米即1毫米。所以kg/m²直接等于毫米水柱数值不变变的只是单位含义。但有一个地方容易出错GLDAS有些变量的单位是kg m-2 s-1属于通量比如蒸散、径流这个可不能直接当水储量加。只看不带s-1的变量比如土壤湿度、雪水当量、冠层截留水才能这么换算。4.2 土壤分层不一致时的组合策略GLDAS-2.1的NOAH模型土壤分四层0到200厘米。但不同版本、不同模型的土壤深度定义完全不同。GLDAS-2.2的CLM模型分得更细最深可到几十米变量名也从逐层命名换成了表层/根区/深层剖面的逻辑。我的建议是把要累加的分量做成配置文件而不是在代码里硬编码TWS_STACK { GLDAS2.1_noah: [ SoilMoi00_10cm_inst, SoilMoi10_40cm_inst, SoilMoi40_100cm_inst, SoilMoi100_200cm_inst, SWE_inst, CanopInt_inst, ], GLDAS2.2_clm: [ SoilMoist_RZ_inst, SWE_inst, CanopInt_inst, ], }另一个容易被忽略的点是雪水当量SWE。在高纬度地区和青藏高原冬季SWE在整个TWS里占比相当可观。有一年处理数据发现某个站点的GRACE冬季TWS异常比GLDAS大很多查到最后是读取GLDAS时SWE变量被无效值掩膜给遮掉了不是物理差是空值处理不当导致的假差异。4.3 经纬度定义顺序和重采样细节GLDAS的lat数组从北到南递减而GRACE球谐合成出的网格通常从南到北排列。如果不统一顺序插值会错位画图也会上下颠倒。用xarray的话一行sortby(lat)解决。从0.25度降到1度有两种常见方案面积平均coarsen(lat4, lon4, boundarytrim).mean()。这个方法把每4×4个0.25度格点合成为1个1度格点物理意义更接近面积平均。双线性插值.interp(lat..., lon...)。速度快适合需要和GRACE格网逐点对齐的场景但会引入平滑效应。我一般区分场景用做泄漏校正或者面积对比时用coarsen面积平均做逐点时间序列对比时用双线性插值。别从头到尾只用一种不然在某些信号梯度大的区域两种方法能差出好几毫米水柱。5. 把read_gldas放进整个解算流水线工程化建议5.1 加日志、加缓存、加断点续算GRACEGLDAS解算不是跑一次就完的事。你会反复调整滤波参数、修改研究区范围、换时间跨度。如果read_gldas每次都要重新解析NetCDF几百个文件跑下来光是I/O就能耗掉大把时间。我的做法是增加一层中间缓存。把每个月算出的TWS_Gldas网格存成.npy或者zarr格式文件名带上GLDAS版本、处理日期、变量清单的哈希值。下次再跑时先查缓存命中就直接读取import numpy as np def load_tws_with_cache(file_path, cache_dircache): cache_path f{cache_dir}/{Path(file_path).stem}_tws.npy if Path(cache_path).exists(): return np.load(cache_path) tws read_gldas_monthly_twsa(file_path) Path(cache_dir).mkdir(exist_okTrue) np.save(cache_path, tws.values) return tws.values日志也要记录完整。每次读取文件至少把文件路径、处理时间、变量列表、数组shape写进日志。这个习惯帮我省了很多排查时间——某一个月的TWSA异常时我能很快确认是原始文件问题还是累计变量被改过。5.2 版本管理别学“IWant!IWant”命名回到标题那个文件夹名。我特别理解叫“read_gldas_GLDAS_IWant!IWant_use”这种名字的心情——同一个脚本改了无数遍想表达“这次真的要用了”。但以我的教训来说把情绪写进文件名的成本很高。三个月后你再看到这个名字根本分不清哪个版本是能跑的哪个版本是废弃的。真正有效的做法是把脚本纳入Git管理每次改动提交commit message写清楚改了什么。脚本命名用read_gldas.py这种稳定名称版本靠Git记录而不是靠文件名后面的“_v2_final_use”。另外强烈建议把整个解算流程的配置抽到一个config.yaml里GLDAS版本、GRACE数据处理机构CSR/JPL/GFZ、RL06还是RL05、滤波类型、高斯半径、研究区范围、时间区间。这样即使半年后回来看也能立刻知道跑出来的是什么参数组合的结果。我在这上面吃过亏曾经有一版结果和另一版结果对不上查了三天才发现是GLDAS版本混用了——一个月的文件来自本文还有配套的精品资源点击获取 GRACE水储量解算中GLDAS数据读取的工程实践与常见错误排查
简介本资源是一套面向地球物理、水文遥感及GRACE重力卫星数据应用研究者的MATLAB工具集聚焦于结合GLDAS陆面模型数据解算区域总水储量变化解决水文质量迁移引起的重力扰动建模与反演难题。压缩包共16个文件包含9个核心MATLAB脚本如main.m主流程、readPotentialCoefficients.m读取球谐系数、gravityDisturbance_fast.m高效重力扰动计算、totalWaterStorage_fast.m水储量快速反演、2个GRACE球谐系数gfc文件ITG-Grace2010系列、2份球谐函数理论PDF讲义、1个海岸线dat数据及辅助脚本整体仅1.52MB轻量实用。已有796人学习下载提供从GLDAS数据读取、勒让德函数计算、Love数加载、高斯滤波到大地水准面修正的完整处理链代码模块清晰、注释充分特别适合初涉GRACE水储量反演的研究者快速上手并理解各环节物理含义与数值实现细节。 看到“GRACE水储量解算_read_gldas_GLDAS_IWant!IWant_use”这个工程目录名我第一反应是亲切——这种命名方式十有八九是同一个脚本改了十几版、被各种read错误折磨到凌晨、最后终于“能用了”的产物。我刚做GRACE水储量解算那会儿也建过类似名字的文件夹。GRACE水储量解算的核心产物是月尺度的陆地水储量异常TWSA而要完成这套解算几乎绕不开GLDAS全球陆面数据同化系统这套陆面再分析数据。GLDAS提供了土壤水、雪水当量、冠层水等分量无论你拿它做地下水储量分离还是做泄漏误差校正第一步都得先把GLDAS数据准确、完整地读进来。这篇就把我从写read_gldas脚本到把它真正跑进GRACE解算流水线的过程捋一遍重点说清楚为什么read这一步最容易拖后腿以及那些年我遇到过的read系列报错到底该怎么治。1. GRACE水储量解算流水线里read_gldas为何卡在最前面1.1 卫星重力怎么变成“水有多深”GRACE和后续的GRACE-FO是两颗相距大约220公里的卫星星间用K波段微波测距持续测量距离变化。当地下水资源增加比如流域遭遇大范围降水入渗或者大型水库蓄水卫星下方的重力就会增强双星距离会跟着发生极其微小的变化。把这种变化按月份累积起来做球谐展开就能得到全球重力场模型的月解通常以球谐系数形式发布比如CSR、JPL、GFZ三家机构的RL06产品。但球谐系数不是水文人员直接能用的“水深”。重力场的变化是大气、海洋、陆水、固体地球内部质量变化的共同结果数据预处理时大气和海洋非潮汐部分已经被一并用去混叠产品扣掉了剩下的信号里陆地水储量变化占主导。要把球谐系数变成空间格网的等效水高需要代入负荷勒夫数公式做一次球面合成输出的单位是毫米水柱这才是TWSA月场。这个过程在数学上不复杂但实际处理里有几个绕不开的坎高阶球谐系数噪声大必须先做去相关滤波常用P4M6再做300公里左右的高斯平滑或者直接用DDK滤波否则TWSA图上全是南北向条带根本没法用。这也是为什么后面做泄漏误差校正时需要一个独立的陆面模型数据来配合——GLDAS就登场了。1.2 GLDAS在解算里的三个硬性用途GLDAS不是随便拿来对比的“参考数据”它在GRACE解算流程里有三个非常具体的用途。第一地下水储量分离。GRACE的TWS是整根土柱的水储量变化包括土壤水、地下水、雪水、地表水、植被冠层截留水。如果研究目标是地下水就必须把土壤水和雪水从GRACE里扣掉GLDAS恰恰提供了这些分量的月尺度估计。经典做法是GWS ≈ TWS_GRACE − TWS_GLDAS其中TWS_GLDAS取到根系深度或200厘米深度的土壤水加上雪水当量加冠层截留。这个式子看着简单但GLDAS分层定义、单位、时间分辨率都会直接影响结果。第二泄漏误差校正。GRACE格网数据经过滤波平滑后信号会从强水储量区“糊”到周边区域比如长江流域汛期的强信号会糊到海里这叫泄漏误差。恢复真实振幅的常用办法是用GLDAS或者水文模型生成一个先验TWS场对它做与GRACE完全相同的滤波和平滑然后比较原始场和平滑场的振幅比得到一个尺度因子再乘回GRACE数据。这个尺度因子的质量直接取决于GLDAS场是不是合理。第三独立交叉验证。GRACE反演出的季节振幅和相位和GLDAS模拟结果能否对得上是判断GRACE解算参数是否合理的第一道检查。如果某流域GRACE夏季TWS异常和GLDAS差得离谱且不是泄漏误差可以解释的那多半是滤波参数或者数据版本出了问题。所以read_gldas在整条流水线里的地位是“地基中的地基”——读错一个变量后面所有分析全白搭。2. GLDAS数据结构和read_gldas脚本的编写思路2.1 数据版本、变量命名和坐标体系我平时用得最多的是GLDAS-2.1陆面模式是NOAH空间分辨率0.25度时间范围从2000年延伸到当前提供3小时步长和月平均两种产品。月尺度水储量解算直接用月平均文件最省事省掉自己聚合几百个3小时文件的麻烦。GLDAS-2.1的变量命名有一个特点看着直观但一个不留神就会写错。土壤水是按深度分层命名的SoilMoi00_10cm_inst代表0到10厘米土壤湿度SoilMoi10_40cm_inst代表10到40厘米后面还有40到100厘米、100到200厘米单位都是kg/m²_inst表示瞬时值。雪水当量是SWE_inst冠层拦截水是CanopInt_inst。坐标体系上0.25度网格的lat数组是从北到南排列的也就是89.875、89.625……一路降到−89.875这个顺序问题很多人第一次画图时才意识到图怎么上下颠倒了。还需要注意一点GLDAS-2.2的CLM模型处理方式不同变量命名变成SoilMoist_SFC_inst表层、SoilMoist_RZ_inst根区、SoilMoist_PROF_inst深层剖面并非NOAH模型逐层命名的风格。你的脚本如果在2.1上跑通了直接拿去读2.2大概率会报变量找不到。所以写脚本一开始就别把变量名焊死在代码里。2.2 用NetCDF4还是xarray两个都要读取NetCDF格式Python生态里最底层的是netCDF4库再上层一点是xarray。我的实际建议是数据读入和重采样用xarray文件健康校验用netCDF4。xarray的好处在于标签索引不用手动管理纬度顺序和掩膜。比如我要把0.25度的GLDAS插值到1度网格直接用.interp(lat..., lon...)就行代码量很少。但xarray在碰到损坏文件时报错信息往往比较绕。有一次我批量处理时看到一个异常xarray报了一长串内部坐标错误后来用netCDF4.Dataset(file)单独一读才发现文件本身就是个只有0KB的坏文件。所以两个库配合用各管一摊效率最高。2.3 先定“整块读”还是“逐层读”写read_gldas脚本有个工程决策一次性把整个NetCDF读进内存还是按需读取变量。我一开始图省事直接用xarray.open_dataset(file_path)让它把整个文件载入结果一个月平均文件就把内存吃掉了不少连续处理几十个文件后Python进程越来越慢最后只能重启。后来改成“按需、逐变量、算完即关”的模式只读目标变量每个变量单独提取算完TWS后立刻close()释放文件句柄。对这种长时间序列批量处理来说这个模式最稳。核心读取函数大概长这样import xarray as xr import numpy as np NOAH_MONTHLY_VARS [ SoilMoi00_10cm_inst, SoilMoi10_40cm_inst, SoilMoi40_100cm_inst, SoilMoi100_200cm_inst, SWE_inst, CanopInt_inst, ] def read_gldas_monthly_twsa(file_path): 读取单个GLDAS-2.1月平均文件返回总水储量项(mm)。 直接使用xarray变量缺失时打印警告避免整个流程中断。 ds xr.open_dataset(file_path, enginenetcdf4) ds ds.sortby(lat) # 统一纬度升序避免和GRACE网格方向不一致 tws None for var in NOAH_MONTHLY_VARS: if var not in ds.variables: print(f[WARN] {var} not found in {file_path}) continue da ds[var].squeeze(dimtime, dropTrue).astype(float64) # 处理填充值GLDAS的无效值不能直接参与累加 fill_value ds[var].attrs.get(_FillValue, -9999.0) da da.where(da ! fill_value) tws da if tws is None else tws da ds.close() return tws.rename(tws_gldas)这个函数里有一个很多人容易忽略的点_FillValue。GLDAS某些格点是无效值值为−9999或者NaN直接累加会把无效污染到整块区域。用.where(da ! fill_value)把无效值置成NaN后面累加时NaN会自然传递你画图或统计时再统一处理缺失区域。3. 下载GLDAS数据read系列报错的根源与排查3.1 retrying连环刷屏连接中断与wget参数调优写批处理脚本下载GLDAS时最烦的就是终端刷屏warning: retrying (retry(total2, connectnone, readnone, redirectnone, statusnone))看着像是无害的警告其实它背后是服务器端连接被重置或者限流。retry(total2)意味着默认只重试两次两次失败后文件就下载失败了。我做2003到2016年月平均文件清单时差不多有五分之一的文件就是因为这个原因缺的。后来我把wget参数调成了这样wget -c -t 10 --waitretry30 --random-wait --timeout60 \ --load-cookies ~/.urs_cookies \ --save-cookies ~/.urs_cookies \ --auth-no-challenge \ -i download_list.txt几个关键参数解释一下-c断点续传。read中断后下次运行能从已下载的部分继续不用重头来。-t 10重试次数从默认2次提到10次。--waitretry30每次重试间隔30秒避免连续请求把限流触发得更厉害。--random-wait让请求间隔随机化降低整体请求频率特征。--timeout60防止连接长期卡在半死状态。另外NASA GES DISC的数据下载需要Earthdata账号认证。如果认证没配好服务器返回的其实是一个HTML登录页面你却把它存成了.nc4这种“文件”一读就是unrecognized format。所以下载之后第一件事不是急着读而是先检查文件大小和文件头。3.2 curl 18 RPC failed大文件传输中断的另一种典型场景如果你是从Git仓库同步GRACE/GLDAS配套脚本或者测试数据可能碰到这样的报错error: rpc failed; curl 18 transfer closed with outstanding read data remain这个错常见于git clone一个携带大二进制文件的仓库。Git传输大文件时默认的HTTP缓冲区相对有限中间如果经过网关或者代理单次请求传输的数据量过大就可能被掐断于是出现curl 18 transfer closed。我试过有效的处理方式git config http.postBuffer 524288000 git config http.version HTTP/1.1第一行把HTTP post缓冲区调大到500MB第二行强制使用HTTP/1.1而不是HTTP/2。HTTP/2的多路复用特性在这种大文件传输下更容易被中间设备判定为异常强制降级到HTTP/1.1反而更稳定。如果这样还是失败那就不要从Git仓库拉二进制数据直接去原始数据服务端下载压缩包绕开Git的传输层问题。3.3 cannot read file. unrecognized file format如何快速定位坏文件读取GLDAS时最直接的报错是cannot read file. unrecognized file format或者拼写略有不同cannot read file. unreconized file format这种报错出现时先别怀疑代码十有八九是文件坏了。GLDAS是标准NetCDF4格式只要文件头完整库就能正常打开。报unrecognized file format通常有三种情况文件是0KB下载彻底失败。文件是HTML登录页认证没通过却保存成了.nc4。文件下载到一半被截断NetCDF头不完整。我写了一个批量健康检查函数每次下载完先跑一遍把坏文件单独拿出来重下from pathlib import Path import netCDF4 as nc def check_gldas_files(folder): bad_files [] nc4_files sorted(Path(folder).glob(*.nc4)) for f in nc4_files: try: with nc.Dataset(f, r) as ds: _ ds.variables[time][:1] except Exception as e: bad_files.append((f.name, str(e))) return bad_files这个函数读每个文件的时间变量第一个值如果连这个都读不出来直接判定为坏文件。实际用下来它能抓出大多数下载残缺的.nc4比我以前一个个手工ncdump效率高太多。4. 从GLDAS到TWSA单位、分层与重采样三大陷阱4.1 单位换算的冷知识kg/m²和mm是同一个数很多教程在读取GLDAS时只说一句“单位换算一下”具体怎么算含糊带过。我第一次读GLDAS时也纠结过SoilMoi00_10cm_inst的单位是kg/m²是不是要乘以一个系数才能变成毫米水柱其实不用。推导一遍1 kg/m²表示1平方米面积上站着1千克水。水的密度是1000 kg/m³1千克水体积是0.001 m³铺在1平方米上厚度就是0.001米即1毫米。所以kg/m²直接等于毫米水柱数值不变变的只是单位含义。但有一个地方容易出错GLDAS有些变量的单位是kg m-2 s-1属于通量比如蒸散、径流这个可不能直接当水储量加。只看不带s-1的变量比如土壤湿度、雪水当量、冠层截留水才能这么换算。4.2 土壤分层不一致时的组合策略GLDAS-2.1的NOAH模型土壤分四层0到200厘米。但不同版本、不同模型的土壤深度定义完全不同。GLDAS-2.2的CLM模型分得更细最深可到几十米变量名也从逐层命名换成了表层/根区/深层剖面的逻辑。我的建议是把要累加的分量做成配置文件而不是在代码里硬编码TWS_STACK { GLDAS2.1_noah: [ SoilMoi00_10cm_inst, SoilMoi10_40cm_inst, SoilMoi40_100cm_inst, SoilMoi100_200cm_inst, SWE_inst, CanopInt_inst, ], GLDAS2.2_clm: [ SoilMoist_RZ_inst, SWE_inst, CanopInt_inst, ], }另一个容易被忽略的点是雪水当量SWE。在高纬度地区和青藏高原冬季SWE在整个TWS里占比相当可观。有一年处理数据发现某个站点的GRACE冬季TWS异常比GLDAS大很多查到最后是读取GLDAS时SWE变量被无效值掩膜给遮掉了不是物理差是空值处理不当导致的假差异。4.3 经纬度定义顺序和重采样细节GLDAS的lat数组从北到南递减而GRACE球谐合成出的网格通常从南到北排列。如果不统一顺序插值会错位画图也会上下颠倒。用xarray的话一行sortby(lat)解决。从0.25度降到1度有两种常见方案面积平均coarsen(lat4, lon4, boundarytrim).mean()。这个方法把每4×4个0.25度格点合成为1个1度格点物理意义更接近面积平均。双线性插值.interp(lat..., lon...)。速度快适合需要和GRACE格网逐点对齐的场景但会引入平滑效应。我一般区分场景用做泄漏校正或者面积对比时用coarsen面积平均做逐点时间序列对比时用双线性插值。别从头到尾只用一种不然在某些信号梯度大的区域两种方法能差出好几毫米水柱。5. 把read_gldas放进整个解算流水线工程化建议5.1 加日志、加缓存、加断点续算GRACEGLDAS解算不是跑一次就完的事。你会反复调整滤波参数、修改研究区范围、换时间跨度。如果read_gldas每次都要重新解析NetCDF几百个文件跑下来光是I/O就能耗掉大把时间。我的做法是增加一层中间缓存。把每个月算出的TWS_Gldas网格存成.npy或者zarr格式文件名带上GLDAS版本、处理日期、变量清单的哈希值。下次再跑时先查缓存命中就直接读取import numpy as np def load_tws_with_cache(file_path, cache_dircache): cache_path f{cache_dir}/{Path(file_path).stem}_tws.npy if Path(cache_path).exists(): return np.load(cache_path) tws read_gldas_monthly_twsa(file_path) Path(cache_dir).mkdir(exist_okTrue) np.save(cache_path, tws.values) return tws.values日志也要记录完整。每次读取文件至少把文件路径、处理时间、变量列表、数组shape写进日志。这个习惯帮我省了很多排查时间——某一个月的TWSA异常时我能很快确认是原始文件问题还是累计变量被改过。5.2 版本管理别学“IWant!IWant”命名回到标题那个文件夹名。我特别理解叫“read_gldas_GLDAS_IWant!IWant_use”这种名字的心情——同一个脚本改了无数遍想表达“这次真的要用了”。但以我的教训来说把情绪写进文件名的成本很高。三个月后你再看到这个名字根本分不清哪个版本是能跑的哪个版本是废弃的。真正有效的做法是把脚本纳入Git管理每次改动提交commit message写清楚改了什么。脚本命名用read_gldas.py这种稳定名称版本靠Git记录而不是靠文件名后面的“_v2_final_use”。另外强烈建议把整个解算流程的配置抽到一个config.yaml里GLDAS版本、GRACE数据处理机构CSR/JPL/GFZ、RL06还是RL05、滤波类型、高斯半径、研究区范围、时间区间。这样即使半年后回来看也能立刻知道跑出来的是什么参数组合的结果。我在这上面吃过亏曾经有一版结果和另一版结果对不上查了三天才发现是GLDAS版本混用了——一个月的文件来自本文还有配套的精品资源点击获取