年の瀬に横浜（一人暮らし新居）と奈良（家族の家）を行ったり来たりのkochizufanです。
引っ越しやら行政手続きやら家財揃えに忙殺されて投稿遅れてすみません。
ネタは、単独画像でない、メッシュ単位に分割された大量のGeoTiff画像等を、多量のGCPから全部ジオリファレンスする際のノウハウについて。
そんなユースケースたいしてないだろうって？いや場合によってはあるんじゃないすか？
例えば...例えば...明治初期のまだ伊能図並み精度でしか測量成果がなかったために位置ズレしてる北海道旧版地形図上のアイヌ地名を、ジオリファレンスで正確な位置に戻すとか...1。
こういうをジオリファレンスする際に、各図郭内に含まれるGCPだけでジオリファレンスすると、各図郭の切れ目で不整合を起こしてしまいます2。
解決法は単純で、gdal_translate等の-gcpオプションは、各図郭のextent外の座標と対応付くgcpでもちゃんと受け付けてくれるので3、全部のgcp4を指定してやれば、隣り合う図郭との不整合も生じずに全てジオリファレンス変換できます。
...と、非常に単純なんですが、いくつか試行錯誤した気をつけるべき点もあるので列挙します。
-
各図郭の外側には透明のバッファ領域を設ける
非線形変形ですので、うにょーんと変形された後の画像は、元の画像のextentからはみ出す場合もあります。
そのような場合、gdalwarpでは、右と下方向5にはみ出した分は勝手に画像が拡大してくれますが、左や上方向にはみ出した分は元のextentで頭切れトンボ？になってしまいます。
これを防ぐために、元画像の外郭にconvertコマンド6等を用いて、透明なバッファ領域を設けましょう。 -
gdalwarpでは変換元の投影系ではなく、GCPの投影系を指定する
元の画像は日本測地系の平面直角だから、-s_srsの指定はそれで...とか考えだすと訳がわからなくなります。
gdal_translateでGCPの座標を埋め込んだ時点で、そのGCPの投影系になっているので、gdalwarpの-s_srsもその定義のprojtext7を指定する必要があります。 -
大量のGCPの指定は、--optfileオプションを使う
GCPの指定が大量であった場合、コマンドラインオプションで指定すると引数が多すぎるというエラーが出て実行できない場合があります。
こういう時は、オプション引数の記述をファイルに書き出して、--optfileオプションでそのファイル名を指定すると、コマンドラインオプションの制限以上のGCP指定を行う事ができます。
以上のような注意点を踏まえ、実際の手順を列挙すると、
-
元画像（GTiffとする）に対し
gdalinfo src.tiff -
convertコマンドで
convert src.tiff -bordercolor "#00000000" -border "1500x1500" mid1.tifで透明バッファ領域を発生させます10。
-
Origin、Pixel Sizeとバッファ幅より、GCPのピクセル座標値（extent外のものも含む）を計算し、
data.gcp-gcp ピクセルX1 ピクセルY1 lat1 lng1 -gcp ピクセルX2 ....といった形でファイルに書き出します。
-
gdal_translateを実行します。
gdal_translate -of "GTiff" -a_srs "+proj=longlat +ellps=WGS84 +datum=WGS84 +units=degrees" --optfile data.gcp tmp1.tif tmp2.tif -
gdalwarpを実行します。
gdalwarp -tps -s_srs "+proj=longlat +ellps=WGS84 +datum=WGS84 +units=degrees" -t_srs "+proj=longlat +ellps=WGS84 +datum=WGS84 +units=degrees" -of "GTiff" -dstnodata 0 tmp2.tif result.tif
これで、複数画像からなる巨大地図でも、切れ目不整合が生じない形でジオリファレンスすることができます。