GRACE/GRACE-FO Harmonic Visualization Program

This (notebook) demonstrates visualizing spherical harmonics from the Gravity Recovery and Climate Experiment (GRACE) and the GRACE Follow-On (GRACE-FO) missions

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 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)

# 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 Visualize Harmonics

These parameters specify corrections and filtering steps

  • GIA Correction

  • Remove Specific Harmonic Fields

  • Redistribute Removed Fields over the Ocean

  • Gaussian Smoothing Radius in kilometers

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

  • Output spherical harmonic units

    1. equivalent water thickness (cm)

    2. geoid height (mm)

# update widgets
widgets.select_corrections(units=['cmwe', 'mmGH'])
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.units])

Convert GRACE/GRACE-FO harmonics to output units

This step extracts the parameters chosen above and then converts the GRACE/GRACE-FO harmonics to the specified units applying the specified corrections and filtering procedures

  • 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

# 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
dfactor = factors.get(gravtk.units.bycode(UNITS))
# units strings for output files and plots
unit_label = ['cm', 'mm']
unit_name = ['Equivalent Water Thickness','Geoid Height']

# converting harmonics to truncated, smoothed coefficients in units
if widgets.destripe.value:
    Ylms = GRACE_Ylms.destripe()
    Ylms.subtract(remove_Ylms.destripe())
else:
    Ylms = GRACE_Ylms.copy()
    Ylms.subtract(remove_Ylms)
# Remove GIA estimate for month
Ylms.subtract(GIA_Ylms)
# smooth harmonics and convert to output units
Ylms.convolve(dfactor*wt)
# create merged masked array
triangle = Ylms.to_masked_array()

Set parameters for creating animation

This step specifies the colormap.

# display widgets for setting GRACE/GRACE-FO regression plot parameters
cmap = gravtk.tools.colormap(vmin=-1, vmax=1)
ipywidgets.VBox([cmap.name,cmap.reverse])

Create animation of GRACE/GRACE-FO months

%matplotlib inline
# plot spherical harmonics for each month
fig, ax1 = plt.subplots(num=1, figsize=(8,4))

# levels and normalization for plot range
cmap.value.set_bad('lightgrey',1.)
# imshow = show image (interpolation nearest for blocks)
im = ax1.imshow(np.ma.zeros((LMAX+1,LMAX+1)), interpolation='nearest',
    cmap=cmap.value, extent=(-LMAX,LMAX,LMAX,0), animated=True)
# Z color limit between -1 and 1
im.set_clim(-1.0,1.0)

# 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 text to label Slm side and Clm side
t1 = ax1.text(0.39, 0.92, '$S_{lm}$', size=24, weight='bold', 
    transform=ax1.transAxes, ha="center", va="center")
t2 = ax1.text(0.61, 0.92, '$C_{lm}$', size=24, weight='bold', 
    transform=ax1.transAxes, ha="center", va="center")
# add x and y labels
ax1.set_ylabel('Degree [l]', fontsize=13)
ax1.set_xlabel('Order [m]', fontsize=13)
ax1.tick_params(axis='both', which='both',
    labelsize=13, direction='in')

# 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='vertical', pad=0.025, shrink=0.85,
    aspect=15, drawedges=False)
# rasterized colorbar to remove lines
cbar.solids.set_rasterized(True)
# Add label to the colorbar
cbar.ax.set_ylabel(unit_name[UNITS-1], labelpad=5, fontsize=13)
cbar.ax.set_xlabel(unit_label[UNITS-1], fontsize=13, rotation=0)
cbar.ax.xaxis.set_label_coords(0.5,1.065)
# ticks lines all the way across
cbar.ax.tick_params(which='both', width=1, length=15, labelsize=13,
    direction='in')
    
# stronger linewidth on frame
[i.set_linewidth(2.0) for i in ax1.spines.values()]
# adjust subplot within figure
fig.patch.set_facecolor('white')
fig.subplots_adjust(left=0.075,right=0.99,bottom=0.07,top=0.99)

# animate frames
def animate_frames(i):
    # set image
    im.set_data(triangle[:,:,i])
    # add date label (year-calendar month e.g. 2002-01)
    year,month = gravtk.time.grace_to_calendar(Ylms.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())