lenapy.utils.harmo

The harmo module provides functions for transforming spherical harmonics datasets into spatial grid representations and vice versa.

This module includes functions to:

  • Convert spherical harmonics datasets into spatial grid DataArray.

  • Convert spatial grid DataArray into spherical harmonics datasets.

  • Compute associated Legendre functions.

  • Calculate mid-month estimates for GRACE data products.

  • Change the associated normalization to a spherical harmonics dataset.

  • Compute scaling factors for unit conversions between spherical harmonics and grid data.

  • Validate spherical harmonics and grid datasets.

lenapy.utils.harmo.sh_to_grid(data: Dataset, unit: Literal['mewh', 'mmgeoid', 'microGal', 'bar', 'mvcu', 'norm'] = 'mewh', errors=False, lmax: int | None = None, mmax: int | None = None, lmin: int | None = None, mmin: int | None = None, used_l: ndarray | None = None, used_m: ndarray | None = None, lonmin: float = -180, lonmax: float = 180, latmin: float = -90, latmax: float = 90, bounds: list = None, dlon: float = 1, dlat: float = 1, longitude: ndarray | None = None, latitude: ndarray | None = None, radians_in: bool = False, force_mass_conservation: bool = False, ellipsoidal_earth: bool = False, include_elastic: bool = True, plm: DataArray = None, normalization_plm: Literal['4pi', 'ortho', 'schmidt'] = '4pi', **kwargs) → DataArray[source]

Transform Spherical Harmonics (SH) dataset into spatial DataArray. With choice for constants, unit, love numbers, degree/order, spatial grid latitude and longitude, Earth hypothesis.

For details on unit transformations, see l_factor_conv() and the associated PDF document.

Parameters:

  • data (xr.Dataset) – xr.Dataset that corresponds to SH data to convert into spatial representation.

  • unit (str, optional) – ‘mewh’, ‘mmgeoid’, ‘microGal’, ‘bar’, ‘mvcu’, or ‘norm’ Unit of the spatial data used in the transformation. Default is ‘mewh’ for meters of Equivalent Water Height. See utils.harmo.l_factor_conv() doc for details on the units.

  • errors (bool, optional) – If True, propagate the errors from eclm and eslm variables to create an error grid. Default is False.

  • lmax (int, optional) – Maximal degree of the SH coefficients to be used.

  • mmax (int, optional) – Maximal order of the SH coefficients to be used.

  • lmin (int, optional) – Minimal degree of the SH coefficients to be used.

  • mmin (int, optional) – Minimal order of the SH coefficients to be used.

  • used_l (np.ndarray, optional) – List of degree to use for the grid computation (if given, lmax and lmin are not considered).

  • used_m (np.ndarray, optional) – List of order to use for the grid computation (if given, mmax and mmin are not considered).

  • lonmin (float, optional) – Minimal longitude of the future grid.

  • lonmax (float, optional) – Maximal longitude of the future grid.

  • latmin (float, optional) – Minimal latitude of the future grid.

  • latmax (float, optional) – Maximal latitude of the future grid.

  • bounds (list, optional) – List of 4 elements with [lonmin, lonmax, latmin, latmax] (if given min/max information are not considered).

  • radians_in (bool, optional) – True if the unit of the given latitude and longitude information is radians. Default is False for degree unit. If radians_in is True and dlat or dlon are given, they are considered as radians.

  • dlon (float, optional) – Spacing of the longitude values.

  • dlat (float, optional) – Spacing of the latitude values.

  • longitude (np.ndarray, optional) – List of longitude to use for the grid computation (if given, others longitude information are not considered).

  • latitude (np.ndarray, optional) – List of latitude to use for the grid computation (if given, others latitude information are not considered).

  • force_mass_conservation (bool, optional) – If True, force that the grid resulting from all coefficients except C0 has a null global mass. Default is False

  • ellipsoidal_earth (bool, optional) – If True, consider the Earth as an ellipsoid following [Ditmar2018]. Default is False for a spherical Earth.

  • include_elastic (bool, optional) – If True, the Earth behavior is elastic. Default is True

  • plm (xr.DataArray, optional) – Precomputed plm values as a xr.DataArray variable. For example with the code : plm = xr.DataArray(compute_plm(lmax, sin_latitude), dims=[‘l’, ‘m’, ‘latitude’], coords={‘l’: data.l, ‘m’: data.m, ‘latitude’: latitude})

  • normalization_plm (str, optional) – If plm need to be computed, choice of the norm corresponding to the SH dataset. Either ‘4pi’, ‘ortho’, or ‘schmidt’ for 4pi normalized, orthonormalized, or Schmidt semi-normalized SH functions, respectively. Default is ‘4pi’.

  • **kwargs – Supplementary parameters used by the function l_factor_conv to modify defaults constants used in the computation for the unit conversion. These parameters include (see l_factor_conv() documentation for more details) : a_earth, gm_earth, f_earth, omega_earth, rho_earth, ds_love

Returns:

xgrid – Spatial representation of the SH Dataset in the chosen unit.

Return type:

xr.DataArray

lenapy.utils.harmo.grid_to_sh(grid: DataArray, lmax: int, unit: Literal['mewh', 'mmgeoid', 'microGal', 'bar', 'mvcu', 'norm'] = 'mewh', mmax: int | None = None, lmin: int = 0, mmin: int = 0, used_l: ndarray | None = None, used_m: ndarray | None = None, ellipsoidal_earth: bool = False, include_elastic: bool = True, plm: DataArray | None = None, normalization_plm: Literal['4pi', 'ortho', 'schmidt'] = '4pi', **kwargs) → Dataset[source]

Transform gravity field spatial representation DataArray into Spherical Harmonics (SH) dataset. With choice for constants, unit of the spatial DataArray, love_numbers, degree/order, Earth hypothesis.

For details on unit transformations, see l_factor_conv() and the associated PDF document.

Parameters:

  • grid (xr.DataArray) – xr.Dataset that corresponds a gravity field spatial representation in a unit to convert into SH.

  • lmax (int) – Maximal degree of the SH coefficients to be computed.

  • unit (str, optional) – ‘mewh’, ‘mmgeoid’, ‘microGal’, ‘bar’, ‘mvcu’, or ‘norm’ Unit of the spatial data used in the transformation. Default is ‘mewh’ for meters of Equivalent Water Height. See constants.l_factor_conv() doc for details on the units.

  • mmax (int, optional) – Maximal order of the SH coefficients to be computed.

  • lmin (int, optional) – Minimal degree of the SH coefficients to be computed. Default is 0.

  • mmin (int, optional) – Minimal order of the SH coefficients to be computed. Default is 0.

  • used_l (np.ndarray, optional) – List of degree to compute for the SH Dataset (if given, lmax and lmin are not considered).

  • used_m (np.ndarray, optional) – List of order to compute for the SH Dataset (if given, mmax and mmin are not considered).

  • ellipsoidal_earth (bool, optional) – If True, consider the Earth as an ellipsoid following [Ditmar2018]. Default is False for a spherical Earth.

  • include_elastic (bool, optional) – If True, the Earth behavior is elastic. Default is True.

  • plm (xr.DataArray, optional) – Precomputed plm values as a xr.DataArray variable. For example with the code : plm = xr.DataArray(compute_plm(lmax, np.sin(np.deg2rad(latitude))), dims=[‘l’, ‘m’, ‘latitude’], coords={‘l’: np.arange(lmax+1), ‘m’: np.arange(lmax+1), ‘latitude’: latitude})

  • normalization_plm (str, optional) – If plm need to be computed, choice of the norm. Either ‘4pi’, ‘ortho’, or ‘schmidt’ for 4pi normalized, orthonormalized, or Schmidt semi-normalized SH functions, respectively. Default is ‘4pi’. Output SH coefficient will be normalized according to this parameter.

  • **kwargs – Supplementary parameters used by the function l_factor_conv to modify defaults constants used in the computation for the unit conversion. These parameters include (see l_factor_conv() documentation for more details) : a_earth, gm_earth, f_earth, omega_earth, rho_earth, ds_love

Returns:

ds_out – SH Dataset computed from the grid with the chosen unit.

Return type:

xr.Dataset

lenapy.utils.harmo.compute_plm(lmax: int, z: ~numpy.ndarray, mmax: int = None, normalization: ~typing.Literal['4pi', 'ortho', 'schmidt'] = '4pi', dtype: complex | float | type[complex] | type[float] = <class 'numpy.longdouble'>) → ndarray[source]

Compute all the associated Legendre functions up to a maximum degree and order using the recursion relation from [Holmes2002] (Adapted from SHTOOLS/pyshtools tools [Wieczorek2018])

Parameters:

  • lmax (int) – Maximum degree of legrendre functions.

  • z (np.ndarray) – Argument of the associated Legendre functions.

  • mmax (int or NoneType, optional) – Maximum order of associated legrendre functions.

  • normalization ({'4pi', 'ortho', 'schmidt'}, optional) – ‘4pi’, ‘ortho’, or ‘schmidt’ for use with geodesy 4pi normalized, orthonormalized, or Schmidt semi-normalized spherical harmonic functions, respectively. Default is ‘4pi’.

  • dtype (dtype, optional) – Data type of the output array. Default is np.float128.

Returns:

plm – Fully-normalized Legendre functions as a 3D array with “l”, “m” and z dimensions.

Return type:

np.ndarray

References

[Holmes2002]

S. A. Holmes and W. E. Featherstone, “A unified approach to the Clenshaw summation and the recursive computation of very high degree and order normalised associated Legendre functions”, Journal of Geodesy, 76, 279–299, (2002). doi: 10.1007/s00190-002-0216-2

[Wieczorek2018]

M. A. Wieczorek and M. Meschede, “SHTools: Tools for working with spherical harmonics”, Geochemistry, Geophysics, Geosystems, 19, 2574–2592, (2018). doi: 10.1029/2018GC007529

lenapy.utils.harmo.mid_month_grace_estimate(begin_time: datetime, end_time: datetime) → datetime[source]

Calculate middle of the month date based on begin_time and end_time for GRACE products. begin_time is rounded to equal to the first day of the month and end_time is rounded to equal to the first day of the month after.

Parameters:

  • begin_time (datetime.datetime) – Date of the beginning of the month.

  • end_time (datetime.datetime) – Date of the end of the month + 1 day.

Returns:

mid_month – Date of the middle of the month.

Return type:

datetime.datetime

lenapy.utils.harmo.change_normalization(ds: Dataset, new_normalization: Literal['4pi', 'ortho', 'schmidt'] = '4pi', old_normalization: Literal['4pi', 'ortho', 'schmidt'] | None = None, apply: bool = False) → Dataset[source]

Spherical Harmonics (SH) dataset are associated with a Legendre polynomial normalization. This function return the dataset updated with the new normalization asked in input (as a deep copy by default). The current normalization can be given in parameters or can be contained in ds.attrs[‘norm’].

Parameters:

  • ds (xr.Dataset) – xr.Dataset that corresponds to SH data associated with a current reference frame (constants whise) to update.

  • new_normalization (str) – New normalization for the SH dataset, either ‘4pi’, ‘ortho’, or ‘schmidt’.

  • old_normalization (str | None, optional) – Current normalization of the SH dataset. Either ‘4pi’, ‘ortho’, ‘schmidt’ or ‘unnorm’ for 4pi normalized, orthonormalized, Schmidt semi-normalized, or unnormalized SH functions, respectively. Default is ‘4pi’. If not provided, uses ds.attrs[‘norm’].

  • apply (bool, optional) – If True, apply the update to the current dataset without making a deep copy. Default is False.

Returns:

ds_out – Updated dataset with the new normalization.

Return type:

xr.Dataset

Examples

>>> ds_new_norm = change_normalization(ds, new_normalization='schmidt')
lenapy.utils.harmo.load_default_love_numbers() → Dataset[source]

Load default Love numbers dataset.

Returns:

ds_love – Love numbers dataset.

Return type:

xr.Dataset

lenapy.utils.harmo.l_factor_conv(l: ndarray, unit: Literal['mewh', 'mmgeoid', 'microGal', 'pascal', 'mvcu', 'norm'] = 'mewh', include_elastic: bool = True, ellipsoidal_earth: bool = False, geocentric_colat: DataArray | None = None, ds_love: Dataset | None = None, a_earth: float | None = None, gm_earth: float | None = None, f_earth: float = 0.003352810681182319, omega_earth: float = 7.292115e-05, rho_earth: float = 5513.407482020534, rho_water: float = 1000, attrs: dict | None = None) → tuple[DataArray, dict][source]

Compute a scale factor for a transformation between spherical harmonics and grid data. Spatial data over the grid are associated with a specific unit. The scale factor is degree-dependent and is computed for the given list of degree l. The scale factor can be estimated using elastic or non-elastic Earth as well as a spherical or ellipsoidal Earth.

Parameters:

  • l (np.ndarray) – Degree for which the scale factor is estimated.

  • unit (str, optional) – ‘mewh’, ‘mmgeoid’, ‘microGal’, ‘pascal’, ‘mvcu’, or ‘norm’ Unit of the spatial data used in the transformation. Default is ‘mewh’ for meters of Equivalent Water Height. ‘mmgeoid’ represents millimeters mmgeoid height, ‘microGal’ represents microGal gravity perturbations, ‘pascal’ represents equivalent surface pressure in pascal and ‘mvcu’ represents meters viscoelastic crustal uplift

  • include_elastic (bool, optional) – If True, the Earth behavior is elastic.

  • ellipsoidal_earth (bool, optional) – If True, consider the Earth as an ellipsoid following [Ditmar2018] and if False as a sphere.

  • geocentric_colat (xr.DataArray | None, optional) – Geocentric colatitude for ellipsoidal Earth radius computation in radians, the dimension is geographic latitude.

  • ds_love (xr.Dataset | None, optional) – Dataset with the l dimension corresponding to degree and with l (and possibly h and k) variables that are Love numbers. Default Love numbers used are from Gegout97.

  • a_earth (float, optional) – Earth semi-major axis [m]. If not provided, uses data.attrs[‘radius’] and if it does not exist, uses LNPY_A_EARTH_GRS80. if it does not exist, uses LNPY_A_EARTH_GRS80.

  • gm_earth (float, optional) – Standard gravitational parameter for Earth [m³.s⁻²]. If not provided, uses data.attrs[‘earth_gravity_constant’] and if it does not exist, uses LNPY_GM_EARTH.

  • f_earth (float, optional) – Earth flattening. Default is LNPY_F_EARTH_GRS80.

  • omega_earth (float, optional) – Earth’s rotation rate [rad.s⁻¹]. Default is LNPY_OMEGA_EARTH.

  • rho_earth (float, optional) – Earth density [kg.m⁻³] for Equivalent Water Height formula. Default is LNPY_RHO_EARTH.

  • rho_water (float, optional) – Water density [kg.m⁻³] for Equivalent Water Height formula. Default is 1000 Kg.m⁻³.

  • attrs (dict | None, optional) – ds.attrs information that might help to estimate l_factor if no parameters are given.

Returns:

  • l_factor (xr.DataArray) – Degree-dependent scale factor.

  • cst (dict) – Attributes for the final dataset

Return type:

tuple[DataArray, dict]

References

[Ditmar2018]

P. Ditmar, “Conversion of time-varying Stokes coefficients into mass anomalies at the Earth’s surface considering the Earth’s oblateness”, Journal of Geodesy, 92, 1401–1412, (2018). doi: 10.1007/s00190-018-1128-0

lenapy.utils.harmo.assert_sh(ds: Dataset) → bool[source]

Verify if the dataset ds has dimensions ‘l’ and ‘m’ as well as variables ‘clm’ and ‘slm’ Raise Assertion error if not.

Parameters:

ds (xr.Dataset) – Dataset to verify.

Returns:

True – Returns True if the dataset ds has dimensions ‘l’ and ‘m’ as well as variables ‘clm’ and ‘slm’.

Return type:

bool

Raises:

AssertionError – This function raise AssertionError is ds is not a xr.Dataset corresponding to spherical harmonics