#!/usr/bin/env python
"""
mascons.py
Written by Tyler Sutterley (07/2026)
Conversion routines for publicly available GRACE/GRACE-FO mascon solutions
PYTHON DEPENDENCIES:
numpy: Scientific Computing Tools For Python (https://numpy.org)
REFERENCES:
grid2mascon.m written by Felix Landerer and David Wiese (JPL)
mascon2grid.m written by Felix Landerer and David Wiese (JPL)
UPDATE HISTORY:
Updated 07/2026: use np.radians to convert from degrees to radians
Updated 03/2023: improve typing for variables in docstrings
Updated 11/2022: use lowercase keyword arguments
Updated 04/2022: updated docstrings to numpy documentation format
Updated 10/2021: publicly released version
Updated 09/2019: check number of latitude points for using regional grids
Updated 08/2018: add GSFC grid and mascon conversion routines
extract number of mascons and number of variables
use lat_bound and lon_bound variables as inputs to JPL conversion
Updated 06/2016: check that mat file exists before attempting conversion
Updated 03/2016: minor clean up of mscdata calculation, update comments
Updated 12/2015: added TRANSPOSE option to output spatial routines
Written 07/2013
"""
import copy
import warnings
import numpy as np
[docs]
def to_gsfc(gdata, lon, lat, lon_center, lat_center, lon_span, lat_span):
"""
Converts an input gridded field to an output GSFC mascon array
:cite:p:`Luthcke:2013ep`
Parameters
----------
gdata: np.ndarray
gridded data
lon: np.ndarray
column vector of defined longitude points
lat: np.ndarray
column vector of defined latitude points
lon_center: np.ndarray
mascon longitudinal center points
lat_center: np.ndarray
mascon latitudinal center points
lon_span: np.ndarray
mascon longitudinal central angles
lat_span: np.ndarray
mascon latitudinal central angles
Returns
-------
data: np.ndarray
row vector of mascons
lat_center: np.ndarray
row vector of latitude values for mascon centers
lon_center: np.ndarray
row vector of longitude values for mascon centers
"""
# number of mascons
nmas = len(lon_center)
# convert mascon centers to -180:180
lon_center = np.where(lon_center > 180, lon_center - 360.0, lon_center)
# remove singleton dimensions
lat = np.squeeze(lat)
lon = np.squeeze(lon)
# for mascons centered on 180: use 0:360
alon = np.copy(lon)
alon = np.where(alon < 0, alon + 360.0, alon)
# loop over each mascon bin and average gdata with the cos-lat weights
# for that bin
mascon_array = {}
mascon_array['data'] = np.zeros((nmas))
mascon_array['lon_center'] = np.zeros((nmas))
mascon_array['lat_center'] = np.zeros((nmas))
for k in range(0, nmas):
# create latitudinal and longitudinal bounds for mascon k
if (lat_center[k] == 90.0) | (lat_center[k] == -90.0):
# NH and SH polar mascons
lon_bound = [0.0, 360.0]
lat_bound = lat_center[k] + np.array([-1.0, 1.0]) * lat_span[k]
else:
# convert from mascon centers to mascon bounds
lon_bound = lon_center[k] + np.array([-0.5, 0.5]) * lon_span[k]
lat_bound = lat_center[k] + np.array([-0.5, 0.5]) * lat_span[k]
# if mascon is centered on +/-180: use 0:360
if (lon_bound[0] <= 180.0) & (lon_bound[1] >= 180.0):
ilon = alon.copy()
elif (lon_bound[0] <= -180.0) & (lon_bound[1] >= -180.0):
lon_bound += 360.0
ilon = alon.copy()
else:
ilon = lon.copy()
# indices for grid points within the mascon
(I,) = np.flatnonzero((lat >= lat_bound[0]) & (lat < lat_bound[1]))
(J,) = np.flatnonzero((ilon >= lon_bound[0]) & (ilon < lon_bound[1]))
I, J = (I[np.newaxis, :], J[:, np.newaxis])
# calculate average data for mascon bin
mascon_array['data'][k] = np.mean(
(np.cos(np.radians(lat[I])) / np.mean(np.cos(np.radians(lat[I]))))
* gdata[I, J]
/ len(I)
)
mascon_array['lat_center'][k] = lat_center[k]
mascon_array['lon_center'][k] = lon_center[k]
# return python dictionary with the mascon array data, lon and lat
return mascon_array
[docs]
def to_jpl(gdata, lon, lat, lon_bound, lat_bound):
"""
Converts an input gridded field to an output JPL mascon array
:cite:p:`Watkins:2015jl`
Parameters
----------
gdata: np.ndarray
gridded data
lon: np.ndarray
column vector of defined longitude points
lat: np.ndarray
column vector of defined latitude points
lon_bound: np.ndarray
mascon longitudinal bounds from coordinate file
lat_bound: np.ndarray
mascon latitudinal bounds from coordinate file
Returns
-------
data: np.ndarray
row vector of mascons
mask: np.ndarray
row vector of mask values showing if mascon has no data
lat: np.ndarray
row vector of latitude values for mascons
lon: np.ndarray
row vector of longitude values for mascons
"""
# mascon dimensions
nmas, nvar = lat_bound.shape
# remove singleton dimensions
lat = np.squeeze(lat)
lon = np.squeeze(lon)
# loop over each mascon bin and average gdata with the cos-lat weights
# for that bin
mascon_array = {}
mascon_array['data'] = np.zeros((nmas))
mascon_array['mask'] = np.zeros((nmas), dtype=bool)
mascon_array['lon'] = np.zeros((nmas))
mascon_array['lat'] = np.zeros((nmas))
for k in range(0, nmas):
# indices for grid points within the mascon
(I,) = np.flatnonzero(
(lat >= lat_bound[k, 1]) & (lat < lat_bound[k, 0])
)
(J,) = np.flatnonzero(
(lon >= lon_bound[k, 0]) & (lon < lon_bound[k, 2])
)
nlt = np.count_nonzero(
(lat >= lat_bound[k, 1]) & (lat < lat_bound[k, 0])
)
I, J = (I[np.newaxis, :], J[:, np.newaxis])
# calculate average data for mascon bin
mascon_array['data'][k] = np.mean(
(np.cos(np.radians(lat[I])) / np.mean(np.cos(np.radians(lat[I]))))
* gdata[I, J]
/ nlt
)
# calculate coordinates of mascon center
mascon_array['lat'][k] = (lat_bound[k, 1] + lat_bound[k, 0]) / 2.0
mascon_array['lon'][k] = (lon_bound[k, 1] + lon_bound[k, 2]) / 2.0
mascon_array['mask'][k] = bool(nlt == 0)
# Do a check at the poles to make the lat/lon equal to +/-90/0
if np.abs(lat_bound[k, 0]) == 90:
mascon_array['lat'][k] = lat_bound[k, 0]
mascon_array['lon'][k] = 0.0
if np.abs(lat_bound[k, 1]) == 90:
mascon_array['lat'][k] = lat_bound[k, 1]
mascon_array['lon'][k] = 0.0
# replace invalid data with 0
mascon_array['data'][mascon_array['mask']] = 0.0
# return python dictionary with the mascon array data, lon and lat
return mascon_array
[docs]
def from_gsfc(
mscdata, grid_spacing, lon_center, lat_center, lon_span, lat_span, **kwargs
):
"""
Converts an input GSFC mascon array to an output gridded field
:cite:p:`Luthcke:2013ep`
Parameters
----------
mscdata: np.ndarray
row vector of mascons
grid_spacing: np.ndarray
spacing of the lat/lon grid
lon_center: float
mascon np.ndarray center points
lat_center: np.ndarray
mascon latitudinal center points
lon_span: np.ndarray
mascon longitudinal central angles
lat_span: np.ndarray
mascon latitudinal central angles
transpose: bool, default False
transpose output matrix (lon/lat)
Returns
-------
mdata: np.ndarray
distributed mass grid
"""
# set default keyword arguments
kwargs.setdefault('transpose', False)
# raise warnings for deprecated keyword arguments
deprecated_keywords = dict(TRANSPOSE='transpose')
for old, new in deprecated_keywords.items():
if old in kwargs.keys():
warnings.warn(
f"""Deprecated keyword argument {old}.
Changed to '{new}'""",
DeprecationWarning,
)
# set renamed argument to not break workflows
kwargs[new] = copy.copy(kwargs[old])
# number of mascons
nmas = len(lon_center)
# convert mascon centers to -180:180
lon_center = np.where(lon_center > 180, lon_center - 360.0, lon_center)
# Define output latitude and longitude grids
lon = np.arange(
-180.0 + grid_spacing / 2.0, 180.0 + grid_spacing / 2.0, grid_spacing
)
lat = np.arange(
90.0 - grid_spacing / 2.0, -90.0 - grid_spacing / 2.0, -grid_spacing
)
nlon, nlat = (len(lon), len(lat))
# for mascons centered on 180: use 0:360
alon = np.copy(lon)
alon = np.where(alon < 0, alon + 360.0, alon)
# loop over each mascon bin and assign value to grid points inside bin:
mdata = np.zeros((nlat, nlon))
for k in range(0, nmas):
# create latitudinal and longitudinal bounds for mascon k
if (lat_center[k] == 90.0) | (lat_center[k] == -90.0):
# NH and SH polar mascons
lon_bound = [0.0, 360.0]
lat_bound = lat_center[k] + np.array([-1.0, 1.0]) * lat_span[k]
else:
# convert from mascon centers to mascon bounds
lon_bound = lon_center[k] + np.array([-0.5, 0.5]) * lon_span[k]
lat_bound = lat_center[k] + np.array([-0.5, 0.5]) * lat_span[k]
# if mascon is centered on +/-180: use 0:360
if (lon_bound[0] <= 180.0) & (lon_bound[1] >= 180.0):
ilon = alon.copy()
elif (lon_bound[0] <= -180.0) & (lon_bound[1] >= -180.0):
lon_bound += 360.0
ilon = alon.copy()
else:
ilon = lon.copy()
# indices for grid points within the mascon
(I,) = np.flatnonzero((lat >= lat_bound[0]) & (lat < lat_bound[1]))
(J,) = np.flatnonzero((ilon >= lon_bound[0]) & (ilon < lon_bound[1]))
I, J = (I[np.newaxis, :], J[:, np.newaxis])
mdata[I, J] = mscdata[k]
# return array
if kwargs['transpose']:
return mdata.T
else:
return mdata
[docs]
def from_jpl(mscdata, grid_spacing, lon_bound, lat_bound, **kwargs):
"""
Converts an input JPL mascon array to an output gridded field
:cite:p:`Watkins:2015jl`
Parameters
----------
mscdata: np.ndarray
row vector of mascons
grid_spacing: np.ndarray
spacing of lat/lon grid
lon_bound: np.ndarray
mascon longitudinal bounds from coordinate file
lat_bound: np.ndarray
mascon latitudinal bounds from coordinate file
transpose: bool, default False
transpose output matrix (lon/lat)
Returns
-------
mdata: np.ndarray
distributed mass grid
"""
# set default keyword arguments
kwargs.setdefault('transpose', False)
# raise warnings for deprecated keyword arguments
deprecated_keywords = dict(TRANSPOSE='transpose')
for old, new in deprecated_keywords.items():
if old in kwargs.keys():
warnings.warn(
f"""Deprecated keyword argument {old}.
Changed to '{new}'""",
DeprecationWarning,
)
# set renamed argument to not break workflows
kwargs[new] = copy.copy(kwargs[old])
# mascon dimensions
nmas, nvar = lat_bound.shape
# Define latitude and longitude grids
# output lon will not include 360
# output lat will not include 90
lon = np.arange(
grid_spacing / 2.0, 360.0 + grid_spacing / 2.0, grid_spacing
)
lat = np.arange(
-90.0 + grid_spacing / 2.0, 90.0 + grid_spacing / 2.0, grid_spacing
)
nlon, nlat = (len(lon), len(lat))
# loop over each mascon bin and assign value to grid points inside bin:
mdata = np.zeros((nlat, nlon))
for k in range(0, nmas):
(I,) = np.flatnonzero(
(lat >= lat_bound[k, 1]) & (lat < lat_bound[k, 0])
)
(J,) = np.flatnonzero(
(lon >= lon_bound[k, 0]) & (lon < lon_bound[k, 2])
)
I, J = (I[np.newaxis, :], J[:, np.newaxis])
mdata[I, J] = mscdata[k]
# return array
if kwargs['transpose']:
return mdata.T
else:
return mdata