※ 本件、別途記事上げますが解決しました。ご協力ありがとうございました。
数日前から@Say_noさんと、米軍の70年前の日本都市地域図、Japan City Plans 1:12,500 U.S. Army Map Service, 1945-1946をタイル地図化するための謀略を行っております。
人手と時間さえかければ、GCP*1を与えまくってのむりくり投影変換=>タイル化、も可能なのですが、人も金もなかりけりなので、投影系の情報を読み解いて一気に変換したいと思っています。
とりあえず、地図上の凡例から、ONE THOUSAND YARD WORLD POLYCONIC GRID BAND IIIn, Zone A(東日本)又はB(西日本)というのが読み取れます。
それをヒントに検索しまくったところ、こちらの記事をみつけまして、
- 当時の米軍は米国内地図で多円錐図法を使っており、基準経線は8度刻みでZone AからGまであったらしい
- 基準緯線は北緯40.5度だったらしい
- その座標系の単位はヤードで、基準経線/緯線での座標値はx:1000000、y:2000000だったらしい
- 測地系はNAD27(準拠楕円体はclrk66)だったらしい
といった情報がありました。
だったら日本の地図作成時も同じ基準流用してるだろ、と思って、全地図見直してみると、
- X座標の1000000は、ZoneAでは東経143度、ZoneBでは東経135度付近にあるっぽい。両者の度間隔は8度
- Y座標の2000000は、北緯40.5度付近にあるっぽい
というのが判って、これは間違いなく流用してるなと。
というわけで、この座標系のproj.4定義は、おそらく
ZoneA:
+proj=poly +lat_0=40.5 +lon_0=143 +x_0=914398.5307444400 +y_0=1828797.0614888800 +ellps=clrk66 +to_meter=0.9143985307444408 +no_defs
ZoneB:
+proj=poly +lat_0=40.5 +lon_0=135 +x_0=914398.5307444400 +y_0=1828797.0614888800 +ellps=clrk66 +to_meter=0.9143985307444408 +no_defs
だろうというのが判りました*2。
が、これでタイル化してやると、領域はほぼ正しく変換されるのですが、微妙に合わない。
場所によって違いますが、だいたい南東方向に数百mずれる。
これはきっと、米軍はあらたに測量し直したわけじゃなくて日本の旧版地形図の測量成果を流用してるでしょうから、旧日本測地系と世界測地系間の楕円体中心ズレがそのまま入ってきてるんじゃね?と思って、試しに旧日本測地系での楕円体ズレパラメータを入れてみました*3。
ZoneA:
+proj=poly +lat_0=40.5 +lon_0=143 +x_0=914398.5307444400 +y_0=1828797.0614888800 +ellps=clrk66 +to_meter=0.9143985307444408 +towgs84=-146.336,506.832,680.254,0,0,0,0 +no_defs
ZoneB:
+proj=poly +lat_0=40.5 +lon_0=135 +x_0=914398.5307444400 +y_0=1828797.0614888800 +ellps=clrk66 +to_meter=0.9143985307444408 +towgs84=-146.336,506.832,680.254,0,0,0,0 +no_defs
すると、数百mもズレるというような事はなくなり、場所によって異なるものの東西はほぼ一致、南北だけ数十mずれる感じにほぼおさまりました。
が、残りの数十mズレをどう抑え込めばいいのか判りません。
旧日本測地系の成果を流用してるのなら、楕円体をbesselにすればいいんじゃね?*4とか、楕円体ズレパラメータもヤードに直すのか?とかいろいろやってみましたが、いっこうに合う様子がない。
そして最終的に、これはもう楕円体が違う時点で、楕円体中心ズレパラメータもだいたい似た領域にいるのは間違いないけど、微妙に修正しないといけないんじゃね?という結論に達しました。
が、その算出方法がわからない。
正確にはパラメータ振って最小二乗法で検定して…とか判るんだけど、使いこなした経験がありません*5。
なので、もしそれが得意な方おられるならば、どなたかにサクッと出してもらえればなー、と思って記事共有します。
やっていただきたい事は下記の通りです。
評価する値を出力するためのproj.4スクリプトは、
cs2cs +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 +no_defs +to +proj=poly +lat_0=40.5 +lon_0=135 +x_0=914398.5307444400 +y_0=1828797.0614888800 +ellps=clrk66 +to_meter=0.9143985307444408 +towgs84=-146.336+δx,506.832+δy,680.254+δz,0,0,0,0 +no_defs -I -f '%10.5f'
です。
入力が米軍地図座標値(ZoneB)、出力がGoogleメルカトル座標値になります*6。
この中のδx,δy,δzの値をいろいろに振ってもらって、下記の米軍地図座標値を入力し、出力として出てくるメルカトル座標値を下記の比較値と比較し、誤差が最小になるδx,δy,δzの値を見つけてください。
パラメータは記事書き始めた時点で6点しかないですが、もうちょっと集めて追って追加します。
GISエキスパートの方、数学の得意な方の、ご協力をいただけるとありがたいです。
データ:
直江津起源データ
入力米軍地図座標値X | 入力米軍地図座標値Y | 比較メルカトルX | 比較メルカトルY |
1316000.25716 | 1601177.52710 | 15390049.84597 | 4463284.59138 |
1313461.53062 | 1600634.89841 | 15387001.91832 | 4462651.74175 |
1317747.73014 | 1603350.73837 | 15392113.70933 | 4465785.62674 |
下松起源データ
684380.42133 | 1215480.00090 | 14680025.19063 | 4028738.91679 |
682682.42843 | 1217635.20624 | 14678142.77804 | 4031064.79392 |
680005.65779 | 1218443.14815 | 14675188.35875 | 4031957.93863 |
*1:Ground Control Point。地図上のここは実際の地物と一致する、と特定できた点のこと
*2:原点の値が1000000,2000000ではなく変な小数値になってるのは、ヤードからの換算をしているためです
*3:すみません、以下一人で大変なので、本当はZoneBの下松、直江津の二ヶ所の地図でしか確認してません。記事作成後、追って追試し、問題があれば修正加えます
*4:よけいエラい事になった
*5:Excelのソルバーレベルなら何度も使ってますが、プログラムした経験は大学の時の演習くらいでしかありません。RとかMATLABとか?みたいな便利なツールもあるんだと思うんですがそれも知りません
*6:経緯度でやらない理由は、経度と緯度で最小二乗法したところで意味がないと思うので