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.value
    # 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)
# inspect GRACE harmonics
GRACE_Ylms

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.0)
# 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('{0:4d}\u2013{1:02d}'.format(year, month))


# set animation
anim = animation.FuncAnimation(fig, animate_frames, frames=nt)
plt.close()
HTML(anim.to_jshtml())