GLDAS 数据产品的下载及使用


全球陆地数据同化系统(Global Land Data Assimilation System, GLDAS) 是美国宇航局(National Aeronautics and Space Administration, NASA) 与美国国家环境预报中心(National Centers for Environmental Prediction, NCEP)国家海洋大气局(National Oceanic and Atmospheric Administration, NOAA)的联合项目。GLDAS 采用先进的数据同化技术将卫星观测数据和地基观测数据整合到一个统一的模型中。目前 GLDAS 旗下共有四种陆地表面模型(Land Surface Model, LSM),分别为 Noah、Mosaic、CLM(Community Land Model)、VIC(Variable Infiltration Capacity)。这些模型以网格数据的形式集成了多种陆地表面场信息,例如全球的降水 (降雨和降雪)、蒸散 (土壤水分蒸发和植被蒸腾)、地表径流、地下径流、土壤湿度、地表温度、地表热流等。模型的空间分辨率有两种,分别为 $1^{\circ}\times1^{\circ}$ 和 $0.25^{\circ}\times0.25^{\circ}$,时间分辨率有三种,分别为 3 个小时、1 天、一个月,该数据产品可从 GES DISC(Goddard Earth Sciences Data and Information Services Center) 进行下载。

GLDAS 数据产品的下载

进入 GES DISC 的搜索页面,输入 gldas

GES DISC
GES DISC 主页

结果如下图所示

GES DISC 搜索结果
GLDAS 搜索结果

根据个人需求(起始时间、结束时间、时间分辨率、空间分辨率)搜索合适的数据。点击 “Subset/Get Data”,进入参数设置画面

GLDAS 数据产品的参数设置

然后设置时间范围、空间范围、变量、数据格式。设置完成后,点击 “Get Data”,服务器会生成一个 url 列表文件 ‘<url.txt>’,有效期为两天。根据该列表文件,这里使用 wget 或 curl 工具下载 GLDAS 数据产品(注意:下载文件前需要注册一个 earthdata 账号),具体步骤如下:

  1. 在主目录下创建 .netrc 文件;
    • cd ~ or cd $HOME
    • touch .netrc
    • echo "machine urs.earthdata.nasa.gov login <uid> password <password>" >> .netrc (这里 ‘<uid>’ 为用户名,'<password>’ 为 Earthdata 登陆密码。)
    • chmod 0600 .netrc (设置权限)
  2. 创建 cookie 文件;
    • cd ~ or cd $HOME
    • touch .urs_cookies
  3. 根据 ‘<url.txt>’ 列表文件,用 wget 或 curl 工具下载数据;
    • wget --load-cookies ~/.urs_cookies --save-cookies ~/.urs_cookies --auth-no-challenge=on --keep-session-cookies --content-disposition -i <url.txt>
    • cat <url.txt> | tr -d '\r' | xargs -n 1 curl -LJO -n -c ~/.urs_cookies -b ~/.urs_cookies

首次下载需要执行上述所有的步骤,之后的下载仅需要执行第三步即可。关于 GLDAS 数据下载的详细说明可参见 Instruction to Downloading Data。GLDAS 数据产品默认为 netcdf 格式,关于该格式的读取和写入可参见 NetCDF4.0 与 HDF5.0 文件的读取与生成

GLDAS 数据产品的处理

陆地水储量变化的估算

使用 GLDAS 数据产品估算一段时期内( 2002 年 4 月到 2016 年 8 月)全球的陆地水储量(Terrestrial Water Storage, TWS)变化。下载 netcdf 格式的 GLDAS 数据产品(空间范围为 [$179.875^{\circ}W \sim 179.875^{\circ}E$, $59.875^{\circ}S \sim 89.875^{\circ}N$],空间分辨率为 $0.25^{\circ}\times0.25^{\circ}$,时间分辨率为 1 个月)。读取所有的 GLDAS 数据产品(共 173 个月文件),然后根据水平衡方程 $$H_{k}=\frac{1}{\rho_{w}} \sum_{i=1}^{k}\left(P_{i}-E T_{i}-R_{i}\right)$$ 计算每个月的陆地水储量变化,其中 $H_k$ 为第 $k$ 个月相对于第 0 个月的陆地水储量变化 (等效水厚);$k$ 的取值为 $1 \sim n$;$n$ 为总月数;$P_i$ 为第 $i$ 个月的降雨量与降雪量之和;$ET_i$ 为第 $i$ 个月的土壤水分蒸发量与植被蒸腾量之和;$R_i$ 为第 $i$ 个月的地表径流量与地下径流量之和。计算 2002 年 4 月到 2016 年 8 月这段时期内的平均陆地水储量变化,然后对每个月的陆地水储量变化进行去平均化处理。由于 GLDAS 数据产品只在陆地覆盖区(不包括南极洲大陆)有数据,而海洋覆盖区均以 $-9999.0$ 进行了填充,因此为了之后的窗口球谐展开,这里已将填充值修改为 0。下面的 Python 小程序估算了 2002 年 4 月到 2016 年 8 月去平均化的全球陆地水储量变化。

陆地水储量变化的单位为 mm(已经将质量转换为等效水厚),其空间范围为 [$179.875^{\circ}W \sim 179.875^{\circ}E$,$59.875^{\circ}S \sim 89.875^{\circ}N$],网格维度为 $600 \times 1440$。为了后面的窗口球谐展开,这里需要调整其空间范围为 [$0.125^{\circ}E \sim 359.875^{\circ}E$,$89.875^{\circ}N \sim 89.875^{\circ}S$],网格维度为 $720 \times 1440$。

陆地水储量变化的球谐展开

SHTOOLS 工具对扩展后的陆地水储量变化进行球谐展开,务必要将网格数据转化为 nxn 或 nx2n 的形式,使数据的第一行对应 $90^{\circ}N$,第一列对应 $0^{\circ}$,因此需要对陆地水储量变化插值到 SHTOOLS 要求的形式。这里用 scipy 中的球面插值函数 RectSphereBivariateSpline 实现对陆地水储量变化的插值操作。

上面的代码实现了对陆地水储量变化的插值,但有几点需要额外注意:

  • RectSphereBivariateSpline 插值函数使用的空间变量为余纬和经度(单位为弧度),其范围分别是 $0\sim\pi$ 和 $0\sim2\pi$ 或 $-\pi\sim\pi$;
  • 插值后的空间范围务必要与插值前的空间范围一致,例如插值前的空间范围为 $0\sim2\pi$,那么插值后的空间范围也必须为 $0\sim2\pi$。若插值前的空间范围为 $-\pi \sim \pi$,而插值后的范围为 $0 \sim 2\pi$,则会得到错误的结果。

对插值后的陆地水储量进行窗口球谐展开前需要建立陆地窗口(不包含南极洲大陆)。这里使用全球地形起伏数据 ETOPO5 来建立陆地窗口:海拔大于零的网格点认为位于陆地并赋值为 1,否则位于海洋并赋值为 0;$60^{\circ}S$ 以南的区域均赋值为 0。关于全球地形起伏数据 ETOPO5 的更多详情可参考用全球地形起伏数据 ETOPO 作图球谐展开与谱分析。下面的小程序完成了陆地窗口的建立,值得注意的是该陆地窗口的阶数为 96,用户可根据需要(例如 60 阶)来自行调整程序的参数。补充:若陆地窗口的阶数为 60,则没必要使用空间分辨率为 $0.25^{\circ}$ 的 GLDAS 数据,使用空间分辨率为 $1^{\circ}$ 的 GLDAS 数据完全满足要求;因为较小的空间分辨率对应较长的计算时间,因此建议 60 阶的用户使用空间分辨率为 $1^{\circ}$ 的 GLDAS 数据。

输出结果为

结果显示 shannon 数为 2467,这表明至少需要 2467 个 slepian 基函数才能较完整地完成陆地水储量变化的窗口球谐展开。根据这 2467 个 slepian 基函数,下面的小程序计算了每个月陆地水储量变化的球谐系数,并将球谐系数保存在文件 TWSCs_interp_land_coeffs.npy 中以供后续的再处理。关于窗口球谐展开的更多详情可参见球谐展开与谱分析

输出结果为

结果显示球谐系数的维度为 (173, 2, 97, 97),其中 173 表示共有 173 个月,2 表示 $\mathcal{C}_{lm}$、$\mathcal{S}_{lm}$,97 为 0…96 阶和次。补充:陆地水储量变化的窗口球谐系数与一般的球谐系数(认为海洋区域与南极洲大陆的填充值 0 是有效值)有些许差别,使用时请注意。


参考

  • Rodell M, Houser P R, Jambor U E A, et al. The global land data assimilation system[J]. Bulletin of the American Meteorological Society, 2004, 85(3): 381-394.

关于 “GLDAS 数据产品的下载及使用” 的 40 个意见

  1. 博士哥哥你好,请问对球谐系数截断到与 CRS 发布的 GRACE 月重力场数据相同的阶数 60阶的话代码需要怎么修改呢

  2. 师兄这篇博文很有用,有一个地方不太懂。为什么 “TWSC 的单位kg/m^2, 除以水的密度(1000kg/m^3)就变成等效水高单位 mm 呢”

      1. 谢谢师兄。比如我想计算 土壤含水量数据 造成的水储量变化,GLDAS中土壤含水量数据分为四层,那么我计算的时候是否需要考虑每层的厚度问题呢?

        1. 把这四层的土壤含水量加起来就行了,但这只能代表表层土壤水。实际上的陆地水储量包括表层土壤水、深层土壤水、地下水、湖泊和洼地的水、冰雪融水等。

      2. 您好,转换等效水高这里有点不明白。1 kg/m^2 / 1000kg/m^3 = 1 mm,可是GLDAS网站给出kg/m^2 转换为mm的方法是:
        [kg m-2] —-> [converting to] —-> [mm]

        kg m^3 1000 mm
        — X —– X ——- = mm
        m^2 1000 kg 1 m

        也就是说 If the assumption is made that the density of the water in the soil is 1000 [kg m-3], then the value in [kg m-2] is identical to the value in [mm]。如果水的密度是1000 [kg m-3],那么kg/m^2与mm大约相同。那么GLDAS中的kg/m^2到底如何转换为等效水高呢?

        1. 单位面积上的质量(kg/m^2)除以水的密度(1000kg/m^3)即为等效水厚;除以冰的密度(917kg/m^3)即为等效冰厚;除以沙的密度(1602kg/m^3)即为等效沙厚,以此类推。[1 kg/m^2] / [1000kg/m^3] = 1 mm 与 GLDAS 网站给出kg/m^2 转换为 mm 的方法一致。

  3. 请问博士大哥,这行代码”precipitation_flux = precipitation_flux[0]“为什么表示降维的意思,能否详细解释一下?

    1. 原始的 data[‘Rainf_f_tavg’][:] 的维度为 (1,nlat,nlon),其中 1 代表时间维,也就是一个月内的降水;为了简化数据,这里 “precipitation_flux[0]” 是将维度为 (1,nlat,nlon) 的数据降为维度为 (nlat,nlon) 的数据。

      1. 谢谢博士大哥,不过我想问的是,这里用的什么”语法“,将三维降为二维?我百度没有查到,请大神给个该语法的链接,多谢。

  4. TWSC_interp = lut.ev(colats_rad,lons_rad).reshape(num_lats,num_lons)
    这里是不是应该是TWSC_interp = lut.ev(colats_rad,lons_rad).reshape(num_lats_gldas,num_lons_gldas)

    1. 确实是 TWSC_interp = lut.ev(colats_rad,lons_rad).reshape(num_lats,num_lons) 而不是 lut.ev(colats_rad,lons_rad).reshape(num_lats_gldas,num_lons_gldas)。num_lats,num_lons 实际上为
      lats = np.linspace(90,-89,180)
      lons = np.linspace(0,359,360)
      num_lats,num_lons = len(lats),len(lons)
      给大家带来的困扰深表歉意,最快下周末我会校准好统一的代码以及详细的注释。

    1. 我的也是,pyshtools一直出现错误“ERROR: Command errored out with exit status 1: python setup.py egg_info Check the logs for full command output.”
      请问你这个问题解决了没呢?感谢!

      1. 根据邓梓锋同学的反馈,win 系统下安装 pyshtools 包(conda 环境下 pip install pyshtools) 显然有较多麻烦,原因可能是由于 pyshtools 是 fortran 程序 和 python 程序的结合体,若系统没有 fortran 语言编译器(gfortran、f95 等),可能会出现错误提示。针对 win 系统下安装 pyshtools 困难的问题,邓同学在经过多次尝试和反复确认之后,得出以下结论:
        – win 系统下安装 Fortran 编译器无法让 pyshtools 成功安装;
        – pyshtools 需要在 Linux 平台才能稳定运行;
        – 推荐使用 VMware Workstation Pro 15.0 的虚拟机进行操作,虚拟机安装的系统推荐使用 Ubuntu 18.04 及以上的版本。该系统版本入门比较简单,VMware Workstation Pro对该系统版本的优化也做得比较好,且能较好地运行pycharm。

  5. 您好,想请教一下我运行该程序时出现mask_new未定义,为什么会出现这个问题呢?Python小白,勿喷

    1. 涉及到 mask_new 的代码行均可以删掉,这并不影响计算的结果。这是早期为了插值 mask 而设置的一个变量,新版本的代码已经将其去掉,近几天内我会对其更新。

  6. 您好,我在官网上下载的ETOPO5数据读取时出现了nicodeDecodeError错误,而从您提供的链接上下载的etopo5.dat文件并没有出现这个错误。您对etopo5.dat进行了其它处理吗?

    1. 全球地形起伏数据 ETOPO5 的空间分辨率为 5 弧分,目前该数据已被分辨率为 1 弧分的 ETOPO1 所取代。官网 NOAA 给出的 ETOPO5.DAT 文件确实是有问题(我这里打开的是乱码,可能是编码问题),从发布时间(2004-02-06)来看,貌似该文件已经不在维护了。链接给出的 etopo5.dat 数据是从 etopo5.nc 文件(记得该文件曾是原来的官方文件之一,可链接已经找不到了)转化而来的,对原始数据没作任何改动。详情可参考用全球地形起伏数据 ETOPO 作图

  7. 请问GLDAS在使用的时候只用了降水、蒸散、地表径流和地下径流,这些的单位应该都不是储量单位kg/m^2,具体换算时用到了month2second和month2threehours,这个具体是多少?
    很多论文都是采用雪水当量、土壤水含量和冠层水含量(也就是GLDAS中三个单位为kg/m^2的项),这两种方法的区别大吗?

    1. 降水、蒸散、地表径流和地下径流均属于过程量,单位是 kg/m^2/3h,因此需要乘以时间;雪水当量、土壤水含量和冠层水含量均属于状态量;已经验证过两种方法的差异,区别很小。

  8. 博士大哥,您好,请问计算陆地水储量变化需要GLDAS中哪些类型数据,0-10cm,10-40cm等等分层数据直接相加,再将所需要的不同类型的数据相加就可以吗?程序中提到了插值,插值是必须进行的步骤吗?

    1. 计算陆地水储量变化需要对 GLDAS 中雪水当量、冠层蓄水量和各层土壤含水量求和;用“降水-蒸发-地表径流-地下径流”计算的结果与用雪水当量、冠层蓄水量和土壤含水量求和的计算结果相当,两者仅有微小的差别;起初以为用“降水-蒸发-地表径流-地下径流”计算的结果能够包含地表水(湖泊)的量,但其实不然,经验证两种方法实际上是等价的,但求和的方法显然简单一些。另外,程序中的插值不是必须的,之所以插值是为了球谐展开方便。

  9. 博士你好!很感谢你的分享,方便了我这个小白学习GLDAS数据处理。
    我想问一下,GLDAS数据是600*1440格网,第一行是不是南纬60度~59.75度,第一列是西经180度~179.75度?就是说南北颠倒了是么?不知道对不对?如果是这样的话,总感觉怪怪的。
    再次感谢!~~

    1. 对于分辨率为 0.25 度的 GLDAS 产品,空间范围为 $179.875^{\circ}W \sim 179.875^{\circ}E$,$-59.875^{\circ}S \sim 89.875^{\circ}N$,网格的维度为 $600\times 1440$。数据的确是从南向北,从西向东开始,你可自行去验证,若发现有任何不对的地方请告知,以便我及时更正。

      1. 嗯嗯 好的。真心感谢你这么及时的回复!以后遇到问题,可能还要过来请教博士呢,谢谢 谢谢!

  10. 博士大哥,您好,很感谢您的分享,我有个小问题想请教一下您,GLDAS数据提取出来的各层土壤含水量,雪水当量,单位是kg/m^2,我看您在评论中的回答需要除以水的密度1000kg/m^3,这样就转化成单位为mm的等效水高,可是我在您写的程序中没有找到除以水密度的过程,而是直接相加,所以没有太懂,希望博士大哥能解决一下我的困惑,感谢!

    1. 假设雪水当量为 10 kg/m^2,将其转换成等效水高的过程为:10 [kg/m^2] / 1000 [kg/m^3] = 0.01 [m], 而 0.01 [m] = 10 [mm],换句话说,10 kg/m^2 的雪水当量等于 10 mm 的等效水高,100 kg/m^2 的土壤含水量等于 100 mm 的等效水高,为了简化程序的计算量,我在数值上省略了除以 1000,再乘以 1000 的步骤,而是直接相加,因而你没有找到除以水密度的过程。

  11. 博士大哥,您好,感谢您的分享,我有一些问题向您请教,我下载的是空间分辨率为1*1,时间分辨率为1个月的GLDAS NOAH模型数据,下载好的数据为nc4文件,我用matlab读取各变量,显示地表径流(Qs_acc)和地下径流(Qsb_acc)的单位为kg/m2,而程序中这两项为过程量,公式中是用降雨和降雪量的和减去土壤蒸发与植被蒸腾之和,再减去地表和地下径流的和计算陆地水储量变化,我看程序中只有降雨量没有降雪量,土壤蒸发和植被蒸腾也只用Evap_tavg一项替代,所以不太明白,向您请教。

    1. 参考 https://hydro1.gesdisc.eosdis.nasa.gov/data/GLDAS/GLDAS_NOAH10_M.2.1/doc/README_GLDAS2.pdf。

      该说明文件的 3.3 Data Interpretation 部分的有以下说明:

    2. Total precipitation rate is the sum of rain and snow precipitation rates. The forcing variable “Rainf_f_tavg” is the total precipitation rate whereas “Rainf_tavg” and “Snowf_tavg” are the liquid precipitation rate and frozen precipitation rate, respectively. 降水量为降雨量与降雪量之和。

      另外,Evap_tavg 是蒸散量,指植被及地面整体向大气输送的水汽总通量,主要包括植被蒸腾、土壤水分蒸发及截留降水或露水的蒸发。

发表评论

电子邮件地址不会被公开。 必填项已用*标注