UTM形式の座標を持つファイルがたくさんあります。座標ごとに、東、北、ゾーンがあります。これをLatLngに変換して、Google Map APIで使用し、地図に情報を表示する必要があります。
これを行うオンライン計算機をいくつか見つけましたが、実際のコードやライブラリはありません。http://trac.osgeo.org/proj4js/はJavascriptのプロジェクションライブラリですが、デモを見るとUTMプロジェクションは含まれていません。
私はまだGISドメイン全体にかなり新鮮なので、私が欲しいのは何かalaです：
(lat,lng) = transform(easting, northing, zone)
私はそれを解決したIBMからのJavaコードを見つけました：http://www.ibm.com/developerworks/java/library/j-coordconvert/index.html
参考までに、必要なメソッドの Python 実装を次に示します。
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)
そして、ここで私はそれが何かのような単純なものだと思ったeasting*x+zone*y.
私が見つけたのは次のサイトです。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パッケージ( 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
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 が既にインストールされていて、作業中の UTM ゾーンのEPSGコードがわかっている場合は、 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 バージョン
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
}
Geography::NationalGridと呼ばれるCPAN経由のperlモジュールがあります。これはeasting/northingをlat/longsに変換できます。それが役立つかもしれません。
あるいは、活字サイトには、緯度/経度と東/北を変換できるスクリプトがたくさんあります。