2011-10-24 16 views
12

私はeasting/northing形式の位置座標を持っていますが、それを適切なlat値に変換して、それをbingマップに配置する必要があります。任意の数式または詳細をeast/northingをlat/lonに変換する方法は?東経度北緯緯度北緯

EDIT:具体的には、私はSVY21が

答えて

13

Eastings北は基点の東と北のそれぞれの距離です。基点は通常緯度と経度であり、eastingsとnorthingsは通常メートルまたはフィートで表されます。しかし、東と北は、それらを正とし、それらが基点の西と南を表現できるようにするために、通常、特定の価値を相殺しています。

一般に、ある座標系から別の座標系への変換は、異なる楕円(地球モデル)とデータムを持つ可能性があるため、単純ではありません。私が理解するように、ある座標系から別の座標系に変換するための式はかなり複雑です。

SVY21は、WGS84とまったく同じデータムと楕円体を使用するため、タスクが簡単になります。 SVY21では、eastingsとnorthingsの基点はBase 7 at Pierce Reservoir,1 degです。 22分。 02.9154秒。北と103度。 49分31.9752秒。 (すなわち、約1.3674765度の緯度と約103.8255487度の経度ですが、よく知られているテキストではそれぞれ1.3666度と103.8333度を使用しています)。東方面のオフセットは28001.642メートル、北向きのオフセットは38744.572メートルです。 EPSGコードは3414です。私はあなたのeastingsとnorthingsがメートルで表されると仮定します。

SVY21はWGS84と同じシステムを使用しているので、あなたがしなければならないすべては、次のとおりです。

  • それぞれのオフセット値によって東距と北距を引きます。
  • 基準点、eastingの絶対値、eastingが正の場合は90度の方位、または270度の方位を指定すると、指定された点の経度を見つけることができます度数が負の場合This linkには関連する式が含まれています。 (この計算では、「開始点からの距離と方位を指定する」の節に示すように、余弦の球法を使用するか、より正確なVincenty's direct formulaを使用することができます。この計算のHaversine公式)
  • 基準点、北向きの絶対値、および北向きが正の場合は0度の方位、または180度の場合の方位を指定して、指定された点の緯度を求めますそれは否定的です。
+0

あなたはHaversineの公式を参照していますか?どのようにeasting/northingからlat/lonへ? – Bahamut

+0

いいえ、私は余弦の球の法則を指していますが、それはVincentyの直接の公式でも働くかもしれません。 –

+0

ありがとう、私はそれを試してみましょう。 – Bahamut

3

異なる座標系の何百もありますWGS84する座標に変換する必要があります - 東距/北距と緯度が/ロング「種類の座標のありますが、それらそれらの座標が得られるシステムを一意的に特定するのに十分ではない。

両方のソースについて、EPSGコード(例:4326,4269,2700,32701)または空間参照系の詳細(データム、投影、素子午線および単位)が必要です選択された宛先フォーマット。 質問タイトルに「GPS」と記載されていますので、必要な緯度/経度はグローバル測位システムで使用されるWGS84データに関連して定義されていると仮定していますが、東西の価値。

あなたが使用される投影の詳細を持ったら、あなたはPROJ.4ライブラリ(http://trac.osgeo.org/proj/)のようなものを使用してコードに変換を行うことができます

+0

私が得た結果はSVY21形式であり、WGS84に変換してビンマップにプロットする必要があります。それはライブラリを使用して可能ですか? – Bahamut

2

perlの中で比較的簡単な解決策があります:

は、すべての最初の、だから、あなたはPerlがインストールされていることを確認してください。あなたはいくつかの方法でこれを行うことができます

ジオ:: HelmertTransform地理:: NationalGrid CAM :: DBF mySociety :: GeoUtil

:次に、 次の4つのモジュールをインストールします。ここで私はそれをやった方法は次のとおりです。

# Geo::HelmertTransform 
wget http://search.cpan.org/CPAN/authors/id/M/MY/MYSOCIETY/Geo-HelmertTransform-1.13.tar.gz 
tar xzf Geo-HelmertTransform-1.13.tar.gz 
perl Makefile.PL 
make 
make install 

# Geography::NationalGrid 
http://search.cpan.org/CPAN/authors/id/P/PK/PKENT/Geography-NationalGrid-1.6.tar.gz 
tar xzf Geography-NationalGrid-1.6.tar.gz 
perl Makefile.PL 
make 
make install 

# CAM::DBF 
wget http://search.cpan.org/CPAN/authors/id/C/CL/CLOTHO/CAM-DBF-1.02.tgz 
tar xzf CAM-DBF-1.02.tgz 
perl Makefile.PL 
make 
make install 

# mySociety::GeoUtil 
# See: http://parlvid.mysociety.org:81/os/ -> https://github.com/mysociety/commonlib/blob/master/perllib/mySociety/GeoUtil.pm 
mkdir -p mySociety 
wget -O mySociety/GeoUtil.pm 'https://raw.githubusercontent.com/mysociety/commonlib/master/perllib/mySociety/GeoUtil.pm' 
  1. はGBのデータを取得します。

をクリックして指示に従うと、英国の「Code-Point®Open」データセットをダウンロードしてください。あなたは次のようにあなたがそれを抽出することができ codepo_gb.zipをダウンロードしたら:

解凍codepo_gb.zip

を解凍したファイルは、あなたが順番に、次のPerlスクリプトを実行することができます 、カレントディレクトリになっていることを仮定しますデータを解析するには、 GB eastings/northingsを抽出し、 緯度/経度に変換します。

use strict; 
use mySociety::GeoUtil qw/national_grid_to_wgs84/; 

while (<>) { 
    my @x=split(/,/); # split csv 
    my ($pc, $east, $north) = ($x[0], $x[10], $x[11]); 
    $pc=~s/\"//g; # remove quotes around postcode 
    my ($lat, $lng) = national_grid_to_wgs84($east, $north, "G"); # "G" means Great Britain 
    print "$pc,$lat,$lng\n"; 
} 

([0]、$ X [10]と$ X $ X、また覚えて...、呼び出しの.plファイルへの最後のコードブロックを保存し、perl script.pl your.csvを呼び出すには、[11 ]は、それぞれ東距と北距、郵便番号の列番号でなければなりません。

に完全な信用http://baroque.posterous.com/uk-postcode-latitudelongitude

+0

は良いでしょうが、私は移植性のために変換の計算や公式を好む。第2に、svy21はシンガポールのため、英国のSVYからWGS84への変換は適用されません。 – Bahamut

+0

私は悪いですが、うまくいけばそれは誰かを助けることができますが、私は解決策を見つけるために急いでこの質問に遭遇し、これは英国の変換のための最も簡単な方法でした。 – rickyduck

+1

だから誰かがこれも便利だと思うかもしれないので、私はそれをdownvoteしなかったのです。 – Bahamut

2

に私は/経度の値を緯度するWGS84のためのT-SQL関数にJavaScript実装を変換しました。気軽異なる座標系が必要な場合は、University of Wisconsinをチェックしてください - 私がソースとして使用したGreen BayのWebページで、更新された定数を取得します。

drop function UF_utm_to_lat 
go 
create function UF_utm_to_lat(@utmz float, @x float, @y float) returns float 
as 
begin 
    --Based on code from this page: http://www.uwgb.edu/dutchs/usefuldata/ConvertUTMNoOZ.HTM 
    declare @latitude float; 
    declare @longitude float; 
    set @latitude = 0.00; 
    set @longitude = 0.00; 

    --Declarations 
    declare @a float; 
    declare @f float; 
    declare @drad float; 
    declare @k0 float; 
    declare @b float; 
    declare @e float; 
    declare @e0 float; 
    declare @esq float; 
    declare @e0sq float; 
    declare @zcm float; 
    declare @e1 float; 
    declare @M float; 
    declare @mu float; 
    declare @phi1 float; 
    declare @C1 float; 
    declare @T1 float; 
    declare @N1 float; 
    declare @R1 float; 
    declare @D float; 
    declare @phi float; 
    declare @lng float; 
    declare @lngd float; 

    --Datum Info here: Name, a, b, f, 1/f 
    --WGS 84 6,378,137.0 6356752.314 0.003352811 298.2572236 

    set @a = 6378137.0; 
    set @b = 6356752.314; 
    set @f = 0.003352811; 
    set @drad = PI()/180.0; 
    set @k0 = 0.9996; --scale on central meridian 

    set @e = SQRT(1.0 - (@b/@a)*(@b/@a)); --Eccentricity 
    --e = Math.sqrt(1 - (b/a)*(b/a));//eccentricity 
    set @e0 = @e/SQRT(1.0 - @e*@e); --Called e prime in reference 
    --e0 = e/Math.sqrt(1 - e*e);//Called e prime in reference 
    set @esq = (1.0 - (@b/@a)*(@b/@a)); --e squared for use in expansions 
    --esq = (1 - (b/a)*(b/a));//e squared for use in expansions 
    set @e0sq = @e*@e/([email protected]*@e); --e0 squared - always even powers 
    --e0sq = e*e/(1-e*e);// e0 squared - always even powers 
    set @zcm = 3.0 + 6.0*(@utmz-1.0) - 180.0; --Central meridian of zone 
    --zcm = 3 + 6*(utmz-1) - 180;//Central meridian of zone 
    set @e1 = (1.0 - SQRT(1.0 - @e*@e))/(1.0 + SQRT(1.0 - @e*@e)); --Called e1 in USGS PP 1395 also 
    --e1 = (1 - Math.sqrt(1 - e*e))/(1 + Math.sqrt(1 - e*e));//Called e1 in USGS PP 1395 also 
    set @M = 0.0 + @y/@k0; --Arc length along standard meridian 
    --M = M0 + y/k0;//Arc length along standard meridian. 
    set @mu = @M/(@a*(1.0 - @esq*(1.0/4.0 + @esq*(3.0/64.0 + 5.0*@esq/256.0)))); 
    --mu = M/(a*(1 - esq*(1/4 + esq*(3/64 + 5*esq/256)))); 
    set @phi1 = @mu + @e1*(3.0/2.0 - 27.0*@e1*@e1/32.0)*SIN(2.0*@mu) + @e1*@e1*(21.0/16.0 - 55.0*@e1*@e1/32.0)*SIN(4.0*@mu); --Footprint Latitude 
    --phi1 = mu + e1*(3/2 - 27*e1*e1/32)*Math.sin(2*mu) + e1*e1*(21/16 -55*e1*e1/32)*Math.sin(4*mu);//Footprint Latitude 
    set @phi1 = @phi1 + @e1*@e1*@e1*(SIN(6.0*@mu)*151.0/96.0 + @e1*SIN(8.0*@mu)*1097.0/512.0); 
    --phi1 = phi1 + e1*e1*e1*(Math.sin(6*mu)*151/96 + e1*Math.sin(8*mu)*1097/512); 
    set @C1 = @e0sq*POWER(COS(@phi1),2.0); 
    --C1 = e0sq*Math.pow(Math.cos(phi1),2); 
    set @T1 = POWER(TAN(@phi1),2.0); 
    --T1 = Math.pow(Math.tan(phi1),2); 
    set @N1 = @a/SQRT(1.0-POWER(@e*SIN(@phi1),2.0)); 
    --N1 = a/Math.sqrt(1-Math.pow(e*Math.sin(phi1),2)); 
    set @R1 = @N1*([email protected]*@e)/(1.0-POWER(@e*SIN(@phi1),2.0)); 
    --R1 = N1*(1-e*e)/(1-Math.pow(e*Math.sin(phi1),2)); 
    set @D = (@x-500000.0)/(@N1*@k0); 
    --D = (x-500000)/(N1*k0); 
    set @phi = (@D*@D)*(1.0/2.0 - @D*@D*(5.0 + 3.0*@T1 + 10.0*@C1 - 4.0*@C1*@C1 - 9.0*@e0sq)/24.0); 
    --phi = (D*D)*(1/2 - D*D*(5 + 3*T1 + 10*C1 - 4*C1*C1 - 9*e0sq)/24); 
    set @phi = @phi + POWER(@D,6.0)*(61.0 + 90.0*@T1 + 298.0*@C1 + 45.0*@T1*@T1 - 252.0*@e0sq - 3.0*@C1*@C1)/720.0; 
    --phi = phi + Math.pow(D,6)*(61 + 90*T1 + 298*C1 + 45*T1*T1 -252*e0sq - 3*C1*C1)/720; 
    set @phi = @phi1 - (@N1*TAN(@phi1)/@R1)*@phi; 
    --phi = phi1 - (N1*Math.tan(phi1)/R1)*phi; 


    set @latitude = FLOOR(1000000.0*@phi/@drad)/1000000.0; 

    set @lng = @D*(1.0 + @D*@D*((-1.0 - 2.0*@T1 - @C1)/6.0 + @D*@D*(5.0 - 2.0*@C1 + 28.0*@T1 - 3.0*@C1*@C1 + 8.0*@e0sq + 24.0*@T1*@T1)/120))/COS(@phi1); 
    set @lngd = @[email protected]/@drad; 
    set @longitude = FLOOR(1000000.0*@lngd)/1000000.0; 


    return @latitude; 
end 
go 
drop function UF_utm_to_long 
go 
create function UF_utm_to_long(@utmz float, @x float, @y float) returns float 
as 
begin 
    --Based on code from this page: http://www.uwgb.edu/dutchs/usefuldata/ConvertUTMNoOZ.HTM 
    declare @latitude float; 
    declare @longitude float; 
    set @latitude = 0.00; 
    set @longitude = 0.00; 

    --Declarations 
    declare @a float; 
    declare @f float; 
    declare @drad float; 
    declare @k0 float; 
    declare @b float; 
    declare @e float; 
    declare @e0 float; 
    declare @esq float; 
    declare @e0sq float; 
    declare @zcm float; 
    declare @e1 float; 
    declare @M float; 
    declare @mu float; 
    declare @phi1 float; 
    declare @C1 float; 
    declare @T1 float; 
    declare @N1 float; 
    declare @R1 float; 
    declare @D float; 
    declare @phi float; 
    declare @lng float; 
    declare @lngd float; 

    --Datum Info here: Name, a, b, f, 1/f 
    --WGS 84 6,378,137.0 6356752.314 0.003352811 298.2572236 

    set @a = 6378137.0; 
    set @b = 6356752.314; 
    set @f = 0.003352811; 
    set @drad = PI()/180.0; 
    set @k0 = 0.9996; --scale on central meridian 

    set @e = SQRT(1.0 - (@b/@a)*(@b/@a)); --Eccentricity 
    --e = Math.sqrt(1 - (b/a)*(b/a));//eccentricity 
    set @e0 = @e/SQRT(1.0 - @e*@e); --Called e prime in reference 
    --e0 = e/Math.sqrt(1 - e*e);//Called e prime in reference 
    set @esq = (1.0 - (@b/@a)*(@b/@a)); --e squared for use in expansions 
    --esq = (1 - (b/a)*(b/a));//e squared for use in expansions 
    set @e0sq = @e*@e/([email protected]*@e); --e0 squared - always even powers 
    --e0sq = e*e/(1-e*e);// e0 squared - always even powers 
    set @zcm = 3.0 + 6.0*(@utmz-1.0) - 180.0; --Central meridian of zone 
    --zcm = 3 + 6*(utmz-1) - 180;//Central meridian of zone 
    set @e1 = (1.0 - SQRT(1.0 - @e*@e))/(1.0 + SQRT(1.0 - @e*@e)); --Called e1 in USGS PP 1395 also 
    --e1 = (1 - Math.sqrt(1 - e*e))/(1 + Math.sqrt(1 - e*e));//Called e1 in USGS PP 1395 also 
    set @M = 0.0 + @y/@k0; --Arc length along standard meridian 
    --M = M0 + y/k0;//Arc length along standard meridian. 
    set @mu = @M/(@a*(1.0 - @esq*(1.0/4.0 + @esq*(3.0/64.0 + 5.0*@esq/256.0)))); 
    --mu = M/(a*(1 - esq*(1/4 + esq*(3/64 + 5*esq/256)))); 
    set @phi1 = @mu + @e1*(3.0/2.0 - 27.0*@e1*@e1/32.0)*SIN(2.0*@mu) + @e1*@e1*(21.0/16.0 - 55.0*@e1*@e1/32.0)*SIN(4.0*@mu); --Footprint Latitude 
    --phi1 = mu + e1*(3/2 - 27*e1*e1/32)*Math.sin(2*mu) + e1*e1*(21/16 -55*e1*e1/32)*Math.sin(4*mu);//Footprint Latitude 
    set @phi1 = @phi1 + @e1*@e1*@e1*(SIN(6.0*@mu)*151.0/96.0 + @e1*SIN(8.0*@mu)*1097.0/512.0); 
    --phi1 = phi1 + e1*e1*e1*(Math.sin(6*mu)*151/96 + e1*Math.sin(8*mu)*1097/512); 
    set @C1 = @e0sq*POWER(COS(@phi1),2.0); 
    --C1 = e0sq*Math.pow(Math.cos(phi1),2); 
    set @T1 = POWER(TAN(@phi1),2.0); 
    --T1 = Math.pow(Math.tan(phi1),2); 
    set @N1 = @a/SQRT(1.0-POWER(@e*SIN(@phi1),2.0)); 
    --N1 = a/Math.sqrt(1-Math.pow(e*Math.sin(phi1),2)); 
    set @R1 = @N1*([email protected]*@e)/(1.0-POWER(@e*SIN(@phi1),2.0)); 
    --R1 = N1*(1-e*e)/(1-Math.pow(e*Math.sin(phi1),2)); 
    set @D = (@x-500000.0)/(@N1*@k0); 
    --D = (x-500000)/(N1*k0); 
    set @phi = (@D*@D)*(1.0/2.0 - @D*@D*(5.0 + 3.0*@T1 + 10.0*@C1 - 4.0*@C1*@C1 - 9.0*@e0sq)/24.0); 
    --phi = (D*D)*(1/2 - D*D*(5 + 3*T1 + 10*C1 - 4*C1*C1 - 9*e0sq)/24); 
    set @phi = @phi + POWER(@D,6.0)*(61.0 + 90.0*@T1 + 298.0*@C1 + 45.0*@T1*@T1 - 252.0*@e0sq - 3.0*@C1*@C1)/720.0; 
    --phi = phi + Math.pow(D,6)*(61 + 90*T1 + 298*C1 + 45*T1*T1 -252*e0sq - 3*C1*C1)/720; 
    set @phi = @phi1 - (@N1*TAN(@phi1)/@R1)*@phi; 
    --phi = phi1 - (N1*Math.tan(phi1)/R1)*phi; 

    set @latitude = FLOOR(1000000.0*@phi/@drad)/1000000.0; 

    set @lng = @D*(1.0 + @D*@D*((-1.0 - 2.0*@T1 - @C1)/6.0 + @D*@D*(5.0 - 2.0*@C1 + 28.0*@T1 - 3.0*@C1*@C1 + 8.0*@e0sq + 24.0*@T1*@T1)/120))/COS(@phi1); 
    set @lngd = @[email protected]/@drad; 
    set @longitude = FLOOR(1000000.0*@lngd)/1000000.0; 


    return @longitude; 
end