Cómo convertir de UTM a LatLng en Python o Javascript
https://stackoverflow.com/questions/343865
-
19-08-2019 - |
italiano
english
français
española
中国
日本の
العربية
Deutsch
한국어
Português
RussianPregunta
Tengo un montón de archivos con coordenadas en forma UTM. Para cada coordenada tengo este, norte y zona. Necesito convertir esto a LatLng para usarlo con la API de Google Map para mostrar la información en un mapa.
He encontrado algunas calculadoras en línea que hacen esto, pero no hay código o bibliotecas reales. http://trac.osgeo.org/proj4js/ es una biblioteca de proyección para Javascript, pero está buscando en la demostración no incluye proyección UTM.
Todavía estoy bastante fresco para todo el dominio SIG, así que lo que quiero es algo:
(lat,lng) = transform(easting, northing, zone)
Solución
Terminé buscando el código Java de IBM que lo resolvió: http://www.ibm.com/developerworks/java/library/j-coordconvert/index.html
Solo como referencia, aquí está mi implementación en Python del método que necesitaba:
import math
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 - math.pow(e, 2) / 4.0 - 3 * math.pow(e, 4) / 64.0 - 5 * math.pow(e, 6) / 256.0))
ei = (1 - math.pow((1 - e * e), (1 / 2.0))) / (1 + math.pow((1 - e * e), (1 / 2.0)))
ca = 3 * ei / 2 - 27 * math.pow(ei, 3) / 32.0
cb = 21 * math.pow(ei, 2) / 16 - 55 * math.pow(ei, 4) / 32
cc = 151 * math.pow(ei, 3) / 96
cd = 1097 * math.pow(ei, 4) / 512
phi1 = mu + ca * math.sin(2 * mu) + cb * math.sin(4 * mu) + cc * math.sin(6 * mu) + cd * math.sin(8 * mu)
n0 = a / math.pow((1 - math.pow((e * math.sin(phi1)), 2)), (1 / 2.0))
r0 = a * (1 - e * e) / math.pow((1 - math.pow((e * math.sin(phi1)), 2)), (3 / 2.0))
fact1 = n0 * math.tan(phi1) / r0
_a1 = 500000 - easting
dd0 = _a1 / (n0 * k0)
fact2 = dd0 * dd0 / 2
t0 = math.pow(math.tan(phi1), 2)
Q0 = e1sq * math.pow(math.cos(phi1), 2)
fact3 = (5 + 3 * t0 + 10 * Q0 - 4 * Q0 * Q0 - 9 * e1sq) * math.pow(dd0, 4) / 24
fact4 = (61 + 90 * t0 + 298 * Q0 + 45 * t0 * t0 - 252 * e1sq - 3 * Q0 * Q0) * math.pow(dd0, 6) / 720
lof1 = _a1 / (n0 * k0)
lof2 = (1 + 2 * t0 + Q0) * math.pow(dd0, 3) / 6.0
lof3 = (5 - 2 * Q0 + 28 * t0 - 3 * math.pow(Q0, 2) + 8 * e1sq + 24 * math.pow(t0, 2)) * math.pow(dd0, 5) / 120
_a2 = (lof1 - lof2 + lof3) / math.cos(phi1)
_a3 = _a2 * 180 / math.pi
latitude = 180 * (phi1 - fact1 * (fact2 + fact3 + fact4)) / math.pi
if not northernHemisphere:
latitude = -latitude
longitude = ((zone > 0) and (6 * zone - 183.0) or 3.0) - _a3
return (latitude, longitude)
Y aquí pensé que era algo simple como facilitar * x + zona * y o algo así.
Otros consejos
Lo que encontré es el siguiente sitio: http: // home. hiwaay.net/~taylorc/toolbox/geography/geoutm.html Tiene un convertidor de JavaScript, debe verificar el algoritmo allí. Desde la página:
Programadores: el código fuente de JavaScript en este documento puede copiarse y reutilizarse sin restricciones.
Según esta página, UTM es compatible con proj4js.
http://trac.osgeo.org/proj4js/wiki/UserGuide#Supportedprojectionclasses
También puede echar un vistazo a GDAL . La biblioteca gdal tiene un excelente soporte para python, aunque puede ser un poco exagerado si solo está haciendo conversión de proyección.
También soy nuevo en esto y he estado estudiando el tema recientemente.
Aquí hay un método que encontré usando el python gdal pacakge (el osr está incluido en gdal). El paquete gdal es bastante poderoso, pero la documentación podría ser mejor.
Esto se deriva de una discusión aquí: 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
Y aquí está el método para convertir de un lat, lon en wgs84 (lo que informa la mayoría de las unidades gps) a 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
También descubrí que si ya tienes instalado django / gdal y conoces el EPSG para la zona UTM en la que está trabajando, puede usar el código Point () método 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"])
Puede usar Proj4js, de la siguiente manera.
Descargue Proj4JS desde GitHub, utilizando este enlace.
El siguiente código convertirá de UTM a latitud de longitud
<html>
<head>
<script src="proj4.js"></script>
<script>
var utm = "+proj=utm +zone=32";
var wgs84 = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs";
console.log(proj4(utm,wgs84,[539884, 4942158]));
</script>
</head>
<body>
</body>
</html>
En este código, la zona UTM es 32, como debería ser obvio. El Este es 539884 y el Norte es 4942158. El resultado es:
[9.502832656648073, 44.631671014204365]
Es decir 44.631671014204365N, 9.502832656648073E. Que tengo verificado es correcto.
Si necesita otras proyecciones, puede encontrar sus cadenas aquí .
Una versión Javascript de la respuesta de 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;
}
////////////////////////////////////////////////////////////////////////////////////////////
//
// ToLL - function to compute Latitude and Longitude given UTM Northing and Easting in meters
//
// Description:
// This member function converts input north and east coordinates
// to the corresponding Northing and Easting values relative to the defined
// UTM zone. Refer to the reference in this file's header.
//
// Parameters:
// north - (i) Northing (meters)
// east - (i) Easting (meters)
// utmZone - (i) UTM Zone of the North and East parameters
// lat - (o) Latitude in degrees
// lon - (o) Longitude in degrees
//
function ToLL(north,east,utmZone)
{
// This is the lambda knot value in the reference
var LngOrigin = DegToRad(utmZone * 6 - 183)
// The following set of class constants define characteristics of the
// ellipsoid, as defined my the WGS84 datum. These values need to be
// changed if a different dataum is used.
var FalseNorth = 0. // South or North?
//if (lat < 0.) FalseNorth = 10000000. // South or North?
//else FalseNorth = 0.
var Ecc = 0.081819190842622 // Eccentricity
var EccSq = Ecc * Ecc
var Ecc2Sq = EccSq / (1. - EccSq)
var Ecc2 = Math.sqrt(Ecc2Sq) // Secondary eccentricity
var E1 = ( 1 - Math.sqrt(1-EccSq) ) / ( 1 + Math.sqrt(1-EccSq) )
var E12 = E1 * E1
var E13 = E12 * E1
var E14 = E13 * E1
var SemiMajor = 6378137.0 // Ellipsoidal semi-major axis (Meters)
var FalseEast = 500000.0 // UTM East bias (Meters)
var ScaleFactor = 0.9996 // Scale at natural origin
// Calculate the Cassini projection parameters
var M1 = (north - FalseNorth) / ScaleFactor
var Mu1 = M1 / ( SemiMajor * (1 - EccSq/4.0 - 3.0*EccSq*EccSq/64.0 -
5.0*EccSq*EccSq*EccSq/256.0) )
var Phi1 = Mu1 + (3.0*E1/2.0 - 27.0*E13/32.0) * Math.sin(2.0*Mu1)
+ (21.0*E12/16.0 - 55.0*E14/32.0) * Math.sin(4.0*Mu1)
+ (151.0*E13/96.0) * Math.sin(6.0*Mu1)
+ (1097.0*E14/512.0) * Math.sin(8.0*Mu1)
var sin2phi1 = Math.sin(Phi1) * Math.sin(Phi1)
var Rho1 = (SemiMajor * (1.0-EccSq) ) / Math.pow(1.0-EccSq*sin2phi1,1.5)
var Nu1 = SemiMajor / Math.sqrt(1.0-EccSq*sin2phi1)
// Compute parameters as defined in the POSC specification. T, C and D
var T1 = Math.tan(Phi1) * Math.tan(Phi1)
var T12 = T1 * T1
var C1 = Ecc2Sq * Math.cos(Phi1) * Math.cos(Phi1)
var C12 = C1 * C1
var D = (east - FalseEast) / (ScaleFactor * Nu1)
var D2 = D * D
var D3 = D2 * D
var D4 = D3 * D
var D5 = D4 * D
var D6 = D5 * D
// Compute the Latitude and Longitude and convert to degrees
var lat = Phi1 - Nu1*Math.tan(Phi1)/Rho1 *
( D2/2.0 - (5.0 + 3.0*T1 + 10.0*C1 - 4.0*C12 - 9.0*Ecc2Sq)*D4/24.0
+ (61.0 + 90.0*T1 + 298.0*C1 + 45.0*T12 - 252.0*Ecc2Sq - 3.0*C12)*D6/720.0 )
lat = RadToDeg(lat)
var lon = LngOrigin +
( D - (1.0 + 2.0*T1 + C1)*D3/6.0
+ (5.0 - 2.0*C1 + 28.0*T1 - 3.0*C12 + 8.0*Ecc2Sq + 24.0*T12)*D5/120.0) / Math.cos(Phi1)
lon = RadToDeg(lon)
// Create a object to store the calculated Latitude and Longitude values
var sendLatLon = new PC_LatLon(lat,lon)
// Returns a PC_LatLon object
return sendLatLon
}
////////////////////////////////////////////////////////////////////////////////////////////
//
// RadToDeg - function that inputs a value in radians and returns a value in degrees
//
function RadToDeg(value)
{
return ( value * 180.0 / Math.PI )
}
////////////////////////////////////////////////////////////////////////////////////////////
//
// PC_LatLon - this psuedo class is used to store lat/lon values computed by the ToLL
// function.
//
function PC_LatLon(inLat,inLon)
{
this.lat = inLat // Store Latitude in decimal degrees
this.lon = inLon // Store Longitude in decimal degrees
}
Hay un módulo perl a través de CPAN llamado Geography :: NationalGrid que puede convertir este / norte a lat / largo. Eso puede ayudar.
Alternativamente, hay muchos scripts en el sitio de tipo móvil que te permiten convertir lat / long y easting / northings.
Un problema que tuve con el uso de proj4js fue que necesitaba la zona exacta como señala @Richard. Encontré un gran recurso aquí que puede convertir WGS a UTM y escribí un contenedor más limpio en JavaScript:
La respuesta de Staale me funcionó con una pequeña modificación: El módulo matemático no puede manejar la serie Pandas, por lo que reemplacé todas las funciones matemáticas con numpy.
Sin embargo, al verificar QGIS, veo una diferencia de aproximadamente 4 m entre las coordenadas UTM y LAT / LON.
Código a continuación:
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)
De esa manera puedo hacer esto directamente:
df['LAT'], df['LON']=utmToLatLng(31, df['X'], df['Y'], northernHemisphere=True)
