GRACE/GRACE-FO Spatial Map Program
This (notebook) demonstrates calculating spatial maps from Gravity Recovery and Climate Experiment (GRACE) and the GRACE Follow-On (GRACE-FO) Level-2 spherical harmonic products.
The primary and secondary instrumentation onboard the GRACE/GRACE-FO satellites are the ranging instrument (GRACE has a microwave ranging instrument, GRACE-FO has both a microwave ranging instrument and a laser interferometer), the global positioning system (GPS), the accelerometers and the star cameras. Data from these instruments are combined to estimate the distance between the two satellites, the positions of the satellites in space, the pointing vector of the satellites and any non-gravitational accelerations the satellites experience. These measurements combined with background gravity field estimates, atmospheric and oceanic variability estimates and tidal estimates are used to create the Level-2 spherical harmonic product of GRACE and GRACE-FO.
There are three main processing centers that create the Level-2 spherical harmonic data as part of the GRACE/GRACE-FO Science Data System (SDS): the University of Texas Center for Space Research (CSR), the German Research Centre for Geosciences (GeoForschungsZentrum, GFZ) and the Jet Propulsion Laboratory (JPL).
The data is freely available in the US from the NASA Physical Oceanography Distributed Active Archive Center (PO.DAAC) and internationally from the GFZ Information System and Data Center.
There are programs within this repository that can sync with both of these data archives podaac_cumulus.py and gfz_isdc_grace_ftp.py.
Some of the text included within this notebook about the data processing details comes from Sutterley and Velicogna (2019) and Sutterley et al. (2020).
Load necessary modules for running the notebook
import numpy as np
import matplotlib
matplotlib.rcParams['mathtext.default'] = 'regular'
matplotlib.rcParams["animation.html"] = "jshtml"
matplotlib.rcParams["animation.embed_limit"] = 50
import matplotlib.pyplot as plt
import matplotlib.animation as animation
import cartopy.crs as ccrs
import ipywidgets
from IPython.display import HTML, Latex
import gravity_toolkit as gravtk
Set the GRACE/GRACE-FO Data Directory
Directory should contain:
Technical notes with SLR coefficients
Subdirectory with geocenter coefficients
Subdirectories for each processing center
# set the directory with GRACE/GRACE-FO data
# update local data with PO.DAAC https servers
widgets = gravtk.tools.widgets()
ipywidgets.VBox([
widgets.directory,
widgets.update,
widgets.endpoint
])
Update Data in Directory
# if updating the local data
if widgets.update.value:
# run podaac sync program to get latest data
!podaac_cumulus.py --directory=$widgets.base_directory --endpoint=$widgets.endpoint
# run GRACE date program to verify months
!run_grace_date.py --directory=$widgets.base_directory --verbose
# get geocenter data from Sutterley and Velicogna (2019)
gravtk.utilities.from_figshare(widgets.base_directory)
Set GRACE/GRACE-FO Parameters
These parameters describe the specific GRACE/GRACE-FO product and the months of data to read
GRACE/GRACE-FO Processing Center
CSR: University of Texas Center for Space Research
GFZ: German Research Centre for Geosciences (GeoForschungsZentrum)
JPL: Jet Propulsion Laboratory
CNES: French Centre National D’Etudes Spatiales
GRACE/GRACE-FO Data Release
GRACE/GRACE-FO Data Product
GAA: non-tidal atmospheric correction
GAB: non-tidal oceanic correction
GAC: combined non-tidal atmospheric and oceanic correction
GAD: GRACE/GRACE-FO ocean bottom pressure product
GSM: corrected monthly GRACE/GRACE-FO static field product
GRACE/GRACE-FO Date Range
# update widgets
widgets.select_product()
# display widgets for setting GRACE/GRACE-FO parameters
ipywidgets.VBox([
widgets.center,
widgets.release,
widgets.product,
widgets.months
])
Set Parameters for Reading GRACE/GRACE-FO Data
These parameters describe processing steps and corrections to be applied when reading the GRACE/GRACE-FO data
Maximum Degree and Order
Geocenter product (Degree 1)
Oblateness product (C20)
Figure axis product (C21 and S21)
Azimuthal dependence product (C22 and S22)
Low Degree Zonal products (C30, C40 and C50)
Pole Tide Correction from Wahr et al. (2015)
Atmospheric Correction as described in Fagiolini et al. (2015)
Geocenter
Measurements of time-variable gravity from the Gravity Recovery and Climate Experiment (GRACE) and the GRACE Follow-On (GRACE-FO) missions are set in a center of mass (CM) reference frame, in which the total degree one variations are inherently zero. The individual contributions to degree one variations in the CM reference frame, such as from oceanic processes or terrestrial water storage change, are not necessarily zero (Wahr et al., 1998). Applications set in a center of figure (CF) reference frame, such as the recovery of mass variations of the oceans, hydrosphere and cryosphere, require the inclusion of degree one terms to be fully accurate (Swenson et al., 2008). The geocenter options above can select the degree one product to include with the GRACE/GRACE-FO derived harmonics of degree two and higher. There are options for using measurements from satellite laser ranging (Cheng, 2013) and calculations from time-variable gravity and ocean model outputs (Swenson et al., 2008).
Low Degree Harmonics
For both GRACE and GRACE-FO, there have been operational issues that have affected the quality of the time-variable gravity fields. During the late stages of the GRACE mission, procedures were enacted to preserve the battery life of the GRACE satellites and extend the mission lifetime. This included turning off the accelerometer and microwave ranging instrument (MWI) to reduce the battery load on the spacecraft during periods of low β angle when solar input to the spacecraft is lowest. In late 2016, the accelerometer onboard GRACE-B was permanently powered down to help maintain the operation of the ranging instrument. For the Follow-On mission, an anomaly shortly after the launch of GRACE-FO caused a malfunction and the accelerometer onboard GRACE-FO 2 (GF2) has been operating in a less optimal Large-Range-Mode. For both of these cases, the GRACE/GRACE-FO processing centers have developed independent methods to spatiotemporally transplant the accelerometer data retrieved from GRACE-A to GRACE-B and from GF1 to GF2 (Bandikova et al., 2019). These single-accelerometer months for both GRACE and GRACE-FO contain significantly more noise, particularly the low-degree zonal harmonics (predominantly C20 and C30 but possibly C40 and C50). C20 has also been difficult for GRACE and GRACE-FO to independently measure throughout both missions. The figure axis harmonics (C21 and S21) may also be contaminated by noise during the single-accelerometer months in the GFZ products Dahle et al, 2019. Measurements from satellite laser ranging (SLR) can provide an independent assessment for some low degree and order spherical harmonics (Cheng et al., 2011; Cheng and Ries, 2018; Loomis et al., 2019; Loomis et al., 2020). The C20, C21, C22, C30, C40, and C50 options below can select the SLR low degree product to use as a replacement for the GRACE/GRACE-FO derived estimates.
Corrections
Prior to GRACE/GRACE-FO Release-6, corrections needed to be applied to compensate for long-period signals in the pole tide that were contaminating the C21 and S21 harmonics (Wahr et al., 2015), as well as for discontinuities in the atmospheric de-aliasing product that were introduced with upgrades in the ECMWF weather prediction model (Fagiolini et al., 2015). The Pole Tide and Atmospheric corrections do not need to be applied to the Release-6 data.
# update widgets
widgets.select_options()
# display widgets for setting GRACE/GRACE-FO read parameters
ipywidgets.VBox([
widgets.lmax,
widgets.mmax,
widgets.geocenter,
widgets.C20,
widgets.CS21,
widgets.CS22,
widgets.C30,
widgets.C40,
widgets.C50,
widgets.pole_tide,
widgets.atm
])
Read GRACE/GRACE-FO data
This step extracts the parameters chosen above and then reads the GRACE/GRACE-FO data applying the specified procedures
# extract values from widgets
PROC = widgets.center.value
DREL = widgets.release.value
DSET = widgets.product.value
months = [int(m) for m in widgets.months.value]
LMAX = widgets.lmax.value
MMAX = widgets.mmax.value
DEG1 = widgets.geocenter.value
SLR_C20 = widgets.C20.value
SLR_21 = widgets.CS21.value
SLR_22 = widgets.CS22.value
SLR_C30 = widgets.C30.value
SLR_C40 = widgets.C40.value
SLR_C50 = widgets.C50.value
POLE_TIDE = widgets.pole_tide.value
ATM = widgets.atm.value
# read GRACE/GRACE-FO data for parameters
start_mon = np.min(months)
end_mon = np.max(months)
missing = sorted(set(np.arange(start_mon,end_mon+1)) - set(months))
Ylms = gravtk.grace_input_months(widgets.base_directory, PROC, DREL, DSET,
LMAX, start_mon, end_mon, missing, SLR_C20, DEG1, MMAX=MMAX,
SLR_21=SLR_21, SLR_22=SLR_22, SLR_C30=SLR_C30, SLR_C40=SLR_C40,
SLR_C50=SLR_C50, POLE_TIDE=POLE_TIDE, ATM=ATM)
# create harmonics object and remove mean
GRACE_Ylms = gravtk.harmonics().from_dict(Ylms)
GRACE_Ylms.mean(apply=True)
# directory of specific GRACE/GRACE-FO product
GRACE_Ylms.directory = Ylms['directory']
# string denoting specific corrections and data used
GRACE_Ylms.title = Ylms['title']
# number of time steps
nt = len(months)
# flag for spherical harmonic order
order_str = f'M{MMAX:d}' if (MMAX != LMAX) else ''
Set Parameters to Convert to Spatial Maps
These parameters specify corrections and filtering steps for converting to the spatial domain at a specified grid spacing
GIA Correction
Remove Specific Harmonic Fields
Redistribute Removed Fields over the Ocean
Gaussian Smoothing Radius in kilometers
Filter (destripe) harmonics (Swenson and Wahr, 2006)
Spatial degree spacing
Spatial degree interval
(-180:180,90:-90)
(degree spacing)/2
Output spatial units
equivalent water thickness (cm)
geoid height (mm)
elastic crustal deformation (mm)
gravitational perturbation (μGal)
equivalent surface pressure (mbar)
Output spatial format
netCDF4
HDF5
Geophysical Leakage
Gravity measurements from GRACE and GRACE-FO are global, near-monthly and are directly related to changes in mass. Several mass transport processes can occur concurrently for a given region, which means that the total time-dependent geopotential from GRACE/GRACE-FO can relate to multiple time-varying components (Wahr et al., 1998). These mass transport processes include but are not limited to terrestrial water storage, glacier and ice sheet mass, atmospheric and oceanic circulation and geodynamic processes. In order to isolate the mass change of a single process, each of the other processes needs to be independently estimated and removed from the GRACE/GRACE-FO data. Uncertainties in the components removed from the GRACE/GRACE-FO data will directly impact the precision of the final mass balance estimate.
Filtering
The GRACE/GRACE-FO coefficients are impacted by random spherical harmonic errors that increase as a function of spherical harmonic degree (Wahr et al., 1998). The impact of these errors can be reduced using Gaussian averaging functions as described in Jekeli, (1981). GRACE/GRACE-FO coefficients are also impacted by correlated north/south “striping” errors, which can be spectrally filtered following Swenson and Wahr (2006).
Units
Spatial fields of surface mass density can be estimated from sets of spherical harmonics if we assume that the mass redistributions are concentrated within a thin layer (thickness ≪ horizontal resolution) after correcting for glacial isostatic adjustment (Wahr et al., 1998). However, we additionally need to compensate for the Earth’s elastic yielding to surface load changes, which induce density anomalies at depth (Wahr et al., 1998). This program accounts for the elastic deformation of the solid Earth using load Love numbers calculated by Han and Wahr (1995) with parameters from the Preliminary reference Earth model (PREM).
This program can output monthly spatial fields in a few different units (mass, deformation, and gravity).
Equivalent water thickness
Geoid height
Vertical elastic crustal deformation
\(a\): average radius of the Earth
\(\rho_e\): average density of the Earth
\(\rho_w\): density of water
\(h_l\), \(k_l\): Load Love numbers of degree \(l\)
\(P_{lm}\): fully-normalized Legendre polynomials of degree \(l\) and order \(m\)
\(C_{lm}\), \(S_{lm}\): cosine and sine spherical harmonics of degree \(l\) and order \(m\)
\(\theta\), \(\phi\): colatitude and longitude in radians
# update widgets
widgets.select_corrections()
widgets.select_output()
# display widgets for setting GRACE/GRACE-FO corrections parameters
ipywidgets.VBox([
widgets.GIA_file,
widgets.GIA,
widgets.remove_file,
widgets.remove_format,
widgets.redistribute_removed,
widgets.mask,
widgets.gaussian,
widgets.destripe,
widgets.spacing,
widgets.interval,
widgets.units,
widgets.output_format])
Convert GRACE/GRACE-FO harmonics to the spatial domain
This step extracts the parameters chosen above and then converts the GRACE/GRACE-FO harmonics to the spatial domain applying the specified corrections and filtering procedures
Set output grid domain
Calculate Fully-Normalized Legendre Polynomials
Read GIA model for correcting GRACE/GRACE-FO data
Read harmonics to be removed from the GRACE/GRACE-FO data
Calculate coefficients for converting to the output units
Convert from the spherical harmonic domain into the spatial domain
Output the monthly spatial maps to netCDF4 or HDF5
# Output spatial data
grid = gravtk.spatial()
grid.time = np.copy(GRACE_Ylms.time)
grid.month = np.copy(GRACE_Ylms.month)
# Output degree spacing
dlon = widgets.spacing.value
dlat = widgets.spacing.value
# Output Degree Interval
INTERVAL = widgets.interval.index + 1
if (INTERVAL == 1):
# (-180:180,90:-90)
nlon = np.int64((360.0/dlon)+1.0)
nlat = np.int64((180.0/dlat)+1.0)
grid.lon = -180 + dlon*np.arange(0,nlon)
grid.lat = 90.0 - dlat*np.arange(0,nlat)
elif (INTERVAL == 2):
# (Degree spacing)/2
grid.lon = np.arange(-180+dlon/2.0,180+dlon/2.0,dlon)
grid.lat = np.arange(90.0-dlat/2.0,-90.0-dlat/2.0,-dlat)
nlon = len(grid.lon)
nlat = len(grid.lat)
# Computing plms for converting to spatial domain
theta = (90.0-grid.lat)*np.pi/180.0
PLM, dPLM = gravtk.plm_holmes(LMAX, np.cos(theta))
# read load love numbers file
# PREM outputs from Han and Wahr (1995)
# https://doi.org/10.1111/j.1365-246X.1995.tb01819.x
love_numbers_file = gravtk.utilities.get_data_path(['data','love_numbers'])
header = 2
columns = ['l','hl','kl','ll']
# LMAX of load love numbers from Han and Wahr (1995) is 696.
# from Wahr (2007) linearly interpolating kl works
# however, as we are linearly extrapolating out, do not make
# LMAX too much larger than 696
# read arrays of kl, hl, and ll Love Numbers
hl,kl,ll = gravtk.read_love_numbers(love_numbers_file, LMAX=LMAX,
HEADER=header, COLUMNS=columns, REFERENCE='CF', FORMAT='tuple')
# read GIA data
GIA = widgets.GIA.value
GIA_Ylms_rate = gravtk.gia(lmax=LMAX).from_GIA(widgets.GIA_model,
GIA=GIA, mmax=MMAX)
gia_str = '' if (GIA == '[None]') else f'_{GIA_Ylms_rate.title}'
# calculate the monthly mass change from GIA
# monthly GIA calculated by gia_rate*time elapsed
# finding change in GIA each month
GIA_Ylms = GIA_Ylms_rate.drift(GRACE_Ylms.time, epoch=2003.3)
GIA_Ylms.month[:] = np.copy(GRACE_Ylms.month)
# if redistributing removed mass over the ocean
if widgets.redistribute_removed.value:
# read Land-Sea Mask and convert to spherical harmonics
ocean_Ylms = gravtk.ocean_stokes(widgets.landmask, LMAX,
MMAX=MMAX, LOVE=(hl,kl,ll))
# read data to be removed from GRACE/GRACE-FO monthly harmonics
remove_Ylms = GRACE_Ylms.zeros_like()
remove_Ylms.time[:] = np.copy(GRACE_Ylms.time)
remove_Ylms.month[:] = np.copy(GRACE_Ylms.time)
# If there are files to be removed from the GRACE/GRACE-FO data
# for each file separated by commas
for f in widgets.remove_files:
if (widgets.remove_format.value == 'netCDF4'):
# read netCDF4 file
Ylms = gravtk.harmonics().from_netCDF4(f)
elif (widgets.remove_format.value == 'HDF5'):
# read HDF5 file
Ylms = gravtk.harmonics().from_HDF5(f)
elif (widgets.remove_format.value == 'index (ascii)'):
# read index of ascii files
Ylms = gravtk.harmonics().from_index(f,format='ascii')
elif (widgets.remove_format.value == 'index (netCDF4)'):
# read index of netCDF4 files
Ylms = gravtk.harmonics().from_index(f,format='netCDF4')
elif (widgets.remove_format.value == 'index (HDF5)'):
# read index of HDF5 files
Ylms = gravtk.harmonics().from_index(f,format='HDF5')
# reduce to months of interest and truncate to range
Ylms = Ylms.subset(months).truncate(LMAX,mmax=MMAX)
# redistribute removed mass over the ocean
if widgets.redistribute_removed.value:
# calculate ratio between total removed mass and
# a uniformly distributed cm of water over the ocean
ratio = Ylms.clm[0,0,:]/ocean_Ylms.clm[0,0]
# for each spherical harmonic
for m in range(0,MMAX+1):
for l in range(m,LMAX+1):
# remove the ratio*ocean Ylms from Ylms
Ylms.clm[l,m,:]-=ratio*ocean_Ylms.clm[l,m]
Ylms.slm[l,m,:]-=ratio*ocean_Ylms.slm[l,m]
# add the harmonics to be removed to the total
remove_Ylms.add(Ylms)
# gaussian smoothing radius in km (Jekeli, 1981)
RAD = widgets.gaussian.value
if (RAD != 0):
wt = 2.0*np.pi*gravtk.gauss_weights(RAD,LMAX)
gw_str = f'_r{RAD:0.0f}km'
else:
# else = 1
wt = np.ones((LMAX+1))
gw_str = ''
# destriping the GRACE/GRACE-FO harmonics
ds_str = '_FL' if widgets.destripe.value else ''
# Setting units factor for output
UNITS = widgets.unit_index
# dfactor is the degree dependent coefficients
# for specific spherical harmonic output units
factors = gravtk.units(lmax=LMAX).harmonic(hl,kl,ll)
# 1: cmwe, centimeters water equivalent
# 2: mmGH, millimeters geoid height
# 3: mmCU, millimeters elastic crustal deformation
# 4: micGal, microGal gravity perturbations
# 5: mbar, millibars equivalent surface pressure
dfactor = factors.get(gravtk.units.bycode(UNITS))
# units strings for output files and plots
unit_label = ['cm', 'mm', 'mm', u'\u03BCGal', 'mb']
unit_name = ['Equivalent Water Thickness', 'Geoid Height',
'Elastic Crustal Uplift', 'Gravitational Undulation',
'Equivalent Surface Pressure']
# converting harmonics to truncated, smoothed coefficients in units
# combining harmonics to calculate output spatial fields
# output spatial grid
grid.data = np.zeros((nlat, nlon, nt))
grid.mask = np.zeros((nlat, nlon, nt), dtype=bool)
for i,grace_month in enumerate(GRACE_Ylms.month):
# GRACE/GRACE-FO harmonics for time t
# and monthly files to be removed
if widgets.destripe.value:
Ylms = GRACE_Ylms.index(i).destripe()
Ylms.subtract(remove_Ylms.index(i).destripe())
else:
Ylms = GRACE_Ylms.index(i)
Ylms.subtract(remove_Ylms.index(i))
# Remove GIA rate for time
Ylms.subtract(GIA_Ylms.index(i))
# smooth harmonics and convert to output units
Ylms.convolve(dfactor*wt)
# convert spherical harmonics to output spatial grid
grid.data[:,:,i] = gravtk.harmonic_summation(Ylms.clm, Ylms.slm,
grid.lon, grid.lat, LMAX=LMAX, MMAX=MMAX, PLM=PLM).T
Output spatial data to file
# output to netCDF4 or HDF5
suffix = dict(netCDF4='nc',HDF5='H5')
file_format = '{0}_{1}_{2}{3}{4}_{5}_L{6:d}{7}{8}{9}_{10:03d}-{11:03d}.{12}'
if widgets.format in ('netCDF4','HDF5'):
FILE = file_format.format(PROC,DREL,DSET,gia_str,GRACE_Ylms.title,
widgets.units.value,LMAX,order_str,gw_str,ds_str,
months[0],months[-1],suffix[widgets.format])
grid.to_file(GRACE_Ylms.directory.joinpath(FILE),
format=widgets.format, varname='z',
units=widgets.units.value, longname=unit_name[UNITS-1],
title='GRACE/GRACE-FO Spatial Data', date=True)
Create animation of GRACE/GRACE-FO months
# slider for the plot min and max for normalization
vmin = np.min(grid.data).astype(np.int64)
vmax = np.ceil(np.max(grid.data)).astype(np.int64)
cmap1 = gravtk.tools.colormap(vmin=vmin, vmax=vmax)
# display widgets for setting GRACE/GRACE-FO regression plot parameters
ipywidgets.VBox([cmap1.range,cmap1.step,cmap1.name,cmap1.reverse])
%matplotlib inline
fig, ax1 = plt.subplots(num=1, nrows=1, ncols=1, figsize=(10.375,6.625),
subplot_kw=dict(projection=ccrs.PlateCarree()))
# levels and normalization for plot range
im = ax1.imshow(np.zeros((nlat, nlon)), interpolation='nearest',
norm=cmap1.norm, cmap=cmap1.value, transform=ccrs.PlateCarree(),
extent=grid.extent, origin='upper', animated=True)
ax1.coastlines('50m')
# add date label (year-calendar month e.g. 2002-01)
time_text = ax1.text(0.025, 0.025, '', transform=fig.transFigure,
color='k', size=24, weight='bold', ha='left', va='baseline')
# Add horizontal colorbar and adjust size
# extend = add extension triangles to upper and lower bounds
# options: neither, both, min, max
# pad = distance from main plot axis
# shrink = percent size of colorbar
# aspect = lengthXwidth aspect of colorbar
cbar = plt.colorbar(im, ax=ax1, extend='both', extendfrac=0.0375,
orientation='horizontal', pad=0.025, shrink=0.85,
aspect=22, drawedges=False)
# rasterized colorbar to remove lines
cbar.solids.set_rasterized(True)
# Add label to the colorbar
cbar.ax.set_xlabel(unit_name[UNITS-1], labelpad=10, fontsize=24)
cbar.ax.set_ylabel(unit_label[UNITS-1], fontsize=24, rotation=0)
cbar.ax.yaxis.set_label_coords(1.045, 0.1)
# Set the tick levels for the colorbar
cbar.set_ticks(cmap1.levels)
cbar.set_ticklabels(cmap1.label)
# ticks lines all the way across
cbar.ax.tick_params(which='both', width=1, length=26, labelsize=24,
direction='in')
# stronger linewidth on frame
ax1.spines['geo'].set_linewidth(2.0)
ax1.spines['geo'].set_capstyle('projecting')
# adjust subplot within figure
fig.patch.set_facecolor('white')
fig.subplots_adjust(left=0.02,right=0.98,bottom=0.05,top=0.98)
# animate frames
def animate_frames(i):
# set image
im.set_data(grid.data[:,:,i])
# add date label (year-calendar month e.g. 2002-01)
year,month = gravtk.time.grace_to_calendar(grid.month[i])
time_text.set_text(u'{0:4d}\u2013{1:02d}'.format(year,month))
# set animation
anim = animation.FuncAnimation(fig, animate_frames, frames=nt)
plt.close()
HTML(anim.to_jshtml())
Time Series Analysis
Using ordinary least-squares to fit a synthetic signal to the GRACE/GRACE-FO derived spatial fields
Can fit both polynomial and cyclical (seasonal) components
Can be used to analyze both the seasonal and long-term change in regional mass distributions
style = {'description_width': 'initial'}
orderText = ipywidgets.BoundedIntText(
value=1,
min=0,
max=4,
step=1,
description='Polynomial Order:',
disabled=False,
style=style,
)
# cyclical options
cyclicLabel = ipywidgets.Label('Cyclical Terms:')
cyclicCheckbox = {}
for key in ['Annual','Semi-Annual']:
cyclicCheckbox[key] = ipywidgets.Checkbox(
value=True,
description=key,
disabled=False,
)
cyclic = ipywidgets.HBox([cyclicLabel,*cyclicCheckbox.values()])
# custom fit terms
termsLabel = ipywidgets.Label('Fit Terms:')
termsCheckbox = ipywidgets.Checkbox(
value=True,
description='S2 Tide',
disabled=False,
)
terms = ipywidgets.HBox([termsLabel,termsCheckbox])
# display widgets for setting GRACE/GRACE-FO regression parameters
ipywidgets.VBox([orderText,cyclic,terms])
# build list of regression fit components
ORDER = orderText.value
PHASES = {'Annual':1.0,'Semi-Annual':0.5}
CYCLES = [v for k,v in PHASES.items() if cyclicCheckbox[k].value]
TERMS = []
if termsCheckbox.value:
TERMS.extend(gravtk.time_series.aliasing_terms(grid.time))
# total number of fit terms
ncomp = (ORDER + 1) + 2*len(CYCLES) + len(TERMS)
# Allocating memory for output variables
out = gravtk.spatial(spacing=grid.spacing, nlon=nlon, nlat=nlat,
extent=grid.extent, fill_value=grid.fill_value)
out.data = np.zeros((nlat, nlon, ncomp))
# update mask and dimensions
out.update_mask()
# calculate the regression coefficients
for i in range(nlat):
for j in range(nlon):
# Calculating the regression coefficients
tsbeta = gravtk.time_series.regress(grid.time, grid.data[i,j,:],
ORDER=ORDER, CYCLES=CYCLES, TERMS=TERMS)
# save regression components
for k in range(0, ncomp):
out.data[i,j,k] = tsbeta['beta'][k]
Set parameters for creating regression plot
This step specifies the regression variable to plot, the colormap, the normalization for the contour plot colors and the step for the color bar.
# strings for polynomial terms
if (ORDER == 0):# Mean
variable_longname = ['Mean']
elif (ORDER == 1):# Trend
variable_longname = ['Constant','Trend']
elif (ORDER == 2):# Quadratic
variable_longname = ['Constant','Linear','Quadratic']
unit_suffix = [' yr$^{{{0:d}}}$'.format(-o) if o else '' for o in range(ORDER+1)]
# strings for cyclical terms
cyclic_longname = {}
cyclic_longname['Annual'] = ['Annual Sine', 'Annual Cosine']
cyclic_longname['Semi-Annual'] = ['Semi-Annual Sine', 'Semi-Annual Cosine']
# strings for custom fit terms
terms_longname = {}
terms_longname['S2 Tide (GRACE)'] = ['S2 Tidal Alias Sine (GRACE)', 'S2 Tidal Alias Cosine (GRACE)']
terms_longname['S2 Tide (GRACE-FO)'] = ['S2 Tidal Alias Sine (GRACE-FO)', 'S2 Tidal Alias Cosine (GRACE-FO)']
# combined strings for all components
variable_longname.extend([i for k,v in cyclic_longname.items() for i in v if cyclicCheckbox[k].value])
unit_suffix.extend(['' for k,v in cyclic_longname.items() for i in v if cyclicCheckbox[k].value])
variable_longname.extend([i for k,v in terms_longname.items() for i in v if termsCheckbox.value])
unit_suffix.extend(['' for k,v in terms_longname.items() for i in v if termsCheckbox.value])
# variable of interest
variableDropdown = ipywidgets.Dropdown(
options=variable_longname,
value=variable_longname[ORDER],
description='Variable:',
disabled=False,
)
# slider for the plot min and max for normalization
i = variableDropdown.index
vmin = np.min(out.data[:,:,i]).astype(np.int64)
vmax = np.ceil(np.max(out.data[:,:,i])).astype(np.int64)
cmap2 = gravtk.tools.colormap(vmin=vmin, vmax=vmax)
# set range and step size for variable
def set_range_and_step(sender):
i = variableDropdown.index
cmin = np.min(out.data[:,:,i]).astype(np.int64)
cmax = np.ceil(np.max(out.data[:,:,i])).astype(np.int64)
cmap2.range.min = cmin
cmap2.range.max = cmax
cmap2.range.value = [cmin,cmax]
cmap2.step.max = cmax - cmin
# watch variable widget for changes
variableDropdown.observe(set_range_and_step)
# display widgets for setting GRACE/GRACE-FO regression plot parameters
ipywidgets.VBox([variableDropdown,cmap2.range,cmap2.step,cmap2.name,cmap2.reverse])
Create plot of regression outputs from GRACE/GRACE-FO data
fig, ax2 = plt.subplots(num=2, nrows=1, ncols=1, figsize=(10.375,6.625),
subplot_kw=dict(projection=ccrs.PlateCarree()))
# levels and normalization for plot range
i = variableDropdown.index
im = ax2.imshow(out.data[:,:,i], interpolation='nearest',
norm=cmap2.norm, cmap=cmap2.value, transform=ccrs.PlateCarree(),
extent=grid.extent, origin='upper')
ax2.coastlines('50m')
# Add horizontal colorbar and adjust size
# extend = add extension triangles to upper and lower bounds
# options: neither, both, min, max
# pad = distance from main plot axis
# shrink = percent size of colorbar
# aspect = lengthXwidth aspect of colorbar
cbar = plt.colorbar(im, ax=ax2, extend='both', extendfrac=0.0375,
orientation='horizontal', pad=0.025, shrink=0.85,
aspect=22, drawedges=False)
# rasterized colorbar to remove lines
cbar.solids.set_rasterized(True)
# Add label to the colorbar
lbl = f'{unit_name[UNITS-1]} [{unit_label[UNITS-1]}{unit_suffix[i]}]'
cbar.ax.set_xlabel(lbl, labelpad=10, fontsize=24)
# Set the tick levels for the colorbar
cbar.set_ticks(cmap2.levels)
cbar.set_ticklabels(cmap2.label)
# ticks lines all the way across
cbar.ax.tick_params(which='both', width=1, length=26, labelsize=24,
direction='in')
# stronger linewidth on frame
ax2.spines['geo'].set_linewidth(2.0)
ax2.spines['geo'].set_capstyle('projecting')
# adjust subplot within figure
fig.patch.set_facecolor('white')
fig.subplots_adjust(left=0.02,right=0.98,bottom=0.05,top=0.98)
plt.show()