とらりもんHOME  Index  Search  Changes  Login

モザイク処理

筑波大学 奈佐原顕郎

モザイク処理

前回のようにGRASSを起動せよ。

地理情報が複数の図幅に分割されている場合、それらを統合して、1枚の大きな図幅にしたいことがよくある。そのような処理を、モザイク処理という。ここでは、GTOPO30の数値地形モデルについて、モザイク処理を行ってみよう。

まず、g.regionコマンドで、regionを、範囲は100E-160E, 20N-60N, 解像度は00:00:30(つまり0.00833333度)とする (日本列島を含む領域の、約1km解像度)。

   g.region w=100 e=160 s=20 n=60 res=0.00833333

次に、GTOPO30のサイトから、 E100N40、 E100N90、 E140N40、 E140N90の4枚の図幅をダウンロードして解凍/展開し、それぞれをインポート。

wget -c ftp://edcftp.cr.usgs.gov/data/gtopo30/global/e100n40.tar.gz
wget -c ftp://edcftp.cr.usgs.gov/data/gtopo30/global/e100n90.tar.gz
wget -c ftp://edcftp.cr.usgs.gov/data/gtopo30/global/e140n40.tar.gz
wget -c ftp://edcftp.cr.usgs.gov/data/gtopo30/global/e140n90.tar.gz
tar zxvf e100n40.tar.gz
tar zxvf e100n90.tar.gz
tar zxvf e140n40.tar.gz
tar zxvf e140n90.tar.gz
r.in.gdal input=E100N40.DEM output=GTOPO_E100N40
r.in.gdal input=E100N90.DEM output=GTOPO_E100N90
r.in.gdal input=E140N40.DEM output=GTOPO_E140N40
r.in.gdal input=E140N90.DEM output=GTOPO_E140N90
r.null map=GTOPO_E100N40 setnull=55537
r.null map=GTOPO_E100N90 setnull=55537
r.null map=GTOPO_E140N40 setnull=55537
r.null map=GTOPO_E140N90 setnull=55537
r.colors map=GTOPO_E100N40 rules=srtm
r.colors map=GTOPO_E100N90 rules=srtm
r.colors map=GTOPO_E140N40 rules=srtm
r.colors map=GTOPO_E140N90 rules=srtm(色づかいを変更し、標高値に応じたグラデーションにする。)

とりあえず、4枚を一緒に表示するだけなら、d.rastコマンドに-oオプション(overlay)をつけることで可能:(追記:GRASS7の場合-oは不要)

d.mon x0
d.erase
d.rast GTOPO_E100N40
d.rast GTOPO_E100N90 -o
d.rast GTOPO_E140N40 -o
d.rast GTOPO_E140N90 -o
GTOPO_4.png

しかしこれでは、何らかの画像処理や演算(傾斜や日射量の算出など)をするのに、いちいち3つの図幅を個別に処理せねばならないし、継目の部分ではうまく処理ができない可能性がある。そこで、この4枚を1つの図幅に統合しよう(モザイク)。それには4枚の図幅に対して演算処理を施す必要があるが、NULLが入っているとまずいので、いったん、NULLを-1のような値に変換しよう: (追記:GRASS7の場合=の前後にスペースが必要です。)

r.mapcalc "dummy1=if(isnull(GTOPO_E100N40),-1,GTOPO_E100N40)"
r.mapcalc "dummy2=if(isnull(GTOPO_E100N90),-1,GTOPO_E100N90)"
r.mapcalc "dummy3=if(isnull(GTOPO_E140N40),-1,GTOPO_E140N40)"
r.mapcalc "dummy4=if(isnull(GTOPO_E140N90),-1,GTOPO_E140N90)"

※"dummy1 = if"のようにスペースを開けないと、上手く動かない可能性があります。(追記:2018/02/08 菊島)

その上で、4枚の図幅から、各ピクセルの最大標高値を抜き出すという演算を施せば、いずれの場所でも有効な標高値となるはずである。

r.mapcalc "GTOPO_Japan=max(dummy1,dummy2,dummy3,dummy4)"

最後に、-1としたNULLを元に戻そう:

r.null map=GTOPO_Japan setnull=-1

ついでにカラーテーブルもいいかんじにする:

r.colors map=GTOPO_Japan color=srtm

※"r.colors map=GTOPO_Japan rules=srtm"(追記:2018/02/08 菊島)

すると、新しい図幅として、GTOPO_Japanというラスターデータが生成される。これが4枚をあわせたモザイクである。

d.erase
d.rast GTOPO_Japan
GRASS_mosaic_Japan.jpg

最後に, ファイルを消去しよう。大学の計算機端末などではディスクがすぐに一杯になってしまうので。

rm E??N??.???
rm e??n??.tar.gz

シェルスクリプトにまとめる

たくさんのコマンドを打ち込んだが、これらはもちろん、以下のようにシェルスクリプトとしてまとめることができる。だれか1人がこのようなシェルスクリプトを作って配布すれば、他の人はこの処理をコマンドひとつでやってしまうことができるのである。

#!/bin/sh
g.region w=100 e=160 s=20 n=60 res=0.00833333
for i in e100n40 e100n90 e140n40 e140n90; do
    wget -c ftp://edcftp.cr.usgs.gov/data/gtopo30/global/$i.tar.gz
    tar zxvf $i.tar.gz
done
for i in E100N40 E100N90 E140N40 E140N90; do
    r.in.gdal input=$i.DEM output=GTOPO_${i}
    r.null map=GTOPO_${i} setnull=55537
    r.mapcalc "dummy${i}=if(isnull(GTOPO_${i}),-1,GTOPO_${i})"
done
r.mapcalc "dummy=max(dummyE100N40,dummyE100N90,dummyE140N40,dummyE140N90)"	
r.mapcalc "GTOPO_Japan=if(dummy==-1,null(),dummy)"
r.colors map=GTOPO_Japan color=srtm
d.rast GTOPO_Japan
for i in E100N40 E100N90 E140N40 E140N90; do
    g.remove rast=dummy${i}
done
rm E??N??.???
rm e??n??.tar.gz

実習課題

上でモザイクしたDEMを用いて、関東平野を対象にして、前回やったように、標高12m未満の地域を抽出して表示してみよ。

GIS入門に戻る。

Last modified:2018/02/08 19:44:17
Keyword(s):
References:[GIS入門]