6. Reference Systems

Cartesian coordinate systems are represented by three perpendicular axes that are defined with respect to an origin, and with directions defined for each axis [101]. Locations of planetary bodies and satellites can be determined in an Earth-centered Earth-Fixed coordinate system (ECEF) [59, 66]. ECEF is a Cartesian coordinate system with \(x\), \(y\), and \(z\) defined with respect to the Earth’s center of mass. The \(z\) axis is aligned with the Earth’s rotation axis, the \(x\) axis is aligned with the intersection of the prime meridian and the equator, and the \(y\) axis is aligned with 90° east longitude and the equator. The \(xy\) plane is also called the equatorial plane.

Figure 6.1: Geocentric Cartesian Coordinate System

As the shape of the Earth can be well approximated by an ellipsoid, positions on the Earth’s surface are often described in terms of geodetic coordinates. Geodetic coordinates (longitude \(\lambda\), latitude \(\varphi\), and height \(h\)) are used to describe the position of a point on the Earth with respect to a defined ellipsoid. Changing the terrestrial reference system can involve both translations and rotations of the reference system [101]. One method of transformation involves converting from a geographic coordinate system into a Cartesian coordinate system, and then performing matrix transformations [Equation 6.1].

(6.1)\[\begin{split}(\lambda, \varphi, h) &\rightarrow (x, y, z) \\ (x, y, z) &\rightarrow (x', y', z') \\ (x', y', z') &\rightarrow (\lambda', \varphi', h')\end{split}\]

The transformation from ellipsoidal coordinates of a point in space to Cartesian coordinates is calculated by:

(6.2)\[\begin{split}x &= (N + h)\cos{\varphi}\cos{\lambda}\\ y &= (N + h)\cos{\varphi}\sin{\lambda}\\ z &= \left((1 - e^2) N + h \right) \sin{\varphi}\end{split}\]

(6.3)\[N = \frac{a}{\sqrt{1 - e^2 \sin^2{\varphi}}}\]

(6.4)\[e^2 = 2f - f^2\]

where \(N\) is the radius of curvature in the prime vertical [Equation 6.3], \(e\) is the ellipsoidal eccentricity, with the squared eccentricity \(e^2\) given by [Equation 6.4], \(a\) is the semi-major axis of the ellipsoid, and \(f\) is the ellipsoidal flattening [40, 101]. Ellipsoid definitions typically specify the semi-major axis (\(a\)) and flattening (\(f\)), and datum definitions additionally include the coordinate system origin [40, 101]. In general, objects move with respect to the coordinate system, and the coordinate system itself moves and rotates in space [101]. Coordinate system definitions, such as the International Terrestrial Reference Frame (ITRF), will often include a time component to account for these changes.

6.1. Geocentric Coordinates

Similar to ECEF cartesian coordinates, geocentric coordinates are defined with respect to the center of the Earth [40, 92]. This is in contrast to geodetic coordinates, which are defined to have latitudes normal to the surface of the Earth.

Figure 6.3: Geodetic and Geocentric Coordinates

Geocentric coordinates are used to estimate Spherical Harmonics coefficients, and for performing coordinate system rotations. Geocentric longitudes are identical to geodetic longitudes, but geocentric latitudes can differ from geodetic latitudes by approximately 0.2 degrees. For applications using a spherical Earth model, the geocentric and geodetic latitudes are identical.

Figure 6.5: Difference between Geodetic and Geocentric Latitude

6.2. Geoid Height

Compared to the reference ellipsoid, a better representation of the figure of the Earth can be defined based on the Earth’s gravitational field [100]. The geoid is an equipotential surface set to coincide with an idealized global mean sea level (i.e. if the oceans were at rest) [40]. It serves as the reference surface to describe any topographic heights above (or below) mean sea level. As with Solid Earth Tides, the geoid has a permanent tide component due to the Earth being in the presence of the Sun and Moon [54, 100].

The distances between the geoid and the reference ellipsoid are called geoid height or geoidal undulations (\(N\)), and the distances between the geoid and points on the Earth’s surface are called orthometric heights (\(H\)) [40].

Figure 6.7: Relationship between the reference ellipsoid, geoid, and surface topography [100, 113]

6.3. Celestial Reference Systems

Celestial reference systems are used to describe the positions of celestial bodies in the sky. Transforming between celestial (\(\mathbf{x}_{CRS}\)) and terrestrial (\(\mathbf{x}_{TRS}\)) reference systems involves a set of transformation matrices for frame bias (\(\mathbf{B}\)), precession (\(\mathbf{P}\)), nutation (\(\mathbf{N}\)), Earth’s rotation (\(\mathbf{T}\)), and polar motion (\(\mathbf{W}\)) [11, 12, 101].

(6.5)\[\mathbf{x}_{CRS} = \mathbf{B}\ \mathbf{P}\ \mathbf{N}\ \mathbf{T}\ \mathbf{W}\ \mathbf{x}_{TRS}\]

In pyTMD, these transformations are used to convert planetary ephemerides from a celestial reference frame to a terrestrial reference frame.

6.4. Barycentric Coordinates

Some tide models are computed using the finite-element method and defined on unstructured meshes. These finite element meshes can use first-order (linear) or higher-order elements. Linear triangular elements have their three nodes at the vertices (\(N_1\), \(N_2\), \(N_3\)). Quadratic triangular elements have three additional nodes at the edge midpoints. The position of any point \(P\) within a triangle can be expressed in terms areal weights known as barycentric coordinates \((\xi, \eta, \lambda)\):

(6.6)\[P = \xi \, P_1 + \eta \, P_2 + \lambda \, P_3\]

Geometrically, these weights are derived from the ratios of sub-triangle areas (\(A_1\), \(A_2\), \(A_3\)) with respect to the total triangle area (\(A\)):

(6.7)\[\xi = \frac{A_1}{A}, \qquad \eta = \frac{A_2}{A}, \qquad \lambda = \frac{A_3}{A} = 1 - \xi - \eta\]

Here, \(A_1\) is the area of the sub-triangle opposite \(N_1\), \(A_2\) is the area of the sub-triangle opposite \(N_2\), and \(A_3\) is the area of the sub-triangle opposite \(N_3\).

Figure 6.8: Linear and quadratic triangular finite elements

The shape functions of finite elements have two important properties:

• Kronecker delta property: the shape functions are equal to one at their respective nodes and equal to zero at all other nodes

• Partition of unity property: the sum of the shape functions is equal to 1 for all points within the element

Because of these properties, barycentric coordinates can be used to help determine if a point is inside any given triangle:

• If all three barycentric coordinates are within \([0, 1]\), then the point is inside the element.

• If any coordinate is negative or beyond that range, then the point lies outside that element.

Important

Knowing the order of the nodes for each element is crucial for calculating the triangle areas and interpolation shape functions.

The standard order for nodes in linear triangular elements is counterclockwise. For a global finite element mesh, an element with a clockwise order might straddle the model boundary. pyTMD checks the winding number of each element to determine the polygon orientation (clockwise or counterclockwise). Quadratic elements can be more complicated and have different node orderings between the vertices and midside edges. Some standard node orderings for quadratic elements are:

• counterclockwise (\(V_1\rightarrow\) \(E_{12}\rightarrow\) \(V_2\rightarrow\) \(E_{23}\rightarrow\) \(V_3\rightarrow\) \(E_{31}\))

• vertices-to-midpoints (\(V_1\rightarrow\) \(V_2\rightarrow\) \(V_3\rightarrow\) \(E_{12}\rightarrow\) \(E_{23}\rightarrow\) \(E_{31}\))

• vertices-to-opposite-midpoints (\(V_1\rightarrow\) \(V_2\rightarrow\) \(V_3\rightarrow\) \(E_{23}\rightarrow\) \(E_{31}\rightarrow\) \(E_{12}\))

Note

For quadratic elements, pyTMD presently only supports the counterclockwise node order for the unstructured (“native”) FES models [52].