🌏

Geotiffを回転させる

2021/12/02に公開

gis

Gazebo

geotiff

tech

何でそんなことしたいの？

Geotiffというのはジオリファレンシングされたtiffなので、それを回転させる事などにさしたる意味はなく、本来はそんなことしなくてもいいはずです。

ところが、GazeboというソフトウェアではDEMの入ったGeotiffを読み込ませると、Geotiffの南北方向とGazebo表示上の南北方向が一致しません（私はGazeboを使うソフトウェアの内部までいじっておらず、表層の部分で目的を達成しようとしているため、この辺の理解はあやふやです）。具体的にはGeotiffの北方向がGazebo表示上ではY軸方向となるのに対して、Gazebo上（周辺ライブラリ？）の北方向はGazebo表示上のX軸方向となります。要は、90度ずれているわけです。

ついでに言うと、GazeboはGeotiffをDEMデータとして読み込ませようとすると正方形として扱って読もうとするため、格納されているDEMデータのサイズが正方形になっていないと中心がずれて表示されたりします。しかし、Geotiffの位置情報を使っていないのかと言えばそうでは無いようで、ジオリファレンシングされていないtiffをDEMデータとして読もうとしてもエラーになります。

Gazebo上でGeotiffを正しく（ズレなく正確なスケールで）表示しようとするととても辛いものがあります（簡単な方法知っている方は是非ご教授願いたい・・・）。

脱線しましたが、目的はGazeboで正しくDEMデータの入ったGeotiffを表示するために、Geotiffを回転させたいわけです。

DEM / Ground model rotation - Gazebo: Q&A Forum

Geotiffとは何なのだろうか？？？

Geotiffというのはジオリファレンシングされたtiffですが、tiffに格納された画像のピクセル毎に位置情報が記録されているわけではありません。tiffは画像データの格納形式として様々なタグを設定できますが、そのタグを使用して画像の4頂点+中心のピクセルについて特定の座標系においての座標値を記録しているのがGeotiffです（多分）。

したがって、Geotiffを表示するソフトウェアではGeotiffのタグから読み取った4頂点+中心ピクセルに対応する位置座標とその座標系から、例えば表示したい地図における座標系の座標値に変換して、特定のピクセルの座標値を既知の5点の座標値から補完して計算した上で地図上にオーバーレイさせるわけです。

つまりは、Geotiffに格納されているデータとは画像とジオリファレンス情報の2つで、この2つは特に密接に関連しているわけではありません。Geotiff処理においては、画像の部分は画像（あるいはそれに類するもの）として処理した上で、その処理結果に対して座標値を適切に指定してやれば良いはずです。

（ここに書いてあることは憶測です、正確なソースを参照したものではありません）

方法

GIMPなどで画像（DEMデータ）の部分を回転させる
- この時、ジオリファレンス情報は失われる
できたtiffに対してgdal_translateによって元のGeotiffの座標値をGCPとして復元する
GCPを与えたtiffに対してgdalwarpによって座標系を指定し幾何補正してGeotiffを得る
gdal_translateによって短辺に合わせてリサイズし正方形にする

GDALを使用するので予めインストールしておきましょう。

1. 画像（DEMデータ）の部分を回転させる

DEM in Geotiffを開いて編集することのできる画像処理ソフトウェアにて、Geotiffの画像部分を回転します。これによってまず、回転させたDEMデータを得ることができます。ただし、多くの場合これをtiffとして出力するとジオリファレンス情報は失われます。

DEMデータは1ピクセルあたり32bit浮動小数点数1つからなるデータなので、正しく開いて情報を保持したまま編集できるソフトはあまりないかもしれません。私はGIMP(2.10)を使用して時計回りに90度回転させました（ただし、時計回りに90度回転をそのままやると回転せず、180度回転させてようやく意図通りになるという謎の現象に出会いました。これがバグなのか仕様なのかはよくわかりません・・・）。

2. gdal_translateによって元のGeotiffの座標値をGCPとして復元する

回転画像を保存する際にGeotiffのGeoの部分が失われているはずなので、回転を考慮した上で元の情報を復帰させます。

まずgdalinfoによって元のGeotiffから情報を取得します。

gdalinfo input.tif

これを実行すると例えば次のように情報が得られます。

Driver: GTiff/GeoTIFF
Files: input.tif
Size is 161, 161
Coordinate System is:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4326]]
Data axis to CRS axis mapping: 2,1
Origin = (137.194833332999991,35.183111111000002)
Pixel Size = (0.000067632850932,-0.000055555552795)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  ( 137.xxxxxxx,  35.yyyyyyy) (137dxx'xx.xx"E, 35dyy'yy.yy"N)
Lower Left  ( 137.xxxxxxx,  35.wwwwwww) (137dxx'xx.xx"E, 35dww'ww.ww"N)
Upper Right ( 137.zzzzzzz,  35.yyyyyyy) (137dzz'zz.zz"E, 35dyy'yy.yy"N)
Lower Right ( 137.zzzzzzz,  35.wwwwwww) (137dzz'zz.zz"E, 35dww'ww.ww"N)
Center      ( 137.aaaaaaa,  35.bbbbbbb) (137daa'aa.aa"E, 35dbb'bb.bb"N)
Band 1 Block=161x12 Type=Float32, ColorInterp=Gray
  Min=122.170 Max=183.190 
  Minimum=122.170, Maximum=183.190, Mean=153.122, StdDev=12.462
  NoData Value=-9999
  Metadata:
    STATISTICS_MAXIMUM=183.19000244141
    STATISTICS_MEAN=153.12228540338
    STATISTICS_MINIMUM=122.16999816895
    STATISTICS_STDDEV=12.461520655668
    STATISTICS_VALID_PERCENT=100

なんか色々出てますが見るべきは次の2つです

Coordinate System isのID
- この例だとID["EPSG",4326]
Corner Coordinatesの5点

このステップで使用するのはCorner Coordinatesの5点の座標値です。上から左上・左下・右上・右下・中央のピクセルに対する座標値を指定しています。この座標値を90度回転させた上で回転画像に対して指定してやります。

指定はgdal_translateによって次のように行います

gdal_translate -gcp 0 0 137.xxxxxxx 35.wwwwwww -gcp 161 0 137.xxxxxxx 35.yyyyyyy -gcp 161 161 137.zzzzzzz 35.yyyyyyy -gcp 0 161 137.zzzzzzz 35.wwwwwww -gcp 80.5 80.5 137.aaaaaaa 35.bbbbbbb input.tif output.tif

オプションの-gcp X Y latitude longitudeによって一つのピクセル位置に対する座標を指定しており、X Yがピクセル座標値、latitude longitudeは緯度経度をそれぞれ指定します。これを画像4隅+画像中心の5点指定します。

この例では入力Geotiffは正方形になっているので中心ピクセルは何も考えなくても大丈夫ですが、長方形画像（普通はそうなるはず）に対しては中心ピクセル座標も回転させなければならないことに注意です（入れ替えるだけですが）。

3. gdalwarpによって座標系を指定し幾何補正する

前段で得られたtiffはまだGeotiffとなってはいません。地理座標をGCPとして指定はしても楕円体や測地系を指定しておらず、ジオリファレンスは完全ではありません。

次は、gdalwarpによって座標系を復元します。

gdalwarp -s_srs EPSG:4326 -r cubic input.tif output.tif

-s_srsオプションに元画像のCoordinate System isのIDからEPSGコードを拾って指定します。これを実行するとtiffがジオリファレンスされ、その座標系と位置情報に基づいて画像がリサイズされます。-r cubicはリサイズに伴う補間方法を指定します。

4. gdal_translateによってリサイズし正方形にする

正方形のGeotiffが欲しいのでなければこのステップは必要ないのですが、今回は画像サイズが正方形になっていて欲しいのでgdal_translateによって短辺に合わせてリサイズしてやります。

gdal_translate -outsize min_len min_len -r cubic input.tif output.tif

min_lenは前段で得られたGeotiffの画像サイズの長辺と短辺の短い方を指定します（別にどっちでもいいですが）。-r cubicは先程同様補間方法です。前段及びこのリサイズ時に適切に補間方法を指定しないとデフォルトのバイリニア補間が行われ、出力にブロック状のノイズが出ることがあります。

ここで得られたtiffが適切？に回転されたはずのGeotiffになります。色々怪しい部分があるので本当に適切かはわかりませんが、一応目的を達成することはできたのでこれで良いことにします・・・