動機

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

用いるデータ

~~本当にやりたいことは別にあるのだが今回はサンプルデータとして、~~ JAXAのALOS World 3D 30m DSMを標高データとして使う。
このDSM、そのままでは解像度が高すぎるので25%にリサイズして利用している。

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

まず取り敢えずどちらのレイヤーも、QGIS経由でRendered imageとして保存しておく。ラジオボタンがRendered imageになっていることに注意。

合成処理

先述の通り、画像の合成法には幾つもの種類があり得手不得手がある。
具体的な演算についてはGIMPのdocumentationに良い記事があったので、それを参考にした。
Layer1.tifはそのままに、globで探したLayer2~以降を重ねていくという実装にした。

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で同じ方法で合成して表示させたのと同じ色になる。めでたしめでたし。

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

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