A package to handle the GRACE and GRACEFO GSM data and GLDAS grid data
Project description
Welcome to the GRACE & GLDAS Tools
The ggtools package is an archive of scientific routines that can be used to handle the GRACE(Gravity Recovery and Climate Experiment) and GRACEFO(Followon) GSM data(RL06 Level2 monthly solutions) and GLDAS grid data.
How to Install
GGTOOLS can be installed with the following two steps
conda install cartopy netcdf4 h5py
pip install ggtools
How to Upgrade
GGTOOLS can be upgraded to the latest version with pip install ggtools upgrade
.
How to use
GRACE
Download static gravity models
Download static gravity models GGM05C and EIGEN6C4 from ICGEM(International Centre for Global Earth Models). This is not a necessary step unless you want to prepare for removing the background(reference) gravity field from the GRACE GSM data later.
from ggtools.gg import static_download for static_model in ['GGM05C','EIGEN6C4']: gravity_model = static_download(static_model) print(gravity_model)
static_models/GGM05C.gfc
static_models/EIGEN6C4.gfc
For more details, please refer to static_download?
.
Check the time coverage of the GRACE GSM data
GRACE GSM data is expressed as a form of Spherical Harmonic Coefficients(SHCs) or Stokes coefficients. Before downloading the GRACE GSM data from ISDC(Information System and Data Center), it's suggested to view the time intervals for all feasible GRACE GSM data so far published by CSR(University of Texas Center for Space Research), GFZ(German Research Centre for Geosciences), and JPL(Jet Propulsion Laboratory). This is not a necessary step, but it can help to understand the outline of the data.
from ggtools.gg import print_gsm_date_coverage for source in ['CSR','GFZ','JPL']: print_gsm_date_coverage(source,96) print_gsm_date_coverage(source,60)
GSM data for CSR_RL06_96 is available from 20020405 to 20191031
GSM data for CSR_RL06_60 is available from 20020405 to 20191031
GSM data for GFZ_RL06_96 is available from 20020405 to 20191031
GSM data for GFZ_RL06_60 is available from 20020405 to 20191031
GSM data for JPL_RL06_96 is available from 20020404 to 20191031
GSM data for JPL_RL06_60 is available from 20020404 to 20191031
For more details, please refer to print_gsm_date_coverage?
.
Download GRACE GSM data
Download all feasible GRACE GSM data published by CSR, GFZ, and JPL so far.
from ggtools.gg import gsm_download for source in ['CSR','GFZ','JPL']: gsm_download(source,96) gsm_download(source,60)
Downloading ... GSM2_20192742019304_GRFO_UTCSR_BB01_0600.gz ... 226 Transfer complete Downloading ... GSM2_20192742019304_GRFO_UTCSR_BA01_0600.gz ... 226 Transfer complete Downloading ... GSM2_20192742019304_GRFO_JPLEM_BB01_0600.gz ... 226 Transfer complete Downloading ... GSM2_20192742019304_GRFO_JPLEM_BA01_0600.gz ... 226 Transfer complete
For more details, please refer to gsm_download?
.
Download SLR C20 data
Download SLR C20 data(RL06 monthly solutions) published by CSR.
from ggtools.gg import slr_c20_download slr_c20_download()
Downloading ... TN11_C20_SLR_RL06.txt ... 226 Transfer complete
For more details, please refer to slr_c20_download?
.
Read GRACE GSM data
Read all 96thdegree GRACE GSM data downloaded previously.
from ggtools.gg import read_gsm csr_gsm = read_gsm('CSR',96) gfz_gsm = read_gsm('GFZ',96) jpl_gsm = read_gsm('JPL',96) # basic information on GRACE GSM data print(csr_gsm,'\n') print(gfz_gsm,'\n') print(jpl_gsm)
title = GRACE & GRACEFO Geopotential Coefficients CSR RL06
max_degree = 96
max_order = 96
degree_order = 96
normalization = fully normalized
institution = UTAUSTIN/CSR
processing_level = 2
product_version = RL06
time_coverage_start = 200204
time_coverage_end = 201910
solution_counts = 178
total_month_counts = 211
missing_month_counts = 33
title = GRACE & GRACEFO Geopotential GSM Coefficients GFZ RL06
max_degree = 96
max_order = 96
degree_order = 96
normalization = fully normalized
institution = GFZ German Research Centre for Geosciences
processing_level = 2
product_version = 6.0
time_coverage_start = 200204
time_coverage_end = 201910
solution_counts = 178
total_month_counts = 211
missing_month_counts = 33
title = GRACE & GRACEFO Geopotential Coefficients JPL RL06
max_degree = 96
max_order = 96
degree_order = 96
normalization = fully normalized
institution = NASA/JPL
processing_level = 2
product_version = 6.0
time_coverage_start = 200204
time_coverage_end = 201910
solution_counts = 169
total_month_counts = 211
missing_month_counts = 42
For more details, please refer to read_gsm?
and csr_gsm?
.
print(csr_gsm.summary) print(csr_gsm.background_gravity) print(csr_gsm.earth_gravity_param) print(csr_gsm.permanent_tide) print(csr_gsm.mean_equator_radius) print(csr_gsm.missing_month,'\n') print(gfz_gsm.background_gravity) print(gfz_gsm.permanent_tide) print(gfz_gsm.mean_equator_radius,'\n') print(jpl_gsm.missing_month) print(jpl_gsm.shc[10],'\n') # geopotential coefficients for the 10th monthly solution print(jpl_gsm.shc_std[10]) # standard deviation in geopotential coefficients
Spherical harmonic coefficients representing an estimate of the mean gravity field of Earth during the specified timespan derived from GRACE & GRACEFO mission measurements. These coefficients represent the full magnitude of land hydrology, ice, and solid Earth processes. Further, they represent atmospheric and oceanic processes not captured in the accompanying GAC product. The 0th and 1st degree terms are excluded from CSR level2.
GGM05C
3.9860044150E+14 m3/s2
inclusive
6.3781363000E+06 m
['200206' '200207' '200306' '201101' '201106' '201205' '201210'
'201303' '201308' '201309' '201402' '201407' '201412' '201507'
'201510' '201511' '201604' '201609' '201610' '201702' '201707'
'201708' '201709' '201710' '201711' '201712' '201801' '201802'
'201803' '201804' '201805' '201808' '201809']
EIGEN6C4
exclusive
6.3781364600E+06 m
['200206' '200207' '200306' '200407' '200408' '200409' '200410'
'201101' '201106' '201204' '201205' '201206' '201207' '201210'
'201303' '201308' '201309' '201402' '201407' '201412' '201501'
'201502' '201507' '201510' '201511' '201604' '201609' '201610'
'201702' '201707' '201708' '201709' '201710' '201711' '201712'
'201801' '201802' '201803' '201804' '201805' '201808' '201809']
[[[ 1.00000000e+00 0.00000000e+00 0.00000000e+00 ... 0.00000000e+00
0.00000000e+00 0.00000000e+00]
[ 0.00000000e+00 0.00000000e+00 0.00000000e+00 ... 0.00000000e+00
0.00000000e+00 0.00000000e+00]
[4.84169287e04 2.29490037e10 2.43933401e06 ... 0.00000000e+00
0.00000000e+00 0.00000000e+00]
...
[ 9.37921024e11 8.38607258e11 2.33486569e10 ... 2.26383594e09
0.00000000e+00 0.00000000e+00]
[ 5.33299728e10 1.81153117e09 2.78305221e10 ... 2.61977837e09
2.59548687e09 0.00000000e+00]
[1.25276141e09 5.11732945e10 1.26614370e09 ... 1.14454762e09
1.62194507e09 2.14363128e09]]
[[ 0.00000000e+00 0.00000000e+00 0.00000000e+00 ... 0.00000000e+00
0.00000000e+00 0.00000000e+00]
[ 0.00000000e+00 0.00000000e+00 0.00000000e+00 ... 0.00000000e+00
0.00000000e+00 0.00000000e+00]
[ 0.00000000e+00 1.38511555e09 1.40030144e06 ... 0.00000000e+00
0.00000000e+00 0.00000000e+00]
...
[ 0.00000000e+00 2.29543385e09 1.35606111e10 ... 1.57903333e09
0.00000000e+00 0.00000000e+00]
[ 0.00000000e+00 6.34755345e10 2.53978387e10 ... 1.20141984e09
3.45902449e10 0.00000000e+00]
[ 0.00000000e+00 1.63243506e09 7.04276590e11 ... 2.73715185e09
4.34383583e10 1.63871799e09]]]
[[[0.0000e+00 0.0000e+00 0.0000e+00 ... 0.0000e+00 0.0000e+00 0.0000e+00]
[0.0000e+00 0.0000e+00 0.0000e+00 ... 0.0000e+00 0.0000e+00 0.0000e+00]
[6.4545e12 3.4606e12 1.3888e12 ... 0.0000e+00 0.0000e+00 0.0000e+00]
...
[1.5314e11 1.5326e11 1.5305e11 ... 4.5916e11 0.0000e+00 0.0000e+00]
[1.4790e11 1.4741e11 1.4724e11 ... 4.1441e11 4.6175e11 0.0000e+00]
[1.5805e11 1.5750e11 1.5728e11 ... 1.1087e10 3.9720e11 4.0593e11]]
[[0.0000e+00 0.0000e+00 0.0000e+00 ... 0.0000e+00 0.0000e+00 0.0000e+00]
[0.0000e+00 0.0000e+00 0.0000e+00 ... 0.0000e+00 0.0000e+00 0.0000e+00]
[0.0000e+00 3.0349e12 1.4492e12 ... 0.0000e+00 0.0000e+00 0.0000e+00]
...
[0.0000e+00 1.5277e11 1.5313e11 ... 4.5258e11 0.0000e+00 0.0000e+00]
[0.0000e+00 1.4695e11 1.4742e11 ... 4.1784e11 4.6037e11 0.0000e+00]
[0.0000e+00 1.5703e11 1.5747e11 ... 1.1097e10 3.9901e11 4.0941e11]]]
Note that CSR and GFZ use different background(reference) gravity models. One is GGM05C and the other is EIGEN6C4. These two gravity models define the same Earth gravity constant, but slightly different Earth's mean equator radius. In addition, both GGM05C and the monthly solutions released by CSR include the permanent tide. In contrast, neither EIGEN6C4 nor the monthly solutions released by GFZ include this item. Therefore, it is necessary to remove the background field or the average field from the monthly solutions.
Deaverage GRACE GSM data
Remove the average field from the monthly solutions.
csr_gsm_deaverage = csr_gsm.deaverage() gfz_gsm_deaverage = gfz_gsm.deaverage() jpl_gsm_deaverage = jpl_gsm.deaverage() print(csr_gsm_deaverage,'\n') print(csr_gsm_deaverage.background_gravity,'\n') print(csr_gsm_deaverage.shc[10],'\n') print(csr_gsm_deaverage.shc_std[10])
title = Deaveraged GRACE & GRACEFO Geopotential Coefficients CSR RL06
max_degree = 96
max_order = 96
degree_order = 96
normalization = fully normalized
institution = UTAUSTIN/CSR
processing_level = 2
product_version = RL06
time_coverage_start = 200204
time_coverage_end = 201910
solution_counts = 178
total_month_counts = 211
missing_month_counts = 33
Average of monthly solutions
[[[ 0.00000000e+00 0.00000000e+00 0.00000000e+00 ... 0.00000000e+00
0.00000000e+00 0.00000000e+00]
[ 0.00000000e+00 0.00000000e+00 0.00000000e+00 ... 0.00000000e+00
0.00000000e+00 0.00000000e+00]
[ 1.84422663e10 9.25679556e11 3.83500305e11 ... 0.00000000e+00
0.00000000e+00 0.00000000e+00]
...
[ 6.10660859e12 2.02340281e13 7.35534442e12 ... 4.30518176e11
0.00000000e+00 0.00000000e+00]
[3.59387039e12 6.02718201e12 1.68845289e12 ... 2.47010217e11
9.95521155e13 0.00000000e+00]
[2.09598368e12 4.38205266e12 3.08504132e12 ... 1.06180191e10
4.22540682e12 9.45004334e12]]
[[ 0.00000000e+00 0.00000000e+00 0.00000000e+00 ... 0.00000000e+00
0.00000000e+00 0.00000000e+00]
[ 0.00000000e+00 0.00000000e+00 0.00000000e+00 ... 0.00000000e+00
0.00000000e+00 0.00000000e+00]
[ 0.00000000e+00 6.95504729e11 4.66008820e12 ... 0.00000000e+00
0.00000000e+00 0.00000000e+00]
...
[ 0.00000000e+00 2.02779108e12 1.28132376e12 ... 2.70660778e11
0.00000000e+00 0.00000000e+00]
[ 0.00000000e+00 6.34766059e12 1.40266107e12 ... 1.63210355e11
1.74411386e12 0.00000000e+00]
[ 0.00000000e+00 2.98368541e12 1.09980733e11 ... 7.87717175e11
7.85024577e12 4.73324005e12]]]
[[[0.000e+00 0.000e+00 0.000e+00 ... 0.000e+00 0.000e+00 0.000e+00]
[0.000e+00 0.000e+00 0.000e+00 ... 0.000e+00 0.000e+00 0.000e+00]
[3.374e13 5.154e13 6.335e13 ... 0.000e+00 0.000e+00 0.000e+00]
...
[1.033e11 1.034e11 1.032e11 ... 1.717e11 0.000e+00 0.000e+00]
[9.930e12 9.886e12 9.858e12 ... 6.068e12 5.054e12 0.000e+00]
[1.056e11 1.054e11 1.053e11 ... 4.096e11 6.989e12 7.420e12]]
[[0.000e+00 0.000e+00 0.000e+00 ... 0.000e+00 0.000e+00 0.000e+00]
[0.000e+00 0.000e+00 0.000e+00 ... 0.000e+00 0.000e+00 0.000e+00]
[0.000e+00 4.225e13 7.317e13 ... 0.000e+00 0.000e+00 0.000e+00]
...
[0.000e+00 1.028e11 1.035e11 ... 1.713e11 0.000e+00 0.000e+00]
[0.000e+00 9.790e12 9.923e12 ... 6.232e12 5.132e12 0.000e+00]
[0.000e+00 1.047e11 1.060e11 ... 4.089e11 6.981e12 7.459e12]]]
For more details, please refer to csr_gsm.deaverage?
.
Debackground GRACE GSM data
Remove the background field from the monthly solutions. This step is not necessary, but it can be used as a way to verify the deaveraged monthly solutions.
csr_gsm_debackground = csr_gsm.debackground() print(csr_gsm_debackground,'\n') print(csr_gsm_debackground.background_gravity,'\n') print(csr_gsm_debackground.shc[10])
title = Debackgrounded GRACE & GRACEFO Geopotential Coefficients CSR RL06
max_degree = 96
max_order = 96
degree_order = 96
normalization = fully normalized
institution = UTAUSTIN/CSR
processing_level = 2
product_version = RL06
time_coverage_start = 200204
time_coverage_end = 201910
solution_counts = 178
total_month_counts = 211
missing_month_counts = 33
GGM05C
[[[ 0.00000000e+00 0.00000000e+00 0.00000000e+00 ... 0.00000000e+00
0.00000000e+00 0.00000000e+00]
[ 0.00000000e+00 0.00000000e+00 0.00000000e+00 ... 0.00000000e+00
0.00000000e+00 0.00000000e+00]
[ 1.74518000e10 5.89105312e11 3.82188300e11 ... 0.00000000e+00
0.00000000e+00 0.00000000e+00]
...
[ 6.74305609e12 1.42961946e12 6.05156366e12 ... 3.27285439e11
0.00000000e+00 0.00000000e+00]
[8.92741659e12 6.19449953e12 2.38098621e12 ... 1.06105863e10
6.68478535e12 0.00000000e+00]
[7.51046521e13 6.47684322e12 1.89032517e13 ... 1.01848744e10
5.33312012e11 6.95728099e12]]
[[ 0.00000000e+00 0.00000000e+00 0.00000000e+00 ... 0.00000000e+00
0.00000000e+00 0.00000000e+00]
[ 0.00000000e+00 0.00000000e+00 0.00000000e+00 ... 0.00000000e+00
0.00000000e+00 0.00000000e+00]
[ 0.00000000e+00 2.56846039e11 1.36937040e11 ... 0.00000000e+00
0.00000000e+00 0.00000000e+00]
...
[ 0.00000000e+00 1.45661890e12 1.10232503e13 ... 1.36114572e11
0.00000000e+00 0.00000000e+00]
[ 0.00000000e+00 8.23449794e12 5.41604372e13 ... 1.59533178e10
6.78363971e12 0.00000000e+00]
[ 0.00000000e+00 4.55778346e12 9.57638370e12 ... 6.04274164e11
2.73714269e11 5.54947649e12]]]
For more details, please refer to csr_gsm.debackground?
.
Read SLR C20 data
Read all feasible SLR C20 RL06 monthly solutions so far.
from ggtools.gg import read_slr_c20 slr_c20 = read_slr_c20() print(slr_c20)
title = Monthly estimates of C20 from 5 SLR satellites based on GRACE RL06 models.
normalization = fully normalized
institution = Center for Space Research, The University of Texas at Austin
product_version = RL06
time_coverage_start = 200204
time_coverage_end = 201910
solution_counts = 178
total_month_counts = 211
missing_month_counts = 33
For more details, please refer to read_slr_c20?
and slr_c20?
.
print(slr_c20.summary,'\n') print(slr_c20.missing_month,'\n') print(slr_c20.date_issued,'\n') print(slr_c20.background_gravity) print(slr_c20.mean_c20,'\n') # c20 in GGM05C print(slr_c20.c20,'\n') print(slr_c20.c20_std,'\n') print(slr_c20.c20_demean) # c20  c20 in GGM05C
As a convenience to users who wish to use a replacement value for C20, a monthly C20 estimate time series is provided. These estimates are obtained from the analysis of Satellite Laser Ranging (SLR) data to five geodetic satellites: LAGEOS1 and 2, Starlette, Stella and Ajisai. The background gravity satellites model used in the SLR analysis is consistent with the GRACE Release06 processing, including the use of the same AtmosphereOcean Dealiasing product.
['200206' '200207' '200306' '201101' '201106' '201205' '201210'
'201303' '201308' '201309' '201402' '201407' '201412' '201506'
'201510' '201511' '201604' '201609' '201610' '201702' '201707'
'201708' '201709' '201710' '201711' '201712' '201801' '201802'
'201803' '201804' '201805' '201808' '201809']
Created November 13, 2019  last month reported is October 2019.
GGM05C
0.00048416945732
[0.00048417 0.00048417 0.00048417 0.00048417 0.00048417 0.00048417
0.00048417 0.00048417 0.00048417 0.00048417 0.00048417 0.00048417
...
0.00048417 0.00048417 0.00048417 0.00048417 0.00048417 0.00048417
0.00048417 0.00048417 0.00048417 0.00048417]
[3.789e11 3.141e11 3.164e11 3.628e11 3.292e11 3.440e11 3.556e11
3.660e11 3.124e11 2.755e11 3.646e11 3.141e11 3.368e11 2.771e11
...
4.315e11 4.219e11 3.312e11 3.169e11 3.659e11 4.319e11 4.160e11
5.326e11 3.543e11 3.717e11]
[ 1.1088e10 7.3380e11 9.6680e11 7.9420e11 6.4740e11 3.8280e11
6.1700e12 4.6160e11 3.5060e11 7.6130e11 5.6550e11 4.5380e11
...
1.0959e10 7.4410e11 8.7660e11 3.0400e11 4.4480e11 9.2330e11
2.3736e10 2.3425e10 2.5003e10 2.5976e10]
Deaverage the SLR C20 monthly solutions
Remove the average C20 from the SLR C20 monthly solutions.
slr_c20_deaverage = slr_c20.deaverage() print(slr_c20_deaverage,'\n') print(slr_c20_deaverage.c20,'\n') print(slr_c20_deaverage.mean_c20) # average of c20
title = Deaveraged Monthly estimates of C20 from 5 SLR satellites based on GRACE RL06 models.
normalization = fully normalized
institution = Center for Space Research, The University of Texas at Austin
product_version = RL06
time_coverage_start = 200204
time_coverage_end = 201910
solution_counts = 178
total_month_counts = 211
missing_month_counts = 33
[ 1.60614513e10 1.23116713e10 4.69406871e11 2.96770871e11
1.50012871e11 1.14603129e11 5.59060129e11 9.58994129e11
...
5.25801288e12 4.25913871e11 1.87620087e10 1.84511087e10
2.00289787e10 2.10022787e10]
0.0004841695070590129
For more details, please refer to slr_c20.deaverage?
.
Combine GRACE GSM data
You may calculate the arithmetic mean of the deaveraged monthly solutions from CSR, GFZ, and JPL, or, you may first calculate the arithmetic mean of the monthly solutions from CSR, GFZ, and JPL to get the combined solutions, then deaverage the combined solutions.
Method 1: deaverage and then calculate arithmetic mean
from ggtools.gg import gsm_average gsm_deaverage_comb = gsm_average([csr_gsm_deaverage,gfz_gsm_deaverage,jpl_gsm_deaverage]) print(gsm_deaverage_comb)
title = Combined Deaveraged GRACE & GRACEFO Geopotential Coefficients RL06
max_degree = 96
max_order = 96
degree_order = 96
normalization = fully normalized
institution = UTAUSTIN/CSR, GFZ German Research Centre for Geosciences, NASA/JPL
processing_level = 2
product_version = RL06
time_coverage_start = 200204
time_coverage_end = 201910
solution_counts = 178
total_month_counts = 211
missing_month_counts = 33
For more details, please refer to gsm_average?
.
Method 2: calculate arithmetic mean and then deaverage
gsm_comb = gsm_average([csr_gsm,gfz_gsm,jpl_gsm]) gsm_comb_deaverage = gsm_comb.deaverage() print(gsm_comb_deaverage)
title = Deaveraged Combined GRACE & GRACEFO Geopotential Coefficients RL06
max_degree = 96
max_order = 96
degree_order = 96
normalization = fully normalized
institution = UTAUSTIN/CSR, GFZ German Research Centre for Geosciences, NASA/JPL
processing_level = 2
product_version = RL06
time_coverage_start = 200204
time_coverage_end = 201910
solution_counts = 178
total_month_counts = 211
missing_month_counts = 33
This method is only applicable to the case where CSR, GFZ, and JPL have a same number of monthly solutions. Generally, CSR, GFZ, and JPL have a different number of monthly solutions, and it is recommended to use the first method to obtain the combined monthly solutions.
Replace C20 in GRACE GSM data with SLR C20
The time interval of the SLR C20 data should be consistent with that of GRACE GSM data. You may replace C20 in deaveraged GRACE GSM data with the deaveraged SLR C20, or, you may first replace C20 in GRACE GSM data with SLR C20, then deaverage the GRACE GSM data.
Reconfirm the time coverage of the SLR C20 data and the GRACE GSM data.
print('SLR C20 solutions start month: ',slr_c20.time_coverage_start) print('SLR C20 solutions end month: ',slr_c20.time_coverage_end) print('Combined deaveraged GRACE GSM solutions start month: ',gsm_deaverage_comb.time_coverage_start) print('Combined deaveraged GRACE GSM solutions end month: ',gsm_deaverage_comb.time_coverage_end)
SLR C20 solutions start month: 200204
SLR C20 solutions end month: 201910
Combined deaveraged GRACE GSM solutions start month: 200204
Combined deaveraged GRACE GSM solutions end month: 201910
Generally, the latest GRACE GSM monthly solution is released one month later than the SLR C20 solution, so the number of monthly solutions for SLR C20 is almost always one more than that for GRACE GSM.
Method 1: deaverage SLR C20 then replace C20
gsm_deaverage_comb_r = gsm_deaverage_comb.replace_slr_c20(slr_c20_deaverage) print(gsm_deaverage_comb_r,'\n') print(gsm_deaverage_comb_r.shc[:,0,2,0]) # C20
title = Combined Deaveraged GRACE & GRACEFO Geopotential Coefficients RL06 with C20 replaced by the SLR measurements
max_degree = 96
max_order = 96
degree_order = 96
normalization = fully normalized
institution = UTAUSTIN/CSR, GFZ German Research Centre for Geosciences, NASA/JPL
processing_level = 2
product_version = RL06
time_coverage_start = 200204
time_coverage_end = 201910
solution_counts = 178
total_month_counts = 211
missing_month_counts = 33
[ 1.60614513e10 1.23116713e10 4.69406871e11 2.96770871e11
1.50012871e11 1.14603129e11 5.59060129e11 9.58994129e11
...
5.25801288e12 4.25913871e11 1.87620087e10 1.84511087e10
2.00289787e10 2.10022787e10]
For more details, please refer to gsm_deaverage_comb.replace_slr_c20?
.
Method 2: replace C20 then deaverage
gsm_comb_r_deaverage = gsm_comb.replace_slr_c20(slr_c20).deaverage() print(gsm_comb_r_deaverage,'\n') print(gsm_comb_r_deaverage.shc[:,0,2,0]) # C20
title = Deaveraged Combined GRACE & GRACEFO Geopotential Coefficients RL06 with C20 replaced by the SLR measurements
max_degree = 96
max_order = 96
degree_order = 96
normalization = fully normalized
institution = UTAUSTIN/CSR, GFZ German Research Centre for Geosciences, NASA/JPL
processing_level = 2
product_version = RL06
time_coverage_start = 200204
time_coverage_end = 201910
solution_counts = 178
total_month_counts = 211
missing_month_counts = 33
[ 1.60614513e10 1.23116713e10 4.69406874e11 2.96770874e11
1.50012874e11 1.14603126e11 5.59060126e11 9.58994125e11
...
5.25801256e12 4.25913874e11 1.87620087e10 1.84511087e10
2.00289787e10 2.10022787e10]
DDK filtering
As a nonisotropic filter, DDK filter is designed to attenuate noise described as striping patterns in GRACE GSM data. There are eight kinds of DDK filters to choose from according to the smoothing strength. From DDK1 to DDK8, the smoothing strength gradually weakens.
Smooth noise using DDK5 filter.
gsm_r_ddk5 = gsm_deaverage_comb_r.filter_ddk('DDK5') print(gsm_r_ddk5) print(gsm_r_ddk5.summary)
title = DDK5 filtered Combined Deaveraged GRACE & GRACEFO Geopotential Coefficients RL06 with C20 replaced by the SLR measurements
max_degree = 96
max_order = 96
degree_order = 96
normalization = fully normalized
institution = UTAUSTIN/CSR, GFZ German Research Centre for Geosciences, NASA/JPL
processing_level = 2
product_version = RL06
time_coverage_start = 200204
time_coverage_end = 201910
solution_counts = 178
total_month_counts = 211
missing_month_counts = 33
Spherical harmonic coefficients representing an estimate of the mean gravity field of Earth during the specified timespan derived from GRACE & GRACEFO mission measurements. These coefficients represent the full magnitude of land hydrology, ice, and solid Earth processes. Further, they represent atmospheric and oceanic processes not captured in the accompanying GAC product. Note that the 2nddegree terms have been replaced with the values from SLR C20. Also note that C20 from SLR also experienced the DDK5 filtering.
For more details, please refer to gsm_deaverage_comb_r.filter_ddk?
.
Note that C20 from SLR also experienced the DDK5 filtering. Alternatively, you may replace C20 after performing the DDK filtering.
gsm_ddk5_r = gsm_deaverage_comb.filter_ddk('DDK5').replace_slr_c20(slr_c20_deaverage)
Substract the raw C20 from DDK5 filtered C20.
print(gsm_r_ddk5.shc[:,0,2,0]  gsm_ddk5_r.shc[:,0,2,0])
[3.07271480e16 5.33381926e17 8.07660043e17 3.21390821e16
4.92343516e17 3.54627589e17 3.12772822e17 4.76970460e17
...
5.22234499e17 5.10455550e17 4.61722884e17 3.88862298e17
4.15231980e17 2.04760123e17 6.78578147e17 4.17457051e17
8.27795363e17 8.68120470e17]
The difference between the two is much smaller than the standard deviation of C20(magnitude of $10^{11}$), so it doesn't matter in using either method.
Gaussian filtering
As an isotropic filter, the Gaussian filter is used to attenuate the striping noise in GRACE GSM data. The smoothing effect can be changed by adjusting the filter radius. Larger filter radius can significantly suppress highdegree noise, but at the same time attenuates the signals, especially the highdegree signal with greater uncertainty. If the DDK filtering has been applied, you don't need to use Gaussian filtering to process it again.
Smooth noise using Gaussian filter with 160km filter radius.
gsm_r_gau = gsm_deaverage_comb_r.filter_gaussian(160) print(gsm_r_gau) print(gsm_r_gau.summary)
title = Gaussian filtered Combined Deaveraged GRACE & GRACEFO Geopotential Coefficients RL06 with C20 replaced by the SLR measurements
max_degree = 96
max_order = 96
degree_order = 96
normalization = fully normalized
institution = UTAUSTIN/CSR, GFZ German Research Centre for Geosciences, NASA/JPL
processing_level = 2
product_version = RL06
time_coverage_start = 200204
time_coverage_end = 201910
solution_counts = 178
total_month_counts = 211
missing_month_counts = 33
Spherical harmonic coefficients representing an estimate of the mean gravity field of Earth during the specified timespan derived from GRACE & GRACEFO mission measurements. These coefficients represent the full magnitude of land hydrology, ice, and solid Earth processes. Further, they represent atmospheric and oceanic processes not captured in the accompanying GAC product. Note that the 2nddegree terms have been replaced with the values from SLR C20. Also note that C20 from SLR also experienced the Gaussian filtering.
For more details, please refer to gsm_deaverage_comb_r.filter_gaussian?
.
Note that C20 from SLR also experienced the Gaussian filtering. Alternatively, you may replace C20 after performing the Gaussian filtering.
gsm_gau_r = gsm_deaverage_comb.filter_gaussian(160).replace_slr_c20(slr_c20_deaverage)
Substract the raw C20 from Gaussian filtered C20.
print(gsm_r_gau.shc[:,0,2,0]  gsm_gau_r.shc[:,0,2,0])
[2.18616534e13 1.67577316e13 6.38921734e14 4.03942446e14
2.04186367e14 1.55989259e14 7.60951083e14 1.30531151e13
...
8.14599500e14 3.35849413e14 5.16146839e14 2.63274544e14
7.15681622e15 5.79722296e14 2.55374514e13 2.51142774e13
2.72619567e13 2.85867403e13]
The difference between the two is still smaller than the standard deviation of C20(magnitude of $10^{11}$), so using either method is acceptable. It is recommended to perform Gaussian filtering after replacing C20.
Convert geopotential coefficients to surface mass anomaly
Three ways to express surface mass anomaly(SMA) are provided. They are equivalent water thickness(EWT), equivalent ice thickness(EIT), and equivalent sand thickness(EST). If no equivalent substance is specified, the default is the equivalent water thickness in [mm w.e.].
gsm_r_ddk5_sma = gsm_r_ddk5.sma() print(gsm_r_ddk5_sma) print(gsm_r_ddk5_sma.summary)
title = Stokes coefficients for Surface Mass Anomaly(SMA) in Equivalent Water Thickness(EWT) derived from the DDK5 filtered Combined Deaveraged GRACE & GRACEFO Geopotential Coefficients RL06 with C20 replaced by the SLR measurements
max_degree = 96
max_order = 96
degree_order = 96
normalization = fully normalized
institution = UTAUSTIN/CSR, GFZ German Research Centre for Geosciences, NASA/JPL
processing_level = 2
product_version = RL06
time_coverage_start = 200204
time_coverage_end = 201910
solution_counts = 178
total_month_counts = 211
missing_month_counts = 33
Spherical harmonic coefficients representing an estimate of the Surface Mass Anomaly(SMA) expressed in terms of Equivalent Water[1000kg/m3] Thickness(EWT) with unit of [mm w.e.] during the specified timespan derived from GRACE & GRACEFO mission measurements. These coefficients represent the full magnitude of land hydrology, ice, and solid Earth processes. Further, they represent atmospheric and oceanic processes not captured in the accompanying GAC product. Note that the 2nddegree terms have been replaced with the values from SLR C20. Also note that C20 from SLR also experienced the DDK5 filtering.
For more details, please refer to gsm_r_ddk5.sma?
.
Note that the DDK filtering is performed before the conversion. However, there is no problem in performing the conversion first and then the filtering, because the difference between the two is orders of magnitude smaller than the uncertainty, although the conversion and filtering do not meet the exchange law.
gsm_r_sma_ddk5 = gsm_deaverage_comb_r.sma().filter_ddk('DDK5') print(gsm_r_sma_ddk5)
title = DDK5 filtered Stokes coefficients for Surface Mass Anomaly(SMA) in Equivalent Water Thickness(EWT) derived from the Combined Deaveraged GRACE & GRACEFO Geopotential Coefficients RL06 with C20 replaced by the SLR measurements
max_degree = 96
max_order = 96
degree_order = 96
normalization = fully normalized
institution = UTAUSTIN/CSR, GFZ German Research Centre for Geosciences, NASA/JPL
processing_level = 2
product_version = RL06
time_coverage_start = 200204
time_coverage_end = 201910
solution_counts = 178
total_month_counts = 211
missing_month_counts = 33
Convert SMA to geopotential coefficients
This is not commonly used.
print(gsm_r_ddk5_sma.gsm())
title = DDK5 filtered Combined Deaveraged GRACE & GRACEFO Geopotential Coefficients RL06 with C20 replaced by the SLR measurements
max_degree = 96
max_order = 96
degree_order = 96
normalization = fully normalized
institution = UTAUSTIN/CSR, GFZ German Research Centre for Geosciences, NASA/JPL
processing_level = 2
product_version = RL06
time_coverage_start = 200204
time_coverage_end = 201910
solution_counts = 178
total_month_counts = 211
missing_month_counts = 33
For more details, please refer to gsm_r_ddk5_sma.gsm?
.
Set up the study area
Two ways are listed to set the study area. one is to provide a boundary file, and the other is to create an ellipse.
Method 1: through a boundary file
Make sure the boundary file exists in the current directory. The first and second columns of the boundary file are latitude and longitude, respectively.
import numpy as np study_area = np.loadtxt('boundaries.txt') print(study_area)
[[32.75 93.75]
[32.25 93.75]
[32.25 92.75]
...
[32.25 97.25]
[32.75 97.25]
[32.75 93.75]]
Method 2: through creating an ellipse
An ellipse is defined by four parameters:
 Center of ellipse
 Rotation angle of the semimajor axis with respect to local east
 Semimajor axis and semiminor axis in degrees
Create an ellipse centered at (30N, 95E), with a semimajor axis of 4.2 degrees and a semiminor axis of 3.2 degrees.
from pyshtools.utils import MakeEllipseCoord lat0,lon0,dec0 = 30,95,0 semi_minor0, semi_major0 = 3.2,4.2 study_area = MakeEllipseCoord(lat0, lon0, dec0, semi_minor0, semi_major0) print(study_area)
[[33.2 95. ]
[33.19970027 95.0667118 ]
[33.19880084 95.13342754]
...
[33.19730097 94.79984886]
[33.19880084 94.86657246]
[33.19970027 94.9332882 ]]
Calculate the ratio of the solid angle formed by the ellipse to the global solid angle(4$\pi$). Note that $area = solid angle \times R^2$, where R is the Earth's mean radius.
from ggtools.gg import solid_angle_ratio area_ratio = solid_angle_ratio(semi_minor0,semi_major0) print(area_ratio)
0.0010421908639105157
For more details, please refer to solid_angle_ratio?
.
Plot the rate of SMA
# set the extent of the drawing area region = [82,106,21,39] # [left lon, right lon, lower lat, upper lat]
Method 1: estimate the rate of SHCs then convert it to grids
gsm_r_ddk5_sma_rate = gsm_r_ddk5_sma.rate() gsm_r_ddk5_sma_rate_grid = gsm_r_ddk5_sma_rate.grid(region) fig_name1 = 'sma_rate_grid_block.png' fig_name2 = 'sma_rate_grid.png' ylabel = 'SMA [mm w.e./yr]' gsm_r_ddk5_sma_rate_grid.plot(fig_name1,ylabel,'block',study_area) gsm_r_ddk5_sma_rate_grid.plot(fig_name2,ylabel,polygons=study_area)
For more details, please refer to gsm_r_ddk5_sma.rate?
, gsm_r_ddk5_sma_rate.grid?
, gsm_r_ddk5_sma_rate_grid?
, and gsm_r_ddk5_sma_rate_grid.plot?
.
Method 2: convert SHCs to grids then estimate the rate of grids
gsm_r_ddk5_sma_grid = gsm_r_ddk5_sma.grid(region) gsm_r_ddk5_sma_grid_rate = gsm_r_ddk5_sma_grid.rate() fig_name1 = 'sma_grid_rate_block.png' fig_name2 = 'sma_grid_rate.png' ylabel = 'SMA [mm w.e./yr]' gsm_r_ddk5_sma_grid_rate.plot(fig_name1,ylabel,'block',study_area) gsm_r_ddk5_sma_grid_rate.plot(fig_name2,ylabel,polygons=study_area)
The calculation will take a few minutes, please be patient.
Since the second method takes a lot of time to estimate the uncertainty of the grid data, the first method is recommended.
Calculate the time series of SMA over the study area and estimate its rate
Method 1: from SHCs to series and rate
Calculate the time series of SMA over the study area.
gsm_r_ddk5_sma_series = gsm_r_ddk5_sma.study_area(study_area) print(gsm_r_ddk5_sma_series,'\n') # time series of SMA print(gsm_r_ddk5_sma_series.qs,'\n') # standard deviation print(gsm_r_ddk5_sma_series.qs_std)
title = Integral(over the study area) of Stokes coefficients for Surface Mass Anomaly(SMA) in Equivalent Water Thickness(EWT) derived from the DDK5 filtered Combined Deaveraged GRACE & GRACEFO Geopotential Coefficients RL06 with C20 replaced by the SLR measurements
time_coverage_start = 200204
time_coverage_end = 201910
solution_counts = 178
total_month_counts = 211
missing_month_counts = 33
equi_material = Water
[ 7.80083803 58.62587324 73.76908473 37.28658384 53.38751612
30.71881195 5.74150633 32.72088735 1.42416991 18.75145636
...
47.60496348 73.31059842 73.00481244 92.60124409 78.85514767
53.20232491 39.92594104 51.5172863 42.71680212 8.95193687
22.16207019 24.10597124 40.20608753]
[ 6.30457663 5.46950829 3.75575696 7.60148364 3.6918224 3.67420728
4.16828886 5.49457631 3.99793319 3.27580461 2.89207024 3.71135591
...
1.82726744 2.03576616 1.80579235 1.81882144 1.75091745 1.79941898
1.75461899 1.83153935 1.98898853 1.77737887]
For more details, please refer to gsm_r_ddk5_sma.study_area?
and gsm_r_ddk5_sma_series?
.
Estimate the rate of the times series of SMA using ILSQM(Iterative Least Square Method).
gsm_r_ddk5_sma_series_rate = gsm_r_ddk5_sma_series.rate() print('rate: ',gsm_r_ddk5_sma_series_rate.qs,' ± ',gsm_r_ddk5_sma_series_rate.qs_std) print('area: ',gsm_r_ddk5_sma_series_rate.area,' km2')
rate: [6.70810897] ± [0.48489388]
area: 508282.0114287184 km2
For more details, please refer to gsm_r_ddk5_sma_series.rate?
.
Plot the time series of SMA.
fig_name = 'sma_series.png' ylabel = 'SMA [Gt]' gsm_r_ddk5_sma_series.plot(fig_name,ylabel)
For more details, please refer to gsm_r_ddk5_sma_series.plot?
.
Estimate the rate of the times series of SMA using IWLSQM(Iterative Weighted Least Square Method).
gsm_r_ddk5_sma_series_rate = gsm_r_ddk5_sma_series.rate('IWLSQM') print('rate: ',gsm_r_ddk5_sma_series_rate.qs,' ± ',gsm_r_ddk5_sma_series_rate.qs_std) print('area: ',gsm_r_ddk5_sma_series_rate.area,' km2')
rate: [6.59799566] ± [0.03844455]
area: 508282.0114287184 km2
Note that the rate of the time series through IWLSQM is close to that through IWLSQM, but the uncertainty of the rate differs by an order of magnitude.
Estimate the rate of SMA and its uncertainty using LSQM and WLSQM. Note that these two methods will not remove any outliers caused by abnormal monthly solutions.
gsm_r_ddk5_sma_series_rate = gsm_r_ddk5_sma_series.rate('LSQM') print('rate: ',gsm_r_ddk5_sma_series_rate.qs,' ± ',gsm_r_ddk5_sma_series_rate.qs_std) print('area: ',gsm_r_ddk5_sma_series_rate.area,' km2','\n') gsm_r_ddk5_sma_series_rate = gsm_r_ddk5_sma_series.rate('WLSQM') print('rate: ',gsm_r_ddk5_sma_series_rate.qs,' ± ',gsm_r_ddk5_sma_series_rate.qs_std) print('area: ',gsm_r_ddk5_sma_series_rate.area,' km2')
rate: [6.59361023] ± [0.49616184]
area: 508282.0114287184 km2
rate: [6.59776226] ± [0.03844441]
area: 508282.0114287184 km2
Method 2: from SHCs rate to series rate
gsm_r_ddk5_sma_rate_region = gsm_r_ddk5_sma_rate.study_area(study_area) print(gsm_r_ddk5_sma_rate_region) print('rate: ',gsm_r_ddk5_sma_rate_region.qs,' ± ',gsm_r_ddk5_sma_rate_region.qs_std) print('area: ',gsm_r_ddk5_sma_rate_region.area,' km2')
title = Integral(over the study area) of Annual change rate of Stokes coefficients for Surface Mass Anomaly(SMA) in Equivalent Water Thickness(EWT) derived from the DDK5 filtered Combined Deaveraged GRACE & GRACEFO Geopotential Coefficients RL06 with C20 replaced by the SLR measurements
time_coverage_start = 200204
time_coverage_end = 201910
solution_counts = 178
total_month_counts = 211
missing_month_counts = 33
equi_material = Water
rate: [6.86790848] ± [0.33702872]
area: 508282.0114287184 km2
Method 3: from grids to series and rate
gsm_r_ddk5_sma_grid_series = gsm_r_ddk5_sma_grid.study_area(study_area) gsm_r_ddk5_sma_grid_series_rate = gsm_r_ddk5_sma_grid_series.rate() print('rate: ',gsm_r_ddk5_sma_grid_series_rate.qs,' ± ',gsm_r_ddk5_sma_grid_series_rate.qs_std) print('area: ',gsm_r_ddk5_sma_grid_series_rate.area,' km2')
rate: [6.70774072] ± [0.48485163]
area: 508317.01700056973 km2
Method 4.1: from grids rate to series rate
The grids rate is obtained by estimating the rate of the grid.
gsm_r_ddk5_sma_grid_rate_region = gsm_r_ddk5_sma_grid_rate.study_area(study_area) print(gsm_r_ddk5_sma_grid_rate_region) print('rate: ',gsm_r_ddk5_sma_grid_rate_region.qs,' ± ',gsm_r_ddk5_sma_grid_rate_region.qs_std) print('area: ',gsm_r_ddk5_sma_grid_rate_region.area,' km2')
title = Annual change rate of Sum(over the study area) of grids expanded from Stokes coefficients for Surface Mass Anomaly(SMA) in Equivalent Water Thickness(EWT) derived from the DDK5 filtered Combined Deaveraged GRACE & GRACEFO Geopotential Coefficients RL06 with C20 replaced by the SLR measurements
time_coverage_start = 200204
time_coverage_end = 201910
solution_counts = 178
total_month_counts = 211
missing_month_counts = 33
equi_material = Water
rate: [6.74723381] ± [0.08916294]
area: 508317.01700056973 km2
Method 4.2: from grids rate to series rate
The grids rate is obtained by converting the SHCs rate.
gsm_r_ddk5_sma_rate_grid_region = gsm_r_ddk5_sma_rate_grid.study_area(study_area) print(gsm_r_ddk5_sma_rate_grid_region) print('rate: ',gsm_r_ddk5_sma_rate_grid_region.qs,' ± ',gsm_r_ddk5_sma_rate_grid_region.qs_std) print('area: ',gsm_r_ddk5_sma_rate_grid_region.area,' km2')
title = Sum(over the study area) of grids expanded from Annual change rate of Stokes coefficients for Surface Mass Anomaly(SMA) in Equivalent Water Thickness(EWT) derived from the DDK5 filtered Combined Deaveraged GRACE & GRACEFO Geopotential Coefficients RL06 with C20 replaced by the SLR measurements
time_coverage_start = 200204
time_coverage_end = 201910
solution_counts = 178
total_month_counts = 211
missing_month_counts = 33
equi_material = Water
rate: [6.8675455] ± [0.05935396]
area: 508317.01700056973 km2
In theory, the four methods mentioned above are equivalent. Although the rates of the time series estimated by these methods are approximately equal, the standard deviations obtained by the first three methods are several times larger than those of the latter method. The first and the third method avoid the rate of SHCs or the rate of grids data, that is, the rate of the time series is obtained by an LSQ fitting directly. Therefore, these two methods are recommended. Compared with the third method, the first method is less computationally expensive, so the first method takes less running time and is recommended.
GLDAS
Download GLDAS data
The GLDAS(Global Land Data Assimilation System) grid data is downloaded from GES DISC(Goddard Earth Sciences Data and Information Services Center) by default. Before downloading the data, make sure you have an EARTHDATA account. If not, please go to the official EARTHDATA website to register one.
Enter your username and password for EARTHDATA.
username,password = 'your_username','your_password'
Download GLDAS monthly grid data with a spatial resolution of 1 degree from 200204 to 201910.
from ggtools.gg import gldas_download start_date,end_date = '200204','201910' gldas_download(username,password,start_date,end_date)
For more details, please refer to gldas_download?
.
Read GLDAS data
Read the GLDAS data from 200204 to 201910. Note that the number of the GLDAS monthly grids should be consistent with the number of the GRACE GSM monthly solutions, if you want to compare the two types of data.
from ggtools.gg import read_gldas gldas = read_gldas('200204','201910') print(gldas,'\n') print(gldas.lats,'\n') print(gldas.lons)
title = GLDAS2.1 LIS land surface model output monthly mean
resolution = 1deg
institution = NASA GSFC
time_coverage_start = 200204
time_coverage_end = 201910
solution_counts = 211
total_month_counts = 211
missing_month_counts = 0
<xarray.DataArray 'lat' (lat: 150)>
array([59.5, 58.5, 57.5, ..., 87.5, 88.5, 89.5], dtype=float32)
Coordinates:
* lat (lat) float32 59.5 58.5 57.5 56.5 55.5 ... 86.5 87.5 88.5 89.5
Attributes:
units: degrees_north
standard_name: latitude
long_name: latitude
vmin: 59.5
vmax: 89.5
<xarray.DataArray 'lon' (lon: 360)>
array([179.5, 178.5, 177.5, ..., 177.5, 178.5, 179.5], dtype=float32)
Coordinates:
* lon (lon) float32 179.5 178.5 177.5 176.5 ... 177.5 178.5 179.5
Attributes:
units: degrees_east
standard_name: longitude
long_name: longitude
vmin: 179.5
vmax: 179.5
For more details, please refer to read_gldas?
and gldas?
.
Calculate Terrestrial Water Storage Change
Two ways to estimate TWSC are provided. The first one is the classic(traditional) technique, i.e. $TWSC = SoilMoi + AccumSnow + Canopy$. The second is through the water balance equation, i.e. $TWSC = Precipitation  WaterEvaporation  SurfaceRunoff  SubsurfaceRunoff$. The difference in results obtained by these two methods is minimal. Note: The TWSC has been deaveraged. The unit of TWSC is [kg/m2] or [mm w.e.].
twsc_grid = gldas.twsc_grid(region) print(twsc_grid) print(twsc_grid.summary)
title = Terrestrial Water Storage Change(TWSC) derived from the GLDAS2.1 LIS land surface model output monthly mean
time_coverage_start = 200204
time_coverage_end = 201910
solution_counts = 211
total_month_counts = 211
missing_month_counts = 0
equi_material = Water
region = [82, 106, 21, 39]
TWSC is estimated by the formulation: [TWSC = SoilMoi(0200cm) + Accum_Snow + Canopy], and it has been converted to Equivalent Water Thickness in mm w.e.
For more details, please refer to gldas.twsc_grid?
.
Estimate the rate of TWSC
twsc_grid_rate = twsc_grid.rate() print(twsc_grid_rate)
title = Annual change rate of Terrestrial Water Storage Change(TWSC) derived from the GLDAS2.1 LIS land surface model output monthly mean
time_coverage_start = 200204
time_coverage_end = 201910
solution_counts = 211
total_month_counts = 211
missing_month_counts = 0
equi_material = Water
region = [82, 106, 21, 39]
Plot the rate of TWSC
fig_name1 = 'twsc_grid_rate_block.png' fig_name2 = 'twsc_grid_rate.png' ylabel = 'TWSC [mm w.e./yr]' twsc_grid_rate.plot(fig_name1,ylabel,'block',study_area) twsc_grid_rate.plot(fig_name2,ylabel,polygons=study_area)
Calculate the time series of TWSC and estimate its rate
Calculate the time series of TWSC.
twsc_grid_series = twsc_grid.study_area(study_area)
Estimate the rate of the times series.
twsc_grid_series_rate = twsc_grid_series.rate() print(twsc_grid_series_rate,'\n') print('rate: ',twsc_grid_series_rate.qs,' ± ',twsc_grid_series_rate.qs_std) print('area: ',twsc_grid_series_rate.area,' km2')
title = Annual change rate of Terrestrial Water Storage Change(TWSC) derived from the GLDAS2.1 LIS land surface model output monthly mean
time_coverage_start = 200204
time_coverage_end = 201910
solution_counts = 211
total_month_counts = 211
missing_month_counts = 0
equi_material = Water
rate: [0.33623102] ± [0.36229108]
area: 504185.85874313745 km2
Plot the time series of TWSC
fig_name = 'twsc_series.png' ylabel = 'TWSC [Gt]' twsc_grid_series.plot(fig_name,ylabel)
Signal leakage correction
Method 1: scale factor
The scale factor(or gain factor) is sensitive to the linear trend of the time series and not to the annual term. Besides, it has local characteristics, i.e., different regions correspond different scale factors. This method is only applicable when the SMA distribution derived from GRACE in the study area is similar to the TWSC distribution derived from GLDAS.
from ggtools.gg import scale_factor twsc_coeffs = gldas.twsc_shc(96) # 96thdegree SHCs for TWSC print(twsc_coeffs,'\n') k = scale_factor(twsc_coeffs,study_area,160) print(k)
title = Terrestrial Water Storage Change(TWSC) derived from the GLDAS2.1 LIS land surface model output monthly mean
max_degree = 89
max_order = 89
degree_order = 96
normalization = fully normalized
institution = NASA GSFC
processing_level = CF1.6
product_version = Noah_v3.3 forced with GDASAGRMETGPCP
time_coverage_start = 200204
time_coverage_end = 201910
solution_counts = 211
total_month_counts = 211
missing_month_counts = 0
0.9867288021260247
For more details, please refer to gldas.twsc_shc?
and scale_factor?
.
gsm_r_ddk5_sma_series_rate_leakage = gsm_r_ddk5_sma_series_rate.leakage_correction('scale_factor',k) print(gsm_r_ddk5_sma_series_rate_leakage,'\n') print('rate: ',gsm_r_ddk5_sma_series_rate_leakage.qs,' ± ',gsm_r_ddk5_sma_series_rate_leakage.qs_std) print('area: ',gsm_r_ddk5_sma_series_rate_leakage.area,' km2')
title = Signal leakage corrected(Scale factor method) Annual change rate of Integral(over the study area) of Stokes coefficients for Surface Mass Anomaly(SMA) in Equivalent Water Thickness(EWT) derived from the DDK5 filtered Combined Deaveraged GRACE & GRACEFO Geopotential Coefficients RL06 with C20 replaced by the SLR measurements
time_coverage_start = 200204
time_coverage_end = 201910
solution_counts = 178
total_month_counts = 211
missing_month_counts = 33
equi_material = Water
rate: [6.61908433] ± [0.47845876]
area: 508282.0114287184 km2
For more details, please refer to gsm_r_ddk5_sma_series_rate.leakage_correction?
.
Method 2 : forward modeling
# calculate the global grid for SMA gsm_r_ddk5_sma_grid = gsm_r_ddk5_sma.grid() gsm_r_ddk5_sma_grid_leakage = gsm_r_ddk5_sma_grid.leakage_correction('forward_model',160) gsm_r_ddk5_sma_grid_leakage_series = gsm_r_ddk5_sma_grid_leakage.study_area(study_area) gsm_r_ddk5_sma_grid_leakage_series_rate = gsm_r_ddk5_sma_grid_leakage_series.rate() print('rate: ',gsm_r_ddk5_sma_grid_leakage_series_rate.qs,' ± ',gsm_r_ddk5_sma_grid_leakage_series_rate.qs_std) print('area: ',gsm_r_ddk5_sma_grid_leakage_series_rate.area,' km2')
rate: [8.37598258] ± [0.50297716]
area: 508317.01700056973 km2
For more details, please refer to gsm_r_ddk5_sma_grid.leakage_correction?
.
Method 3 : space domain
This method requires setting the locations of the mascons in advance. Take the glacier catalog with an area of not less than 10 $km^2$ in Southeast Tibet as an example.
import numpy as np from ggtools.gg import generate_nodes glaciers = np.loadtxt('glaciers.txt') print(glaciers)
[[30.650817 99.443696]
[29.91388 99.635285]
[28.492921 98.660107]
...
[28.985472 97.526066]
[30.292734 90.38734 ]
[30.415653 90.54994 ]]
Generate the mascon nodes(centers) based on the glacier catalog.
nodes = generate_nodes(gsm_r_ddk5_sma_rate_grid,glaciers,study_area) gsm_r_ddk5_sma_rate_grid_leakage = gsm_r_ddk5_sma_rate_grid.leakage_correction('space_domain',160,nodes,study_area) gsm_r_ddk5_sma_rate_grid_leakage_series = gsm_r_ddk5_sma_rate_grid_leakage.study_area(study_area) print('rate: ',gsm_r_ddk5_sma_rate_grid_leakage_series.qs,' ± ',gsm_r_ddk5_sma_rate_grid_leakage_series.qs_std) print('area: ',gsm_r_ddk5_sma_rate_grid_leakage_series.area,' km2')
rate: [8.55902161] ± [0.]
area: 508317.01700056973 km2
Because it involves the uncertainties of estimated parameters in Tikhonov regularization, the standard deviation for the rate is temporarily set to 0. This will be improved in the next version of ggtools.
Plot the rate of mascons.
fig_name1 = 'SMA_SpaceD_block.png' fig_name2 = 'SMA_SpaceD.png' ylabel = 'SMA [mm w.e./yr]' gsm_r_ddk5_sma_rate_grid_leakage.plot(fig_name1,ylabel,'block',study_area,nodes.nodes) gsm_r_ddk5_sma_rate_grid_leakage.plot(fig_name2,ylabel,polygons=study_area,circles=nodes.nodes)
Method 4 : inverse filtering
This method is equivalent to the forward modeling. The algorithm is extremely simple and avoids the iterative process in the forward modeling.
gsm_r_ddk5_sma_leakage = gsm_r_ddk5_sma.leakage_correction('filter_inverse',160) gsm_r_ddk5_sma_leakage_series = gsm_r_ddk5_sma_leakage.study_area(study_area) gsm_r_ddk5_sma_leakage_series_rate = gsm_r_ddk5_sma_leakage_series.rate() print('rate: ',gsm_r_ddk5_sma_leakage_series_rate.qs,' ± ',gsm_r_ddk5_sma_leakage_series_rate.qs_std) print('area: ',gsm_r_ddk5_sma_leakage_series_rate.area,' km2')
rate: [8.37974852] ± [0.50322993]
area: 508282.0114287184 km2
Method 5 : spectral domain
This method requires setting the locations of the mascons in advance. Still take the glacier catalog with an area of not less than 10 $km^2$ in Southeast Tibet as an example.
# global grids gsm_r_ddk5_sma_rate_grid = gsm_r_ddk5_sma_rate.grid() nodes = generate_nodes(gsm_r_ddk5_sma_rate_grid,glaciers,study_area) gsm_r_ddk5_sma_rate_leakage_grid = gsm_r_ddk5_sma_rate.leakage_correction('spectral_domain',160,nodes,study_area,'windows',area_ratio) gsm_r_ddk5_sma_rate_leakage_grid_series = gsm_r_ddk5_sma_rate_leakage_grid.study_area(study_area) print('rate: ',gsm_r_ddk5_sma_rate_leakage_grid_series.qs,' ± ',gsm_r_ddk5_sma_rate_leakage_grid_series.qs_std) print('area: ',gsm_r_ddk5_sma_rate_leakage_grid_series.area,' km2')
rate: [7.33355992] ± [0.]
area: 508317.01700056973 km2
Because it involves the uncertainties of estimated parameters in Tikhonov regularization, the standard deviation for the rate is temporarily set to 0. This will be improved in the next version of ggtools. Note that the 2x area ratio multiplication is just to use more tapers in windowed spherical harmonics.
gsm_r_ddk5_sma_rate_leakage_grid_crop = gsm_r_ddk5_sma_rate_leakage_grid.set_region(region) fig_name1 = 'SMA_SpectralD_block.png' fig_name2 = 'SMA_SpectralD.png' ylabel = 'SMA [mm w.e./yr]' gsm_r_ddk5_sma_rate_leakage_grid_crop.plot(fig_name1,ylabel,'block',study_area,nodes.nodes) gsm_r_ddk5_sma_rate_leakage_grid_crop.plot(fig_name2,ylabel,polygons=study_area,circles=nodes.nodes)
Method 6 : deconvolution
In Forward Modeling, the purpose of the signal leakage correction is to find a solution for the equation $\mathcal F(X) = Y$, where $X$ is the true mass distribution to be determined, $Y$ is the apparent mass distribution, and $\mathcal F$ represents the combination of spherical harmonic truncations and filter smoothing. One way to solve this equation is by iteration, which happens to be what the Forward Modeling does. Another possible approach is to find the expression of $\mathcal F$. Fortunately, $\mathcal F(X)$ can be expressed as a convolution, i.e., $\mathcal F(X)\equiv X\otimes PSF$, where PSF is a convolution kernel, i.e., the point spread function. According to $X = \mathcal F^{1}(Y)$, we can deduce $X = Y \otimes/\otimes PSF$, where $\otimes/\otimes$ denotes the deconvolution operator. The unsupervised WienerHunt algorithm can be employed to implement the deconvolution, and it has been one of the modules among the scikitimage package.
'Equivalent' Gaussian filter radius for DDK filtering
Place a unit mass with an equivalent water thickness of 1000mm at the North Pole and expand it to a specific degree, and perform DDK filtering and Gaussian filtering, respectively. Change the Gaussian filtering radius to maximize the correlation coefficient of the filtered spectrum. At this time, the Gaussian filtering radius is the 'equivalent' radius corresponding to the DDK filtering. Note: This method is suitable for DDK4 ~ DDK8 filtering. For 'equivalent radius' of DDK1 ~ DDK3, please refer to Kusche 2009.
from ggtools.gg import ddk_gaussian
for i in range(1,9):
ddk_gaussian('DDK'+str(i),96)
ddk_gaussian('DDK5',96,'visible')
Correlation: 0.9776
Approximate equivalent Gaussian filter radius for DDK1: 345
Correlation: 0.9671
Approximate equivalent Gaussian filter radius for DDK2: 255
Correlation: 0.9593
Approximate equivalent Gaussian filter radius for DDK3: 200
Correlation: 0.9578
Approximate equivalent Gaussian filter radius for DDK4: 190
Correlation: 0.9570
Approximate equivalent Gaussian filter radius for DDK5: 160
Correlation: 0.9583
Approximate equivalent Gaussian filter radius for DDK6: 150
Correlation: 0.9660
Approximate equivalent Gaussian filter radius for DDK7: 125
Correlation: 0.9714
Approximate equivalent Gaussian filter radius for DDK8: 110
Correlation: 0.9570
Approximate equivalent Gaussian filter radius for DDK5: 160
(a) Unit mass after expanding and truncating up to 96th degree. (b) after DDK5 filtering (c) after Gaussian filtering with a radius of 160km
<center>Cross section for the unit mass</center> For more details, please refer to `ddk_gaussian?`.
GRACE minus TWSC
# for signal leakage correction with inverse filtering
gsm_twsc_series = gsm_r_ddk5_sma_leakage_series  twsc_grid_series
gsm_twsc_series_rate = gsm_twsc_series.rate()
fig_name = 'GRACE_TWSC_series.png'
ylabel = 'GRACETWSC [Gt]'
gsm_twsc_series.plot(fig_name,ylabel,kernel='rbf')
Change log

1.1.7 — Apr 12, 2020
 Fixed the issue that it fails to download GLDAS data under windows platform.

1.1.6 — Mar 31, 2020
 Fixed the issue that the boundary file(defines the study area) goes across the prime meridian. From now on, boundary files that pass through the prime meridian are supported.
Next release

Complete the help documentation

Improve the code structure to make it easier to read

Add the degree and order (d/o) one(Geocenter) correction

Add outliers elimination in Gaussian Process Regression(GPR)

Add a module to handle Glacial Isostatic Adjustment(GIA) effects

Add other destriping filters, such as P4M6

Add other map projections, such as AlbersEqualArea

Find a way to quickly calculate the uncertainty of the grid data

Estimate the uncertainty of parameters in Tikhonov regularization
Acknowledgments
Thank the ISDC for sharing the GRACE & GRACEFO data, GES DISC for the GLDAS land surface model data, and the Cold and Arid Regions Environmental and Engineering Research Institute(CAREERI) for providing the Second Glacier Inventory Dataset of China(SGIDC). Many appreciations to the contributors of SHTOOLS and GRACEfilter.
Project details
Download files
Download the file for your platform. If you're not sure which to choose, learn more about installing packages.