Rendered imageなGeoTIFF画像を色々なモードで合成する

##動機
毎回QGISで複数枚のラスターを合成させていると読み出しの動作が遅いので、予めRenderedなGeoTIFFを使うようにしようと思い立った。しかしQGISはプリントコンポーザを介した方法でしか複数枚のラスターを描画して出力する方法を提供していない(多分)。
そこで**pythonのgdalモジュールを使ったスクリプトによって、各レイヤーごとに出力したRendered imageを合成(blend)することにした**。状況によって使いたい合成方法が変わるので、よく使う乗算(multiply), オーバーレイ(overlay), スクリーン(screen), ソフトライ(Soft light), ハードライト(Hard light)など実装した。

##用いるデータ
~~本当にやりたいことは別にあるのだが今回はサンプルデータとして、~~ JAXAの[ALOS World 3D 30m DSM](http://www.eorc.jaxa.jp/ALOS/en/aw3d30/)を標高データとして使う。
このDSM、そのままでは解像度が高すぎるので25%にリサイズして利用している。
![Fig1.jpeg](https://qiita-image-store.s3.amazonaws.com/0/172235/996e1941-e64f-8835-ec18-85124b1dd4f9.jpeg)

今回はTopoUSM図と段彩図(Elevation tints)を、それぞれLayer1.tifとLayer2.tifとして合成することにする。
【宣伝】TopoUSMに関しては： [DEM可視化ツール TopoUSM - Qiita](http://qiita.com/geoign/items/1126591f2f4de8075c44) をどうぞ。

まず取り敢えずどちらのレイヤーも、QGIS経由でRendered imageとして保存しておく。ラジオボタンがRendered imageになっていることに注意。
<img width="616" alt="fig2.jpg" src="https://qiita-image-store.s3.amazonaws.com/0/172235/53ae2476-5655-dff6-38f4-80e5bb51bfd1.jpeg">

##合成処理
先述の通り、画像の合成法には幾つもの種類があり得手不得手がある。
具体的な演算については[GIMPのdocumentation](https://docs.gimp.org/2.6/ja/gimp-concepts-layer-modes.html)に良い記事があったので、それを参考にした。
Layer1.tifはそのままに、globで探したLayer2~以降を重ねていくという実装にした。

```py3:GTiff_RenderStack.py
import glob, gdal
infiles = glob.glob('Layer?.tif')
mode = 'multiply'
outfile = 'tmp.tif'

ds0 = gdal.Open(infiles[0])
R,G,B = [ds0.GetRasterBand(i).ReadAsArray() for i in range(1,4)]
for infile in infiles[1:]:
    ds = gdal.Open(infile)
    R1,G1,B1 = [ds.GetRasterBand(i).ReadAsArray() for i in range(1,4)]
    print('%s:\t%s\t%s\t%s' % (infile, R1.shape, G1.shape, B1.shape))
    if mode=='multiply':
        R,G,B = R*(R1/255), G*(G1/255), B*(B1/255)
    elif mode=='overlay':
        R = R/255 * (R+R1/255*2*(255-R))
        G = G/255 * (G+G1/255*2*(255-G))
        B = B/255 * (B+B1/255*2*(255-B))
    elif mode=='screen':
        R = 255 - (255-R)/255*(255-R1)
        G = 255 - (255-G)/255*(255-G1)
        B = 255 - (255-B)/255*(255-B1)
    elif mode=='hardlight':
        R = (255 - (255-2*(R1-128))/256*(255-R))*(R1>128) + R1/128*R*(R1<=128)
        G = (255 - (255-2*(G1-128))/256*(255-G))*(G1>128) + G1/128*G*(G1<=128)
        B = (255 - (255-2*(B1-128))/256*(255-B))*(B1>128) + B1/128*B*(B1<=128)
    elif mode=='softlight':
        R = ((255 - (255-R)/255*(255-R1) + (255-R)/255*R1)*R)/255
        G = ((255 - (255-G)/255*(255-G1) + (255-G)/255*G1)*G)/255
        B = ((255 - (255-B)/255*(255-B1) + (255-B)/255*B1)*B)/255
R,G,B = R.astype(np.byte), G.astype(np.byte), B.astype(np.byte)

ds1 = gdal.GetDriverByName('GTiff').CreateCopy(outfile, ds0)
ds1.GetRasterBand(1).WriteArray(R)
ds1.GetRasterBand(2).WriteArray(G)
ds1.GetRasterBand(3).WriteArray(B)
ds1 = None
```

##出力結果
当然だがQGISで同じ方法で合成して表示させたのと同じ色になる。めでたしめでたし。
![overlay-hardlight2.jpg](https://qiita-image-store.s3.amazonaws.com/0/172235/a9776a04-65e5-9bc6-0fa2-7c9f8ebb50e7.jpeg)
![screen-multiply.jpg](https://qiita-image-store.s3.amazonaws.com/0/172235/30c4e70c-20d7-f4db-5d16-0853007a867b.jpeg)
![softlight.jpg](https://qiita-image-store.s3.amazonaws.com/0/172235/7a815bc8-ffe2-4968-dda2-e245bf5501ee.jpeg)

今回背景に使うTopoUSMはベースがグレー(階調が128/255)なので、オーバーレイやハードライトは鮮やかに色が着く。一方スクリーンは明るくなりすぎ、乗算では暗くなりすぎる。そして割りとお気に入りのソフトライトでは、ハードライトよりやや落ち着いた感じの色乗りになる。

因みにもし例えばTopoUSMではなくベースが白(階調が255/255)になるような傾斜図に重ねる場合、やはり乗算が鉄板になる。また暗い図ではスクリーンが良い感じになることもある。