Convert from various non-Magnetic Coordinate Systems
Usually non-magnetic coordinates for space sciences are provided in geodetic
(geographic) coordinates, where the latitude is the angle between the
equatorial plane and the normal to the ellipsoid surface at the desired
location. There are several different reference ellipsoids that may be used,
but the most common (and the one used by Apex
is WGS84, the
World Geodetic System 1984 (as described in, for example, Snay and Soler [1]).
However, it frequently makes sense for a particular instrument to use a
different reference ellipsoid or another type of coordinate system.
Different types of geodetic or geocentric coordinate conversions are not supported within apexpy, because they are not a magnetic coordinate transformations. We recommend first converting your other non-magnetic coordinates to geodetic WGS84 coordinates and then performing the desired conversion to apex or quasi-dipole coordinates.
Geocentric Example
One commonly encountered alternative to geodetic latitude is geocentric latitude. Geocentric latitude is the angle between the equatorial plane and a line joining the center of the Earth and the desired location. If you have data in geocentric coordinates you can convert them to geodetic using a simple equation layed out in equation 3-28 of Snyder [2].
import numpy as np
def gc_to_gd(gc_lat, e_sq):
"""Convert from geocentric to geodetic
Parameters
----------
gc_lat : array-like
Geocentric latitude in degrees
e_sq : float
First eccentricity squared, unitless
Returns
-------
gd_lat : array-like
Geodetic latitude in degrees
"""
gd_lat = np.arctan(np.tan(np.radians(gc_lat)) / (1.0 - e_sq))
return np.degrees(gd_lat)
The function above requires the first eccentricity of the reference ellipsoid.
This example uses the pyproj
library [3] to get the WGS84 ellipsoid data,
but the function shown will take any float. This lets you decide the level of
precision you need in your coordinate transformation.
import pyproj
# Initialize the WGS84 reference ellipsoid
wgs84_geod = pyproj.crs.GeographicCRS(name='WGS84').get_geod()
print("{:.12f}".format(wgs84_geod.es))
0.006694379990
Now, try converting from geocentric to quasi-dipole coordinates. In this example you need to supply height, because the coordinate transformation from geodetic to quasi-dipole takes place by converting from geodetic to apex and then from apex to quasi-dipole.
import apexpy
# Define the starting values
year = 2015.3
gc_lat = 45.0
glon = 0.0
height = 0.0
# Get the quasi-dipole coordiantes
apex_out = apexpy.Apex(date=year)
qlat, qlon = apex_out.geo2qd(gc_to_gd(gc_lat, wgs84_geod.es), glon, height)
print("{:.12f}, {:.12f}".format(qlat, qlon))
39.852313995361, 76.711242675781