Wie von UTM zu LatLng in Python oder Javascript konvertieren
https://stackoverflow.com/questions/343865
-
19-08-2019 - |
italiano
english
français
española
中国
日本の
العربية
Deutsch
한국어
Português
RussianFrage
Ich habe eine Reihe von Dateien mit Koordinaten in UTM-Form. jede Koordinate Ich habe nämlich Ost, Nord und Zone. Ich brauche dies für die Verwendung mit Google Map API LatLng zu konvertieren, die Informationen in einer Karte zu zeigen.
Ich habe einig Online-Rechner gefunden, der dies tut, aber keinen eigentlichen Code oder Bibliotheken. http://trac.osgeo.org/proj4js/ ist eine Projektion Bibliothek für Javascript, aber auf der Suche an der Demo spielt es keine UTM-Projektion ist.
Ich bin immer noch ziemlich frisch auf die gesamte GIS-Domäne so, was ich will, ist etwas ala:
(lat,lng) = transform(easting, northing, zone)
Lösung
Ich landete Java-Code von IBM zu finden, die es zu lösen: http://www.ibm.com/developerworks/java/library/j-coordconvert/index.html
Nur als Referenz, hier ist mein Python-Implementierung der Methode, die ich benötigen:
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)
Und ich dachte, es war etwas Einfaches wie easting*x+zone*y oder so etwas.
Andere Tipps
Was ich fand, ist die folgende Website: http: // nach Hause. hiwaay.net/~taylorc/toolbox/geography/geoutm.html Es hat einen Javascript-Konverter, sollten Sie den Algorithmus dort überprüfen. Von der Seite:
Programmierer. Der JavaScript-Quellcode in diesem Dokument kopiert und ohne Einschränkung wieder verwendet werden können
Nach dieser Seite UTM wird von proj4js unterstützt.
http://trac.osgeo.org/proj4js/wiki/UserGuide#Supportedprojectionclasses
Sie können auch einen Blick auf GDAL zu nehmen. Die gdal Bibliothek verfügt über ausgezeichnete Python Unterstützung, obwohl es ein bisschen übertrieben sein, wenn Sie nur Projektion Konvertierung.
Ich bin neu in dieser so gut und wurde vor kurzem zum Thema Studium auf.
Hier ist eine Methode, die ich mit dem Python gdal pacakge (die osr gefunden Paket wird in gdal enthalten). Das gdal Paket ist ziemlich mächtig, aber die Dokumentation könnte besser sein.
Dies ist aus einer Diskussion hier abgeleitet: 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
Und hier ist das Verfahren von einem lat zum Umwandeln lon in WGS84 (was die meisten GPS-Geräte Bericht) 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
Ich fand auch, dass, wenn Sie bereits django / gdal installiert habe und Sie kennen die EPSG Code für die UTM-Zone, die Sie gerade arbeiten, können Sie einfach die Point() verwenden Transformation () Methode.
from django.contrib.gis.geos import Point
utm2epsg = {"54N": 3185, ...}
p = Point(lon, lat, srid=4326) # 4326 = WGS84 epsg code
p.transform(utm2epsg["54N"])
Sie könnten Proj4js wie folgt verwenden.
Download Proj4JS von GitHub, mit diesem Link.
Der folgende Code wird von UTM Länge Breite konvertieren
<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>
In diesem Code ist die UTM-Zone 32, wie offensichtlich sein sollte. Der Easting ist 539.884 und das Northing ist 4942158. Das Ergebnis ist:
[9.502832656648073, 44.631671014204365]
Was ist 44.631671014204365N zu sagen, 9.502832656648073E. Was ich habe verifiziert korrekt ist.
Wenn Sie andere Projektionen benötigen, können Sie ihre Saiten finden hier .
Eine Javascript-Version Staale Antwort
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
}
Es ist ein Perl-Modul über CPAN Geographie :: NationalGrid genannt, die easting / northing umwandeln / longs lat. Das kann helfen.
Alternativ gibt es viele Skripte auf dem beweglichen Lettern Website , dass Sie bei lat / long und easting / Koordinaten X konvertieren.
Ein Problem, das ich mit der Verwendung von proj4js hatte war, dass es die genaue Zone benötigt als @ Richard weist darauf hin. Ich fand eine große Ressource hier die WGS zu UTM umwandeln und schrieb eine sauberere Wrapper in JavaScript:
Die Antwort von Staale arbeitete für mich mit einer kleinen Modifikation - Das Mathematikmodul nicht verarbeiten kann Pandas Serie, so ersetzte ich alle mathematischen Funktionen mit numpy.
Doch in QGIS Überprüfung, sehe ich etwa 4m Differenz zwischen den UTM und LAT / LON-Koordinaten.
-Code unten:
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)
So kann es dies direkt tun können:
df['LAT'], df['LON']=utmToLatLng(31, df['X'], df['Y'], northernHemisphere=True)
