日本測地系から世界測地系への座標変換

日本測地系の緯度・経度情報から世界測地系の緯度・経度情報に変換するAppleScriptです。

AppleScriptに35種類の関数計算機能を付加する「calcLibAS」の計算精度を確認するための実例でもあります。だいたい期待の精度は出せている感じです。

より精度を求められるような用途には、AppleScriptの数値変数自体が浮動小数点演算で10桁の計算精度しかないので、数値計算自体を外部にすべて出してしまうか、より精度の高い変数型を作ってしまうか（文字列で扱うようにすれば問題なさそう）というところでしょうか。

「日本測地系から世界測地系への座標変換（近似式）」のほうはObjective-Cから、「日本測地系から世界測地系への直行座標変換」のほうはRubyからAppleScriptに書き換えてみたものです。

Rubyのプログラムが割と単調だったので、読んでいて疲れてしまいました。変数名の付け方が途中で変わっているのは、「とりあえず動けばいい」ぐらいで組んでいたのと、変数名にいちいち「Local」とか付けるのに疲れたからです。

Piyomaru Script Assistantで変数名のみ一括置換してもいいわけですが、それはそれ、、、

AppleScript名：日本測地系から世界測地系への座標変換（近似式）

— Created 2018-07-13 by Takaaki Naganoya
— 2018 Piyomaru Software
use AppleScript version "2.4"
use scripting additions
use framework "Foundation"

set jpLat to 42.16593046887
set jpLng to 142.771877437246
set {wdLat, wdLng, wdH} to tky2wgsLinea(jpLat, jpLng, 0) of me –日本測地系から世界測地系に換算
–>　　{42.168515890674, 142.768119997121, 0}

set {latJ, lonJ, heightJ} to wgs2tkyLinea(wdLat, wdLng, wdH) of me –世界測地系から日本測地系に換算
–>　　{42.16593052955, 142.77187743437, 0}

–日本測地系から世界測地系に換算
–http://d.hatena.ne.jp/iRSS/20111112/1321113232
on tky2wgsLinea(latJ, lonJ, heightJ)
　　set heightW to heightJ
　　set lonW to lonJ – latJ * 4.6038E-5 – lonJ * 8.3043E-5 + 0.01004
　　set latW to latJ – latJ * 1.0695E-4 + lonJ * 1.7464E-5 + 0.0046017
　　return {latW, lonW, heightW}
end tky2wgsLinea

–世界測地系から日本測地系に換算
–https://altarf.net/computer/技術的なポエム/3332
on wgs2tkyLinea(latW, lonW, heightW)
　　set heightJ to heightW
　　set latJ to (latW * 1.000106961) – (lonW * 1.7467E-5) – 0.004602017
　　set lonJ to (lonW * 1.000083049) + (latW * 4.6047E-5) – 0.010041046
　　return {latJ, lonJ, heightJ}
end wgs2tkyLinea

★Click Here to Open This Script　

AppleScript名：日本測地系から世界測地系への直行座標変換

— Created 2018-07-13 by Takaaki Naganoya
— 2018 Piyomaru Software
use AppleScript version "2.4"
use scripting additions
use framework "Foundation"
use math : script "calcLibAS" –http://piyocast.com/as/archives/3523

set jpLat to 42.16593046887
set jpLng to 142.771877437246

set {wdLat, wdLng, wdH} to JapanGeodeticSystem’s convertIntoWgs(jpLat, jpLng)
–>　　{42.168504277889, 142.768075650137, 61.914452206343}

–https://altarf.net/computer/ruby/3347
–日本測地系から世界測地系に換算
script JapanGeodeticSystem
　　–ラジアン(度)
　　property rdNum : pi / 180
　　
　　–日本測地系の定数(ベッセル楕円体)
　　property rJP : 6.377397155E+6 –赤道半径
　　property fJP : 1 / 299.1528128 –扁平率
　　property e2JP : (2 * fJP) – (fJP * fJP) –第一離心率
　　
　　–世界測地系の定数(WGS84)
　　property rWS : 6.378137E+6 –赤道半径
　　property fWS : 1 / 298.257223563 –扁平率
　　property e2WS : (2 * fWS) – (fWS * fWS) –第一離心率
　　
　　–並行移動量(m)
　　property dxNum : -148
　　property dyNum : 507.0
　　property dzNum : 681.0
　　
　　–楕円体高
　　property heightNum : 0
　　
　　
　　–日本測地系から世界測地系に変換する
　　on convertIntoWgs(aLat, aLong)
　　　　set {x, y, z} to llh2xyz(aLat, aLong, heightNum, rJP, e2JP, rdNum) of me
　　　　
　　　　set x to x + dxNum
　　　　set y to y + dyNum
　　　　set z to z + dzNum
　　　　
　　　　set {lat1, lng1, h1} to xyz2llh(x, y, z, rWS, e2WS, rdNum) of me
　　　　return {lat1, lng1, h1}
　　end convertIntoWgs
　　
　　
　　–座標系の変換(緯度経度 -> xyz)
　　on llh2xyz(latLocal, lngLocal, hLocal, aLocal, e2Local, rdLocal)
　　　　set latLocal to latLocal * rdLocal
　　　　set lngLocal to lngLocal * rdLocal
　　　　
　　　　set sbNum to calcSin(latLocal) of math
　　　　set cbNum to calcCos(latLocal) of math
　　　　
　　　　set rnLocal to aLocal / (calcSqrt(1 – e2Local * sbNum * sbNum) of math)
　　　　
　　　　set xLocal to (rnLocal + hLocal) * cbNum * (calcCos(lngLocal) of math)
　　　　set yLocal to (rnLocal + hLocal) * cbNum * (calcSin(lngLocal) of math)
　　　　set zLocal to (rnLocal * (1 – e2Local) + hLocal) * sbNum
　　　　
　　　　return {xLocal, yLocal, zLocal}
　　end llh2xyz
　　
　　
　　–座標系の変換(xyz -> 緯度経度)
　　on xyz2llh(x, y, z, a, e2, rdLocal)
　　　　set bda to calcSqrt(1 – e2) of math
　　　　set p to calcSqrt(x * x + y * y) of math
　　　　set t to calcAtan2(z, p * bda) of math
　　　　set stNum to calcSin(t) of math
　　　　set ctNum to calcCos(t) of math
　　　　set b to calcAtan2(z + e2 * a / bda * stNum * stNum * stNum, p – e2 * a * ctNum * ctNum * ctNum) of math
　　　　
　　　　set l to calcAtan2(y, x) of math
　　　　set sb to calcSin(b) of math
　　　　set rn to a / (calcSqrt(1 – e2 * sb * sb) of math)
　　　　set h to p / (calcCos(b) of math) – rn
　　　　set l1 to b / rdLocal
　　　　set l2 to l / rdLocal
　　　　return {l1, l2, h}
　　end xyz2llh
　　
end script