こんにちは、さくらインターネット研究所の菊地です。
いきなりですが、例えば、
- 複数の自律したドローンが撮影した空撮画像を組み合わせたい
- 地図上に、分析のために別の画像を貼り付けたい
といった地図上に画像を貼り付けたいとしたときに、どのように作業しますか？ 特に、位置や大きさを正確に扱いたいときにはどうしたらいいでしょう。今回、そういった地図画像処理時のノウハウを若干ですがまとめましたので公開いたします。
さくらインターネットではビッグデータ活用のビジネス展開の一環として、衛星データ活用プラットフォームの構築・サービスに向けた取り組みなどを実施しております（参考）が、今回研究所ブログでは、自前の画像を地図と組み合わせて扱う場合の、位置情報の埋め込みや画像の加工に関する基本的な知識をまとめてみました。
今回のノウハウでできるようになること
まず、イメージしやすいように、結果画像を用いて説明してみます。
左の図は、さくらインターネットのオフィスのある位置にさくらインターネットのロゴマークを貼り付けてみた図です。これはマウス操作で画像貼り付け位置を指定したのではなく、緯度経度の数値を画像に埋め込み、その画像を地図ソフトで読み込ませたものです。右の図は元画像をクリッピングし、回転させて貼り付けたものです。
これらのような、位置の指定や画像の加工などに関するノウハウになります。
前提となる知識
前提と書いていますが、むしろ、本記事記載のコマンドを試すときにわからないことがあったら参照する、という感じで使ってください。なお菊地の理解に基づいて書いておりますので、不足や誤っている可能性もあります。併せて表記したリンク先などの参考資料も読んでみてください。
地図画像データフォーマット
GeoTIFF
地図に画像を貼り込むためには、画像データそのものと、その画像が地図上のどの位置に貼り込まれるのかという情報が必要です。これらを扱うには複数の方法がありますが、今回は、GeoTIFFという形式の画像ファイルを使います。GeoTIFFは名前の通りTIFF画像に位置情報を埋め込んだものです。
地図に関する知識
測地系
地球上での位置を数値（緯度経度）として扱う際に、その座標をどのように定義しているかの方式を測地系といいます。測地系は開発の経緯などから複数のものがあります（米国の測地系WGS84や日本測地系2000、日本測地系2011など）が、現在において実用上は、WGS84を使う、と覚えておけば十分でしょう。WGS84はGPSで使われている測地系です。
（余談）なお、かつてアップルがiOS6発表時に独自マップを作ったところその出来がひどかったという事件がありました。このときの問題の一つに、地図上に表示されている建物などの位置が実際の場所とずれているという現象がありましたが、これは測地系を正しく扱っていなかったことによるものと言われています（（旧）日本測地系による位置情報でWGS84の地図システム上に表示したため）。しかし現在は、日本測地系（日本測地系2000/日本測地系2011）はWGS84とほぼ同じとなっているため、仮に日本測地系で記述されている位置情報を利用する場合でも、このような問題は発生しません。
緯度経度の度分秒表記と10進数表記
位置を緯度経度で表すときには、一般的に日本語では「北緯35度41分39秒800、東経193度41分44秒00」といった具合に、度・分・秒という表記を使います。数字記号で表現する場合には「35° 41’39.800″N, 139°41’44.00″E」などと書きます。このときの分と秒の数え方は60進法です。しかしこれはコンピュータ上で扱う際に面倒なので、分・秒の部分を10進数に変換した形式というものが用いられています。先程の緯度経度を10進表記で表すと「35.69438888,139.69555555」となります。ごく簡単に言えば、分を1/60し、秒を1/3600して足します。詳しくは、下記のサイトなどを参照してください。
（余談）「°」を非日本語環境でキーボードから入力するのは至難の業です（私はいつもコピペで対応しています…）。これが大変なのは各国共通なのか「°」の代わりに「d」を使う実装もあります（本記事で紹介するgdalなど）。この他、上述の度分秒10進表記のほか、東経西経北緯南緯の（EWNS）の表し方、その順番などが、コマンド・ツールによってまちまちで、これを間違うとコマンドが動かなかったり、思ったとおりの場所ではなく見当外れな場所を指し示してしまったりといった事が起こります。こういった細かい所の違いが「地図情報を扱うのは大変」という印象を作り出している一因になっていると感じています。
平面直角座標系
上の測地系と言葉が似てますが全く違うものですので混同しないようにしましょう。緯度経度は地球上の一点の位置を表すものですが、二点間の距離を表す場合などには不便です（同じ水平距離の1度が緯度によって異なる長さになる、など）。このため、地面を平面とみなして直角座標で扱えるようにしたものが平面直角座標系です。ただしこの平面直角座標系は球形の地球を平面にして扱うことで歪みが発生するため、一つの座標系で広い範囲をカバーできません。そのため日本においては19の平面直角座標系が定義されています（同様のものが各国にあります）。このため、地図に画像データを貼り付ける際には、どの地点に画像データを貼り付けたいかによって19の平面直角座標系を使い分ける必要があります。
平面直角座標系の必要性や定義などについては、以下の国土地理院のサイトにその名の通り非常にわかりやすく書いてありますので、そちらを参照してください。変換ツールもあります。Wikipediaの記述も参考になります。
ツール・コマンド
実際に画像ファイルに位置情報を埋め込んだり、それを地図上に表示して結果を確認するのには以下のツールを使いました。
gdalコマンド群
Geospatial Data Abstraction Library – GDALは、地図情報を取り扱うための様々なコマンド・ライブラリを集めたツール集です。統合型のGISアプリケーション上ではなく、コマンドラインやプログラムから地図データなどを扱いたいときには必ずお世話になるツールということのようです。
Google Earth プロ
GeoTiff画像を読み込んで地図上に表示するおなじみのツール。Google Earthは最近はブラウザ上でも動きますが、GeoTiff画像のインポートには「プロ」というスタンドアローン版をインストールする必要があります。
名前からして有料アプリなのかと思いがちですが、現在はどうやら無料なようです。Google Earth（無印）がブラウザ向け、プロがスタンドアローンということのようです。
環境の構築
gdalコマンドをインストールするには以下の２つの方法があるようです。なお、以降の動作はMacOS10.13環境下で確認しています。
ソースコードを持ってきて自前でコンパイル・インストールする
ソースコードはここにあるので、最新版をダウンロードしてきます。記事執筆時点では、バージョン2.3.2です。その後、展開してconfigure:make:make installでいけます。
PC103429:gdal s-kikuchi$ tar zxvf gdal-2.3.2.tar.gz PC103429:gdal s-kikuchi$ ls gdal-2.3.2/ gdal-2.3.2.tar.gz PC103429:gdal s-kikuchi$ cd gdal-2.3.2 PC103429:gdal-2.3.2 s-kikuchi$ ./configure (...) PC103429:gdal-2.3.2 s-kikuchi$ make (...長い時間かかります) PC103429:gdal-2.3.2 s-kikuchi$ make install (...) PC103429:gdal-2.3.2 s-kikuchi$
止まること無くコンパイル・インストールできますが、時間がかかります（MBP13 2017の環境で10分くらい）。
Docker環境でLinux用パッケージを利用する
docker環境がすでにあるなら次のdockerで環境を作るほうが簡単です。ただしこの場合、用意されているパッケージが古いという問題があります。gdal_translateのバージョンは1.11.3でした。しかし本稿のコマンドを試す範囲内ではこのバージョンでも問題ありません。
適当な作業ディレクトリに以下のようなDockerfileを用意します。
FROM ubuntu RUN apt-get update && apt-get -y upgrade \ && apt-get install -y python-gdal \ && apt-get install -y gdal-bin
以下のようにコンテナをビルドして、実行させます。
PC103429:gdal s-kikuchi$ PC103429:gdal s-kikuchi$ ls Dockerfile PC103429:gdal s-kikuchi$ docker build -t 'my/ubuntu' . Sending build context to Docker daemon 1.128GB Step 1/2 : FROM ubuntu ---> 6a2f32de169d Step 2/2 : RUN apt-get update && apt-get -y upgrade && apt-get install -y python-gdal && apt-get install -y gdal-bin ---> Using cache ---> ba64cb2771fb Successfully built ba64cb2771fb Successfully tagged my/ubuntu:latest PC103429:gdal s-kikuchi$ mkdir work PC103429:gdal s-kikuchi$ ls Dockerfile work/ PC103429:gdal s-kikuchi$ docker run -v $(pwd)/work:/work -it my/ubuntu bash root@9dd686378313:/#
ホストとコンテナ内で作業ディレクトリを共有しておくと、環境を隔離しつつ画像データファイルなどはホスト側環境ですぐ用意・確認できるので大変便利です。
画像の処理の作業は、dockerコンテナ内で/workディレクトリ下で実行します。以下の各実行結果は基本的このコンテナ環境下で実施しているものです。
位置情報埋め込み画像(GeoTiff)を作る方法
gdal_translateコマンドを利用します。
緯度経度で指定する
まずは、最も基本だと多くの人が思うであろう、緯度経度による指定方法です。
root@9dd686378313:/work# root@9dd686378313:/work# gdal_translate -of "GTiff" -co "COMPRESS=LZW" -a_srs WGS84 -a_ullr 139.69555555 35.69438888 139.69951631666666 35.691382775 sakura_mark.tiff sakura_mark_shinjuku.tif Input file size is 358, 334 0...10...20...30...40...50...60...70...80...90...100 - done. root@9dd686378313:/work# ls sakura_mark.tiff sakura_mark_shinjuku.tif root@9dd686378313:/work#
コマンドリファレンスは以下となります。
$ gdal_translate -of "GTiff" -co "COMPRESS=LZW" -a_srs WGS84 -a_ullr 貼り付け先左上X座標 貼り付け先左上Y座標 貼り付け先右下X座標 貼り付け先右下Y座標 入力ファイル 出力ファイル
-of “GTiff”と-co “COMPRESS=LZW”はそれぞれ出力フォーマットと画像圧縮アルゴリズムの指定です。
-a_srs WGS84は、測地系の指定でWGS84、-a_ullrは入力画像を貼り付ける左上と右下の座標で、それぞれ経度緯度、10進数表記で指定します。
ここで経度緯度の順番に気をつける必要があります。地図を北を上にして見たときに、東西軸＝経度がX座標に相当し、南北軸＝緯度がY座標に相当します。経度＝X座標は東に向かうほど値が大きくなり、緯度=Y座標は北に向かうほど値が大きくなります。緯度については画像を貼り付けるときは右下Y座標は左上Y座標より値が小さくなることに注意してください。
変換（座標情報の埋め込み）ができたら、gdalinfoコマンドで確認してみましょう。
root@9dd686378313:/work# gdalinfo sakura_mark_shinjuku.tif
Driver: GTiff/GeoTIFF
Files: sakura_mark_shinjuku.tif
Size is 358, 334
Coordinate System is:
GEOGCS["WGS 84",
DATUM["WGS_1984",
SPHEROID["WGS 84",6378137,298.257223563,
AUTHORITY["EPSG","7030"]],
AUTHORITY["EPSG","6326"]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433],
AUTHORITY["EPSG","4326"]]
Origin = (139.695555549999995,35.694388879999998)
Pixel Size = (0.000011063594041,-0.000009000314371)
Metadata:
AREA_OR_POINT=Area
TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)
TIFFTAG_XRESOLUTION=72
TIFFTAG_YRESOLUTION=72
Image Structure Metadata:
COMPRESSION=LZW
INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left ( 139.6955555, 35.6943889) (139d41'44.00"E, 35d41'39.80"N)
Lower Left ( 139.6955555, 35.6913828) (139d41'44.00"E, 35d41'28.98"N)
Upper Right ( 139.6995163, 35.6943889) (139d41'58.26"E, 35d41'39.80"N)
Lower Right ( 139.6995163, 35.6913828) (139d41'58.26"E, 35d41'28.98"N)
Center ( 139.6975359, 35.6928858) (139d41'51.13"E, 35d41'34.39"N)
Band 1 Block=358x7 Type=Byte, ColorInterp=Red
Band 2 Block=358x7 Type=Byte, ColorInterp=Green
Band 3 Block=358x7 Type=Byte, ColorInterp=Blue
root@9dd686378313:/work#
Coordinate SystemとしてWGS 84が使用されていることと、Corner Coordinatesに四隅の座標が緯度経度で入っていることがわかると思います。
座標の値がそれっぽい値になっていたら、Google Earthプロに読み込ませてみましょう。「ファイル」メニュー内に「インポート」という項目がありますので、そこで該当ファイルを読み込ませます。
結果はこちら。
さくらのロゴの左上部分が、さくらインターネットの新宿のオフィスの上に来ているのがわかりますでしょうか。
なお本筋から少し外れますが、Google Erathに貼り付けると、視点を変えて俯瞰で見るなども自由自在です。さらにこの図では、画像を海抜50mの位置に置いて、建物の3D表示もONにしています。ビルが画像を突き抜けていてなかなか不思議な景色です。（画像自体は同じ位置に貼り付いているのですが、建物の3D表示による差異のため、画面の雰囲気はかなり変わります。）
インポート画像のパラメータ調整ダイアログ。標高もこのように調整できます。
ところで、上のように座標を埋め込むにあたって、その座標の値をどのように取得したらいいのでしょうか。貼り付け先左上座標は実際に貼り付けたい位置を緯度経度で取得すればよいのでそれほど難しくないでしょう。問題は貼り付け先右下座標です。今回のさくらのマーク画像のような（ほぼ）正方形の画像を例えば地図上で100m四方にして貼り付けたいとしたときに、左上座標の緯度経度から東に100m、南に100m行った地点の緯度経度をどうやって取得したらいいのでしょうか。さらに言えばこのさくらのマークは実は正確な正方形ではなく、358x334pixelの画像データです。このアスペクト比を保ったままで、東に358m、南に334m移動した右下座標の緯度経度をどのように取得すればよいのでしょうか。言うまでもなくもし右下座標の位置を間違えると、期待した大きさで貼り付けられなかったり、アスペクト比が狂ったりしてしまうのです。
緯度経度による指定は、実は思ったほど簡単ではないのです。そしてこのことが、平面直角座標系が存在する理由なのです。
ある基準となる地点から指定の距離だけ離れた別の場所の緯度経度を取得する一つのやり方として、平面直角座標系を経由して変換する方法があります。
- まず最初の基準地点（貼り付け先左上座標）の緯度経度を求めます。私はこの国土地理院の「平面直角座標への換算」ページを使いました。地図で任意のポイントを指定することができます。
- その基準地点の平面直角座標上の位置を取得します。同ページを使うと自動的に値が計算され表示されています。
- 基準地点（貼り付け先左上座標）から指定の距離だけ離れた地点を平面直角座標で計算します。東に358m、南に334m移動した地点としたければ、基準地点の平面直角座標のX座標に358足し、同Y座標から344引きます。平面直角座標は単位がm（メートル）です。またY座標は北に行くほど値が大きくなるので（WGS84と同じ）、北を上にして画像を貼る場合には引き算することに注意しましょう。
- ステップ3で求めた、貼り付け先右下座標の平面直角座標での値を緯度経度に（逆）変換します。同じく国土地理院のこのページが便利です。
計算にあたっては、加減算の別だけでなく、緯度経度、X座標Y座標の組み合わせを取り違えないように注意しましょう。１箇所でも間違えると地図画像が思ったとおりになりません。
基準点からの相対距離（平面直角座標系）で指定する
緯度経度を使わず平面直角座標系で指定する場合には以下のようになります。
root@9dd686378313:/work# root@9dd686378313:/work# gdal_translate -of "GTiff" -co "COMPRESS=LZW" -a_srs EPSG:2451 -a_ullr -12469.0624 -33897.3015 -12111.0624 -34231.3015 sakura_mark.tiff sakura_mark_shinjuku_heimen.tif Input file size is 358, 334 0...10...20...30...40...50...60...70...80...90...100 - done. root@9dd686378313:/work# ls sakura_mark.tiff sakura_mark_shinjuku_heimen.tif root@9dd686378313:/work#
コマンドリファレンスは以下となります。
$ gdal_translate -of "GTiff" -co "COMPRESS=LZW" -a_srs 平面直角座標系番号 -a_ullr 貼り付け先左上X座標 貼り付け先左上Y座標 貼り付け先右下X座標 貼り付け先右下Y座標 入力ファイル 出力ファイル
-a_srs EPSG:2451は、平面直角座標系で関東付近をカバーする9系を使用する、という意味になります。-a_ullrは入力画像を貼り付ける左上と右下の座標で、指定した平面直角座標系での原点からの位置（単位はメートル）で指定します。
ここで指定する平面直角座標系の番号はEPSGコードで、EPSGコードとはGISで使用される様々な（各国で独自に定義される）座標系などを統一した番号付けで扱うために作られたコード体系です。日本の平面直角座標系とそのEPSGコードとの関係について図に示します。画像を貼り付ける先の地点に応じて適切なコードを指定する必要があります。例えば関東なら平面直角座標系9はEPSG:2451であり、大阪なら平面直角座標系6=EPSG:2448です。
画像を貼り付けたい座標（左上座標）が平面直角座標系でどのような値になっているかは、この国土地理院の「平面直角座標への換算」ページを使って調べることができます。また、貼り付け先右下座標は、平面直角座標系では単純にX,Yの移動量を加減算するだけで良いので簡単に算出できます。（繰り返しになりますが）東に358m、南に334m移動した地点としたければ、左上座標の平面直角座標のX座標に358足し、同Y座標から344引きます。
画像ファイルができたら、内容を確認してみましょう。
root@9dd686378313:/work# gdalinfo sakura_mark_shinjuku_heimen.tif
Driver: GTiff/GeoTIFF
Files: sakura_mark_shinjuku_heimen.tif
Size is 358, 334
Coordinate System is:
PROJCS["JGD2000 / Japan Plane Rectangular CS IX",
GEOGCS["JGD2000",
DATUM["Japanese_Geodetic_Datum_2000",
SPHEROID["GRS 1980",6378137,298.2572221010002,
AUTHORITY["EPSG","7019"]],
TOWGS84[0,0,0,0,0,0,0],
AUTHORITY["EPSG","6612"]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433],
AUTHORITY["EPSG","4612"]],
PROJECTION["Transverse_Mercator"],
PARAMETER["latitude_of_origin",36],
PARAMETER["central_meridian",139.8333333333333],
PARAMETER["scale_factor",0.9999],
PARAMETER["false_easting",0],
PARAMETER["false_northing",0],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]],
AUTHORITY["EPSG","2451"]]
Origin = (-12469.062400000000707,-33897.301500000001397)
Pixel Size = (1.000000000000000,-1.000000000000000)
Metadata:
AREA_OR_POINT=Area
TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)
TIFFTAG_XRESOLUTION=72
TIFFTAG_YRESOLUTION=72
Image Structure Metadata:
COMPRESSION=LZW
INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left ( -12469.062, -33897.302) (139d41'44.00"E, 35d41'39.80"N)
Lower Left ( -12469.062, -34231.302) (139d41'44.02"E, 35d41'28.96"N)
Upper Right ( -12111.062, -33897.302) (139d41'58.24"E, 35d41'39.82"N)
Lower Right ( -12111.062, -34231.302) (139d41'58.26"E, 35d41'28.98"N)
Center ( -12290.062, -34064.302) (139d41'51.13"E, 35d41'34.39"N)
Band 1 Block=358x7 Type=Byte, ColorInterp=Red
Band 2 Block=358x7 Type=Byte, ColorInterp=Green
Band 3 Block=358x7 Type=Byte, ColorInterp=Blue
root@9dd686378313:/work#
Coordinate SystemとしてJGD2000 / Japan Plane Rectangular CS IXが使用されていることと、Corner Coordinatesに四隅の座標が平面直角座標系9系での原点からの相対位置（メートル）で入っていることが確認できると思います。
Google Earthプロに読み込ませてみると、緯度経度で指定したときと全く同じ位置に貼り付けられているのがわかります。
画像の一部分だけを抽出する方法
元画像に余計な部分が含まれる場合など画像の一部分だけ使いたいときは、gdal_translateコマンドで指定領域だけをクリッピングすることができます。
$ gdal_translate -of "GTiff" -co "COMPRESS=LZW" -srcwin 元ファイルクリッピング領域左上X座標 元ファイルクリッピング領域左上Y座標 クリッピング領域サイズX クリッピング領域サイズY -a_srs 平面直角座標系番号 -a_ullr 貼り付け先左上X座標 貼り付け先左上Y座標 貼り付け先右下X座標 貼り付け先右下Y座標 入力ファイル 出力ファイル
-srcwin 指定で、元ファイルのピクセル単位で領域を指定します。貼り付け先座標の指定は任意の座標系を使用することができます。
試しにさくらインターネットのマークの左上1/4だけクリッピングしてみます。
root@9dd686378313:/work# gdal_translate -of "GTiff" -co "COMPRESS=LZW" -srcwin 0 0 179 167 -a_srs EPSG:2451 -a_ullr -12469.0624 -33897.3015 -12290.0624 -34064.3015 sakura_mark.tiff sakura_mark_shinjuku_quarter1.tif
Input file size is 358, 334
0...10...20...30...40...50...60...70...80...90...100 - done.
root@9dd686378313:/work# gdalinfo sakura_mark_shinjuku_quarter1.tif
Driver: GTiff/GeoTIFF
Files: sakura_mark_shinjuku_quarter1.tif
Size is 179, 167
Coordinate System is:
PROJCS["JGD2000 / Japan Plane Rectangular CS IX",
GEOGCS["JGD2000",
DATUM["Japanese_Geodetic_Datum_2000",
SPHEROID["GRS 1980",6378137,298.2572221010002,
AUTHORITY["EPSG","7019"]],
TOWGS84[0,0,0,0,0,0,0],
AUTHORITY["EPSG","6612"]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433],
AUTHORITY["EPSG","4612"]],
PROJECTION["Transverse_Mercator"],
PARAMETER["latitude_of_origin",36],
PARAMETER["central_meridian",139.8333333333333],
PARAMETER["scale_factor",0.9999],
PARAMETER["false_easting",0],
PARAMETER["false_northing",0],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]],
AUTHORITY["EPSG","2451"]]
Origin = (-12469.062400000000707,-33897.301500000001397)
Pixel Size = (1.000000000000000,-1.000000000000000)
Metadata:
AREA_OR_POINT=Area
TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)
TIFFTAG_XRESOLUTION=72
TIFFTAG_YRESOLUTION=72
Image Structure Metadata:
COMPRESSION=LZW
INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left ( -12469.062, -33897.302) (139d41'44.00"E, 35d41'39.80"N)
Lower Left ( -12469.062, -34064.302) (139d41'44.01"E, 35d41'34.38"N)
Upper Right ( -12290.062, -33897.302) (139d41'51.12"E, 35d41'39.81"N)
Lower Right ( -12290.062, -34064.302) (139d41'51.13"E, 35d41'34.39"N)
Center ( -12379.562, -33980.802) (139d41'47.56"E, 35d41'37.09"N)
Band 1 Block=179x15 Type=Byte, ColorInterp=Red
Band 2 Block=179x15 Type=Byte, ColorInterp=Green
Band 3 Block=179x15 Type=Byte, ColorInterp=Blue
root@9dd686378313:/work#
こんな感じです。
この他、-srcwinの代わりに-projwinパラメータを指定することで、元ファイルのピクセル単位ではなく、座標系での指定でクリッピングすることも可能なようです（ここでは未試用のため割愛）。
画像を回転させて貼り付ける方法
貼り付け画像を地形や道路の碁盤目に合わせたい場合などで回転させたい場合があります。画像の回転（や移動）にはアフィン変換という行列演算を用いますが、原理はさておき画像を回転させるだけであれば以下の関係式を理解しておけば十分です。
r x cos(θ) + r y sin(θ) + ⊿x = ax + by + c
r x sin(θ) – r y cos(θ) + ⊿y = dx + ey + f
rはスケール、θは回転角、⊿x、⊿yはそれぞれXY方向の移動量です。それぞれのパラメータに操作したい希望の値を代入し、abcdefの各値を求めます。例えば、30度回転させて（100, 200）移動させたいとすると、cos(30°)=0.866025、sin(30°)=0.5ですので、
1 * x * 0.866025 + 1 * y * 0.5 + 100
1 * x * 0.5 – 1 * y * 0.866025 + 200
という式になり、a = 0.866025, b = 0.5, c = 100, d = 0.5, e = -0.866025, f = 200となります。eの値はマイナスですので注意してください。
回転など、この6パラメータを指定する場合はこれまでのコマンドライン引数でやる方法ではできず、ワールドファイルと呼ばれる別ファイルでパラメータを指定します。ワールドファイルは、貼り付けたい元画像と同じファイル名で拡張子を.tfwとしたもの（tiffの場合）で、テキストエディタなどで内容を簡単に作成することができます。
フォーマットは以下の通りで、abcdefの各パラメータをadbecfの順に改行で区切った単純なテキストファイルになります。
a d b e c f
実際にやってみた例は以下のようになります。30度回転した画像を新宿に貼る場合です。
root@9dd686378313:/work# cat sakura_mark.tfw 0.866025403784439 0.5 0.5 -0.866025403784439 -12445.080947277414581 -33980.8015 root@9dd686378313:/work# gdal_translate -of "GTiff" -co "COMPRESS=LZW" -srcwin 179 0 179 167 -a_srs EPSG:2451 sakura_mark.tiff sakura_mark_shinjuku_quarter2_kaiten1.tif Input file size is 358, 334 0...10...20...30...40...50...60...70...80...90...100 - done. root@9dd686378313:/work#
30度回転と貼り付け先場所の指定は.tfwファイルで、元画像の右上1/4クリップはコマンドライン引数で指定しています。
貼り付け結果です。
なおワールドファイル（.tfw）を使う方法では、回転させない場合でも貼り付け先左上座標（c,f）とスケール（a,e）を指定すれば（b,dは0）、貼り付け先右下座標の計算は自動でやってくれるので実はコマンドライン引数で指定するより楽でおすすめです。
複数画像を組み合わせる方法
作成したGeoTiff画像を複数組み合わせて一つの画像にするには、gdal_merge.pyというコマンドを利用します。
root@9dd686378313:/work# gdal_merge.py -co "COMPRESS=LZW" -o 出力ファイル名 入力ファイル1 入力ファイル2...
-o で出力先ファイル名を指定します。複数の入力ファイルで領域が重なる場合、後から指定したファイルの内容で上書きされます。これで困る場合は、事前にgdal_translateコマンドを使って位置を移動させておく必要があります。
gdal_translateコマンドで1/4づつ切り出した画像を組み合わせた例です。わざと隙間を開けてみたのですが、残念ながらデータが存在しないところを透過にすることはできないようです。-init 引数を指定することで初期値をRGB変えて色をRGB指定で変更することはできますが、透過（アルファチャンネル）の設定はできないようです。
まとめ
gdal_translateコマンドとgdal_merge.pyコマンドを使うことで、画像を切り取ったり、組み合わせて任意の場所に貼り付ける事ができるようになります。座標指定の意味がわかってくるまでは大変ですが、そこを超えてしまえばだいたい思ったとおりの操作ができるようになります。本メモがその一助になれば幸いです。
（おまけ）
全長70kmの宇宙船が東京の上空1000mに浮いている場合を想像して、画像を合成してみました。
艦首（？）部分を新宿に置いてみたのですが、こんなに大きいとは思いませんでした…。艦尾が九十九里浜に達しそうです。
ま、こんな楽しみ方もあるということで。