座標変換モジュール
地球は曲面ですが、地球上の位置を示す際、狭い範囲の測量であれば誤差も少ないことから、「平面」として計算したほうが便利です。
平面直角座標系とは、平面として計算を行えるように定められた座標系であり、比較的狭い範囲を扱う場合に適しています。
平面直角座標系は、日本を19のゾーンに分割して横メルカトル図法で投影し、各ゾーンに座標原点を設けて、
その原点を通る子午線をX軸、これに直交する方向をY軸としたものです。
各座標系（1系〜19系）の原点の値は、それぞれX=Y=0メートルとなります。
直交座標系と緯度経度の相互変換式は複雑難解なものとなっています。
また、日本測地系と世界測地系では、楕円体の長半径及び逆扁平率に差があり、計算結果が異なります。
下記は平面直角座標系と緯度経度を相互変換するプロシージャ(平面直角座標変換.bas)の一部です。
これをインポートすれば簡単に直交座標の変換が可能となります。
どうしても計算式の元データを知りたい方は平面直角座標の計算方法を参考ください。
平面直角座標x,yから緯度求めるプロシージャにφ(系図番号、Ｘ値、Ｙ値、測地系モード(0:日本測地系、1:世界測地系))入力
経度求めるプロシージャも同様にλ(系図番号、Ｘ値、Ｙ値、測地系モード(0:日本測地系、1:世界測地系))入力
'平面直角座標x,yから緯度を求める計算
Public Function φ(Keizu As Integer, X As Double, Y As Double, Mode As Integer) As Double
Dim φ1 As Double, N1 As Double, n12 As Double, t1 As Double
NewMode = Mode
φ1 = φn1(Keizu, X)
n12 = (e2 / (1 - e2)) * Cos(φ1) ^ 2
N1 = a / Sqr(1 - e2 * Sin(φ1) ^ 2)
t1 = Tan(φ1)
φ = φ1 - t1 / 2 / (N1 ^ 2) * (1 + n12) * (Y / m0) ^ 2 + _
t1 / 24 / (N1 ^ 4) * (5 + 3 * t1 ^ 2 + 6 * n12 - 6 * t1 ^ 2 * n12 - 3 * n12 ^ 2 - 9 * t1 ^ 2 * n12 ^ 2) * (Y / m0) ^ 4 _
- t1 / 720 / (N1 ^ 6) * (61 + 90 * t1 ^ 2 + 45 * t1 ^ 4 + 107 * n12 - 162 * t1 ^ 2 * n12 - 45 * t1 ^ 4 * n12) * (Y / m0) ^ 6 _
+ t1 / 40320 / (N1 ^ 8) * (1385 + 3633 * t1 ^ 2 + 4095 * t1 ^ 4 + 1575 * t1 ^ 6) * (Y / m0) ^ 8
φ = φ * 180 / PI
End Function
'平面直角座標x,yから緯度を求める計算
Public Function λ(Keizu As Integer, X As Double, Y As Double, Mode As Integer) As Double
Dim φ1 As Double, N1 As Double, n12 As Double, t1 As Double ', λ0 As Double
NewMode = Mode
φ1 = φn1(Keizu, X)
n12 = (e2 / (1 - e2)) * Cos(φ1) ^ 2
N1 = a / Sqr(1 - e2 * Sin(φ1) ^ 2)
t1 = Tan(φ1)
λ = GXo(Keizu) * PI / 180 + 1 / N1 / Cos(φ1) * (Y / m0) - 1 / 6 / N1 ^ 3 / Cos(φ1) * (1 + 2 * t1 ^ 2 + n12) * (Y / m0) ^ 3 _
+ 1 / 120 / N1 ^ 5 / Cos(φ1) * (5 + 28 * t1 ^ 2 + 24 * t1 ^ 4 + 6 * n12 + 8 * t1 ^ 2 * n12) * (Y / m0) ^ 5 _
- 1 / 5040 / N1 ^ 7 / Cos(φ1) * (61 + 662 * t1 ^ 2 + 1320 * t1 ^ 4 + 720 * t1 ^ 6 * n12) * (Y / m0) ^ 7
λ = λ * 180 / PI
End Function
このプロシージャだけを見ると簡単ですが、φn1()プロシージャにデータ処理で初期値がない場合は各変数を計算します。
緯度・緯度から平面直角座標x,yを求めるプロシージャは下記のようになります。
'緯度・経度から平面直角座標xを求める
Public Function X座標(Keizu As Integer, E As Double, ns As Double, Mode As Integer) As Double
Dim N As Double, T As Double, dE As Double, Er As Double, Nr As Double, n2 As Double
If a = 0 Then Call 初期値
dE = (E - GXo(Keizu)) * PI / 180
Er = E * PI / 180: Nr = ns * PI / 180
N = a / Sqr(1 - e2 * Sin(Nr) ^ 2)
n2 = (e2 / (1 - e2)) * Cos(Nr) ^ 2
T = Tan(Nr)
X座標 = ((S(Nr) - S(GYo(Keizu) * PI / 180)) + N / 2 * Cos(Nr) ^ 2 * T * dE ^ 2 _
+ N / 24 * Cos(Nr) ^ 4 * T * (5 - T ^ 2 + 9 * n2 + 4 * n2 ^ 2) * dE ^ 4 _
- N / 720 * Cos(Nr) ^ 6 * T * (-61 + 58 * T ^ 2 - T ^ 4 - 270 * n2 + 330 * T ^ 2 * n2) * dE ^ 6 _
- N / 40320 * Cos(Nr) ^ 8 * T * (-1385 + 3111 * T ^ 2 - 543 * T ^ 4 + T ^ 6) * dE ^ 8) * m0
End Function
'緯度・経度から平面直角座標Yを求める
Public Function Y座標(Keizu As Integer, E As Double, ns As Double, Mode As Integer) As Double
Dim N As Double, T As Double, dE As Double, Er As Double, Nr As Double, n2 As Double
NewMode = Mode
If a = 0 Or NewMode <> OldMode Then Call 初期値
dE = (E - GXo(Keizu)) * PI / 180
Er = E * PI / 180: Nr = ns * PI / 180
N = a / Sqr(1 - e2 * Sin(Nr) ^ 2)
n2 = (e2 / (1 - e2)) * Cos(Nr) ^ 2
T = Tan(Nr)
Y座標 = (N * Cos(Nr) * dE _
- N / 6 * Cos(Nr) ^ 3 * (-1 + T ^ 2 - n2) * dE ^ 3 _
- N / 120 * Cos(Nr) ^ 5 * (-5 + 18 * T ^ 2 - T ^ 4 + 58 * T ^ 2 * n2) * dE ^ 5 _
- N / 5040 * Cos(Nr) ^ 7 * (-61 + 479 * T ^ 2 - 179 * T ^ 4 + T ^ 6) * dE ^ 7) * m0
End Function
一般に同一系図Noで、同一測地系を連続して計算することが多いため、変数を維持し計算速度が上がるように作っています。
日本測地系と世界測地系の変換
2000年以前は日本測地系を使っていましたが2001年からは世界測地系が標準で使われるようになりました。
世界測地系はＧＰＳで用いられている緯度経度ですが、日本測地系は地球の中心位置と扁平率が違うため、精度の高い変換はかなり厄介になります。
下記のような簡易式で計算することもできますが、数ｍの誤差が発生します。
Public Function 世界座標簡易変換Y(Y As Double, X As Double) As Double
世界座標簡易変換Y = Y - 0.00010695 * Y + 0.000017464 * X + 0.0046017
End FunctionPublic Function 世界座標簡易変換X(Y As Double, X As Double) As Double
世界座標簡易変換X = X - 0.000046038 * Y - 0.000083043 * X + 0.01004
End FunctionPublic Function 日本測地簡易変換Y(Y As Double, X As Double) As Double
日本測地簡易変換Y = Y + Y * 0.00010696 - X * 0.000017467 - 0.004602
End FunctionPublic Function 日本測地簡易変換X(Y As Double, X As Double) As Double
日本測地簡易変換X = X + Y * 0.000046047 + X * 0.000083049 - 0.010041
End Function
楕円形地球体と原点の違いだけを正しく修正するだけならモロデンスキー法がありこれで計算すればいいのですが、2000年以前の測量には計測誤差がかなりあり、計測誤差も含めて1kmメッシュ毎に修正値を測量結果が測地成果2000として公開されています。
このデータをダウンロードし、北海道付近の値をバイナリファイル化(TKJ2_600.dat)し高速で参照できるようにしたのがTKJ2000.basです。インポートしたExcelファイルと同一フォルダにTKJ2_600.datを保存する。
また、測地成果2000は地上部分だけのデータしかないため、海上部分はモロデンスキー法により求めています。
変換プロシージャは下記になり、Ｘ：緯度、Ｙ：経度を入力すると変換結果が得られます。
実際の処理はかなり多いのですが地上部分であればかなり高速で処理します。
Public Function 世界座標変換Y(Y As Double, X As Double) As Double Dim D As Single
D = TYJ2_Dxy(X, Y, 0)
If D = 0 Then '海の上
世界座標変換Y = Molodensky(Y, X, 0, 0, 0)
Else
世界座標変換Y = D + Y
End IfEnd Function
Public Function 世界座標変換X(Y As Double, X As Double) As Double Dim D As Single
D = TYJ2_Dxy(X, Y, 1)
If D = 0 Then '海の上
世界座標変換X = Molodensky(Y, X, 0, 1, 0)
Else
世界座標変換X = D + X
End IfEnd Function
Public Function 日本測地変換Y(Y As Double, X As Double) As Double 日本測地変換Y = Y - TYJ2_Dxy(X, Y, 0)
If 日本測地変換Y = Y Then '海の上
日本測地変換Y = Molodensky(Y, X, 0, 0, 1)
End IfEnd Function
Public Function 日本測地変換X(Y As Double, X As Double) As Double 日本測地変換X = X - TYJ2_Dxy(X, Y, 1)
If 日本測地変換X = X Then '海の上
日本測地変換X = Molodensky(Y, X, 0, 1, 1)
End IfEnd Function
主要計算処理部分
'新図葉からTYJ_Dx
Public Function TYJ_Dx(NZUYo As Long, ByRef dy As Single) As SingleDim I As Long, Ip As Long, Ig As Long, Id As Long
If TYJ_N = 0 Then 'まだデータを読み込んでいない
Open ActiveWorkbook.Path & "\TKJ2_600.dat" For Binary As #56
Get #56, , TYJ_N '
ReDim TYJ(TYJ_N), TYJ_xy(TYJ_N, 1) ', FedNo(HiDtn)
Get #56, , TYJ '
Get #56, , TYJ_xy
Close #56
End If
Id = -1
Ip = Int(TYJ_N / 2)
For I = 0 To 16
If TYJ(Ig + Ip) < NZUYo Then
Ig = Ig + Ip
ElseIf TYJ(Ig + Ip) = NZUYo Then
Id = Ig + Ip: Exit For
End If
Ip = Int(Ip / 2)
If Ip + Ig > TYJ_N Then Exit For
Next
If Id = -1 Then
For I = 1 To 16
If Ig + I > TYJ_N Then Exit For
If TYJ(Ig + I) = NZUYo Then Id = Ig + I: Exit For
Next
End If
If Id = -1 Then
TYJ_Dx = 0: dy = 0: Exit Function
End If
TYJ_Dx = TYJ_xy(Id, 0): dy = TYJ_xy(Id, 1)End Function
'座標からTYJ2
Public Function TYJ2_Dxy(X As Double, Y As Double, SY As Long) As SingleDim dy As Single
TYJ2_Dxy = TYJ_Dx(NZUYo(X, Y), dy)
If SY = 1 Then TYJ2_Dxy = dyEnd Function
おまけ jpgファイルにあるGPS情報の取得
JPGファイルの規格はこちら見てもよくわかりません。
JPGファイルの記述にはBig EndianとLittle Endianがあり、かなり複雑ですが下記のプロシージャにディレクトリーとファイル名を入れると座標が得られます。
'JPGファイルのGPS情報取り出し SW=0:緯度 SW=1:経度 を返す
Public Function JPGGPS(FName As String, SW As Integer) As Variant
省略
End Function
'経度情報はPublic変数で保存されているため緯度計算後であればこれで取得もできる
Public Function GPS_Ny() As Double
GPS_Ny = Gny
End Function
JPGPS.basをインポートするだけです。