Contour Map on SAGA GIS

SAGA GIS等高線図を作成しました。いつものASTER GDEM¹を使います。

SAGA GISでのコンターマップの作成

等高線を引くためのモジュールは、Shapes – Grid / Contour Lines from Gridです。

設定項目は、

• 標高データのグリッド
• 最少標高
• 最高標高

を指定します。最高標高と最低標高の指定は、特定の高さの範囲の等高線を引きたいときや、海底地形を含むグリッドデータで、海底のみのコンターを引きたいときとかに便利です。今回はASTER GDEMの画像一枚すべてを作成してみました。
できました。

ズームしてみるとこんな感じです。

国土地理院のウォッちずで、正解の地図を確認しようとしたところ、

平素より、地図閲覧サービス（ウォッちず）をご利用頂きありがとうございます。
ウォッちずは平成26年4月1日以降地図を更新しておらず、平成26年10月1日をめどにサイトを停止する予定です。
なお、最新の地図は「地理院地図」でご覧いただけます。

とのことです。
地理院地図は、ウォッちずで提供されていた1:25,000の地図のほか、基準点や土壌図、断層、火山情報等のさまざまな情報をオーバーレイできるように整備したサイトのようです。しかも、webサイトへの埋め込み機能があるようです。
上の地図は北岳周辺なので、地理院地図と比較してみました。
1:25,000の地図はさすがに地表の細かい崖まで書かれていますが、20万分の1程度の地図であれば、そこそこの再現率ではないでしょうか。

---

1 ASTER GDEMはMETI/NASAに帰属します。^↩

SAGA GISで陰影図の作成

DEMから陰影図を作成します。DEMはASTER GDEM^*1を利用しました。日本で一番ダウンロードされているだろうと思われる、富士山を含む地域のDEMです。

DEMの読み込み

まずはDEMをSAGA GISに読み込みます。

workspaceのDataタブにDEMデータをドラッグアンドドロップします。ASTEER GDEMをダウンロードした際に同梱されるnumファイルは、標高データを作成するために何枚の衛星画像を用いたかや、欠損値や異常値の補完に用いた画像に関するデータが書かれているため、今回は利用しません。
データツリーに追加されたDEMファイルをダブルクリックすると、ワークスペースのMapsタブにDEMが追加され、地図が表示されます。

陰影図の作成

読み込んだDEMにTerrarian Analysis – Lightning, Visibility 内のAnalytical Hillshadingモジュールを利用します。

入力としてDEMを指定す他はデフォルトのパラメータで、実行します。(上記画像参照)。
ワークスペースのデータタブにAnalytical Hillshadingのグリッドデータが追加されるため、ダブルクリックしてDEMと同じ地図に追加します。

色づけ

色づけしたいレイヤをMapsタブで選択します。まずはDEMから作業をするので、Analytical Hillshadingを、ダブルクリックして非表示にします。非表示レイヤはmapsタブ内では[ ]でくくられます。
DEMのレイヤを選択すると、右側のウィンドウ(Object Propery Window)のsettingタブでDEMの色を設定できるようになります。colorsをクリックして出てくる色設定画面から、Presetを選択、defoltを指定します。

Objet Property Windowに戻り、Applyを押すと、DEMに色が付きます。

Analytical Hillshadingを再度ダブルクリックして表示させ、colorsのTypeでShadeを選択し、Applyします。

色合いはあんまりよくないけれど、完成です。
*1:ASTER GDEMはMETI/NASAに帰属する

R コマンドメモ

Rで使った、今後もちょくちょく使いそうなコマンドをメモ。自分で使うときに見返すためのメモのため、使い方の中のコマンドには、事前にxやiに適当な変数が代入されていたり、for loop内での使用や直前にグラフの描写がされていることをを想定していたりします。

カテゴリ	説明	コマンド	使い方
作業準備	ワークスペースのオブジェクトをリストに出力	ls()	ls()
	ワークスペース内のオブジェクトを削除	rm(list=ls())	rm(list=ls())
	グラフを消す	graphics.off()	graphics.off()
	作業ディレクトリの表示	getwd()	getwd()
	作業ディレクトリの変更	setwd()	getwd()
データの読み込み	テキストファイルのデータを読み込む	read.table()	read.table("data\\rain.tsv", row.names=1, header=TRUE)
テキスト操作	変数に入っている文字列、変数等をつなぎ合わせる	paste()	paste("Type", x, sep= "-")
	各列に同じ操作を繰りかえし行う.	apply()	apply(x, 2, sd) または apply(x, 2, function(i){         runs.test(as.factor(i < mean(i)))     } )
オブジェクトの確認	代入と同時に値を表示	()	(x <- 3 +5)
	データの最初の部分を確認	head()	head(dataframe)
分析	主成分分析	prcomp()	pc.df <- prcomp(df, scale=TRUE)
	主成分分析結果の表示	summary()	summary(pc.df)
	クラスター分析	hclust()	x <- hclust(x, method="ward") plot(x)
	クラスター分析のデンドログラムを資格で囲む	rect.hclust()	rect.hclust(x, k=i, border=rainbow(i))
グラフ描写	描写画面の分割	par()	par(mfcol=c(2,2))
	レーダーチャートの描写	radarchart()	radarchart(rbind(max=a, min=b), x, maxmin=TRUE)
出力	画像としてグラフを保存	dev.copy()	dev.copy(png,file=paste("out/", i, ".png", sep="")) dev.off()
	テキストファイルに出力	print() sink()	sink(paste("out/", dir, "Summary_", case,".txt", sep="")) print(summary(x)) sink()
	テキストファイルにデータフレームの出力	write.table()	write.table(x, paste("out/", i,".txt", sep=""), append=FALSE, quote=FALSE,row.names=FALSE, sep="\t")
その他	ヘルプを表示	help()	help()

複数のデータを一つのグラフに重ねて描写する場合、軸は繰り返し描写されないよう、最後のデータ描写時に描写する。

for (i in 1:n-1){
    plot(x[i], type="l", col=rainbow(n)[i], ylim=c(a, b), xlab="", ylab="", axes=FALSE)
    par(new="TRUE")
}
plot(x[n], type="l", col=rainbow(n)[n], ylim=c(a, b), xlab="day", ylab="price")

解析を行う場合、その解析が適応できるデータ分布となっているか事前に判定すること。(正規分布を仮定している解析ではないか、独立標本を仮定していないか、など)

SAGA GIS

QGISをインストールしたら、SAGA GISも同時にインストールされていました。以前はなかった気がするので、スタンドアロンのインストーラではなく、ネットワークインストーラを利用したからでしょうか。

大学ではGRASSを教えてもらったのですが、GRASSはマップセットの設定とかがうまく理解できていないため、GUIで何とかなりそうなSAGA GISでまずはいろいろやってみようかと思います。

QGISはUIが多言語化されているし、書籍やネット上の情報もほかのフリーのGISソフトと比べて入手しやすいですが、地形解析のプラグインや3D viewの機能があまりよくありません。プラグインは2.0になってから対応するものが減ったのでしょうか。

QGISのメニューの中に、GRASSのマップセットを開くメニューがあるので、QGISの中でGRASSと連携して解析できるのだろうけれど、ひとまずこちらで。

とりあえず、以下を目標に調べてみます。

• 鳥瞰図の作成
• 水系図の作成
• 特定の地点でのCatchment Areaの計算
• 降水量を与えた時の流量の推定

河川流量は、どういった流出モデルやパラメータを使うかによって全然違うだろうけれど、どうせ正確な雨量データなんて入手できないので、桁が違わなければいいかなと考えています。