Contents

Coordinate systems

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 showing a reference ellipsoid along with an observation point "p". It also shows the geodetic coordinates of this observation point: its longitude, latitude and geodetic height.

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)

Ellipsoid(name='WGS84', semimajor_axis=6378137, flattening=0.0033528106647474805, geocentric_grav_const=398600441800000.0, angular_velocity=7.292115e-05, long_name='World Geodetic System 1984', reference='Hofmann-Wellenhof, B., & Moritz, H. (2006). Physical Geodesy (2nd, corr. ed. 2006 edition ed.). Wien\u202f; 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 showing an observation point "p" defined in a spherical coordinate system.

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: [[6371008.77141506 6371008.77141506 6371008.77141506 6371008.77141506
  6371008.77141506 6371008.77141506]
 [6371008.77141506 6371008.77141506 6371008.77141506 6371008.77141506
  6371008.77141506 6371008.77141506]
 [6371008.77141506 6371008.77141506 6371008.77141506 6371008.77141506
  6371008.77141506 6371008.77141506]
 [6371008.77141506 6371008.77141506 6371008.77141506 6371008.77141506
  6371008.77141506 6371008.77141506]
 [6371008.77141506 6371008.77141506 6371008.77141506 6371008.77141506
  6371008.77141506 6371008.77141506]
 [6371008.77141506 6371008.77141506 6371008.77141506 6371008.77141506
  6371008.77141506 6371008.77141506]]

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.18102866 -35.18102866 -35.18102866 -35.18102866 -35.18102866
  -35.18102866]
 [-34.17864829 -34.17864829 -34.17864829 -34.17864829 -34.17864829
  -34.17864829]
 [-33.17604904 -33.17604904 -33.17604904 -33.17604904 -33.17604904
  -33.17604904]
 [-32.17323399 -32.17323399 -32.17323399 -32.17323399 -32.17323399
  -32.17323399]
 [-31.17020649 -31.17020649 -31.17020649 -31.17020649 -31.17020649
  -31.17020649]
 [-30.16697016 -30.16697016 -30.16697016 -30.16697016 -30.16697016
  -30.16697016]]
height: [[  -69.07752705   -69.07752705   -69.07752705   -69.07752705
    -69.07752705   -69.07752705]
 [ -418.12875198  -418.12875198  -418.12875198  -418.12875198
   -418.12875198  -418.12875198]
 [ -762.34749847  -762.34749847  -762.34749847  -762.34749847
   -762.34749847  -762.34749847]
 [-1101.31193327 -1101.31193327 -1101.31193327 -1101.31193327
  -1101.31193327 -1101.31193327]
 [-1434.60646097 -1434.60646097 -1434.60646097 -1434.60646097
  -1434.60646097 -1434.60646097]
 [-1761.8222431  -1761.8222431  -1761.8222431  -1761.8222431
  -1761.8222431  -1761.8222431 ]]

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: [[6371008.77141506 6371008.77141506 6371008.77141506 6371008.77141506
  6371008.77141506 6371008.77141506]
 [6371008.77141506 6371008.77141506 6371008.77141506 6371008.77141506
  6371008.77141506 6371008.77141506]
 [6371008.77141506 6371008.77141506 6371008.77141506 6371008.77141506
  6371008.77141506 6371008.77141506]
 [6371008.77141506 6371008.77141506 6371008.77141506 6371008.77141506
  6371008.77141506 6371008.77141506]
 [6371008.77141506 6371008.77141506 6371008.77141506 6371008.77141506
  6371008.77141506 6371008.77141506]
 [6371008.77141506 6371008.77141506 6371008.77141506 6371008.77141506
  6371008.77141506 6371008.77141506]]

previous

Plot a prism layer in 3D

next

Forward modelling