GRACE/GRACE-FO Error Visualization Program

This (notebook) demonstrates visualizing spherical harmonic errors from the Gravity Recovery and Climate Experiment (GRACE) and the GRACE Follow-On (GRACE-FO) Level-2 products.

Load necessary modules for running the notebook

import numpy as np
import matplotlib
matplotlib.rcParams['mathtext.default'] = 'regular'
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import ipywidgets
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)

# 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)
# number of time steps
nt = len(months)

Set Parameters to Convert to Spatial Error Maps

These parameters specify corrections and filtering steps for converting to the spatial domain at a specified grid spacing

  • Gaussian Smoothing Radius in kilometers

  • Filter (destripe) harmonics (Swenson and Wahr, 2006)

  • Spatial degree spacing

  • Spatial degree interval

    1. (-180:180,90:-90)

    2. (degree spacing)/2

# update widgets
widgets.select_corrections()
# display widgets for setting GRACE/GRACE-FO corrections parameters
ipywidgets.VBox([
    widgets.gaussian,
    widgets.destripe,
    widgets.spacing,
    widgets.interval])

Estimate GRACE/GRACE-FO errors and convert to spatial domain

  • Set output grid domain

  • Calculate Fully-Normalized Legendre Polynomials

  • Calculate coefficients for converting to the output units

  • Convert errors from the spherical harmonic domain into the spatial domain

# Output spatial data
grid = gravtk.spatial()
# 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))
# square of legendre polynomials truncated to order MMAX
mm = np.arange(0,MMAX+1)
PLM2 = PLM[:,mm,:]**2
# Calculating cos(m*phi)^2 and sin(m*phi)^2
phi = grid.lon[np.newaxis,:]*np.pi/180.0
ccos = np.cos(np.dot(mm[:,np.newaxis],phi))**2
ssin = np.sin(np.dot(mm[:,np.newaxis],phi))**2
    
# 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')

# gaussian smoothing radius in km (Jekeli, 1981)
RAD = widgets.gaussian.value
if (RAD != 0):
    wt = 2.0*np.pi*gravtk.gauss_weights(RAD,LMAX)
else:
    # else = 1
    wt = np.ones((LMAX+1))

# destriping the GRACE/GRACE-FO harmonics
if widgets.destripe.value:
    Ylms = GRACE_Ylms.destripe()
else:
    Ylms = GRACE_Ylms.copy()

# dfactor is the degree dependent coefficients
# for converting to spherical harmonic output units
factors = gravtk.units(lmax=LMAX).harmonic(hl,kl,ll).mmwe
# mmwe, millimeters water equivalent
dfactor = factors.get('mmwe')
# units strings for output plots
unit_label = 'mm'
unit_name = 'Equivalent Water Thickness'

# Delta coefficients of GRACE time series (Error components)
delta_Ylms = gravtk.harmonics(lmax=LMAX,mmax=MMAX)
delta_Ylms.clm = np.zeros((LMAX+1, MMAX+1))
delta_Ylms.slm = np.zeros((LMAX+1, MMAX+1))
# Smoothing Half-Width (CNES is a 10-day solution)
# All other solutions are monthly solutions (HFWTH for annual = 6)
if ((PROC == 'CNES') and (DREL in ('RL01','RL02'))):
    HFWTH = 19
else:
    HFWTH = 6
# Equal to the noise of the smoothed time-series
# for each spherical harmonic order
for m in range(0,MMAX+1):# MMAX+1 to include MMAX
    # for each spherical harmonic degree
    for l in range(m,LMAX+1):# LMAX+1 to include LMAX
        # Delta coefficients of GRACE time series
        for cs,csharm in enumerate(['clm','slm']):
            # calculate GRACE Error (Noise of smoothed time-series)
            # With Annual and Semi-Annual Terms
            val1 = getattr(Ylms, csharm)
            smth = gravtk.time_series.smooth(Ylms.time, val1[l,m,:],
                HFWTH=HFWTH)
            # number of smoothed points
            nsmth = len(smth['data'])
            tsmth = np.mean(smth['time'])
            # GRACE delta Ylms
            # variance of data-(smoothed+annual+semi)
            val2 = getattr(delta_Ylms, csharm)
            val2[l,m] = np.sqrt(np.sum(smth['noise']**2)/nsmth)
                
# convolve delta harmonics with degree dependent factors
delta_Ylms = delta_Ylms.convolve(dfactor*wt)
# smooth harmonics and convert to output units
YLM2 = delta_Ylms.power(2.0).scale(1.0/nsmth)
# Calculate fourier coefficients
d_cos = np.zeros((MMAX+1,nlat))# [m,th]
d_sin = np.zeros((MMAX+1,nlat))# [m,th]
# Calculating delta spatial values
for k in range(0,nlat):
    # summation over all spherical harmonic degrees
    d_cos[:,k] = np.sum(PLM2[:,:,k]*YLM2.clm, axis=0)
    d_sin[:,k] = np.sum(PLM2[:,:,k]*YLM2.slm, axis=0)

# Multiplying by c/s(phi#m) to get spatial maps (lon,lat)
grid.data = np.sqrt(np.dot(ccos.T,d_cos) + np.dot(ssin.T,d_sin)).T
grid.mask = np.zeros_like(grid.data, dtype=bool)

Create plot of GRACE/GRACE-FO degree spectra

fig, ax1 = plt.subplots(num=1, nrows=1, ncols=1)
ax1.plot(delta_Ylms.l[1:], delta_Ylms.amplitude[1:], color='red', linewidth=2)
ax1.set_xlabel('Degree [l]', fontsize=13)
ax1.set_ylabel(f'{unit_name} [{unit_label}]', fontsize=13)
ax1.set_yscale('log')
ax1.set_xlim(0, LMAX)
ax1.grid(True, which='both', linestyle='-', color='grey')
ax1.set_title('GRACE/GRACE-FO Error Degree Amplitude')
plt.show()

Create GRACE/GRACE-FO error map

# slider for the plot min and max for normalization
vmax = np.ceil(np.max(grid.data)).astype(np.int64)
cmap = gravtk.tools.colormap(vmin=0, vmax=vmax)
# display widgets for setting GRACE/GRACE-FO plot parameters
ipywidgets.VBox([cmap.range,cmap.step,cmap.name,cmap.reverse])

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
im = ax2.imshow(grid.data, interpolation='nearest',
    norm=cmap.norm, cmap=cmap.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
cbar.ax.set_xlabel(f'{unit_name} [{unit_label}]',
    labelpad=10, fontsize=24)
# Set the tick levels for the colorbar
cbar.set_ticks(cmap.levels)
cbar.set_ticklabels(cmap.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()