ちずぶらりHackers

オープンデータ化も進んでいる古地図、絵地図を扱うiOSアプリ、ちずぶらりをハックし倒します。iOS / Androidでのクローンアプリの開発も予定しています。

『Japan City Plans』古地図の投影系パラメータ推定にご協力お願いします

※ 本件、別途記事上げますが解決しました。ご協力ありがとうございました。

数日前から@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:経緯度でやらない理由は、経度と緯度で最小二乗法したところで意味がないと思うので

© TileMapJp/歴史国土/地図タイル工法協会