.. _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: - :ref:Cartesian coordinates  - :ref:Geodetic coordinates  - :ref:Spherical coordinates  .. _cartesian_coordinates: 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: .. jupyter-execute:: 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: .. jupyter-execute:: # 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) 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: .. jupyter-execute:: # 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]) .. _geodetic_coordinates: Geodetic coordinates -------------------- Geodetic (or geographic) coordinates use a reference ellipsoid for defining the *longitude* (:math:\lambda), *latitude* (:math:\varphi) and *height* (:math:h) coordinates of points (see :ref:geodetic coordinates figure). The reference ellipsoid is a mathematical surface that approximates the figure of the Earth and it's used to define point coordinates. .. figure:: ../_static/figures/geodetic-coordinate-system.svg :name: geodetic coordinates figure :width: 90% :alt: 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 (:math:X, :math:Y, :math:Z). Where :math:a and :math:b are the semimajor and semiminor axes of the ellipsoid, while the :math:\lambda, :math:\varphi and :math:h represent the geodetic coordinates of the point **p** in this geodetic coordinate system, where :math:\lambda is the *longitude*, :math:\varphi the *latitude* and :math:h the *height*. The :math:\phi is the *spherical latitude* of point **p** (see :ref: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 :mod:verde. .. jupyter-execute:: 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) Some processes need to know the reference ellipsoid used to define the geodetic coordinates of points. :mod:boule offers several ellipsoids that are commonly used on geophysical applications. Harmonica will ask for a :class:boule.Ellipsoid instance as argument if it needs the reference ellipsoid. Lets define the WGS84 ellipsoide using :mod:boule: .. jupyter-execute:: import boule as bl ellipsoid = bl.WGS84 print(ellipsoid) 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 :mod:pyproj. As an example, lets project the *longitude* and *latitude* coordinates of the previously generated grid using a Mercator projection: .. jupyter-execute:: 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) 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 --------------------- 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* (:math:\lambda), *spherical latitude* (:math:\phi) and *radius* (:math:r) (see :ref:spherical coordinates figure). .. 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:: ../_static/figures/spherical-coordinate-system.svg :name: spherical coordinates figure :width: 50% :alt: Figure showing an observation point "p" defined in a spherical coordinate system. Figure: Spherical coordinates Point **p** defined in a spherical coordinate system, where :math:\lambda is the *longitude*, :math:\phi the *latitude* and :math:r the *radius*. The spherical coordinates are defined upon a geocentric Cartesian system (:math:X, :math:Y, :math: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*. .. jupyter-execute:: 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) We can convert spherical coordinates to geodetic ones through :meth:boule.Ellipsoid.spherical_to_geodetic: .. jupyter-execute:: coordinates_geodetic = ellipsoid.spherical_to_geodetic(*coordinates) longitude, latitude, height = coordinates_geodetic[:] print("longitude:", longitude) print("latitude:", latitude) print("height:", height) While the conversion of spherical coordinates into geodetic ones can be carried out through :meth:boule.Ellipsoid.geodetic_to_spherical: .. jupyter-execute:: 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)