#!/usr/bin/env python
u"""
C20.py
Written by Tyler Sutterley (05/2023)
Reads in C20 spherical harmonic coefficients derived from SLR measurements
Dataset distributed by NASA PO.DAAC
TN-05_C20_SLR.txt
TN-07_C20_SLR.txt
TN-11_C20_SLR.txt
TN-14_C30_C30_GSFC_SLR.txt
Dataset distributed by UTCSR
C20_RL05.txt
Datasets distributed by GFZ
GFZ_RL06_C20_SLR.dat
GRAVIS-2B_GFZOP_GRACE+SLR_LOW_DEGREES_0002.dat
CALLING SEQUENCE:
SLR_C20 = gravity_toolkit.SLR.C20(SLR_file)
INPUTS:
SLR_file:
RL04: TN-05_C20_SLR.txt
RL05: TN-07_C20_SLR.txt
RL06: TN-11_C20_SLR.txt
CSR: C20_RL05.txt
GFZ: GFZ_RL06_C20_SLR.dat
GFZ (combined): GRAVIS-2B_GFZOP_GRACE+SLR_LOW_DEGREES_0002.dat
OUTPUTS:
data: SLR degree 2 order 0 cosine stokes coefficients (C20)
error: SLR degree 2 order 0 cosine stokes coefficient error (eC20)
month: GRACE/GRACE-FO month of measurement (April 2002 = 004)
date: date of SLR measurement
OPTIONS:
AOD: remove background De-aliasing product from the SLR solution (for CSR)
HEADER: file contains header text to be skipped (default: True)
PYTHON DEPENDENCIES:
numpy: Scientific Computing Tools For Python
https://numpy.org
https://numpy.org/doc/stable/user/numpy-for-matlab-users.html
dateutil: powerful extensions to datetime
https://dateutil.readthedocs.io/en/stable/
PROGRAM DEPENDENCIES:
time.py: utilities for calculating time operations
REFERENCES:
Cheng, M. and Tapley, B. D., "Variations in the Earth's oblateness during
the past 28 years", Journal of Geophysical Research: Solid Earth,
109(B9), B09402, 2004. 10.1029/2004JB003028
Loomis, B. D., Rachlin, K. E., and Luthcke, S. B., "Improved Earth
Oblateness Rate Reveals Increased Ice Sheet Losses and Mass-Driven Sea
Level Rise", Geophysical Research Letters, 46(12), 6910-6917, 2019.
https://doi.org/10.1029/2019GL082929
Koenig, R., Schreiner, P, and Dahle, C. "Monthly estimates of C(2,0)
generated by GFZ from SLR satellites based on GFZ GRACE/GRACE-FO
RL06 background models." V. 1.0. GFZ Data Services, (2019).
https://doi.org/10.5880/GFZ.GRAVIS_06_C20_SLR
UPDATE HISTORY:
Updated 05/2023: use pathlib to define and operate on paths
Updated 01/2023: refactored satellite laser ranging read functions
Updated 04/2022: updated docstrings to numpy documentation format
include utf-8 encoding in reads to be windows compliant
Updated 09/2021: use functions for converting to and from GRACE months
Updated 05/2021: added GFZ SLR and GravIS oblateness solutions
define int/float precision to prevent deprecation warning
Updated 02/2021: use adjust_months function to fix special months cases
replaced numpy bool to prevent deprecation warning
Updated 12/2020: using utilities from time module
Updated 08/2020: flake8 compatible regular expression strings
Updated 07/2020: added function docstrings
Updated 08/2019: add catch to verify input SLR file exists
Updated 07/2019: added tilde-expansion of input SLR file
Updated 06/2019: added new GRACE-FO special month (October 2018)
Updated 11/2018: new TN-11 files only list GRACE months available
Updated 06/2016: added option HEADER for files that do not have header text
Updated 05/2016: added option AOD to not remove the AOD correction
Updated 03/2016: minor update to read PO.DAAC
Updated 05/2015: minor change to file determination (only regular expressions)
Updated 02/2015: updated UT/CSR portion and comments
Updated 09/2014: rewrite of the TN-07 read program
using regular expressions and convert_calendar_decimal
Updated 01/2014: updated to use UT/CSR monthly time-series
as an alternative to PO.DAAC as it is updated more regularly
Updated 05/2013: adapted for python
Updated 09/2012: Changed month scheme to output.
Used to remove the GRACE missing months in this program by feeding in the GRACE months
BUT, as the new SLR files start with an earlier date, decided to parallel
the degree-1 read program, and remove the missing months in the read_grace program
Updated 06/2012: OVERHAUL of dating and modification for 'special' GRACE months
Initiated from an incorrect date tag in the SLR data file
New dating will convert from the MJD file into date fraction
Some GRACE 'months' have the accelerometer turned off
for half the month to preserve battery power
These months use half of the prior month in the GRACE global gravity solution
For these months the SLR file has a second dataline for the modified period
Will use these marked (*) data to replace the GRACE C2,0
ALSO converted the mon and slrdate inputs into options
Updated 01/2012: Updated to feed in SLR file from outside
Will accommodate upcoming GRACE RL05, which will use different SLR files
Written 12/2011
"""
import re
import pathlib
import numpy as np
import gravity_toolkit.time
# PURPOSE: read oblateness data from Satellite Laser Ranging (SLR)
[docs]
def C20(SLR_file, AOD=True, HEADER=True):
r"""
Reads *C*\ :sub:`20` spherical harmonic coefficients from SLR measurements
Parameters
----------
SLR_file: str
Satellite Laser Ranging file
AOD: bool, default True
Remove background De-aliasing product from the
Center for Space Research (CSR) SLR solutions
HEADER: bool, default True
File contains header text to be skipped
Returns
-------
data: float
SLR degree 2 order 0 cosine stokes coefficients
error: float
SLR degree 2 order 0 cosine stokes coefficient error
month: int
GRACE/GRACE-FO month of measurement
date: float
date of SLR measurement
"""
# check that SLR file exists
SLR_file = pathlib.Path(SLR_file).expanduser().absolute()
if not SLR_file.exists():
raise FileNotFoundError(f'{str(SLR_file)} not found in file system')
# output dictionary with data variables
dinput = {}
# determine if imported file is from PO.DAAC or CSR
if bool(re.search(r'C20_RL\d+', SLR_file.name ,re.I)):
# SLR C20 file from CSR
# Just for checking new months when TN series isn't up to date as the
# SLR estimates always use the full set of days in each calendar month.
# format of the input file (note 64 bit floating point for C20)
# Column 1: Approximate mid-point of monthly solution (years)
# Column 2: C20 from SLR (normalized)
# Column 3: Delta C20 relative to a mean value (1E-10)
# Column 4: Solution sigma (1E-10)
# Column 5: Mean value of Atmosphere-Ocean De-aliasing model (1E-10)
# Columns 6-7: Start and end dates of data used in solution
dtype = {}
dtype['names'] = ('time','C20','delta','sigma','AOD','start','end')
dtype['formats'] = ('f','f8','f','f','f','f','f')
# header text is commented and won't be read
file_input = np.loadtxt(SLR_file, dtype=dtype)
# date and GRACE/GRACE-FO month
dinput['time'] = file_input['time']
dinput['month'] = gravity_toolkit.time.calendar_to_grace(dinput['time'])
# monthly spherical harmonic replacement solutions
dinput['data'] = file_input['C20'].copy()
# monthly spherical harmonic formal standard deviations
dinput['error'] = file_input['sigma']*1e-10
# Background gravity model includes solid earth and ocean tides, solid
# earth and ocean pole tides, and the Atmosphere-Ocean De-aliasing
# product. The monthly mean of the AOD model has been restored.
if AOD:
# Removing AOD product that was restored in the solution
dinput['data'] -= file_input['AOD']*1e-10
elif bool(re.search(r'GFZ_(RL\d+)_C20_SLR', SLR_file.name, re.I)):
# SLR C20 file from GFZ
# Column 1: MJD of BEGINNING of solution span
# Column 2: Year and fraction of year of BEGINNING of solution span
# Column 3: Replacement C(2,0)
# Column 4: Replacement C(2,0) - mean C(2,0) (1.0E-10)
# Column 5: C(2,0) formal error (1.0E-10)
with SLR_file.open(mode='r', encoding='utf8') as f:
file_contents = f.read().splitlines()
# number of lines contained in the file
file_lines = len(file_contents)
# counts the number of lines in the header
count = 0
# Reading over header text
while HEADER:
# file line at count
line = file_contents[count]
# find PRODUCT: within line to set HEADER flag to False when found
HEADER = not bool(re.match(r'PRODUCT:+',line))
# add 1 to counter
count += 1
# number of months within the file
n_mon = file_lines - count
# date and GRACE/GRACE-FO month
dinput['time'] = np.zeros((n_mon))
dinput['month'] = np.zeros((n_mon),dtype=np.int64)
# monthly spherical harmonic replacement solutions
dinput['data'] = np.zeros((n_mon))
# monthly spherical harmonic formal standard deviations
dinput['error'] = np.zeros((n_mon))
# time count
t = 0
# for every other line:
for line in file_contents[count:]:
# find numerical instances in line including exponents,
# decimal points and negatives
line_contents = re.findall(r'[-+]?\d*\.\d*(?:[eE][-+]?\d+)?',line)
# check if line has G* or Gm flags
if bool(re.search(r'(G\*|Gm)',line)):
# reading decimal year for start of span
dinput['time'][t] = np.float64(line_contents[1])
# Spherical Harmonic data for line
dinput['data'][t] = np.float64(line_contents[2])
dinput['error'][t] = np.float64(line_contents[4])*1e-10
# GRACE/GRACE-FO month of SLR solutions
dinput['month'][t] = gravity_toolkit.time.calendar_to_grace(
dinput['time'][t], around=np.round)
# add to t count
t += 1
# truncate variables if necessary
for key,val in dinput.items():
dinput[key] = val[:t]
elif bool(re.search(r'GRAVIS-2B_GFZOP', SLR_file.name, re.I)):
# Combined GRACE/SLR solution file produced by GFZ
# Column 1: MJD of BEGINNING of solution data span
# Column 2: Year and fraction of year of BEGINNING of solution span
# Column 3: Replacement C(2,0)
# Column 4: Replacement C(2,0) - mean C(2,0) (1.0E-10)
# Column 5: C(2,0) formal standard deviation (1.0E-12)
with SLR_file.open(mode='r', encoding='utf8') as f:
file_contents = f.read().splitlines()
# number of lines contained in the file
file_lines = len(file_contents)
# counts the number of lines in the header
count = 0
# Reading over header text
while HEADER:
# file line at count
line = file_contents[count]
# find PRODUCT: within line to set HEADER flag to False when found
HEADER = not bool(re.match(r'PRODUCT:+',line))
# add 1 to counter
count += 1
# number of months within the file
n_mon = file_lines - count
# date and GRACE/GRACE-FO month
dinput['time'] = np.zeros((n_mon))
dinput['month'] = np.zeros((n_mon), dtype=int)
# monthly spherical harmonic replacement solutions
dinput['data'] = np.zeros((n_mon))
# monthly spherical harmonic formal standard deviations
dinput['error'] = np.zeros((n_mon))
# time count
t = 0
# for every other line:
for line in file_contents[count:]:
# find numerical instances in line including exponents,
# decimal points and negatives
line_contents = re.findall(r'[-+]?\d*\.\d*(?:[eE][-+]?\d+)?',line)
count = len(line_contents)
# check for empty lines
if (count > 0):
# reading decimal year for start of span
dinput['time'][t] = np.float64(line_contents[1])
# Spherical Harmonic data for line
dinput['data'][t] = np.float64(line_contents[2])
dinput['error'][t] = np.float64(line_contents[4])*1e-10
# GRACE/GRACE-FO month of SLR solutions
dinput['month'][t] = gravity_toolkit.time.calendar_to_grace(
dinput['time'][t], around=np.round)
# add to t count
t += 1
# truncate variables if necessary
for key,val in dinput.items():
dinput[key] = val[:t]
elif bool(re.search(r'TN-(11|14)', SLR_file.name, re.I)):
# SLR C20 RL06 file from PO.DAAC
with SLR_file.open(mode='r', encoding='utf8') as f:
file_contents = f.read().splitlines()
# number of lines contained in the file
file_lines = len(file_contents)
# counts the number of lines in the header
count = 0
# Reading over header text
while HEADER:
# file line at count
line = file_contents[count]
# find PRODUCT: within line to set HEADER flag to False when found
HEADER = not bool(re.match(r'PRODUCT:+',line,re.IGNORECASE))
# add 1 to counter
count += 1
# number of months within the file
n_mon = file_lines - count
# date and GRACE/GRACE-FO month
dinput['time'] = np.zeros((n_mon))
dinput['month'] = np.zeros((n_mon),dtype=np.int64)
# monthly spherical harmonic replacement solutions
dinput['data'] = np.zeros((n_mon))
# monthly spherical harmonic formal standard deviations
dinput['error'] = np.zeros((n_mon))
# time count
t = 0
# for every other line:
for line in file_contents[count:]:
# find numerical instances in line including exponents,
# decimal points and negatives
line_contents = re.findall(r'[-+]?\d*\.\d*(?:[eE][-+]?\d+)?',line)
# check for empty lines as there are
# slight differences in RL04 TN-05_C20_SLR.txt
# with blanks between the PRODUCT: line and the data
count = len(line_contents)
# if count is greater than 0
if (count > 0):
# modified julian date for line
MJD = np.float64(line_contents[0])
# converting from MJD into month, day and year
YY,MM,DD,hh,mm,ss = gravity_toolkit.time.convert_julian(
MJD+2400000.5, format='tuple')
# converting from month, day, year into decimal year
dinput['time'][t] = gravity_toolkit.time.convert_calendar_decimal(
YY, MM, day=DD, hour=hh)
# Spherical Harmonic data for line
dinput['data'][t] = np.float64(line_contents[2])
dinput['error'][t] = np.float64(line_contents[4])*1e-10
# GRACE/GRACE-FO month of SLR solutions
dinput['month'][t] = gravity_toolkit.time.calendar_to_grace(
dinput['time'][t], around=np.round)
# add to t count
t += 1
# truncate variables if necessary
for key,val in dinput.items():
dinput[key] = val[:t]
else:
# SLR C20 file from PO.DAAC
with SLR_file.open(mode='r', encoding='utf8') as f:
file_contents = f.read().splitlines()
# number of lines contained in the file
file_lines = len(file_contents)
# counts the number of lines in the header
count = 0
# Reading over header text
while HEADER:
# file line at count
line = file_contents[count]
# find PRODUCT: within line to set HEADER flag to False when found
HEADER = not bool(re.match(r'PRODUCT:+',line))
# add 1 to counter
count += 1
# number of months within the file
n_mon = file_lines - count
# GRACE/GRACE-FO dates
date_conv = np.zeros((n_mon))
# monthly spherical harmonic replacement solutions
C20_input = np.zeros((n_mon))
# monthly spherical harmonic formal standard deviations
eC20_input = np.zeros((n_mon))
# flag denoting if replacement solution
slr_flag = np.zeros((n_mon),dtype=bool)
# time count
t = 0
# for every other line:
for line in file_contents[count:]:
# find numerical instances in line including exponents,
# decimal points and negatives
line_contents = re.findall(r'[-+]?\d*\.\d*(?:[eE][-+]?\d+)?',line)
# check for empty lines as there are
# slight differences in RL04 TN-05_C20_SLR.txt
# with blanks between the PRODUCT: line and the data
count = len(line_contents)
# if count is greater than 0
if (count > 0):
# modified julian date for line
MJD = np.float64(line_contents[0])
# converting from MJD into month, day and year
YY,MM,DD,hh,mm,ss = gravity_toolkit.time.convert_julian(
MJD+2400000.5, format='tuple')
# converting from month, day, year into decimal year
date_conv[t] = gravity_toolkit.time.convert_calendar_decimal(
YY, MM, day=DD, hour=hh)
# Spherical Harmonic data for line
C20_input[t] = np.float64(line_contents[2])
eC20_input[t] = np.float64(line_contents[4])*1e-10
# line has * flag
if bool(re.search(r'\*',line)):
slr_flag[t] = True
# add to t count
t += 1
# truncate for RL04 if necessary
date_conv = date_conv[:t]
C20_input = C20_input[:t]
eC20_input = eC20_input[:t]
slr_flag = slr_flag[:t]
# GRACE/GRACE-FO month of SLR solutions
mon = gravity_toolkit.time.calendar_to_grace(date_conv,around=np.round)
# number of unique months
dinput['month'] = np.unique(mon)
n_uniq = len(dinput['month'])
# Removing overlapping months to use the data for
# months with limited GRACE accelerometer use
dinput['time'] = np.zeros((n_uniq))
dinput['data'] = np.zeros((n_uniq))
dinput['error'] = np.zeros((n_uniq))
# New SLR datasets have * flags for the modified GRACE periods
# these GRACE months use half of a prior month in their solution
# this will find these months (marked above with slr_flag)
for t in range(n_uniq):
count = np.count_nonzero(mon == dinput['month'][t])
# there is only one solution for the month
if (count == 1):
i = np.nonzero(mon == dinput['month'][t])
dinput['time'][t] = date_conv[i]
dinput['data'][t] = C20_input[i]
dinput['error'][t] = eC20_input[i]
# there is a special solution for the month
# will the solution flagged with slr_flag
elif (count == 2):
i = np.nonzero((mon == dinput['month'][t]) & slr_flag)
dinput['time'][t] = date_conv[i]
dinput['data'][t] = C20_input[i]
dinput['error'][t] = eC20_input[i]
# The 'Special Months' (Nov 2011, Dec 2011 and April 2012) with
# Accelerometer shutoffs make the relation between month number
# and date more complicated as days from other months are used
# For CSR and GFZ: Nov 2011 (119) is centered in Oct 2011 (118)
# For JPL: Dec 2011 (120) is centered in Jan 2012 (121)
# For all: May 2015 (161) is centered in Apr 2015 (160)
# For GSFC: Oct 2018 (202) is centered in Nov 2018 (203)
dinput['month'] = gravity_toolkit.time.adjust_months(dinput['month'])
# return the SLR-derived oblateness solutions
return dinput