いろいろ調べると判ってくるもんだなー。
Proj4Textが大分暗号ではなくなってきた。
場所によって異なるGoogle Mercator定義
Google Mercator (EPSG:3857)のProj.4定義ですが、あるところでは、
+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs
と定義されてるのが、別のところでは
+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +a=6378137 +b=6378137 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
と書いてあったりして、定義が違う。
nadgrids設定の有無で、変換結果が異なる
大きな差は、nadgridsという定義(wktextもあるけど、これは今回調べなかったのでおいておいて)。
これがあるかないかで結果が違う。
$ cs2cs -f %.10f +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +to +proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +a=6378137 +b=6378137 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
135 35 0
15028131.2570919339 4139363.7974272855 -6995.7693383824
$ cs2cs -f %.10f +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +to +proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +a=6378137 +b=6378137 +towgs84=0,0,0,0,0,0,0 +nadgrids=@NULL +units=m +no_defs
135 35 0
15028131.2570919376 4163881.1440642942 0.0000000000
その差、20000m以上…
え？メルカトル投影なので、赤道上でなければ誇張されているとはいえ、10kmレベルでズレるって事？
nadgrids=@nullは準拠楕円体の違いによる経緯度変換を抑制する設定らしい
なんぞこれ？と思って調べてみると、
Changing Ellipsoid / Why can't I convert from WGS84 to Virtual Earth Mercator?
to convert between WGS84 and mercator on the sphere there will be substantial shifts in the Y mercator coordinates. This is because internally cs2cs is having to adjust the lat/long coordinates from being on the sphere to being on the WGS84 datum which has a quite differently shaped ellipsoid.
In this case, and many other cases using spherical projections, the desired approach is to actually treat the lat/long locations on the sphere as if they were on WGS84 without any adjustments when using them for converting to other coordinate systems. The solution is to "trick" PROJ.4 into applying no change to the lat/long values when going to (and through) WGS84. This can be accomplished by asking PROJ to use a null grid shift file for switching from your spherical lat/long coordinates to WGS84.
地図投影系としては楕円体近似に球を採用している座標系に対し、楕円体として測量された経緯度をそのまま持ってきて、投影だけ球とみなして投影する、そのための設定がnadgrids=@nullらしい。
考えてみりゃ確かにそうだ。
これまでの連続記事で見てきたように、towgs84で地球中心にバインディングしたら、準拠楕円体が違う時点で経緯度がまず変わる。
違った値になった経緯度で投影したら、ズレるのは当たり前。
towgs84、或いはdatum無指定で経緯度ベースで結びつけるのと似ている
経緯度の値を変換させないと言う事は、towgs84或いはdatum無指定で、GNDをふらふらさせて変換するのと同じような話、という事だ。
試してみる。
WGS84側のdatum指定を外して、
$ cs2cs +proj=longlat +ellps=WGS84 +no_defs +to +proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +a=6378137 +b=6378137 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
135 35 0
15028131.2570919339 4163881.1440642942 0.0000000000
ほらほらキタワァ。
nadgrids指定の場合と完全一致。
それならば、towgs84指定しなけりゃいいんじゃないの=>そうでもない、違いもある
それなら、towgs84を指定しなきゃいいんじゃないの、と思うかもしれないがそうでもない。
使う経緯度がWGS84だけなら大した差はないのだけれども、日本測地系なんかが入ってくると話が変わる。
下がその例。
$ cs2cs -f %.10f +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +to +proj=longlat +ellps=bessel +towgs84=-146.336,506.832,680.254,0,0,0,0 +no_defs
135 35 0
135.0027926037 34.9968021017 -50.4782586880
これで得られた日本測地系経緯度を、Google Mercatorに変換してWGS84からの直接変換と近似の値が出れば、正しい座標系定義となるだろうが、
$ cs2cs -f %.10f +proj=longlat +ellps=bessel +towgs84=-146.336,506.832,680.254,0,0,0,0 +no_defs +to +proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +a=6378137 +b=6378137 +towgs84=0,0,0,0,0,0,0 +nadgrids=@NULL +units=m +no_defs
135.0027926037 34.9968021017 -50.4782586880
15028131.2570958808 4163881.1440635058 0.0000000019
towgs84とnadgridsを指定すると、WGS84から直接指定とほぼ等しい結果が出ているが、
$ cs2cs -f %.10f +proj=longlat +ellps=bessel +towgs84=-146.336,506.832,680.254,0,0,0,0 +no_defs +to +proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +a=6378137 +b=6378137 +units=m +no_defs
135.0027926037 34.9968021017 -50.4782586880
15028442.1283138040 4163446.5709505822 -50.4782586880
towgs84とnadgridsを指定しないと、大分ズレた値が出ている。
要するに、towgs84を指定しないと、経緯度側がWGS84だろうがTokyoだろうが、そのまま経緯度として無変換で持ってきてしまうのに対し、
towgs84は指定するけどnadgrids=@NULLも指定する、というのは、WGS84以外の経緯度からはWGS84には一旦変換するけど、WGS84からはそのまま持ってくる、という意味になり、WGS84以外の測地系が混ざれば全然異なる結果になる。
よってより正しいGoogle Mercatorの定義は
よって、これまでの議論より、冒頭に挙げたGoogle Mercatorの定義にはどちらにも問題がある。
一方はtowgs84定義がなく、一方はnadgrids定義がない。
wktextの意味とかまだ判らないので、これが絶対正しいというわけではないが、よりベターなGoogle Mercatorの定義は、
+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0.0 +k=1.0 +units=m +towgs84=0,0,0,0,0,0,0 +nadgrids=@null +wktext +no_defs
になるだろう。
nadgrids, wktextについてもっと調べてみる
nadgrids=@nullの意味は判ったが、nadgridsには他にもいろんな設定があるようだし、wktextの意味は判らないので、引き続き調べてみよう。
[追記]
Proj.4のWikiの方のFAQに、wktextの説明もあった。
When using this coordinate system with GDAL/OGR it is helpful to include the +wktext so the exact proj.4 string will be preserved in the WKT representation (otherwise key parameters like +nadgrids=@null will be dropped)
どうも、WKTでの投影系定義との相互運用の際に、必要なフラグらしい。
WKTに詳しくないので重要さが判らないし、またproj.4自体の動作にも影響なさそうだが、nadgrids=@null が影響ある設定として挙げられているので、付けておいた方がよさそう。