前記事でお願いしたご協力ですが、自分でも効率的なアルゴリズムなんかは判らないものの、とりあえずパラメータ振ればいいんだろ、と半自動半手作業でパラメータ推定し、少なくとも局所最適の解も見つけたのですが、適用した結果はまだ納得のいくものではありませんでした。
局所最適に捕まって実際の解に辿り着かなかったのか、そもそもやり方が間違ってたのかは謎です。
ただ、得たパラメータでの結果はとんでもない結果ではなく、方角は一致して経緯度に対し斜め方向に平行移動した結果とかだったので、方法的には間違ってなくて局所最適解をみつけてしまったのだと思います。
で、いかんともし難いので、ちょっと考えてみました。
なぜこんな、独自投影系パラメータを推定しないといけないような投影系ができてしまったのか。
時は70年前に遡ります。
以下、取材等を行ったわけではなく、この地図に起きていた事から考えた想像図です。
戦争に勝ち日本を占領した米軍。
統治のために早く地図を整備する必要があったけれど、新しく全土を測量し直す余裕もない。
困った挙げ句やっちゃったのが、日本が昔から自分たちの測地系（旧日本測地系）で測量してきた成果を、自分たちの測地系（NAD27）に等しいと考えて同じ経緯度値を採用し、その結果で多円錐図法の投影を行っちゃった。
すばやく統治のために必要な地図も整備できて、めでたしめでたし…じゃねーよ！
お前らがサボったから今俺が苦労してるんだよ！
と、ここまで妄想歴史を想像して突っ込んだところで、気付いてしまいました。
なんで一段階変換にこだわってたんだろう？
歴史上起きた事象の逆をなぞればええ話やん！
旧日本測地系の地図をNAD27とみなしたために地図投影系がおかしな事になってるなら、一度この地図をNAD27 /w 多円錐図法パラメータからNAD27経緯度パラメータでの投影系に変換し、「その結果」を日本測地系の地図とみなして、他の測地系、投影系に変換すればええやん！
で、やってみると…備後でした！
Courtesy of the University of Texas Libraries, The University of Texas at Austin.
OpenStreetMap
いやどこからどうみても備中です本当にありが(ry*1
こんなしょうもないネタのためにサクッと変換するのも苦にならないほど簡単に、かつこれまでより桁違いに正確に、変換できるようになりました。
なんか問題にぶち当たった時は、難しく考えるより基礎に立ち戻る事が重要と改めて感じました。
で、結論としてJapan City Plansの地図群をタイル地図変換するには、以下の手順で可能な事がわかりました。
地図上のピクセル座標点と、地図座標系座標値の対応付け点を4点ほど取得。
4点もいるかどうかは判りませんが、4隅に近いところで調べておくと確実かと。
地図座標系座標値は、下の図のように、地図上で座標罫線っぽいの（これは経緯度線ではなく地図座標線です）が交差しているところで読み取ります。
この例だと1399、1208と書かれていますが、下3桁が省略されていますので実際の座標値は1399000,1208000です。
この座標値と、左上原点でのピクセル座標値との対照付けを行います。
いわゆるジオリファレンス作業ですが、地形の特徴点とか探す必要もないのでサクサク進むはず。
QGISを使える人ですと、ジオリファレンサーというプラグインで画像表示、ピクセル特定と座標の入力、ファイルへの書き出しまで出来ますが、それで標準化すると使える人しかできない作業になってしまいますので、ペイントツールや画像ビューアでも手分け作業できるよう、QGIS必須にしない方がいいのかなと思います。
調べた対応点リストは、下記のフォーマットでファイルにしておいておきます。
例：並びは地図座標X,地図座標Y,ピクセル座標X,ピクセル座標Yの順
889000.000000000000000 1296000.000000000000000 416.950351068112354 499.489423234597439 889000.000000000000000 1290000.000000000000000 432.044043767982259 3212.580686036285897 895000.000000000000000 1287000.000000000000000 3162.744614719519177 4543.341259074852132 897000.000000000000000 1296000.000000000000000 4053.272484011868983 452.950537409990659
ワールドファイルを生成し、元jpg画像を位置情報付き画像化
こちらから、対応座標データから画像用位置情報付与メタファイル（通称ワールドファイル）を生成するスクリプトを取得します。
実行するとmatfuncライブラリがないと言われる場合は、同じ場所にあるので取得できますが、pythonのstatlibモジュールを導入済みの場合は、ソースコードの from matfunc を from statlib.matfunc に変えるだけでも動きます。
取得できれば、下記を実行。
python gcp2wld.py -i [対応点情報ファイル名] > [地図画像ファイルの拡張子除いた名前].jpw
これで、対応ライブラリでは地図画像ファイルと連動して動作し、位置メタ情報を提供するファイルが出来ました。
米軍がやった投影変換の逆変換作業
米軍の使った地図投影系は、今日の日本全体で3つありますが、沖縄は今回の地図画像中には含まれていない*2ので、実際には東経143度を中心としたZONE A、135度を中心としたZONE Bの2投影系があります。
どっちが使われているかは、地図の凡例領域を見れば書いてあります。
それぞれの座標系の変換パラメータ*3は、
ZONE A
+proj=poly +lat_0=40.5 +lon_0=143 +x_0=914398.5307444408 +y_0=1828797.0614888816 \ +ellps=clrk66 +to_meter=0.9143985307444408 +no_defs
ZONE B
+proj=poly +lat_0=40.5 +lon_0=135 +x_0=914398.5307444408 +y_0=1828797.0614888816 \ +ellps=clrk66 +to_meter=0.9143985307444408 +no_defs
のような感じになります。
このパラメータを使って、下記のコマンドを実行します（\部分の改行は含みません）。
gdalwarp -srcnodata 0 -dstalpha -s_srs '地図投影系に従った変換パラメータ' \ -t_srs '+proj=longlat +ellps=clrk66 +datum=NAD27 +no_defs' 地図画像名.jpg 中間生成物名.tif
これで、米軍が元にした日本の元地図に等しい経緯度座標地図画像ができました。
所望の座標系、フォーマットのデータに加工
先のコマンドで出来た中間生成物は、日本測地系*4地図画像です。
作業の都合上、GeoTiffの中には北米のNAD27測地系*5の情報が埋め込まれてしまってますので、パラメータ無指定で処理するとそれとして扱われてしまいますが、日本測地系のパラメータを与えてやればちゃんとそのように動作します。
よって、例えばgdal2tilesでタイル画像にしたいなら、
python gdal2tiles.py -s '+proj=longlat +ellps=bessel +towgs84=-146.336,506.832,680.254' \ -a 0 中間生成物名.tif
みたいにしてやれば、ぴったり合うタイル地図ができますし、世界測地系GeoTiffに直したいのならばそれに適したコマンドを打てばOKです。
タイル化地図提供サービスを検討中
変換情報に関してはという感じですが、今Tilemap.JP周りでこの地図のタイル提供を検討しています。
既にテキサス大の了承は得ており、今のところCC0あるいはCC-BYで出す予定です。
手が回っていない事もありまだ公開まで少しかかるかもしれませんが、気長にお待ちいただければ幸いです。
