動機
毎回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~以降を重ねていくという実装にした。
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)になるような傾斜図に重ねる場合、やはり乗算が鉄板になる。また暗い図ではスクリーンが良い感じになることもある。