[追記]本稿後段のcs2csでプレーン楕円体同士が無変換になってしまう問題は、次稿で原因が判りました。ご興味があれば併せてお読みください。
10年来の位置情報知識の土台がガラガラ崩れていくのが聞こえる…。
普通にcs2csで測地系/投影系変換叩いているだけなら一生気付く事もなかった問題ですが、emscriptenでproj4のJavaScriptコンパイルに手を出す => だったらちゃんと動いてるか確認のためのテストデータ作らないといけないよね、からずるずると深みに。
projコマンドとcs2csとの違い
projコマンドはWGS84から指定した投影系に変換するコマンドと思っていたので、
$ proj +投影系定義
と
$ cs2cs +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +to +投影系定義
は等価だと思ってたのですが、テストデータ作成のため両方叩いてみると、
$ proj -f %.10f +proj=tmerc +lat_0=0 +lon_0=-62 +k=0.9995000000000001 +x_0=400000 +y_0=0 +ellps=clrk80 +towgs84=725,685,536,0,0,0,0 +units=m +no_defs
-61.35 15.42 0
469737.4006643681 1704564.7371689796 0
$ cs2cs -f %.10f +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +to +proj=tmerc +lat_0=0 +lon_0=-62 +k=0.9995000000000001 +x_0=400000 +y_0=0 +ellps=clrk80 +towgs84=725,685,536,0,0,0,0 +units=m +no_defs
-61.35 15.42 0
468774.3915903920 1704156.5394733255 14.4524246091
あれ。なんだこれ？
等価のはずなのに値が違う…。
いくら調べても原因が判らないので、proj4のメーリングリストに質問投げつつ調査を続けると、towgs84スイッチを持つ定義でだけテスト結果でエラーが出ているのに気付く。
もしかしてprojではtowgs84が効かないの？と思って試してみると、
$ proj -f %.10f +proj=tmerc +lat_0=0 +lon_0=-62 +k=0.9995000000000001 +x_0=400000 +y_0=0 +ellps=clrk80 +units=m +no_defs
-61.35 15.42 0
469737.4006643681 1704564.7371689796 0
$ cs2cs -f %.10f +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +to +proj=tmerc +lat_0=0 +lon_0=-62 +k=0.9995000000000001 +x_0=400000 +y_0=0 +ellps=clrk80 +units=m +no_defs
-61.35 15.42 0
469737.4006643681 1704564.7371689796 0.0000000000
ビンゴ。
同時にML側でも答えてくれた方がいて、やっぱりtowgs84は無視されるとの内容で裏付けが取れた。
ひえー、proj4使い出して10年近くなるけど気付いてなかった！
なるほど、projコマンドはprojectionの名前のとおり、汎用座標変換コマンドじゃなくて、経緯度から地図投影系に投影するコマンドだったわけね。
でも、なんで元経緯度座標系が、WGS84なんだろ？
単に経緯度から地図投影系への変換なら、各地図投影系の準拠楕円体での経緯度からの投影にすればいいのに…。
とそこまで考えたところで、嫌な可能性に思い当たりました。
まさか…まさか？
楕円体中心ズレの定義されない楕円体同士は、cs2csにかけても変換結果は同じ
towgs84定義のないプレーンな楕円体同士の経緯度変換、やってみました。
$ cs2cs -f %.10f +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +to +proj=longlat +ellps=clrk80 +no_defs
-61.35 15.42 0
-61.3500000000 15.4200000000 0.0000000000
$ cs2cs -f %.10f +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +to +proj=longlat +ellps=bessel +no_defs
-61.35 15.42 0
-61.3500000000 15.4200000000 0.0000000000
$ cs2cs -f %.10f +proj=longlat +ellps=clrk80 +no_defs +to +proj=longlat +ellps=bessel +no_defs
-61.35 15.42 0
-61.3500000000 15.4200000000 0.0000000000
同じだ…同じだよ…寸分の差もないよ…。
高度を0以外にして変換しても一緒でした。
その昔、なんか日本測地系とか世界測地系とかいうのがあって系が違うんだ、とか覚えた時、その違いは選んだ楕円体が違うからなんだ、と教えられた気がする。
測量誤差による楕円体中心ズレの差も（あと、もちろんメッシュそのものの非線形歪みも）あるけどそれは従で、飽くまで主は楕円体そのものの選び方なんだ、という印象を受けた気がする。
その後、自分でいろいろ試して、パラメータとしては明らかに中心ズレの方が大きいのは知ってたけど、GIS的な考え方としては楕円体選択が違うから経緯度がズレる、というのが正しいと思い込んでいた。
まさか、中心ズレのないプレーンな楕円体同士で変換すると全く経緯度差がなくて、事実上100%楕円体中心ズレ起因でズレるというのは、今初めて気付かされた。
なんか10年間に身につけた知識の基礎がガラガラと崩れてる気がする…。
まあ極座標的に考えると、正しいのかもしれないが
どの楕円体を採用しようと楕円体表面は楕円体表面に変換されると仮にして、その上で経緯度だけで考えると、経緯度というのは極座標だから、どんな楕円体を選ぼうと中心ズレがなければ一致する、というのは正しいのかもしれない。
が、towgs84パラメータというのは極座標ではなく楕円体中心を原点とした3次元直交座標系ベースのパラメータだし、3次元直交座標系ベースで考えれば、同じ経緯度でも違う楕円体の表面は違うxyz座標値、つまり違う点になるはず。
そして、高度の定義にもよるけれども、高度の方向が楕円体中心を通る線ではなく局地での楕円体表面の法線方向であるならば、楕円体表面からの高度が0でない点は、別の楕円体上に移した際には（楕円体中心ズレがなくとも）経緯度が異なる、ということもあるはず。
準拠楕円体の違いによる経緯度の差というのは、そういう形で発生するとずっと思い込んでいたんだけど、まさか中心ズレがなければ、どの楕円体でも同じ経緯度、同じ高度に変換されるとは、思ってもいませんでした。
いやあ、驚いた。
ほんと。