全球陆地数据同化系统(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 天、一个月。GLDAS 数据可从 GES DISC(Goddard Earth Sciences Data and Information Services Center) 进行下载。

GLDAS 数据的下载

这里提供了两种下载方式,一种为直接在官网下载,另一种为使用 GGTOOLS 程序包进行下载。

方法 1

进入 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>” 为登陆密码。)
    • 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

方法 2

GGTOOLS(GRACE & GLDAS Tools) 程序包为 GLDAS 数据的下载提供了一个便捷通道。下载 2010 年 3 月份到 2015 年 11 月份的空间分辨率为 $1^{\circ}$、时间分辨率为 1 个月的 GLDAS NOAH 产品

下载 2002 年 4 月份到 2019 年 10 月份的空间分辨率为 $0.25^{\circ}$、时间分辨率为 1 个月的 GLDAS NOAH 产品

目前,GGTOOLS 仅支持时间分辨率为 1 个月 NOAH 产品。

GLDAS 数据的处理

GLDAS 数据的空间范围为南纬 $60^{\circ}$ 到北极,海洋区域均以 $-9999.0$ 进行了填充,数据默认为 netcdf 格式。关于 netcdf 格式文件的读取和写入请参见 NetCDF4.0 与 HDF5.0 文件的读取与生成

陆地水储量变化的估算

计算陆地水储量变化(Changes of Terrestrial Water Storage, TWSC) 的方式有两种,一种为对土壤水(0~200cm)、积雪融水、植被冠层水求和的经典方式,另一种是直接利用水平衡方程 $$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$ 个月的地表径流量与地下径流量之和。

利用空间分辨率为 $1^{\circ}$ 的 GLDAS 数据计算 2002 年 4 月到 2019 年 10 月的藏东南地区的陆地水储量变化。在读取 GLDAS 数据前,务必保证已经下载了数据,空间分辨率为 $1^{\circ}$ 的数据默认存放于目录 GLDAS/NOAH10/中,分辨率为 $0.25^{\circ}$ 的数据默认存放于目录 GLDAS/NOAH025/ 中。

输出结果为

设置区域的范围为 [82$^{\circ}$E, 106$^{\circ}$E, 21$^{\circ}$N, 39$^{\circ}$N],然后用经典的方法计算陆地水储量变化,其单位为 [mm w.e.](已将单位面积上的质量转换为等效水厚)。

输出结果为

若采用水平衡方程的方法计算陆地水储量变化,只需添加 mode = 'wbe' 即可,详情请参考 gldas.twsc_grid?

输出结果为

用经典的方法和水平衡方程法分别计算陆地水储量变化,结果表明两者的差别很小,读者可自行验证。

估计陆地水储量的变化率

默认使用迭代最小二乘法(Iterative Least Square method, ILSQM)拟合陆地水储量的变化率。其他的三种选择分别为最小二乘法(Least Square method, LSQM)、加权最小二乘法(Weighted Least Square method, WLSQM)、迭代加权最小二乘法(Iterative Weighted Least Square method, IWLSQM)。迭代与非迭代的区别在于迭代法会根据 3$\sigma$ 原则剔除异常点。对于 GLDAS 数据而言,没有发现任何异常数据,因此迭代和非迭代得到的结果一样。加权和非加权的区别在于数据的不确定度是否作为权重因子参加到最小二乘参数的拟合中。因为 GLDAS 数据没有提供任何不确定度信息,因此,加权法在这种情况下等价于非加权法。

输出结果为

更多详情请参考 twsc_grid.rate?

对藏东南地区的陆地水储量的变化率作图

输出结果为

陆地水储量变化的时间序列

设定研究区域

设定研究区域的方式有以下两种:

  1. 提供研究区域的边界文件 boundaries.txt
  2. 自定义椭圆作为研究区域的边界
方法 1

边界文件由两列组成,第一列为纬度,第二列为经度。

输出结果为

方法 2

创建以 (30$^{\circ}$N,95$^{\circ}$E) 为中心,半长轴为 4.2$^{\circ}$,半短轴为 3.2$^{\circ}$ 的椭圆

输出结果为

对研究区域作图仅需添加 polygons=study_area 即可

输出结果为

计算研究区域内陆地水储量的时间序列、变化率、研究区域的面积

输出结果为

对陆地水储量的时间序列作图

输出结果为

上图画出了陆地水储量的时间序列(加号)及其变化趋势(红色实线)、线性趋势(虚线),长期趋势(绿色实线)。长期趋势包含了线性趋势和年际变化,即扣除了周年项的变化趋势。

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

球谐展开陆地水储量到 96 阶。这一步是非必要的,除非与同阶的 GRACE 反演的质量变化做比较。

输出结果为

根据 Driscoll and Healy’s (1994) 样本定理,空间分辨率为 1$^{\circ}$ 的 GLDAS 数据实际上最大能展开到 89 阶。

计算尺度因子

若两个时间序列 $y_1$、$y_2$ 之间满足 $y_2(t)=ky_1(t)$,则 k 为尺度因子。在实际应用中,$y_2$ 代表“真实”的陆地水储量(未做任何处理的陆地水储量),$y_1$ 代表“视”陆地水储量(高斯滤波光滑后的陆地水储量)。

输出结果为

结果显示尺度因子接近 1,这说明研究区域选的足够大使其能够容纳扩散后(滤波光滑导致)的质量异常。若研究区域选的很小,或者质量异常位于研究区域的边缘地带,则尺度因子越偏离 1。尺度因子只对时间序列的线性趋势敏感,而对周年项不敏感。除此之外,不同的研究区域对应不同大小的尺度因子。使用尺度因子法对同阶的 GRACE 反演的质量变化进行信号泄漏改正,需要满足研究区域内 GRACE 反演的质量变化分布与 GLDAS 计算的陆地水储量变化的分布相似,否则使用该方法是不准确的。

尺度因子敏感性验证

设两个时间序列为 $y_1 = t + \sin(t) + 0.1$,$y_2 = 3t + 2\sin(t) – 0.1$,计算尺度因子

输出结果为

设两个时间序列为 $y_1 = t + \sin(t) + 0.1$, $y_2 = 0.1t + 2\sin(t) – 0.1$,计算尺度因子

输出结果为

结果表明尺度因子只对线性趋势敏感,而对振幅不敏感。


参考

  • 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.
  • GRACE & GLDAS Tools

84
Leave a Reply

avatar
25 Comment threads
59 Thread replies
2 Followers
 
Most reacted comment
Hottest comment thread
21 Comment authors
Name*lcx366Leafzhongletaomy Recent comment authors
  Subscribe  
newest oldest most voted
Notify of
JOJO
Guest
JOJO

爱分享的博士哥哥,大爱哟!

妮娜
Guest
妮娜

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

徐永明
Guest
徐永明

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

zhangchen
Guest
zhangchen

学长想问一下,windows系统下wget如何获取EARTHDATA的cookies?

丫丫
Guest
丫丫

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

赵广超
Guest
赵广超

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)

2
Guest
2

pyshtools这个库一直装不好呀

Wendy
Guest
Wendy

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

Wendy
Guest
Wendy

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

lucialin
Guest
lucialin

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

王涛
Guest
王涛

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

sugaror
Guest
sugaror

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

Wt
Guest
Wt

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

Wt
Guest
Wt

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

Chloe
Guest
Chloe

lcx366老师你好,非常感谢你的分享,收益很大,但是想请教几个问题,
1. GLDAS数据从nc提取并计算得到降水,径流,陆地水储量等变量后可以直接用于其他分析和计算吗?
2. 后面水储量插值和球谐展开是为了和GRACE计算得到的TWSC结果进行匹配吗?
期待你的回复,再次感谢!

abcat
Guest
abcat

您好,感谢分享!想请教几个问题:
1.使用GLDAS数据前必须进行插值和球谐展开吗?比如要得到流域尺度的月土壤水含量、蒸散 这些是否需要进行如上的处理呢?
2.从GLDAS网格化数据中直接提取的数值是什么数据呢?

Viva La Vida
Guest
Viva La Vida

博士您好,请问ggtools您是在什么编辑器上使用的,ubuntu上的pycharm,使用ggtools,有很多处UnboundLocalError:local variable ‘files’ referenced before assignment. 有很多函数变量不太会申明全局变量或者局部变量,比如用xarray读取的datum = xr.open_dataset(filename), datum我该用矩阵定义吗? 期待您的回复。

Viva La Vida
Guest
Viva La Vida

from ggtools.gg import read_gldas
gldas = read_gldas(‘2002-04′,’2019-10′)
print(gldas)
print(gldas.lats,’\n’)
print(gldas.lons)

Traceback(most recent call last):
File”/home/Li/PycharmProjects/GLDAS”,line2,in
gldas = gldas = read_gldas(‘2002-04′,’2019-10’)
File”/home/Li/.local/lib/python3.6/site-packages/ggtools/gg/gldas_utils.py”,line 153, in read_gldas files = np.sort(files)
UnboundLocalError:local variable ‘files’ referenced before assignment.

我觉得是函数里面的变量python当做了局部变量,就去源代码改了一下,改完之后后面datum也是这个原因
UnboundLocalError:local variable ‘datum’ referenced before assignment.

Viva La Vida
Guest
Viva La Vida

这个下载的是NOAH1.0数据吗

Viva La Vida
Guest
Viva La Vida

试着按照ggtools的程序下载了数据,下载的是GLDAS NOAH2.1版本的数据,在Ubuntu中pycharm确实使用不了,老是出现其他的错误,然后下载了ggtools的包,把源代码放在一个文件夹里。改了源代码引用路径,然后用什么函数,就从py文件调用函数就可以用了。应该是我之前下载的很多包与这个冲突了吧,按照我的水平也实在给不了一个确切的答案。总算是解决了,楼主的代码牛逼!
我按Reply还是会另外成为一楼,一个问题出现了好几层楼,看的不舒服,就把之前的删了吧。😂
再请问一下,我们论文中一般使用的就是TWSC,画趋势图也是直接用的TWSC,变化率代表的是什么意思?最后在画趋势图中的给出的TWSC【GT】代表的是取平均化的每月TWS变化量吗?我如何把他转化为等效水高?期待您的回复感谢博士。

taomy
Guest
taomy

博士,您好!利用您上述代码处理GLDAS数据时,第一步读取GLDAS运行代码出现错误,具体如下:
from ggtools.gg import read_gldas
# 读取 GLDAS 数据
gldas = read_gldas(‘2002-04′,’2017-06′)
print(gldas,’\n’)
print(gldas.lats,’\n’)
print(gldas.lons)

提示错误:
KeyError: [, (‘/home/taomy/GLDAS/NOAH10/GLDAS_NOAH10_M.A200204.021.nc4’,), ‘r’, ((‘clobber’, True), (‘diskless’, False), (‘format’, ‘NETCDF4’), (‘persist’, False))]

During handling of the above exception, another exception occurred:

OSError: [Errno -51] NetCDF: Unknown file format: b’/home/taomy/GLDAS/NOAH10/GLDAS_NOAH10_M.A200204.021.nc4′

自己也搜索了网上的解决方案还是没有解决,希望博士大哥帮助一哈!提前感谢!

Name*
Guest
Name*

博士师兄您好!我也是使用了您的代码,在运行读取数据的时候出现了错误

代码:
from ggtools.gg import read_gldas
# 读取 GLDAS 数据
gldas = read_gldas(‘2002-04′,’2019-10′)
print(gldas)
print(gldas.lats,’\n’)
print(gldas.lons)

错误:
File “C:\ProgramData\Anaconda3\lib\site-packages\ggtools\gg\gldas_utils.py”, line 151, in read_gldas
files = np.sort(files)

UnboundLocalError: local variable ‘files’ referenced before assignment

我也查看了”GLDAS\NOAH025″路径下.nc4文件已下载成功,这个错误该如何解决?自己刚接触这一领域,冒昧向您请教,期待您的回复!非常感谢!

zhongle
Guest
zhongle

博士您好,我用ggtools下载GLDAS数据时用的代码是:
from ggtools.gg import gldas_download
username,password = ‘。。。。。。’,’。。。。。。’
start_date,end_date = ‘2002-04′,’2019-10′
gldas_download(username,password,start_date,end_date)
出现如下情况:
C:\ProgramData\Miniconda3\python.exe E:/GGTOOLS-master/ggtools/tryfive.py
Traceback (most recent call last):
File “E:/GGTOOLS-master/ggtools/tryfive.py”, line 11, in
gldas_download(username,password,start_date,end_date)
File “C:\ProgramData\Miniconda3\lib\site-packages\ggtools\gg\gldas_utils.py”, line 53, in gldas_download
if not path.exists(home+’/.netrc’):
TypeError: unsupported operand type(s) for +: ‘NoneType’ and ‘str’
我是在windows上的pycharm上运行的,出现了这样的情况,不太会处理,想请教一下博士大哥

zhongle
Guest
zhongle

博士您好,在已经做出GRACE和GLDAS空间分布图后,想通过二者相减模拟地下水储量的空间分布图,ggtools可以做到吗?应该添加什么代码呢?

Leaf
Guest
Leaf

老师您好,请问用空间分辨率为0.25度的GLDAS数据和分辨率为0.5度的GRACE数据进行比对,相较于处理空间分辨率为1度的数据,代码应该作何变动呢?