diff --git a/docs/XYZ.md b/docs/XYZ.md index b3737c52ac03dbe1551758a448b849179b04c0eb..497b656ad1ead59ade945c99d8df0b9aa7f1090a 100644 --- a/docs/XYZ.md +++ b/docs/XYZ.md @@ -1,102 +1,180 @@ -#Algorithm Theoretical Basis for "Geomag XYZ" +# Algorithm Theoretical Basis for "Geomag XYZ" E. Joshua Rigler <[erigler@usgs.gov](mailto:erigler@usgs.gov)> -#Summary +# Summary -Mathematical underpinnings and general algorithm considerations are presented for converting geomagnetic observations from so-called HEZ coordinates, used by the USGS geomagnetism program, into XYZ coordinates, used by a growing number of international geomagnetism programs, as well as various academic and commercial entities. Inverse transformations are also provided. +Mathematical underpinnings and general algorithm considerations are presented +for converting geomagnetic observations from so-called HEZ coordinates, used by +the USGS geomagnetism program, into XYZ coordinates, used by a growing number of +international geomagnetism programs, as well as various academic and commercial +entities. Inverse transformations are also provided. -#Background and Motivation +# Background and Motivation -Historically, the most common coordinate system used to specify measured geomagnetic fields has been HDZ, where: +Historically, the most common coordinate system used to specify measured +geomagnetic fields has been HDZ, where: -- H is the magnitude of the geomagnetic field vector tangential to the Earth’s surface; -- D is the declination, or clockwise angle from the vector pointing to the geographic north pole to the H vector; -- Z is the downward component of the geomagnetic field. +- `H` is the magnitude of the geomagnetic field vector tangential to the Earth's surface; +- `D` is the declination, or clockwise angle from the vector pointing to the geographic north pole to the H vector; +- `Z` is the downward component of the geomagnetic field. +This coordinate system is useful for navigation (it is the natural coordinate +system for a magnetic compass), and any scientific analysis conducted in a +geomagnetic field-aligned reference frame, but is somewhat awkward for most +other applications. A more generally useful set of coordinates for scientific +analysis and engineering applications is the XYZ system: -This coordinate system is useful for navigation (it is the natural coordinate system for a magnetic compass), and any scientific analysis conducted in a geomagnetic field-aligned reference frame, but is somewhat awkward for most other applications. A more generally useful set of coordinates for scientific analysis and engineering applications is the XYZ system: +- `X` is the magnitude of the geographic north pole component of the H vector; +- `Y` is the magnitude of the east component of the H vector; +- `Z` is the downward component of the geomagnetic field, same as before. +Conversion between these two coordinate systems involves relatively straight- +forward trigonometry (see Eqs. 1, 2, and 3). However, in practice, a 3-axis +magnetometer necessarily takes on a fixed orientation upon installation. For +USGS observatories, this is aligned with the average magnetic north vector and +downward, with the final axis completes a right-handed 3-dimensional coordinate +system (roughly eastward). This is often referred to as HEZ coordinates, but for +the remainder of this document we will refer to it as heZ, to avoid confusion +with more traditional definitions of H and E(==Y). +The purpose of this document then is to provide a mathematical and algorithmic +description of how one converts data measured in heZ coordinates to true HDZ, +and finally to XYZ. -- X points to the geographic north pole; -- Y points eastward; -- Z, as before, points downward. +# Math and Theory +First, following definitions in the previous section, the conversion from +cylindrical HDZ to Cartesian XYZ is very straight-forward trigonometry: +- <a name="eq1"></a>Equation 1:`X = H cos(D)` +- <a name="eq2"></a>Equation 2:`Y = H sin(D)` +- <a name="eq3"></a>Equation 3:`Z = Z` -Conversion between these two coordinate systems involves relatively straight-forward trigonometry (see Eqs. 1, 2, and 3). However, in practice, a 3-axis magnetometer necessarily takes on a fixed orientation upon installation. For USGS observatories, this is aligned with the average magnetic north vector and downward, with the final axis completes a right-handed 3-dimensional coordinate system (roughly eastward). This is often referred to as HEZ coordinates, but for the remainder of this document we will refer to it as heZ, to avoid confusion with more traditional definitions of H and E(==Y). +However, as noted previously, the USGS aligns its magnetometers with the +magnetic north upon installation at an observatory, meaning raw data is +generated in heZ coordinates, where "h" is is the primary axis in a fixed +reference frame, "e" is the secondary axis in this reference frame, and "Z" is +the tertiary axis, which remains common for all reference frames discussed in +this document. + -The purpose of this document then is to provide a mathematical and algorithmic description of how one converts data measured in heZ coordinates to true HDZ, and finally to XYZ. +The figure above illustrates how the same full magnetic field vector **F**, can +be represented in heZ, HDZ, and XYZ coordinates. Red objects are specific to the +magnetometer's reference frame, while blue objects are specific to the +geographic reference frame. Black is common to all frames considered here, and +dashed lines help define Cartesian grids. +One thing that is not labeled in this figure is the angle d (see Eq. 4](#eq4)), +which is the difference between declination D, and a declination baseline (D0, +or DECBAS). +The equations Eqs [4](#eq4), [5](#eq5), [6](#eq6) describe how to convert the horizontal components of a USGS +magnetometer's raw data element into more standard H and D components. -#Math and Theory +- <a name="eq4"></a>Equation 4:`d = arctan(e/h)` +- <a name="eq5"></a>Equation 5:`D = D0 + d` +- <a name="eq6"></a>Equation 6:`H = `sqrt(h*h + e*e)` = h / cos(d)` -First, following definitions in the previous section, the conversion from cylindrical HDZ to Cartesian XYZ is very straight-forward trigonometry: +To inverse transform from XY to HD: -- X = H cos(D) -- Y = H sin(D) -- Z = Z +- <a name="eq7"></a>Equation 7:`H = sqrt(X*X + Y*Y)` +- <a name="eq8"></a>Equation 8:`D = arctan(Y/X)` -However, as noted previously, the USGS aligns its magnetometers with the magnetic north upon installation at an observatory, meaning raw data is generated in heZ coordinates, where "h" is is the primary axis in a fixed reference frame, "e" is the secondary axis in this reference frame, and "Z" is the tertiary axis, which remains common for all reference frames discussed in this document. +...and from HD to he: - +- <a name="eq9"></a>Equation 9:`d = D - D0` +- <a name="eq10"></a>Equation 10:`h = sqrt(H*H / 1 + tan2(d)) = H cos(d)` +- <a name="eq10"></a>Equation 11:`e = h * tan(d)` -The figure above illustrates how the same full magnetic field vector **F**, can be represented in heZ, HDZ, and XYZ coordinates. Red objects are specific to the magnetometer’s reference frame, while blue objects are specific to the geographic reference frame. Black is common to all frames considered here, and dashed lines help define Cartesian grids. +It is worth noting that there is potential for mathematically undefined results +in several of the preceding equations, where infinite ratios are a possible +argument to the arctan() function. However, Python's Numpy package, and indeed +most modern math libraries, will return reasonable answers in such situations +(hint: arctan(Inf)==pi/2). Internally we use the atan2 function, that handles +this problem. -One thing that is not labeled in this figure is the angle d (see Eq. 4), which is the difference between declination D, and a declination baseline (D0, or DECBAS). +# Algorithm +The geomag algorithm references the 3 reference frames as geo for +geographic/cartesian, obs for observatory and mag for magnetic/cylindrical. In +reference to this article -The next 3 equations describe how to convert the horizontal components of a USGS magnetometer’s raw data element into more standard H and D components. +- geo is XYZ +- obs is heZ +- mag is HDZ -- d = arctan(e/h) -- D = D0 + d -- H = `sqrt(h*h + e*e)` = h / cos(d) +The underlying library provides calculations for both the basic conversions, +such as get_get_y_from_mag, which is based off of Y = H sin(D), and higher level +conversions, such as get_geo_from_mag. (Which converts HD to XY) -To inverse transform from XY to HD: +Upper libraries only provide higher level conversions, ie get_geo_from_mag. This +is the level most users should be accessing. -- H = `sqrt(X*X + Y*Y)` -- D = arctan(Y/X) +## Declination Angular Units -...and from HD to he: +The library internally uses radians for all angles, and factories convert into +this standard. You can optionally reference the +ChannelConverter.get_radians_from_minutes. -- d = D - D0 -- h = sqrt(H*H / 1 + tan2(d)) = H cos(d) +##Declination Baseline -It is worth noting that there is potential for mathematically undefined results in several of the preceding equations, where infinite ratios are a possible argument to the arctan() function. However, Python’s Numpy package, and indeed most modern math libraries, will return reasonable answers in such situations (hint: arctan(Inf)==pi/2). +The library IAGA factory attempts to parse DECBAS from the iaga comments +section. You can optionally link to metadata for where users can find this. +##Data Flags -#Algorithm -The geomag algorithm references the 3 reference frames as geo for geographic/cartesian, obs for observatory and mag for magnetic/cylindrical. In reference to this article -- geo is XYZ -- obs is heZ -- mag is HDZ - -The underlying library provides calculations for both the basic conversions, such as get_get_y_from_mag, which is based off of Y = H sin(D), and higher level conversions, such as get_geo_from_mag. (Which converts HD to XY) -Upper libraries only provide higher level conversions, ie get_geo_from_mag. This is the level most users should be accessing. +The library internally represents data gaps as NaN, and factories convert to +this where possible. +Note: In the algorithm, all channels are uppercase. We use context (ie obs vs. +mag), to differentiate between h,e and HD. This mirrors the various data +formats, (ie IAGA2002, etc). -#Practical Considerations +# Practical Considerations -##Magnetic Intensity Units +## Magnetic Intensity Units -It is understood that all raw data inputs are provided in units of nanoTesla (nT). Of course this is not required for the equations to be valid, but it is incumbent on the programmer to make sure all input data units are the same, and that output units are defined accurately. +It is understood that all raw data inputs are provided in units of nanoTesla +(nT). Of course this is not required for the equations to be valid, but it is +incumbent on the programmer to make sure all input data units are the same, and +that output units are defined accurately. -##Declination Angular Units +## Declination Angular Units -The equations in the preceding section are relatively simple to code up, with the standard caveat that angles must be appropriate for the trigonometric functions (e.g., if sin/cos/tan expect radians, be sure to provide parameters in radians). One thing that can potentially complicate this is that IAGA standards require declination angles to be in minutes of arc. Furthermore, D0 (DECBAS) is not very well-defined by IAGA standards, but is typically reported in tenths of minutes of arc. None of these are difficult to convert, but it is incumbent on the programmer to make sure they know what units are being used for the inputs. +The equations in the preceding section are relatively simple to code up, with +the standard caveat that angles must be appropriate for the trigonometric +functions (e.g., if sin/cos/tan expect radians, be sure to provide parameters in +radians). One thing that can potentially complicate this is that IAGA standards +require declination angles to be in minutes of arc. Furthermore, +D0 (DECBAS) is not very well-defined by IAGA standards, but is +typically reported in tenths of minutes of arc. None of these are difficult to +convert, but it is incumbent on the programmer to make sure they know what units +are being used for the inputs. -##Declination Baseline +## Declination Baseline -Declination baseline is not well-defined by IAGA standards. The typical method used to publish it with actual data is to include it in the metadata. For older IMF formatted files, it is part of the periodic block header. For IAGA2002 formatted file, it may be in the file header, but is not required. To the best of my knowledge, if it is not included, one should assume it is zero, but no corroborating documentation could be found to justify this statement. +Declination baseline is not well-defined by IAGA standards. The typical method +used to publish it with actual data is to include it in the metadata. For older +IMF formatted files, it is part of the periodic block header. For IAGA2002 +formatted file, it may be in the file header, but is not required. To the +best of my knowledge, if it is not included, one should assume it is zero, but +no corroborating documentation could be found to justify this statement. ## Declination in USGS Variations Data -The USGS variations data is actually published in hdZ coordinates. If one wishes to apply equations in the preceding section to USGS variations data, they must first convert "d" back into "e" via Eq. 11. +The USGS variations data is actually published in hdZ coordinates. If one +wishes to apply equations in the preceding section to USGS variations data, they +must first convert "d" back into "e" via Eq. 11. ## Data Flags -It should go without saying that bad data in one coordinate system is bad data in another. However, on occasion, operational USGS Geomagnetism Program code has been discovered where coordinate transformations were applied before checking data flags. This is not an issue if data flags are NaN (not-a-number values), but more typical for Geomag data, these are values like 99999, which can lead to seemingly valid, but erroneous values at times when the raw data were known to be bad. +It should go without saying that bad data in one coordinate system is bad data +in another. However, on occasion, operational USGS Geomagnetism Program code has +been discovered where coordinate transformations were applied +before checking data flags. This is not an issue if data flags are NaN +(not-a-number values), but more typical for Geomag data, these are values like +99999, which can lead to seemingly valid, but erroneous values at times when the +raw data were known to be bad.