GRACE/GRACE-FO Geostrophic Current Map Program
This (notebook) demonstrates using Gravity Recovery and Climate Experiment (GRACE) and the GRACE Follow-On (GRACE-FO) Level-2 spherical harmonic products to calculate geostrophic currents following Wahr et al. (2002).
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 matplotlib.offsetbox as offsetbox
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
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)
# 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 of Geostrophic Current
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)
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 geostrophic current 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). We additionally need to compensate for the Earth’s elastic yielding to surface load changes, which induce density anomalies at depth within the solid Earth (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). Monthly spatial fields of zonal (\(\nu_\mathrm{zonal}\)) and meridional (\(\nu_\mathrm{merid}\)) geostrophic current will be calculated following Wahr et al. (2002). This method is not accurate for latitudes within 10° of the equator where the geostrophic approximation is not valid.
Zonal Geostrophic Currents
Meridional Geostrophic Currents
\(\Omega_e\): average angular velocity of the Earth
\(g_e\): standard gravitational acceleration of the Earth
\(\rho_w\): average density of seawater at depth
\(\rho_e\): average density of the Earth
\(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])
Convert GRACE/GRACE-FO harmonics to geostrophic currents in 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 spatial data
grid = gravtk.spatial(fill_value=np.nan)
grid.time = np.copy(GRACE_Ylms.time)
grid.month = np.copy(GRACE_Ylms.month)
# Read Smoothed Ocean and Land Functions
# will mask out land regions in the final current maps
LANDMASK = gravtk.utilities.get_data_path(['data','land_fcn_300km.nc'])
landsea = gravtk.spatial().from_netCDF4(LANDMASK,
date=False, varname='LSMASK')
# degree spacing and grid dimensions
# will create GRACE spatial fields with same dimensions
dlon,dlat = landsea.spacing
nlat, nlon = landsea.shape
# shift landsea mask to have longitudes -180:180
landsea.mask, landsea.lon = gravtk.tools.shift_grid(180.0 + dlon,
landsea.mask, landsea.lon, CYCLIC=360)
# grid latitude and longitude
grid.lon = np.copy(landsea.lon)
grid.lat = np.copy(landsea.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))
RAD = widgets.gaussian.value
# 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
LOVE = gravtk.read_love_numbers(love_numbers_file, LMAX=LMAX,
HEADER=header, COLUMNS=columns, REFERENCE='CF', FORMAT='class')
# 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=LOVE)
# 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)
# converting harmonics to truncated, smoothed coefficients in units
# combining harmonics to calculate output spatial fields
# output geostrophic current grid
grid.data = np.zeros((nlat, nlon, 2,nt))
grid.mask = np.ones((nlat, nlon, 2,nt), dtype=bool)
# mask equatorial regions due to hydrostrophic inaccuracies
valid, = np.nonzero((np.abs(grid.lat) > 10))
grid.mask[valid,:,:,:] = False
# set land values from land-sea mask to invalid
indy,indx = np.nonzero(np.logical_not(landsea.mask))
grid.mask[indy,indx,:,:] = True
# for each GRACE/GRACE-FO month
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))
# convert spherical harmonics to output spatial grid
currents = gravtk.geostrophic_currents(Ylms.clm, Ylms.slm,
grid.lon, grid.lat[valid], LMAX=LMAX, MMAX=MMAX,
RAD=RAD, LOVE=LOVE, PLM=PLM)
# transpose to outputs to latxlon
grid.data[valid,:,:,i] = currents.transpose(1,0,2)
# update the mask and replace fill values
grid.update_mask();
Create animation of GRACE/GRACE-FO months
# slider for the plot min and max for normalization
vmin = np.nanmin(grid.data).astype(np.int64)
vmax = np.ceil(np.nanmax(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,ax2) = plt.subplots(num=1, nrows=2, ncols=1, figsize=(10.375,11.625),
sharex=True, sharey=True, subplot_kw=dict(projection=ccrs.PlateCarree()))
# levels and normalization for plot range
im1 = ax1.imshow(np.zeros((nlat, nlon)), interpolation='nearest',
norm=cmap1.norm, cmap=cmap1.value, transform=ccrs.PlateCarree(),
extent=grid.extent, origin='upper', animated=True)
im2 = ax2.imshow(np.zeros((nlat, nlon)), interpolation='nearest',
norm=cmap1.norm, cmap=cmap1.value, transform=ccrs.PlateCarree(),
extent=grid.extent, origin='upper', animated=True)
# add date label (year-calendar month e.g. 2002-01)
time_text = ax1.text(0.025, 0.015, '', transform=fig.transFigure,
color='k', size=24, weight='bold', ha='left', va='baseline')
# Add colorbar
# Add an axes at position rect [left, bottom, width, height]
cbar_ax = fig.add_axes([0.095, 0.075, 0.81, 0.03])
# extend = add extension triangles to upper and lower bounds
# options: neither, both, min, max
cbar = fig.colorbar(im1, cax=cbar_ax, extend='both',
extendfrac=0.0375, drawedges=False, orientation='horizontal')
# rasterized colorbar to remove lines
cbar.solids.set_rasterized(True)
# Add label to the colorbar
cbar.ax.set_title('Geostrophic Current', fontsize=18, rotation=0, y=-1.65, va='top')
cbar.ax.set_xlabel('cm/s', fontsize=18, rotation=0, va='center')
cbar.ax.xaxis.set_label_coords(1.085, 0.5)
# 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=25, labelsize=18,
direction='in')
# add labels, coastlines and adjust frames
labels = ['Zonal', 'Meridional']
for i, ax in enumerate([ax1, ax2]):
# add current label
at = offsetbox.AnchoredText(labels[i],
loc=3, pad=0, borderpad=0.25, frameon=True,
prop=dict(size=24, weight='bold', color='k'))
at.patch.set_boxstyle("Square,pad=0.2")
at.patch.set_edgecolor("white")
ax.axes.add_artist(at)
# add coastlines
ax.coastlines('50m')
# stronger linewidth on frame
ax.spines['geo'].set_linewidth(2.0)
ax.spines['geo'].set_zorder(10)
ax.spines['geo'].set_capstyle('projecting')
# adjust subplot within figure
fig.patch.set_facecolor('white')
fig.subplots_adjust(left=0.01, right=0.99, bottom=0.12, top=0.97,
hspace=0.05, wspace=0.05)
# animate frames
def animate_frames(i):
# set image
im1.set_data(grid.data[:,:,0,i])
im2.set_data(grid.data[:,:,1,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())