Преобразование из географических в геомагнитные координаты

Question

Я пытаюсь преобразовать географические и геомагнитные координаты.Я нашел следующий сценарий Пролога, но я не понимаю его достаточно, чтобы преобразовать его сам.Целевым языком является Java, но все, что понятно, прекрасно (C, Python, VB и т. Д.).

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

Если кто-то может помочь с преобразованием этого скрипта или объяснить, чтоименно это и делает (эти операции с массивами меня сбивают с толку), я бы очень признателен.

Спасибо

Answer 1

В зависимости от применения высота может быть важной переменной в этом преобразовании координат, поскольку геомагнитные координаты являются отображением дипольного магнитного поля Земли.

В Python вы можете легко преобразовать географические координаты в геомагнитные (и наоборот).и наоборот) с SpacePy (http://sourceforge.net/projects/spacepy/).

Поскольку вы ищете исходный код для преобразования в Java, SpacePy реализует библиотеку Fortran International Radiation Belt Environment Modeling (IRBEM), источник которой доступен (http://irbem.svn.sourceforge.net/viewvc/irbem/web/index.html)

В Python, если другие ищут быстрое решение:

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

Я сделал это в коде Python и попытался проверить с помощью этого сайта http://wdc.kugi.kyoto -u.ac.jp / igrf / gggm / index.html . Я обнаружил, что

1. Магнитный полюс такой же, как в 1995 году.
2. Даже если я установлю вышеприведенные вычисления для использования значения за 1995 год, я не получу правильного соответствия.

Я использовал значение для Киото, Япония (35N, 135,45 Вт). Рассчитанная веб-страница (25,18, -155.80). Я получил (25.33580652, -155.82724011). Так что я не совсем уверен, может ли это быть полезным ...

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