Как преобразовать от UTM до LatLng в Python или JavaScript

Я нашел следующий сайт: http://home.hiwaay.net/~taylorc/toolbox/geography/geoutm.html У него есть конвертер JavaScript, вы должны проверить алгоритм там. Со страницы:

Программисты: Исходный код JavaScript в этом документе может быть скопирован и использован без ограничений.

Согласно этой странице, UTM поддерживается proj4js.

http://trac.osgeo.org/proj4js/wiki/UserGuide#Supportedprojectionclasses

можно также хотеть смотреть на GDAL. gdal библиотека имеет превосходную поддержку Python, хотя это может быть немного излишества, если Вы только делаете преобразование проекции.

Я также новичок в этом и недавно изучал эту тему.

Вот метод, который я нашел, используя python gdal pacakge (пакет osr включен в gdal). Пакет gdal довольно мощный, но документация может быть лучше.

Это вытекает из обсуждения здесь: http://www.mail-archive.com/gdal-dev@lists.osgeo.org/msg12398.html

import osr

def transform_utm_to_wgs84(easting, northing, zone):
    utm_coordinate_system = osr.SpatialReference()
    utm_coordinate_system.SetWellKnownGeogCS("WGS84") # Set geographic coordinate system to handle lat/lon
    is_northern = northing > 0    
    utm_coordinate_system.SetUTM(zone, is_northern)

    wgs84_coordinate_system = utm_coordinate_system.CloneGeogCS() # Clone ONLY the geographic coordinate system 

    # create transform component
    utm_to_wgs84_transform = osr.CoordinateTransformation(utm_coordinate_system, wgs84_coordinate_system) # (<from>, <to>)
    return utm_to_wgs84_transform.TransformPoint(easting, northing, 0) # returns lon, lat, altitude

А вот метод преобразования из lat, lon в wgs84 (о чем сообщает большинство единиц gps) в utm:

def transform_wgs84_to_utm(lon, lat):    
    def get_utm_zone(longitude):
        return (int(1+(longitude+180.0)/6.0))

    def is_northern(latitude):
        """
        Determines if given latitude is a northern for UTM
        """
        if (latitude < 0.0):
            return 0
        else:
            return 1

    utm_coordinate_system = osr.SpatialReference()
    utm_coordinate_system.SetWellKnownGeogCS("WGS84") # Set geographic coordinate system to handle lat/lon  
    utm_coordinate_system.SetUTM(get_utm_zone(lon), is_northern(lat))

    wgs84_coordinate_system = utm_coordinate_system.CloneGeogCS() # Clone ONLY the geographic coordinate system 

    # create transform component
    wgs84_to_utm_transform = osr.CoordinateTransformation(wgs84_coordinate_system, utm_coordinate_system) # (<from>, <to>)
    return wgs84_to_utm_transform.TransformPoint(lon, lat, 0) # returns easting, northing, altitude    

Я также обнаружил, что если у вас уже установлен django / gdal и вы знаете код EPSG для зоны UTM, над которой вы работаете, вы можете просто использовать метод Point() transform () .

from django.contrib.gis.geos import Point
utm2epsg = {"54N": 3185, ...}
p = Point(lon, lat, srid=4326) # 4326 = WGS84 epsg code
p.transform(utm2epsg["54N"])

Javascript-версия Staale ответа

function utmToLatLng(zone, easting, northing, northernHemisphere){
        if (!northernHemisphere){
            northing = 10000000 - northing;
        }

        var a = 6378137;
        var e = 0.081819191;
        var e1sq = 0.006739497;
        var k0 = 0.9996;

        var arc = northing / k0;
        var mu = arc / (a * (1 - Math.pow(e, 2) / 4.0 - 3 * Math.pow(e, 4) / 64.0 - 5 * Math.pow(e, 6) / 256.0));

        var ei = (1 - Math.pow((1 - e * e), (1 / 2.0))) / (1 + Math.pow((1 - e * e), (1 / 2.0)));

        var ca = 3 * ei / 2 - 27 * Math.pow(ei, 3) / 32.0;

        var cb = 21 * Math.pow(ei, 2) / 16 - 55 * Math.pow(ei, 4) / 32;
        var cc = 151 * Math.pow(ei, 3) / 96;
        var cd = 1097 * Math.pow(ei, 4) / 512;
        var phi1 = mu + ca * Math.sin(2 * mu) + cb * Math.sin(4 * mu) + cc * Math.sin(6 * mu) + cd * Math.sin(8 * mu);

        var n0 = a / Math.pow((1 - Math.pow((e * Math.sin(phi1)), 2)), (1 / 2.0));

        var r0 = a * (1 - e * e) / Math.pow((1 - Math.pow((e * Math.sin(phi1)), 2)), (3 / 2.0));
        var fact1 = n0 * Math.tan(phi1) / r0;

        var _a1 = 500000 - easting;
        var dd0 = _a1 / (n0 * k0);
        var fact2 = dd0 * dd0 / 2;

        var t0 = Math.pow(Math.tan(phi1), 2);
        var Q0 = e1sq * Math.pow(Math.cos(phi1), 2);
        var fact3 = (5 + 3 * t0 + 10 * Q0 - 4 * Q0 * Q0 - 9 * e1sq) * Math.pow(dd0, 4) / 24;

        var fact4 = (61 + 90 * t0 + 298 * Q0 + 45 * t0 * t0 - 252 * e1sq - 3 * Q0 * Q0) * Math.pow(dd0, 6) / 720;

        var lof1 = _a1 / (n0 * k0);
        var lof2 = (1 + 2 * t0 + Q0) * Math.pow(dd0, 3) / 6.0;
        var lof3 = (5 - 2 * Q0 + 28 * t0 - 3 * Math.pow(Q0, 2) + 8 * e1sq + 24 * Math.pow(t0, 2)) * Math.pow(dd0, 5) / 120;
        var _a2 = (lof1 - lof2 + lof3) / Math.cos(phi1);
        var _a3 = _a2 * 180 / Math.PI;

        var latitude = 180 * (phi1 - fact1 * (fact2 + fact3 + fact4)) / Math.PI;

        if (!northernHemisphere){
          latitude = -latitude;
        }

        var longitude = ((zone > 0) && (6 * zone - 183.0) || 3.0) - _a3;

        var obj = {
              latitude : latitude,
              longitude: longitude
        };


        return obj;
      }

Ответ Стейла сработал для меня с небольшой модификацией - математический модуль не может обрабатывать серии Pandas, поэтому я заменил все математические функции на numpy.

Однако, проверяя в QGIS, я вижу разницу примерно в 4 м между координатами UTM и LAT / LON.

Код ниже:

import numpy as np

def utmToLatLng(zone, easting, northing, northernHemisphere=True):
    if not northernHemisphere:
        northing = 10000000 - northing

a = 6378137
e = 0.081819191
e1sq = 0.006739497
k0 = 0.9996

arc = northing / k0
mu = arc / (a * (1 - np.power(e, 2) / 4.0 - 3 * np.power(e, 4) / 64.0 - 5 * np.power(e, 6) / 256.0))

ei = (1 - np.power((1 - e * e), (1 / 2.0))) / (1 + np.power((1 - e * e), (1 / 2.0)))

ca = 3 * ei / 2 - 27 * np.power(ei, 3) / 32.0

cb = 21 * np.power(ei, 2) / 16 - 55 * np.power(ei, 4) / 32
cc = 151 * np.power(ei, 3) / 96
cd = 1097 * np.power(ei, 4) / 512
phi1 = mu + ca * np.sin(2 * mu) + cb * np.sin(4 * mu) + cc * np.sin(6 * mu) + cd * np.sin(8 * mu)

n0 = a / np.power((1 - np.power((e * np.sin(phi1)), 2)), (1 / 2.0))

r0 = a * (1 - e * e) / np.power((1 - np.power((e * np.sin(phi1)), 2)), (3 / 2.0))
fact1 = n0 * np.tan(phi1) / r0

_a1 = 500000 - easting
dd0 = _a1 / (n0 * k0)
fact2 = dd0 * dd0 / 2

t0 = np.power(np.tan(phi1), 2)
Q0 = e1sq * np.power(np.cos(phi1), 2)
fact3 = (5 + 3 * t0 + 10 * Q0 - 4 * Q0 * Q0 - 9 * e1sq) * np.power(dd0, 4) / 24

fact4 = (61 + 90 * t0 + 298 * Q0 + 45 * t0 * t0 - 252 * e1sq - 3 * Q0 * Q0) * np.power(dd0, 6) / 720

lof1 = _a1 / (n0 * k0)
lof2 = (1 + 2 * t0 + Q0) * np.power(dd0, 3) / 6.0
lof3 = (5 - 2 * Q0 + 28 * t0 - 3 * np.power(Q0, 2) + 8 * e1sq + 24 * np.power(t0, 2)) * np.power(dd0, 5) / 120
_a2 = (lof1 - lof2 + lof3) / np.cos(phi1)
_a3 = _a2 * 180 / np.pi

latitude = 180 * (phi1 - fact1 * (fact2 + fact3 + fact4)) / np.pi

if not northernHemisphere:
    latitude = -latitude

longitude = ((zone > 0) and (6 * zone - 183.0) or 3.0) - _a3

return (latitude, longitude)

Таким образом, я могу сделать это напрямую:

df['LAT'], df['LON']=utmToLatLng(31, df['X'], df['Y'], northernHemisphere=True)

Через CPAN есть модуль Perl под названием Geography :: NationalGrid, который может преобразовывать восток / север в широту / долготу. Это может помочь.

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