Coordinate systems#
Gravity and magnetic data are usually observed at points above the Earth surface, which can be referenced through different coordinate systems. When processing these data, applying corrections, performing forward models and inversions, we should take into account which type of coordinate system are we using and if our choice introduces errors into these computations.
Harmonica can handle three different 3D coordinate systems:
Cartesian coordinates#
Cartesian (or plain) coordinates are referenced by a system defined by three orthogonal axis: two for the horizontal directions and one for the vertical one. Harmonica name the two horizontal axis as easting and northing, while the vertical one is referenced as upward, emphasizing its pointing direction. This means that positive values implies points above the zeroth plane, while negative values refer to points below. Harmonica assumes that the easting, northing and upward coordinates are given in meters.
For example, we can define a prism below the horizontal plane at zeroth height:
import verde as vd
# Define boundaries of the rectangular prism (in meters)
west, east, south, north = -20, 20, -20, 20
bottom, top = -40, -20
prism = [west, east, south, north, bottom, top]
Or a regular grid of points at 100 meters above the zeroth plane:
# Define a regular grid of observation points (coordinates in meters)
coordinates = vd.grid_coordinates(
region=(-40, 40, -40, 40), shape=(5, 5), extra_coords=100
)
easting, northing, upward = coordinates[:]
print("easting:", easting)
print("northing:", northing)
print("upward:", upward)
easting: [[-40. -20. 0. 20. 40.]
[-40. -20. 0. 20. 40.]
[-40. -20. 0. 20. 40.]
[-40. -20. 0. 20. 40.]
[-40. -20. 0. 20. 40.]]
northing: [[-40. -40. -40. -40. -40.]
[-20. -20. -20. -20. -20.]
[ 0. 0. 0. 0. 0.]
[ 20. 20. 20. 20. 20.]
[ 40. 40. 40. 40. 40.]]
upward: [[100. 100. 100. 100. 100.]
[100. 100. 100. 100. 100.]
[100. 100. 100. 100. 100.]
[100. 100. 100. 100. 100.]
[100. 100. 100. 100. 100.]]
Because the upward axis points in the upward direction, it’s easy to check if a given point is lower or higher than another one:
# Define two points
point_1 = (30, 20, -67)
point_2 = (30, 20, -58)
print("Point 1 is higher than point 2?", point_1[2] > point_2[2])
Point 1 is higher than point 2? False
Geodetic coordinates#
Geodetic (or geographic) coordinates use a reference ellipsoid for defining the longitude (\(\lambda\)), latitude (\(\varphi\)) and height (\(h\)) coordinates of points (see Figure: Geodetic coordinates). The reference ellipsoid is a mathematical surface that approximates the figure of the Earth and it’s used to define point coordinates.
Figure: Geodetic coordinates#
Reference ellipsoid and a point p along with a geocentric Cartesian system (\(X\), \(Y\), \(Z\)). Where \(a\) and \(b\) are the semimajor and semiminor axes of the ellipsoid, while the \(\lambda\), \(\varphi\) and \(h\) represent the geodetic coordinates of the point p in this geodetic coordinate system, where \(\lambda\) is the longitude, \(\varphi\) the latitude and \(h\) the height. The \(\phi\) is the spherical latitude of point p (see Spherical coordinates). This figure is a modified version of [Oliveira2021].
Harmonica assumes that longitude, latitude are given in decimal degrees and the ellipsoidal height is given in meters. Positive values of height refer to points outside the ellipsoid, while negative values refer to points that live inside it. Spatial data are usually given in geodetic coordinates, along with the reference ellipsoid on which they are defined.
For example, let’s define a regular grid of points (separated by equal
angles) at 2km above the ellipsoid using verde.
coordinates = vd.grid_coordinates(
region=(-70, -65, -35, -30), shape=(6, 6), extra_coords=2e3
)
longitude, latitude, height = coordinates[:]
print("longitude:", longitude)
print("latitude:", latitude)
print("height:", height)
longitude: [[-70. -69. -68. -67. -66. -65.]
[-70. -69. -68. -67. -66. -65.]
[-70. -69. -68. -67. -66. -65.]
[-70. -69. -68. -67. -66. -65.]
[-70. -69. -68. -67. -66. -65.]
[-70. -69. -68. -67. -66. -65.]]
latitude: [[-35. -35. -35. -35. -35. -35.]
[-34. -34. -34. -34. -34. -34.]
[-33. -33. -33. -33. -33. -33.]
[-32. -32. -32. -32. -32. -32.]
[-31. -31. -31. -31. -31. -31.]
[-30. -30. -30. -30. -30. -30.]]
height: [[2000. 2000. 2000. 2000. 2000. 2000.]
[2000. 2000. 2000. 2000. 2000. 2000.]
[2000. 2000. 2000. 2000. 2000. 2000.]
[2000. 2000. 2000. 2000. 2000. 2000.]
[2000. 2000. 2000. 2000. 2000. 2000.]
[2000. 2000. 2000. 2000. 2000. 2000.]]
Some processes need to know the reference ellipsoid used to define the
geodetic coordinates of points. boule offers several ellipsoids that
are commonly used on geophysical applications. Harmonica will ask for
a boule.Ellipsoid instance as argument if it needs the reference
ellipsoid.
Lets define the WGS84 ellipsoide using boule:
import boule as bl
ellipsoid = bl.WGS84
print(ellipsoid)
WGS84 - World Geodetic System (1984)
Oblate ellipsoid:
• Semimajor axis: 6378137 m
• Flattening: 0.0033528106647474805
• GM: 398600441800000.0 m³/s²
• Angular velocity: 7.292115e-05 rad/s
Source:
Hofmann-Wellenhof, B., & Moritz, H. (2006). Physical Geodesy (2nd,
corr. ed. 2006 edition ed.). Wien; New York: Springer.
Some other processes are only designed to work under Cartesian coordinates.
We can easily transform geodetic coordinates to Cartesian by applying
map projections. This can be done through pyproj.
As an example, lets project the longitude and latitude coordinates of the previously generated grid using a Mercator projection:
import pyproj
# Define a Mercator projection through pyproj
projection = pyproj.Proj(proj="merc", ellps="WGS84")
# Project the longitude and latitude coordinates of the grid points
longitude, latitude = coordinates[:2]
easting, northing = projection(longitude, latitude)
print("easting:", easting)
print("northing:", northing)
easting: [[-7792364.35552915 -7681044.86473588 -7569725.3739426 -7458405.88314933
-7347086.39235606 -7235766.90156278]
[-7792364.35552915 -7681044.86473588 -7569725.3739426 -7458405.88314933
-7347086.39235606 -7235766.90156278]
[-7792364.35552915 -7681044.86473588 -7569725.3739426 -7458405.88314933
-7347086.39235606 -7235766.90156278]
[-7792364.35552915 -7681044.86473588 -7569725.3739426 -7458405.88314933
-7347086.39235606 -7235766.90156278]
[-7792364.35552915 -7681044.86473588 -7569725.3739426 -7458405.88314933
-7347086.39235606 -7235766.90156278]
[-7792364.35552915 -7681044.86473588 -7569725.3739426 -7458405.88314933
-7347086.39235606 -7235766.90156278]]
northing: [[-4139372.7622473 -4139372.7622473 -4139372.7622473 -4139372.7622473
-4139372.7622473 -4139372.7622473 ]
[-4004909.10948031 -4004909.10948031 -4004909.10948031 -4004909.10948031
-4004909.10948031 -4004909.10948031]
[-3872033.73289718 -3872033.73289718 -3872033.73289718 -3872033.73289718
-3872033.73289718 -3872033.73289718]
[-3740670.1135821 -3740670.1135821 -3740670.1135821 -3740670.1135821
-3740670.1135821 -3740670.1135821 ]
[-3610745.18533098 -3610745.18533098 -3610745.18533098 -3610745.18533098
-3610745.18533098 -3610745.18533098]
[-3482189.08540862 -3482189.08540862 -3482189.08540862 -3482189.08540862
-3482189.08540862 -3482189.08540862]]
Remember that this process implies projecting the geodetic coordinates onto
a flat surface, what carries projection errors. For small regions, these
errors are small, but for regional and global regions, these can heavily
increase. Projections can also be used to recover geodetic coordinates from
Cartesian ones, by setting the inverse argument to True.
Spherical coordinates#
Spherical coordinates (a.k.a spherical geocentric coordinates) are defined by a coordinate system whose origin is located on the center of the Earth. Each point can be represented by its longitude (\(\lambda\)), spherical latitude (\(\phi\)) and radius (\(r\)) (see Figure: Spherical coordinates).
Important
The longitude coordinates defined in spherical coordinates and in geodetic coordinates are equivalent. Nevertheless, the spherical latitude and the (geodetic) latitude are not. They would be the same if the reference ellipsoid were a sphere.
The longitude and spherical latitude are angles given in decimal degrees, while the radius is the Euclidean distance between the point and the origin of the system (in meters). Although this reference system is rarely used for storing data, it’s used for some non-Cartesian forward models, like tesseroids (spherical prisms).
Figure: Spherical coordinates#
Point p defined in a spherical coordinate system, where \(\lambda\) is the longitude, \(\phi\) the latitude and \(r\) the radius. The spherical coordinates are defined upon a geocentric Cartesian system (\(X\), \(Y\), \(Z\)) whose origin is located in the center of the Earth.
Let’s define a regular grid of points in spherical coordinates, located at the same radius equal to the mean radius of the Earth.
coordinates = vd.grid_coordinates(
region=(-70, -65, -35, -30),
shape=(6, 6),
extra_coords=ellipsoid.mean_radius,
)
longitude, sph_latitude, radius = coordinates[:]
print("longitude:", longitude)
print("spherical latitude:", sph_latitude)
print("radius:", radius)
longitude: [[-70. -69. -68. -67. -66. -65.]
[-70. -69. -68. -67. -66. -65.]
[-70. -69. -68. -67. -66. -65.]
[-70. -69. -68. -67. -66. -65.]
[-70. -69. -68. -67. -66. -65.]
[-70. -69. -68. -67. -66. -65.]]
spherical latitude: [[-35. -35. -35. -35. -35. -35.]
[-34. -34. -34. -34. -34. -34.]
[-33. -33. -33. -33. -33. -33.]
[-32. -32. -32. -32. -32. -32.]
[-31. -31. -31. -31. -31. -31.]
[-30. -30. -30. -30. -30. -30.]]
radius: [[6370994.40182758 6370994.40182758 6370994.40182758 6370994.40182758
6370994.40182758 6370994.40182758]
[6370994.40182758 6370994.40182758 6370994.40182758 6370994.40182758
6370994.40182758 6370994.40182758]
[6370994.40182758 6370994.40182758 6370994.40182758 6370994.40182758
6370994.40182758 6370994.40182758]
[6370994.40182758 6370994.40182758 6370994.40182758 6370994.40182758
6370994.40182758 6370994.40182758]
[6370994.40182758 6370994.40182758 6370994.40182758 6370994.40182758
6370994.40182758 6370994.40182758]
[6370994.40182758 6370994.40182758 6370994.40182758 6370994.40182758
6370994.40182758 6370994.40182758]]
We can convert spherical coordinates to geodetic ones through
boule.Ellipsoid.spherical_to_geodetic:
coordinates_geodetic = ellipsoid.spherical_to_geodetic(*coordinates)
longitude, latitude, height = coordinates_geodetic[:]
print("longitude:", longitude)
print("latitude:", latitude)
print("height:", height)
longitude: [[-70. -69. -68. -67. -66. -65.]
[-70. -69. -68. -67. -66. -65.]
[-70. -69. -68. -67. -66. -65.]
[-70. -69. -68. -67. -66. -65.]
[-70. -69. -68. -67. -66. -65.]
[-70. -69. -68. -67. -66. -65.]]
latitude: [[-35.18102907 -35.18102907 -35.18102907 -35.18102907 -35.18102907
-35.18102907]
[-34.17864869 -34.17864869 -34.17864869 -34.17864869 -34.17864869
-34.17864869]
[-33.17604943 -33.17604943 -33.17604943 -33.17604943 -33.17604943
-33.17604943]
[-32.17323438 -32.17323438 -32.17323438 -32.17323438 -32.17323438
-32.17323438]
[-31.17020688 -31.17020688 -31.17020688 -31.17020688 -31.17020688
-31.17020688]
[-30.16697054 -30.16697054 -30.16697054 -30.16697054 -30.16697054
-30.16697054]]
height: [[ -83.44704281 -83.44704281 -83.44704281 -83.44704281
-83.44704281 -83.44704281]
[ -432.49826962 -432.49826962 -432.49826962 -432.49826962
-432.49826962 -432.49826962]
[ -776.71701812 -776.71701812 -776.71701812 -776.71701812
-776.71701812 -776.71701812]
[-1115.68145507 -1115.68145507 -1115.68145507 -1115.68145507
-1115.68145507 -1115.68145507]
[-1448.97598505 -1448.97598505 -1448.97598505 -1448.97598505
-1448.97598505 -1448.97598505]
[-1776.19176956 -1776.19176956 -1776.19176956 -1776.19176956
-1776.19176956 -1776.19176956]]
While the conversion of spherical coordinates into geodetic ones can be
carried out through boule.Ellipsoid.geodetic_to_spherical:
coordinates_spherical = ellipsoid.geodetic_to_spherical(*coordinates_geodetic)
longitude, sph_latitude, radius = coordinates_spherical[:]
print("longitude:", longitude)
print("spherical latitude:", sph_latitude)
print("radius:", radius)
longitude: [[-70. -69. -68. -67. -66. -65.]
[-70. -69. -68. -67. -66. -65.]
[-70. -69. -68. -67. -66. -65.]
[-70. -69. -68. -67. -66. -65.]
[-70. -69. -68. -67. -66. -65.]
[-70. -69. -68. -67. -66. -65.]]
spherical latitude: [[-35. -35. -35. -35. -35. -35.]
[-34. -34. -34. -34. -34. -34.]
[-33. -33. -33. -33. -33. -33.]
[-32. -32. -32. -32. -32. -32.]
[-31. -31. -31. -31. -31. -31.]
[-30. -30. -30. -30. -30. -30.]]
radius: [[6370994.40182757 6370994.40182757 6370994.40182757 6370994.40182757
6370994.40182757 6370994.40182757]
[6370994.40182757 6370994.40182757 6370994.40182757 6370994.40182757
6370994.40182757 6370994.40182757]
[6370994.40182758 6370994.40182758 6370994.40182758 6370994.40182758
6370994.40182758 6370994.40182758]
[6370994.40182757 6370994.40182757 6370994.40182757 6370994.40182757
6370994.40182757 6370994.40182757]
[6370994.40182758 6370994.40182758 6370994.40182758 6370994.40182758
6370994.40182758 6370994.40182758]
[6370994.40182758 6370994.40182758 6370994.40182758 6370994.40182758
6370994.40182758 6370994.40182758]]