Converting from geographic coordinates to geomagnetic coordinates

Question

I am trying to convert geographic and geomagnetic coordinates. I found the following Prolog script, but I don't understand it enough to convert it myself. The target language is Java, but everything is clear (excellent, C, Python, VB).

http://idlastro.gsfc.nasa.gov/ftp/pro/astro/geo2mag.pro

If someone can either help with the conversion of this script, or explain exactly what he is doing (these array operations are puzzling me), I would really appreciate it.

thanks

Answer 1

Depending on the application, altitude can be an important variable in this coordinate transformation, since the geomagnetic coordinates are a representation of the Earth's dipole magnetic field.

In Python, you can easily convert geographic coordinates to geomagnetic (and vice versa) using SpacePy ( http://sourceforge.net/projects/spacepy/ ).

Since you are looking for source code to convert to Java, SpacePy is introducing the Fortran International Modelation Belt Environment Modeling (IRBEM) library, the source of which is available ( http://irbem.svn.sourceforge.net/viewvc/irbem)./web/index.html )

In Python, in case others are looking for a quick solution:

import spacepy.coordinates as coord from spacepy.time import Ticktock import numpy as np def geotomag(alt,lat,lon): #call with altitude in kilometers and lat/lon in degrees Re=6371.0 #mean Earth radius in kilometers #setup the geographic coordinate object with altitude in earth radii cvals = coord.Coords([np.float(alt+Re)/Re, np.float(lat), np.float(lon)], 'GEO', 'sph',['Re','deg','deg']) #set time epoch for coordinates: cvals.ticks=Ticktock(['2012-01-01T12:00:00'], 'ISO') #return the magnetic coords in the same units as the geographic: return cvals.convert('MAG','sph') 

Answer 2

I did this for python code and tried to check using this website http://wdc.kugi.kyoto-u.ac.jp/igrf/gggm/index.html . I found that

• The magnetic pole is the year 1995.
• Even if I set the above calculation to use the value for 1995, I did not quite understand correctly.

I used the value for Kyoto, Japan (35 N, 135.45 W). Web page calculation (25.18, -155.80). I got (25.33580652, -155.82724011). Therefore, I am not entirely sure that this can be really used ...

 import numpy as np from numpy import pi, cos, sin, arctan2, sqrt, dot def geo2mag(incoord): """geographic coordinate to magnetic coordinate: incoord is numpy array of shape (2,*) array([[glat0,glat1,glat2,...], [glon0,glon1,glon2,...]) where glat, glon are geographic latitude and longitude (or if you have only one point it is [[glat,glon]]) returns array([mlat0,mlat1,...], [mlon0,mlon1,...]]) """ # SOME 'constants'... lon = 288.59 # or 71.41W lat = 79.3 r = 1.0 # convert first to radians lon, lat = [x*pi/180 for x in lon,lat] glat = incoord[0] * pi / 180.0 glon = incoord[1] * pi / 180.0 galt = glat * 0. + r coord = np.vstack([glat,glon,galt]) # convert to rectangular coordinates x = coord[2]*cos(coord[0])*cos(coord[1]) y = coord[2]*cos(coord[0])*sin(coord[1]) z = coord[2]*sin(coord[0]) xyz = np.vstack((x,y,z)) # computer 1st rotation matrix: geo2maglon = np.zeros((3,3), dtype='float64') geo2maglon[0,0] = cos(lon) geo2maglon[0,1] = sin(lon) geo2maglon[1,0] = -sin(lon) geo2maglon[1,1] = cos(lon) geo2maglon[2,2] = 1. out = dot(geo2maglon , xyz) tomaglat = np.zeros((3,3), dtype='float64') tomaglat[0,0] = cos(.5*pi-lat) tomaglat[0,2] = -sin(.5*pi-lat) tomaglat[2,0] = sin(.5*pi-lat) tomaglat[2,2] = cos(.5*pi-lat) tomaglat[1,1] = 1. out = dot(tomaglat , out) mlat = arctan2(out[2], sqrt(out[0]*out[0] + out[1]*out[1])) mlat = mlat * 180 / pi mlon = arctan2(out[1], out[0]) mlon = mlon * 180 / pi outcoord = np.vstack((mlat, mlon)) return outcoord if __name__ == '__main__': mag = geo2mag(np.array([[79.3,288.59]]).T).T print mag # should be [90,0] mag = geo2mag(np.array([[90,0]]).T).T print mag # should be [79.3,*] mag = geo2mag(np.array([ [79.3,288.59], [90,0] ]).T).T print mag # should be [ [90,0]. [79.3,*] ] # kyoto, japan mag = geo2mag(np.array([[35.,135.45]]).T).T print mag # should be [25.18, -155.80], according to # this site using value for 1995 # http://wdc.kugi.kyoto-u.ac.jp/igrf/gggm/index.html