Welcome to pygeoid’s documentation!

pygeoid

Local gravity field modelling with Python.

The main purpose of the package is to realise common algorithms and methods in physical geodesy for the local gravity field modelling.

The project is in the early planning and development stage.

Install

To install pygeoid simply use pip as usual:

pip install pygeoid

Alternatively, you can use pip to install the latest version of the source code from GitHub:

pip install --upgrade git+https://github.com/ioshchepkov/pygeoid

Documentation

Full documentation is available at https://pygeoid.readthedocs.io

Code of Conduct

Please note that PyGeoid is released with a Code of Conduct (see CODE_OF_CONDUCT.md). By participating in this project you agree that you will follow all its terms.

License

This project is licensed under the MIT License - see the LICENSE.txt file for details

Installing

Installing with pip

You can install pygeoid using the pip package manager:

pip install pygeoid

Installing the latest development version

You can use pip to install the latest version of the source code from GitHub:

python -m pip install --upgrade git+https://github.com/ioshchepkov/pygeoid

API Reference

Coordinate Systems and Transformations (pygeoid.coordinates)

Reference Ellipsoid (pygeoid.coordinates.ellipsoid.Ellipsoid)

Geometry of the reference ellipsoid.

class pygeoid.coordinates.ellipsoid.Ellipsoid(ellps: Optional[str] = None, **kwargs)

Class represents an ellipsoid of revolution and its geometry.

This class uses proj.Geod class from pyproj package, so any valid init string for Proj are accepted as arguments. See pyproj.Geod.__new__ documentation (https://pyproj4.github.io/pyproj/stable/api/geod.html) for more information.

Parameters:ellps (str, optional) – Ellipsoid name, most common ellipsoids are accepted. Default is ‘GRS80’.

equatorial_radius

Return semi-major or equatorial axis radius, in metres.

polar_radius

Return semi-minor or polar axis radius, in metres.

flattening

Return flattening of the ellipsoid.

Notes

The flattening of the ellipsoid \(f\) is

\[f = \frac{a - b}{a},\]

where \(a\) and \(b\) – equatorial and polar axis of the ellipsoid respectively.

reciprocal_flattening

Return reciprocal (inverse) flattening.

eccentricity

Return first eccentricity.

Notes

The first eccentricity of the ellipsoid \(e\) is

\[e = \sqrt{\frac{a^2 - b^2}{a^2}},\]

where \(a\) and \(b\) – equatorial and polar axis of the ellipsoid respectively.

eccentricity_squared

Return first eccentricity squared.

second_eccentricity

Return second eccentricity.

Notes

The second eccentricity of the ellipsoid \(e'\) is

\[e' = \sqrt{\frac{a^2 - b^2}{b^2}}\]

where \(a\) and \(b\) – equatorial and polar axis of the ellipsoid respectively.

second_eccentricity_squared

Return second eccentricity squared.

linear_eccentricity

Return linear eccentricity, in metres.

Notes

The linear eccentricity of the ellipsoid \(E\) is

\[E = ae,\]

where \(a\) – equatorial radius of the ellipsoid, \(e\) – (first) eccentricity.

polar_curvature_radius

Return polar radius of curvature, in metres.

Notes

The polar radius of curvature of the ellipsoid \(c\) is

\[c = \frac{a^2}{b},\]

where \(a\) and \(b\) – equatorial and polar axis of the ellipsoid respectively.

quadrant_distance

Return arc of meridian from equator to pole, in metres.

Notes

The arc length of meridian from equator to pole is

\[Q = c\frac{\pi}{2}\left( 1 - \frac{3}{4}e'^2 + \frac{45}{64}e'^4 + \frac{175}{256}e'^6 + \frac{11025}{16384}e'^8\right),\]

where \(c\) – polar radius of curvature, \(e'\) – second eccentricity.

surface_area

Return surface area of the ellipsoid, in squared metres.

Notes

The surface area of the ellipsoid is

\[A = 2\pi a^2 \left[1 + \frac{1 - e^2}{2e} \ln{\left( \frac{1 + e}{1 - e}\right)}\right],\]

where \(a\) – equatorial axis of the ellipsoid, \(e\) – (first) eccentricity.

volume

Return volume of the ellipsoid, in cubical metres.

Notes

The volume of the ellipsoid is

\[V = \frac{4}{3}\pi a^2 b,\]

where \(a\) and \(b\) – equatorial and polar axis of the ellipsoid respectively.

mean_radius(kind: str = 'arithmetic')

Return the radius of a sphere.

Parameters:kind ({'arithmetic', 'same_area', 'same_volume'}, optional) –

Controls what kind of radius is returned.

• ’arithmetic’ returns the arithmetic mean value

  \(R_m\) of the 3 semi-axis of the ellipsoid.

• ’same_area’ returns the authalic radius \(R_A\) of

  the sphere with the same surface area as the ellipsoid.

• ’same_volume’ returns the radius \(R_V\) of

  the sphere with the same volume as the ellipsoid.

Default is ‘arithmetic’.

Returns:Mean radius of the ellipsoid, in metres. Return type:float

Notes

The arithmetic mean radius of the ellipsoid is

\[R_m = \frac{2a + b}{2},\]

where \(a\) and \(b\) are equatorial and polar axis of the ellipsoid respectively.

A sphere with the same surface area as the elliposid has the radius

\[R_A = \sqrt{\frac{A}{4\pi}},\]

where \(A\) is the surface area of the ellipsoid.

A sphere with the same volume as the ellipsoid has the radius

\[R_V = a^2 b.\]

meridian_curvature_radius(lat: Unit("deg")) -> Unit("m")

Return radius of curvature of meridian normal section.

Parameters:lat (Quantity) – Geodetic latitude. Returns:Value of the radius of curvature of meridian normal section. Return type:Quantity

Notes

The radius of curvature of meridian normal section \(M\) is

\[M = \frac{c}{V^3},\]

where \(c\) – polar radius of curvature, \(V\) – auxiliary function which depends on geodetic latitude.

prime_vertical_curvature_radius(lat: Unit("deg")) -> Unit("m")

Return radius of curvature of prime vertical normal section.

Parameters:lat (Quantity) – Geodetic latitude. Returns:Value of the radius of curvature of prime vertical normal section. Return type:Quantity

Notes

The radius of curvature of prime vertical \(N\) is

\[N = \frac{c}{V},\]

where \(c\) – polar radius of curvature, \(V\) – auxiliary function which depends on geodetic latitude.

mean_curvature(lat: Unit("deg")) → <Quantity 1. 1 / m>

Return mean curvature, in inverse metres.

Parameters:lat (Quantity) – Geodetic latitude. Returns:Value of the mean curvature. Return type:Quantity

Notes

The mean curvature is \(1/\sqrt{MN}\), where \(M\) – radius of curvature of meridian normal section, \(N\) – radius of curvature of prime vertical.

gaussian_curvature(lat: Unit("deg")) → <Quantity 1. 1 / m>

Return Gaussian curvature, in inverse metres.

Parameters:lat (Quantity) – Geodetic latitude. Returns:Value of the Gaussian radius of curvature. Return type:Quantity

Notes

The Gaussian curvature is \(1/MN\), where \(M\) – radius of curvature of meridian normal section, \(N\) – radius of curvature of prime vertical.

average_curvature(lat: Unit("deg")) → <Quantity 1. 1 / m>

Return average curvature, in inverse metres.

Parameters:lat (Quantity) – Geodetic latitude. Returns:Value of the average curvature. Return type:Quantity

Notes

The average curvature is

\[\frac{1}{2} \left( \frac{1}{M} + \frac{1}{N} \right),\]

where \(M\) – radius of curvature of meridian normal section, \(N\) – radius of curvature of prime vertical.

meridian_arc_distance(lat1: Unit("deg"), lat2: Unit("deg")) -> Unit("m")

Return the distance between two parallels lat1 and lat2.

Parameters:

• lat1 (Quantity) – Geodetic latitude of the first point.
• lat2 (Quantity) – Geodetic latitude of the second point.

Returns:

The distance between two parallels.

Return type:

Quantity

parallel_arc_distance(lat: Unit("deg"), lon1: Unit("deg"), lon2: Unit("deg"))

Return the distance between two points on a parallel.

Parameters:

• lat (Quantity) – Geodetic latitude of the parallel.
• lon1 (Quantity) – Geodetic longitude of the first point.
• lon2 (Quantity) – Geodetic longitude of the second point.

Returns:

The distance between two meridians along the parallel.

Return type:

Quantity

fwd(lat: Unit("deg"), lon: Unit("deg"), azimuth: Unit("deg"), distance: Unit("m"))

Solve forward geodetic problem.

Returns latitudes, longitudes and back azimuths of terminus points given latitudes lat and longitudes lon of initial points, plus forward `azimuth`s and `distance`s.

This method use pyproj.Geod.fwd as a backend.

Parameters:

• lat (Quantity) – Geodetic latitude of the initial point.
• lon (Quantity) – Longitude of the initial point.
• azimuth (Quantity) – Geodetic azimuth.
• distance (Quantity) – Distance.

Returns:

• lat (~astropy.units.Quantity) – Geodetic latitude of the terminus point.
• lon (~astropy.units.Quantity) – Longitude of the terminus point.
• back_azimuth (~astropy.units.Quantity) – Back geodetic azimuth.

inv(lat1: Unit("deg"), lon1: Unit("deg"), lat2: Unit("deg"), lon2: Unit("deg"))

Solve inverse geodetic problem.

Returns forward and back azimuths, plus distances between initial points (specified by lat1, lon1) and terminus points (specified by lat1, lon2).

This method use pyproj.Geod.inv as a backend.

Parameters:

• lat1 (Quantity) – Geodetic latitude of the initial point.
• lon1 (Quantity) – Longitude of the initial point.
• lat2 (Quantity) – Geodetic latitude of the terminus point.
• lon2 (Quantity) – Longitude of the terminus point.

Returns:

• azimuth (~astropy.units.Quantity) – Geodetic azimuth.
• back_azimuth (~astropy.units.Quantity) – Back geodetic azimuth.
• distance (~astropy.units.Quantity) – Distance, in metres.

npts(lat1: Unit("deg"), lon1: Unit("deg"), lat2: Unit("deg"), lon2: Unit("deg"), npts: int) -> Unit("deg")

Return equaly spaced points along geodesic line.

Given a single initial point and terminus point (specified by lat1, lon1 and lat2, lon2), returns a list of longitude/latitude pairs describing npts equally spaced intermediate points along the geodesic between the initial and terminus points.

This method use pyproj.Geod.npts as a backend.

Parameters:

• lat1 (Quantity) – Geodetic latitude of the initial point.
• lon1 (Quantity) – Longitude of the initial point.
• lat2 (Quantity) – Geodetic latitude of the terminus point.
• lon2 (Quantity) – Longitude of the terminus point.
• npts (int) – Number of intermediate points.

Returns:

points – List of latitudes and longitudes of the intermediate points.

Return type:

Quantity list of tuples

circle_radius(lat: Unit("deg")) -> Unit("m")

Return the radius of the parallel, in metres.

Parameters:lat (Quantity) – Geodetic latitude.

Notes

The radius of the parallel \(\phi\) is

\[r_\phi = N \cos{\phi},\]

where \(N\) – radius of curvature of prime vertical, \(\phi\) – geodetic latitude.

polar_equation(lat: Unit("deg")) -> Unit("m")

Return radius of the ellipsoid with respect to the origin.

Parameters:lat (Quantity) – Geocentric latitude. Returns:Geocentric radius of the parallel. Return type:Quantity

Notes

The polar equation of the ellipsoid is

\[r = \frac{ab}{\sqrt{a^2\sin^2{\vartheta} + b^2\cos^2{\vartheta}}},\]

where \(a\) and \(b\) – equatorial and polar axis of the ellipsoid respectively, \(\vartheta\) – geocentric latitude.

geocentric_latitude(lat: Unit("deg")) -> Unit("deg")

Convert geodetic latitude to geocentric latitude.

Parameters:lat (Quantity) – Geodetic latitude. Returns:Geocentric (spherical) latitude. Return type:Quantity

Notes

The relationship between geodetic \(\phi\) and geocentric \(\vartheta\) latitudes is

\[\vartheta = \tan^{-1}{\left(\left(1 - f\right)^2\tan\phi\right)},\]

where \(f\) – flattening of the ellipsoid.

reduced_latitude(lat: Unit("deg")) -> Unit("deg")

Convert geodetic latitude to reduced (parametric) latitude.

Parameters:lat (Quantity) – Geodetic latitude. Returns:Reduced latitude. Return type:Quantity

Notes

The relationship between geodetic \(\phi\) and reduced \(\beta\) latitudes is

\[\beta = \tan^{-1}{\left(\left(1 - f\right)\tan\phi\right)},\]

where \(f\) – flattening of the ellipsoid.

authalic_latitude(lat: Unit("deg")) -> Unit("deg")

Convert geodetic latitude to authalic latitude.

Authalic latitude will return a geocentric latitude on a sphere having the same surface area as the ellipsoid. It will preserve areas with relative to the ellipsoid. The authalic radius can be calculated from mean_radius(kind=’same_area’) method.

Parameters:lat (Quantity) – Geodetic latitude. Returns:Authalic latitude. Return type:Quantity

Position 3D (pygeoid.coordinates.position)

Coordinate Transformations (pygeoid.coordinates.transform)

Geodetic Integrals (pygeoid.integrals)

Stokes’s Kernel and Integral (pygeoid.integrals.stokes)

Stokes integral and kernel.

class pygeoid.integrals.stokes.StokesKernel

Stokes kernel class.

kernel(spherical_distance: Unit("deg"))

Evaluate Stokes spherical kernel.

This method will calculate the original Stokes’s function.

Parameters:spherical_distance (Quantity) – Spherical distance, in radians.

Notes

In closed form, Stokes’s kernel depends on the spherical distance \(\psi\) by [1]_:

\[S\left(\psi\right) = \dfrac{1}{\sin{(\psi / 2)}} - 6\sin{(\psi/2)} + 1 - 5\cos{\psi} - 3\cos{\psi} \ln{\left[\sin{(\psi/2)} + \sin^2{(\psi/2)}\right]}.\]

References

[1]	Heiskanen WA, Moritz H (1967) Physical geodesy. Freeman, San Francisco

derivative_spherical_distance(spherical_distance)

Evaluate Stokes’s spherical kernel derivative.

The derivative of the Stokes function is the Vening-Meinesz function.

Parameters:spherical_distance (Quantity) – Spherical distance.

Notes

The derivative of Stokes’s kernel is the Vening-Meinesz and it depends on the spherical distance \(\psi\) by [1]_:

\[\dfrac{d S\left(\psi\right)}{d\psi} = - \dfrac{\cos{(\psi / 2)}}{2\sin^2{(\psi / 2)}} + 8\sin{\psi} - 6\cos{(\psi / 2)} - 3\dfrac{1 - \sin{(\psi / 2)}}{\sin{\psi}} + 3\sin{\psi}\ln{\left[\sin{(\psi/2)} + \sin^2{(\psi/2)}\right]}.\]

References

[1]	Heiskanen WA, Moritz H (1967) Physical geodesy. Freeman, San Francisco

name

Return kernel name.

class pygeoid.integrals.stokes.StokesExtendedKernel

name

Return kernel name.

Vening Meinesz Kernel and Integral (pygeoid.integrals.veningmeinesz)

Vening Meinesz kernel.

class pygeoid.integrals.veningmeinesz.VeningMeineszKernel

Vening Meinesz kernel class.

kernel(spherical_distance: Unit("deg"))

Evaluate Vening Meinesz kernel.

Parameters:

• spherical_distance (float or array_like of floats) – Spherical distance, in radians.
• degrees (bool, optional) – If True, the spherical distance is given in degrees, otherwise radians.

Notes

The derivative is the Vening-Meinesz and it depends on the spherical distance \(\psi\) by [1]:

\[V\left(\psi\right) = \dfrac{d S\left(\psi\right)}{d\psi} = - \dfrac{\cos{(\psi / 2)}}{2\sin^2{(\psi / 2)}} + 8\sin{\psi} - 6\cos{(\psi / 2)} - 3\dfrac{1 - \sin{(\psi / 2)}}{\sin{\psi}} + 3\sin{\psi}\ln{\left[\sin{(\psi/2)} + \sin^2{(\psi/2)}\right]},\]

which is the derivative of Stokes’s kernel with respect to \(\psi\).

References

[1]	Heiskanen WA, Moritz H (1967) Physical geodesy. Freeman, San Francisco

name

Return kernel name.

Truncation coefficients (pygeoid.integrals.truncation)

Truncation coefficients for Stokes’s integral.

pygeoid.integrals.truncation.molodensky_truncation_coefficients(spherical_distance: float, degree_n: int, method: str = 'hagiwara', **kwargs) → numpy.ndarray

Evaluate Molodensky’s truncation coefficients Qn.

Compute sequence of Molodensky’s truncation coefficients for all degrees from 0 to degree_n (inclusive).

Parameters:

• spherical_distance (float) – Spherical distance, in degrees.
• degree_n (int) – Maximum degree of the coefficients.
• method ({'hagiwara', 'numerical'}, optional) –

  Controls how coefficients are calculated.

  • ’hagiwara’ calculates coefficients by Hagiwara (1976) recurrence relations.
  • ’numerical’ calculates coefficients by numerical integration.

  Default is ‘hagiwara’.

• **kwargs – Keyword arguments for scipy.intagrate.quad if method is ‘numerical’.

Returns:

Molodensky’s truncation coefficient for all degrees from 0 to degree_n (inclusive).

Return type:

array_like of floats

Notes

The Molodensky’s truncation coefficients \(Q_n\) of degree \(n\) are defined as [1]_:

\[Q_n \left(\psi_0\right) = \int\limits_{\psi_0}^{\pi} S\left(\psi\right) P_n \left(\cos{\psi} \right) \sin{\psi} d\psi,\]

where \(S\left(\psi\right)\) – Stokes function, \(P_n\) – Legendre polynomial, \(\psi\) – spherical distance.

The function calculates this integral by Hagiwara’s [2] method or by the numerical integration with scipy.intagrate.quad.

References

[1]	Molodensky MS, Yeremeyev VF, Yurkina MI (1962) Methods for study of the external gravitational field and figure of the Earth. Translated from Russian, Isreali Programme for Scientific Translations, Jerusalem

[2]	Hagiwara Y (1976) A new formula for evaluating the truncation error coefficient. Bulletin Géodésique 50:131–135.

pygeoid.integrals.truncation.paul_coefficients(spherical_distance: float, n: int, k: int = None, method: str = 'paul', **kwargs) → numpy.ndarray

Return Paul’s coefficients.

In the original article (1973) the Paul’s coefficients are denoted as Rnk, but in many later articles they are denoted as enk.

Parameters:

• spherical_distance (float or array_like of floats) – Spherical distance.
• n,k (int) – Degrees of the coefficients. k by default is None, i.e. it is equal to n.
• method ({'paul', 'numerical'}, optional) –

  Controls how coefficients are calculated.

  • ’paul’ calculate coefficients by Paul (1973) recurrence relations.
  • ’numerical’ calculate coefficients by numerical integration.

  Default is ‘paul’.

• **kwargs – Keyword arguments for scipy.intagrate.quad if method is ‘numerical’.

Notes

The Pauls’s coefficients \(e_{nk}\) of degrees \(n\) and \(k\) are defined as [1]_:

\[e_{nk} \left(\psi_0\right) = \int\limits_{\psi_0}^{\pi} P_n \left(\cos{\psi}\right) P_k \left(\cos{\psi}\right) \sin{\psi} d\psi,\]

where \(P_n\) and \(P_k\) are Legendre polynomial of degrees \(n\) and \(k\) respectively, \(\psi\) is the spherical distance.

Note that in the original article [1]_ the Paul’s coefficients are denoted as \(R_{n,k}\).

References

[1]	Paul MK (1973) A method of evaluating the truncation error coefficients for geoidal height. Bull Géodésique 110:413–425

Gravity Reduction (pygeoid.reduction)

Normal Gravity Field (pygeoid.reduction.normal)

Gravity field and geometry of the level ellipsoid.

class pygeoid.reduction.normal.Centrifugal(omega=<Quantity 7.292115e-05 rad / s>)

Centrifugal potential and its derivatives.

Parameters:omega (float) – Angular rotation rate of the body, in rad/s. Default value is the angular speed of Earth’s rotation 7292115e-11 rad/s

potential(lat: Unit("deg"), radius: Unit("m")) -> Unit("m2 / s2")

Return centrifugal potential.

Parameters:

• lat (Quantity) – Spherical (geocentric) latitude.
• radius (Quantity) – Radius.

r_derivative(lat: Unit("deg"), radius: Unit("m"))

Return radial derivative.

Parameters:

• lat (Quantity) – Spherical (geocentric) latitude.
• radius (Quantity) – Radius.

lat_derivative(lat: Unit("deg"), radius: Unit("m"))

Return latitude derivative.

Parameters:

• lat (Quantity) – Spherical (geocentric) latitude.
• radius (Quantity) – Radius.

gradient(lat: Unit("deg"), radius: Unit("deg")) -> Unit("m / s2")

Return centrifugal force.

Parameters:

• lat (Quantity) – Spherical (geocentric) latitude.
• radius (Quantity) – Radius.

class pygeoid.reduction.normal.LevelEllipsoid(ellps=None, **kwargs)

Class represents the gravity field of the level ellipsoid.

This class intialize Ellipsoid class from pygeoid.coordinates.ellipsoid, so all geometrical methods and parameters are available too.

Parameters:ellps ({'GRS80', 'WGS84', 'PZ90', 'GSK2011'}, optional) – Ellipsoid name. Default is ‘GRS80’.

j2

Return dynamic form factor J2.

gm

Return geocentric gravitational constant.

omega

Return angular velocity, in radians.

m

Auxiliary constant.

Notes

\[m = \frac{{\omega}^2 a^2 b}{GM}.\]

surface_potential

Return normal gravity potential on the ellipsoid.

Value of the normal gravity potential on the ellipsoid, or on the equipotential surface U(x, y, z) = U_0.

gravitational_potential(rlat: Unit("deg"), u_ax: Unit("m")) -> Unit("m2 / s2")

Return normal gravitational potential V.

Calculate normal gravitational potential from the rigorous formula.

Parameters:

• rlat (Quantity) – Reduced latitude.
• u_ax (Quantity) – Polar axis of the ellipsoid passing through the point.

Returns:

Normal gravitational potential.

Return type:

Quantity

gravity_potential(rlat: Unit("deg"), u_ax: Unit("m")) -> Unit("m2 / s2")

Return normal gravity potential U.

Calculate normal gravity potential from the rigorous formula.

Parameters:

• rlat (Quantity) – Reduced latitude.
• u_ax (Quantity) – Polar axis of the ellipsoid passing through the point.

Returns:

Normal gravity potential.

Return type:

Quantity

equatorial_normal_gravity

Return normal gravity at the equator.

polar_normal_gravity

Return normal gravity at the poles.

mean_normal_gravity

Return mean normal gravity over ellipsoid.

gravity_flattening

Return gravity flattening.

f* = (gamma_p - gamma_e) / gamma_e

conventional_gravity_coeffs()

Return coefficients for the conventional gravity formula.

gamma_0 = gamma_e*(1 + beta*sin(lat)**2 - beta1*sin(2*lat)**2)

surface_normal_gravity(lat: Unit("deg")) -> Unit("m / s2")

Return normal gravity on the ellipsoid.

Calculate normal gravity value on the level ellipsoid by the rigorous formula of Somigliana.

Parameters:lat (Quantity) – Geodetic latitude. Returns:Normal gravity on the ellipsoid. Return type:Quantity

surface_vertical_normal_gravity_gradient(lat: Unit("deg"))

Return the vertical gravity gradient on the ellipsoid.

Vertical gradient of the normal gravity at the reference ellipsoid.

Parameters:lat (Quantity) – Geodetic latitude.

height_correction(lat: Unit("deg"), height: Unit("m")) -> Unit("m / s2")

Return height correction.

Second-order approximation formula is used instead of -0.3086*height.

Parameters:

• lat (Quantity) – Geodetic latitude.
• height (Quantity) – Geodetic height.

normal_gravity(rlat: Unit("deg"), u_ax: Unit("m")) -> Unit("m / s2")

Return normal gravity.

Calculate normal gravity value at any arbitrary point by the rigorous closed formula.

Parameters:

• rlat (Quantity) – Reduced latitude.
• u_ax (Quantity) – Polar axis of the ellipsoid passing through the point.

Returns:

Normal gravity.

Return type:

Quantity

j2n(n: int) -> Unit(dimensionless)

Return even zonal coefficients J with a degree of 2*n.

If n = 0, the function will return -1. If n = 1, the function will return J2 (dynamic form factor). If n = 2, the function will return J4. If n = 3, the function will return J6. If n = 4, the function will return J8.

Parameters:n (int) – Degree of the J coefficient.

gravitational_potential_sph(lat: Unit("deg"), radius: Unit("m"), n_max: int = 4) -> Unit("m2 / s2")

Return normal gravitational potential V.

Calculate normal gravitational potential from spherical approximation.

Parameters:

• lat (Quantity) – Spherical (geocentric) latitude.
• radius (Quantity) – Radius, in metres.
• n_max (int) – Maximum degree.

gravity_potential_sph(lat: Unit("deg"), radius: Unit("m"), n_max: int = 4) -> Unit("m2 / s2")

Return normal gravitational potential V.

Calculate normal gravitational potential from spherical approximation.

Parameters:

• lat (Quantity) – Spherical (geocentric) latitude.
• radius (Quantity) – Radius, in metres.
• n_max (int) – Maximum degree.

pygeoid.reduction.normal.surface_normal_gravity_clairaut(lat: Unit("deg"), model: str = None) -> Unit("m / s2")

Return normal gravity from the first Clairaut formula.

Parameters:

• lat (Quantity) – Geodetic latitude.
• model ({'helmert', 'helmert_14mGal', '1930', '1930_14mGal', '1967', '1980'}) – Which gravity formula will be used.

Topography effects (pygeoid.reduction.topography)

Topographic reduction in gravity and geoid modelling

pygeoid.reduction.topography.bouguer_plate(height: Unit("m"), density: Unit("kg / m3") = <Quantity 2670. kg / m3>) -> Unit("mGal")

Return an attraction of an infinite Bouguer plate.

Parameters:

• height (Quantity) – Height above sea level.
• density (Quantity) – Density of the prism. Default is 2670 kg/m**3.

Notes

\[F_B = 2\pi G\delta H,\]

where \(G\) – gravitational constant, \(\delta\) – density, \(H\) – height above sea level.

pygeoid.reduction.topography.spherical_bouguer_cap(height: Unit("m"), density: Unit("kg / m3") = <Quantity 2670. kg / m3>) -> Unit("mGal")

Return spherical Bouguer correction.

Parameters:

• height (Quantity) – Height above sea level.
• density (Quantity) – Density of the prism. Default is 2670 kg/m**3.

Notes

The corected (spherical) Bouguer attraction accounts the curvature of the Earth. It is calclated by the closed-form formula for a spherical cap of radius 166.7 km [1]:

\[F_{SB} = 2\pi G ((1 + \mu) H - \lambda R),\]

where \(G\) – gravitational constant, \(\delta\) – density, \(H\) – height above sea level, \(\lambda\) and \(\mu\) – dimensionless coefficients, \(R = R_e + H\) – sum of the mean radius of the Earth and the height.

References

[1]	LaFehr, T.R., 1991. An exact solution for the gravity

curvature (Bullard B) correction. Geophysics, 56(8), pp.1179-1184.

Atmosphere correction (pygeoid.reduction.atmosphere)

Calculate atmospheric correction for the gravity anomalies.

pygeoid.reduction.atmosphere.ussa76_density(alt_arr: Unit("km") = <Quantity 0. km>) -> Unit("kg / m3")

Return atmospheric density from USSA76 model.

Refer to the following document (2 doc codes) for details of this model:

NOAA-S/T 76-1562 NASA-TM-X-74335

All assumptions are the same as those in the source documents.

Derived from: https://github.com/mattljc/atmosphere-py

Parameters:alt_arr (Quantity) – Altitude above sea level.

pygeoid.reduction.atmosphere.iag_atm_corr_sph(density_function: Callable[[astropy.units.quantity.Quantity], astropy.units.quantity.Quantity], height: Unit("m"), height_max: Unit("m"), samples=10000.0) -> Unit("mGal")

Return atmospheric correction to the gravity anomalies by IAG approach.

This function numerically integrates samples from density function by trapezoidal rule. The spherical layering of the atmosphere is considered.

IAG approach:

g_atm = G*M(r) / r**2

inf /

M(r) = 4*pi*| rho(r) * r**2 dr

/ r

Parameters:

• density_function (callable) – The density_funtion is called for all height samples to calculate density of the atmosphere.
• height (Quantity) – Height above sea level.
• height_max (Quantity) – Maximum height of the atmosphere layer above sea level.
• samples (float) – Number of samples for integration. Default is 1e4.
• Quantity – Atmospheric correction.

pygeoid.reduction.atmosphere.grs80_atm_corr_interp(height: Unit("m"), kind: str = 'linear') -> Unit("mGal")

Return GRS 80 atmospheric correction, in mGal.

Interpolated from the table data [1]_.

Note: If height < 0 m or height > 40000 m, then correction is extrapolated

Parameters:

• height (Quantity) – Height above sea level.
• kind (str or int, optional) – Specifies the kind of interpolation as a string (‘linear’, ‘nearest’, ‘zero’, ‘slinear’, ‘quadratic’, ‘cubic’ where ‘zero’, ‘slinear’, ‘quadratic’ and ‘cubic’ refer to a spline interpolation of zeroth, first, second or third order) or as an integer specifying the order of the spline interpolator to use. Default is ‘linear’.

Returns:

Atmospheric correction.

Return type:

Quantity

References

[1]	Moritz, H. (1980). Geodetic reference system 1980.

Bulletin Géodésique, 54(3), 395-405

pygeoid.reduction.atmosphere.wenzel_atm_corr(height: Unit("m")) -> Unit("mGal")

Return atmospheric correction by Wenzel, in mGal.

Parameters:height (Quantity) – Height above sea level. Returns:Atmospheric correction. Return type:Quantity

References

[1]	Wenzel, H., 1985, Hochauflosende Kugelfunktionsmodelle fur des

Gravitationspotential der Erde [1]: Wissenschaftliche arbeiten der Fachrichtung Vermessungswesen der Universitat Hannover, 137

pygeoid.reduction.atmosphere.pz90_atm_corr(height: Unit("m")) -> Unit("mGal")

Return PZ-90 atmospheric correction, in mGal.

Parameters:height (Quantity) – Height above sea level. Returns:Atmospheric correction. Return type:Quantity

Spherical Harmonics (pygeoid.sharm)

Global Gravity Field Model (pygeoid.sharm.ggm.GlobalGravityFieldModel)

Legendre functions (pygeoid.sharm.legendre)

Legender polynomials and functions

pygeoid.sharm.legendre.lplm_d1(lmax, x)

Associated Legendre functions of the first kind and their derivatives.

Parameters:

• lmax (int) – Maximum degree and order.
• x (float) – Input value.

Returns:

• Plm(x) ((lmax + 1, lmax + 1) array) – Values of Plm(x) for all orders and degrees.
• Plm_d(x) ({(lmax + 1, lmax + 1), None} array) – Values of Plm_d(x) for all orders and degrees. Will return None if x is very close or equal to -1 or 1.

pygeoid.sharm.legendre.lplm(lmax, x)

Associated Legendre functions of the first kind.

Parameters:

• lmax (int) – Maximum degree and order of the Plm(x).
• x (float) – Input value.

Returns:

Plm(x) – Values of Plm(x) for all orders and degrees.

Return type:

(lmax + 1, lmax + 1) array

Simpe Bodies (pygeoid.simple)

Prism (pygeoid.simple.prism)

class pygeoid.simple.prism.Prism(bounds, density=1.0)

External gravitational field of the right rectangular prism.

Parameters:

• bounds ((x1, x2, y1, y2, z1, z2), floats) – West, east, south, north, top and bottom of the prism, in metres.
• density (float) – Density of the prism, in kg/m**3.

Notes

The formulas from Nagy et al.[1]_ are used in this class.

References

[1]	Nagy, D., Papp, G. and Benedek, J., 2000. The gravitational potential

and its derivatives for the prism. Journal of Geodesy, 74(7-8), pp.552-560.

density

Return density value.

gravitation(x, y, z)

Return gravitation value.

invariants(x, y, z)

Return invariants of the gradient tensor.

tensor(x, y, z)

Return gradient tensor.

Tides (pygeoid.tides)

Laplace’s tidal representation (pygeoid.tides.laplace)

class pygeoid.tides.laplace.LaplaceTidalEquation(bodies: list = ['moon', 'sun'], parts: list = ['zonal', 'tesseral', 'sectorial'], love: dict = {'h': <<class 'astropy.constants.constant.Constant'> name='Nominal degree 2 Love number h2' value=0.6078 uncertainty=0 unit='' reference='IERS Conventions(2010), IERS Technical Note 36, Verlagdes Bundesamts für Kartographie und Geodäsie, Frankfurt am Main, Germany.'>, 'k': <<class 'astropy.constants.constant.Constant'> name='Nominal degree 2 Love number k2' value=0.29525 uncertainty=0 unit='' reference='IERS Conventions(2010), IERS Technical Note 36, Verlagdes Bundesamts für Kartographie und Geodäsie, Frankfurt am Main, Germany.'>, 'l': <<class 'astropy.constants.constant.Constant'> name='Nominal degree 2 Shida number l2' value=0.0847 uncertainty=0 unit='' reference='IERS Conventions(2010), IERS Technical Note 36, Verlagdes Bundesamts für Kartographie und Geodäsie, Frankfurt am Main, Germany.'>}, gravimetric_factor: float = 1.1563, diminishing_factor: float = 0.6947, ephemeris: str = None)

Class for Laplace’s tidal representation.

This class is a realision of the Laplace’s tidal representaion for the second-degree tidal potential:

W_2 = D (Z + T + S),

where D is called Doodson’s constant and Z, T and S are zonal, tesseral and sectorial parts of the tidal potential respectively. This representation is often applied for practical calculations of the total tidal potential and its derivatives, instead of more advanced and complicated harmonic decomposition.

If some tidal effect calculated from this class with all parts included (by default) and fot the elastic Earth is applied as a correction (effect taken with the negative sign), then the corrected value will be in the conventional tide free system, because a permanent part of the tide will also be removed. If mean or zero tide systems are desired, then permanent part of the tide should be restored in some way.

bodies

List of solar system bodies, the total tide from which will be taken into account. By default, only the tides from the Moon and the Sun is accounted since they are the largest.

Type:list, optional

parts

Which parts of Laplace’s tidal eqution are involved. There are ‘zonal’, ‘tesseral’ and ‘sectorial’ parts. By default, all parts are included.

Type:list, optional

love

Approximate Love numbers for elastic deformation calculations. This dictionary should contain ‘k’, ‘l’ and ‘h’ numbers of the second degree. Default values are taken from IERS Conventions 2010.

Type:dict, optional

gravimetric_factor

Gravimetric factor (delta). If None it will be calculated from Love numbers as delta = 1 + h - 3/2*k. Default value is approximate gravimetric factor from PREM model.

Type:float, optional

diminishing_factor

Diminishing factor (gamma). If None then it will be calculated from Love numbers as gamma = 1 + k - h. Default value is approximate diminishing factor from PREM model.

Type:float, optional

ephemeris

Ephemeris to use. If not given, use the one set with astropy.coordinates.solar_system_ephemeris.set (which is set to ‘builtin’ by default).

Type:str, optional

Notes

One can check which bodies are covered by a given ephemeris using:

>>> from astropy.coordinates import solar_system_ephemeris
>>> solar_system_ephemeris.bodies

References

[1]	Melchior, P. (1983), The Tides of the Planet Earth. 2nd edition, Pergamon Press.

potential(time: astropy.time.core.Time, lat: Unit("deg"), lon: Unit("deg"), r: Unit("m")) -> Unit("m2 / s2")

Return tidal potential.

Parameters:

• time (Time) – Time of observation.
• lat (Quantity) – Geocentric (spherical) latitude.
• lon (Quantity) – Geocentric (spherical) longitude.
• r (Quantity) – Geocentric radius (radial distance).

Returns:

potential – Tidal potential.

Return type:

Quantity

potential_lat_derivative(time: astropy.time.core.Time, lat: Unit("deg"), lon: Unit("deg"), r: Unit("m")) -> Unit("m2 / s2")

Return tidal potential latitude derivative.

Parameters:

• time (Time) – Time of observation.
• lat (Quantity) – Geocentric (spherical) latitude.
• lon (Quantity) – Geocentric (spherical) longitude.
• r (Quantity) – Geocentric radius (radial distance).

Returns:

potential_dlat – Tidal potential latitude derivative.

Return type:

Quantity

potential_lon_derivative(time: astropy.time.core.Time, lat: Unit("deg"), lon: Unit("deg"), r: Unit("m")) -> Unit("m2 / s2")

Return tidal potential longitude derivative.

Parameters:

• time (Time) – Time of observation.
• lat (Quantity) – Geocentric (spherical) latitude.
• lon (Quantity) – Geocentric (spherical) longitude.
• r (Quantity) – Geocentric radius (radial distance).

Returns:

potential_dlon – Tidal potential longitude derivative.

Return type:

Quantity

potential_r_derivative(time: astropy.time.core.Time, lat: Unit("deg"), lon: Unit("deg"), r: Unit("m")) -> Unit("m / s2")

Return tidal potential radial derivative.

Parameters:

• time (Time) – Time of observation.
• lat (Quantity) – Geocentric (spherical) latitude.
• lon (Quantity) – Geocentric (spherical) longitude.
• r (Quantity) – Geocentric radius (radial distance).

Returns:

potential_dr – Tidal potential radial derivative.

Return type:

Quantity

gradient(time: astropy.time.core.Time, lat: Unit("deg"), lon: Unit("deg"), r: Unit("m")) -> Unit("m / s2")

Return tidal potential gradient (tidal acceleration).

Parameters:

• time (Time) – Time of observation.
• lat (Quantity) – Geocentric (spherical) latitude.
• lon (Quantity) – Geocentric (spherical) longitude.
• r (Quantity) – Geocentric radius (radial distance).

Returns:

grad – Tidal potential radial derivative.

Return type:

Quantity

displacement(time: astropy.time.core.Time, lat: Unit("deg"), lon: Unit("deg"), r: Unit("m"), gravity: Unit("m / s2") = <<class 'astropy.constants.codata2018.CODATA2018'> name='Standard acceleration of gravity' value=9.80665 uncertainty=0.0 unit='m / s2' reference='CODATA 2018'>, elastic: bool = True) -> Unit("m")

Return direct tidal deformation for elastic or equilibrium Earth.

Direct effect is an equipotential surface deformation caused by tidal potential. Equilibrium tidal deformation is caused only by theoretical tide, elastic deformation takes into account the elastic properties of the Earth.

Parameters:

• time (Time) – Time of observation.
• lat (Quantity) – Geocentric (spherical) latitude.
• lon (Quantity) – Geocentric (spherical) longitude.
• r (Quantity) – Geocentric radius (radial distance).
• gravity (Quantity) – Mean global gravity of the Earth.
• elastic (bool, optional) – If True, elastic deformation (with approximate Love numbers) will be returned, otherwise equilibrium deformation. Default is True.

Returns:

• u_lat (~astropy.units.Quantity) – Deformation in the north-south direction.
• u_lon (~astropy.units.Quantity) – Deformation in the east-west direction.
• u_r (~astropy.units.Quantity) – Deformation in the radial direction.

deformation_potential(time: astropy.time.core.Time, lat: Unit("deg"), lon: Unit("deg"), r: Unit("m")) -> Unit("m2 / s2")

Return deformation tidal potential.

An incremental (deforamtion) potential is produced by the re-distribution of masses which is factorized by a number, the second Love number k, on the earth’s surface.

This is also called the indirect effect.

Parameters:

• time (Time) – Time of observation.
• lat (Quantity) – Geocentric (spherical) latitude.
• lon (Quantity) – Geocentric (spherical) longitude.
• r (Quantity) – Geocentric radius (radial distance).

Returns:

delta_potential – Deformation potential.

Return type:

Quantity

gravity(time: astropy.time.core.Time, lat: Unit("deg"), lon: Unit("deg"), r: Unit("m"), elastic: bool = True) -> Unit("m / s2")

Return tidal gravity variation.

Parameters:

• time (Time) – Time of observation.
• lat (Quantity) – Geocentric (spherical) latitude.
• lon (Quantity) – Geocentric (spherical) longitude.
• r (Quantity) – Geocentric radius (radial distance).
• elastic (bool, optional) – If True then the Earth is elastic (deformable) and gravity change is multiplied by the gravimetric factor specified in the class instance.

Returns:

delta_g – Tidal gravity variation.

Return type:

Quantity

tilt(time: astropy.time.core.Time, lat: Unit("deg"), lon: Unit("deg"), r: Unit("m"), azimuth: Unit("deg") = None, gravity: Unit("m / s2") = <<class 'astropy.constants.codata2018.CODATA2018'> name='Standard acceleration of gravity' value=9.80665 uncertainty=0.0 unit='m / s2' reference='CODATA 2018'>, elastic: bool = True) -> Unit(dimensionless)

Return tidal tilt.

Parameters:

• time (Time) – Time of observation.
• lat (Quantity) – Geocentric (spherical) latitude.
• lon (Quantity) – Geocentric (spherical) longitude.
• r (Quantity) – Geocentric radius (radial distance).
• azimuth (Quantity, optional) – If given then tilt will be returned in that direction. If None then full tilt angle will be returned. Default is None.
• elastic (bool, optional) – If True then the Earth is elastic (deformable) and tilt is multiplied by the combination of Love numbers (h - (1 + k)).
• gravity (Quantity, optional) – Mean global gravity of the Earth.

Returns:

tilt – Tidal tilt.

Return type:

Quantity

geoidal_height(time: astropy.time.core.Time, lat: Unit("deg"), lon: Unit("deg"), r: Unit("m"), gravity: Unit("m / s2") = <<class 'astropy.constants.codata2018.CODATA2018'> name='Standard acceleration of gravity' value=9.80665 uncertainty=0.0 unit='m / s2' reference='CODATA 2018'>) -> Unit("m")

Return tidal variation of the geoidal height.

Geoidal heights are affected by direct and indirect effects.

Parameters:

• time (Time) – Time of observation.
• lat (Quantity) – Geocentric (spherical) latitude.
• lon (Quantity) – Geocentric (spherical) longitude.
• r (Quantity) – Geocentric radius (radial distance).
• gravity (Quantity, optional) – Mean global gravity of the Earth.

Returns:

delta_geoid – Geoidal height tidal variation.

Return type:

Quantity

physical_height(time: astropy.time.core.Time, lat: Unit("deg"), lon: Unit("deg"), r: Unit("m"), gravity: Unit("m / s2") = <<class 'astropy.constants.codata2018.CODATA2018'> name='Standard acceleration of gravity' value=9.80665 uncertainty=0.0 unit='m / s2' reference='CODATA 2018'>) -> Unit("m")

Return tidal variation of the (absolute) physical heights.

Tidal variations of physically defined heights like orthometric, normal, or normal-orthometric heights can be derived as the difference between ellipsoidal and geoidal heights.

Parameters:

• time (Time) – Time of observation.
• lat (Quantity) – Geocentric (spherical) latitude.
• lon (Quantity) – Geocentric (spherical) longitude.
• r (Quantity) – Geocentric radius (radial distance).
• gravity (Quantity, optional) – Mean global gravity of the Earth.

Returns:

delta_h – Physical heights tidal variation.

Return type:

Quantity

physical_height_difference(time: astropy.time.core.Time, lat: Unit("deg"), lon: Unit("deg"), r: Unit("m"), azimuth: Unit("deg"), length: Unit("m"), gravity: Unit("m / s2") = <<class 'astropy.constants.codata2018.CODATA2018'> name='Standard acceleration of gravity' value=9.80665 uncertainty=0.0 unit='m / s2' reference='CODATA 2018'>) -> Unit("m")

Return tidal variation of the physical heights difference.

Parameters:

• time (Time) – Time of observation.
• lat (Quantity) – Geocentric (spherical) latitude.
• lon (Quantity) – Geocentric (spherical) longitude.
• r (Quantity) – Geocentric radius (radial distance).
• azimuth (Quantity) – Azimuth of the levelling line.
• length (Quantity) – Length of the levelling line.
• gravity (Quantity, optional) – Mean global gravity of the Earth.

Returns:

delta_h – Physical heights difference tidal variation.

Return type:

Quantity

Indices and tables

• Index
• Module Index
• Search Page