# つくばもくもくかい+CS立体図

## つくばFOSS4Gもくもく会・CS立体図勉強会

日時：2017年1月30日、13時から
場所：つくば市吾妻交流センター

まずはマニュアルをみながらQGISでひとさらいしています.

### １．目的
QGISで、CS立体図を作れるようにする。

いまやっと開始！13:20 -
小橋、到着。18:00 -

### ２．兵庫県さんの作ったマニュアルの確認

- 印刷物を3部、持参します。事前に渡せるなら、渡した方がいいかしら？→三島さんに渡した

作業環境
朝日：Win10 64bit OSGeo4W 2.18.2→問題なくできた
三島：Win8.1 64bit QGIS 2.16.3 → カラーが反転しない。2.18に変更。できた。
勝部：Win7 64bit QGIS 2.8 → SAGAが動かない。2.14.11に変更。できた。
岩崎：
OSGeoLive 10 QGIS 2.14.3 → Processing からSAGAの機能を呼び出すところでエラー。SAGA本体をつかってGausian filterと、General Curvatureを計算する。→エキスポートして、QGIS上で表示。
MacOS X El Capitan QGIS 2.18→問題なく実行できた

その他注意事項

- 整数値DEMだと、おかしくなる
- ５ｍ→６ｍといった感じで、標高がデジタルにかわるため
- DEM自体がなめらかである必要がある。
- Tiffのファイルサイズを小さくするために、整数型にしてあるファイルがあります。このようなファイルでCSを作ると、等高線の様な線が入ります。小数点以下の値を持ったDEMを使用する。（戸田）
- 曲率処理は10MDEMでは必要ない　ガウジアンではなく単純平均(近傍）でスムージングさせる程度
- レーザープロファイラから作成されたDEM等を使う場合はガウジアンを使って前処理
- Curvatureの計算は、degreeがいいのか、radianがいいのか。検討。

　どちらでも可。１ラジアン＝1（80／π）°のため、比例定数がかかるだけ。GISでは表示の際の最大、最小値のレンジが変わるだけなので、見た目には影響しない。（戸田）

- 平らなところを強調したい場合→色を変えてやるといい。
- トライアンドエラーでやるといい。
- 表示のレンジを狭くすれば、平らな場所でも微地形が強調されます。（戸田）
- 細かいDEMの場合、フィルターをかけたほうがいい。
- 重み付け関数を使ってフィルターをかけている。単純な移動平均ではよくない。
- 単純な移動平均だと、ゴーストの様な線が入ります。（戸田）
- Gaussian Filterのときの12とはなに？
- というか、Gaussian Filterの威力が重要。
- 三大地形量（標高、曲率、傾斜）に、色を付けただけ。が、それが重要。
- 12とは検索半径です。QGISの場合はメートルかピクセルかは不明。１ｍメッシュDEMの場合は１２ｍ。実際の地形にして２４ｍの地形規模（半径１２ｍ）を強調したいので、１２と入力。地すべりなど、より大きな地形を判読したい場合は、もっと大きな数字にする。（戸田）
- 曲率と平滑化が何を意味しているのかを、どう書けばいいのか悩んでいる。誰か書いて（戸田）。

### ３．グラフィックモデラーへの実装→４．Pythonでの実装

- モデラーだと配色やレイヤの順序の指定ができないので朝日さんがpyQGIS使って書いている

ベースの色が茶色

もともと谷に青で尾根が茶色という配色
谷を強調するのが最初(湧水地点を検出するため)

- その他、項目があれば、追加してください。

### 5. こけたところ

- QGIS 2.16.3 でカラーマップの反転がなぜかできなかった(win8.1 スタンドアロン版)
- 上にも書いたが、OSGeoLive10ではSAGAをProcessingから読み出せないので、SAGAで直接計算した。
- QGIS 2.14の問題？

### 6．質問

- 事前に質問があれば、こちらにどうぞ！
- 喜多さんから。「作り方の注意点や、アドバイスなどがあれば」
- 静岡のCS立体図はあまり立体に見えません（個人的見解）なんかテカテカしているのは何故？　喜多
- 長野県は、できあがったものに、コントラストと明るさの調整をしている。静岡はしてないので、てらてらしてるのではないか？（戸田）
- お勧めは、コントラスト＝+20、明るさ＝-10。モニターやプリンターによって見え方が違うので、適宜調整してください。（戸田）
- ありがとうございます。コントラストの調整が必要なのですね。（喜多）
- 戸田さんに追加質問：戸田さんの方で全国のCS立体図を作ったと聞きましたが、いつ公開されますか？　（喜多）
- 今実施しているプロジェクトの知財の制約などがあるため、4月以降にはオープンにできます。ESRI、G空間情報センター、Avenzaで公開できるように準備を進めています。公開になりましたら、またご連絡します。（戸田）

- CS立体図で平地の地形、例えば台地とか、谷津とかを表現するには、どの辺を工夫すればいいのか？（いわさき）
- しつもん：近傍平均をとるときっていつもgrassつかっているんだけど QGISの他のツールでさっくりできないかしら
- 曲率処理平滑化は必須ですか？たとえば、35.123mなどmmまである標高データで点間隔が3mmくらいの場合です（山口より）。
- 平滑化ではないのか？平滑化は、必要である。平滑化をどのようにするかによって、強調したい地形が見えてくる。（岩崎代筆）
- ありがとうございます。平滑化です。強調したいのは周囲との差が数ミリ～数センチくらいの凸凹です。なるほどトライアンドエラーですね（山口）。
- はじめまして。戸田です。おっしゃる通りです。私も平滑化についてはトライアンドエラーで自分が見たい地形規模が強調されるフィルターを作りました。
- QGISではガウジアンフィルターの標準偏差と検索範囲を調整すればよいと思います。お試しください。（戸田）
- ありがとうございます。やってみます。（山口）たくさんヒントをありがとうございます。
- 質問：DEMの平滑化をQGISで行うには？（喜多）
- SAGA-RasterFilter-GaussianFilterで行うようです。（戸田）
- 下記の処理（1mDEM）では兵庫県マニュアルにならってStandard Deviation=3, Search Mode=Circle, Search Radius=12で処理しています。（勝部）

- 質問：10mメッシュDEMでも平滑化処理は必要ですか？現状でも自分では特に気にならないです。（喜多）
- 地理院の１０ｍメッシュDEMは、作成時に既にある程度平滑化されているので処理しなくても大丈夫ですが、私は半径３０ｍから５０ｍの移動平均を行います。（戸田）
- 戸田さん、ありがとうございました。大変勉強になりました。平滑化にも挑戦してみます。（喜多）

### 参考

- sagaのソース(傾斜、曲率)

Do_MaximumSlope()のとこでいいのかな？
http://saga.sourcearchive.com/documentation/2.0.7plus-pdfsg-2build2/Morphometry_8cpp_source.html

- esriのcurvature

http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//00q90000000t000000

- gaussian filter

https://sourceforge.net/p/saga-gis/wiki/grid_filter_1/
http://docs.qgis.org/2.6/ja/docs/user_manual/processing_algs/saga/grid_filter/gaussianfilter.html
standard diviation : 3
search radius : 12

- QGISでの実装（簡易版？） by 喜多さん
- 簡易版というか、立体図のみ作った。
- 実はただしいかよくわかってないけど、なんとなく立体に見えるしいいかと思って作ってみた　喜多
- http://koutochas.seesaa.net/article/444171690.html
- leafletでの実装（簡易版？） by 松澤さん
- https://github.com/frogcat/csmap
- http://qiita.com/frogcat/items/b9aaa3df866a89d46ef5
- Pythonでの実装（簡易版？） by 和山さん
- https://github.com/makinux/PyCsmap

==========
勝手にコード書いてます、pythonコンソールから入力ファイル決め打ちで計算まで
ToDo:

ソースの整理はあと
一番上のslopeレイヤの反転と透過ができてない
ランクの数字がでてない
processingのスクリプトにまとめる or プラグインにする

from PyQt4.QtCore import *
from PyQt4.QtGui import *
from osgeo import gdal
import processing

def set_layer_style(layer, color_name, rank, reverse, min=None, max=None):
if min is None or max is None:
data = gdal.Open(layer.dataProvider().dataSourceUri())
band = data.GetRasterBand(1)
(min,max) = band.ComputeRasterMinMax(1)

sa = (max - min)/rank

if color_name == 'WhiteToBlack':
lst = [ QgsColorRampShader.ColorRampItem(min, QColor(255,255,255), str(round(min,3))), QgsColorRampShader.ColorRampItem(max, QColor(0,0,0), str(round(max,3))) ]
else :
lst = []
for i in range(rank):
idx = rank-1-i if reverse else i
lst.append(QgsColorRampShader.ColorRampItem(min+sa*i, QgsColorBrewerPalette.listSchemeColors(color_name, rank)[idx], str(round(min+sa*i,3))))

myRasterShader = QgsRasterShader()
myColorRamp = QgsColorRampShader()

myColorRamp.setColorRampItemList(lst)
myColorRamp.setColorRampType(QgsColorRampShader.INTERPOLATED)
myRasterShader.setRasterShaderFunction(myColorRamp)

myPseudoRenderer = QgsSingleBandPseudoColorRenderer(layer.dataProvider(), layer.type(), myRasterShader)

layer.setRenderer(myPseudoRenderer)
layer.renderer().setOpacity(0.5)

layer.triggerRepaint()
iface.legendInterface().refreshLayerSymbology(layer)

dem = QgsRasterLayer(r'C:\tmp\testdem.tif', "DEM")
QgsMapLayerRegistry.instance().addMapLayer(dem)

dem_result = processing.runalg("saga:slopeaspectcurvature", dem, 6, 0, 0, None, None, None, None, None, None,None, None, None, None, None, None)

gaussian = processing.runalg("saga:gaussianfilter", dem, 3, 1, 12, None)
result = processing.runalg("saga:slopeaspectcurvature", gaussian["RESULT"], 6, 0, 0, None, None, None, None, None, None,None, None, None, None, None, None)

#gaussian_layer = processing.load(gaussian["RESULT"])
curvature_layer = processing.load(result["C_GENE"])
slope_layer = processing.load(dem_result["SLOPE"])
curvature_layer2 = processing.load(result["C_GENE"])
slope_layer2 = processing.load(dem_result["SLOPE"])

set_layer_style(curvature_layer, "Blues", 9, True, -0.2, 0.2)
set_layer_style(curvature_layer2, "RdBu", 9, True, -0.2, 0.2)
set_layer_style(slope_layer, "Oranges", 9, False)
set_layer_style(slope_layer2, "WhiteToBlack", 2, False)

なんかできたもの(3レイヤ版)


国土地理院 基盤地図情報　数値標高モデル(10mメッシュ)　を利用(福島県いわき市付近)

### ７．CS立体図のキモとなる平滑化処理と曲率の計算について。
　～地形学的観点から思うところを述べます。（勝部）

曲率がどんな地形を表現しているかは、下記のモデル図が分かり易いです。
＃出典：「建設技術者のための地形図読図入門〈第3巻〉段丘・丘陵・山地」　鈴木隆介著


SAGAを使うといろんな曲率が計算されますが代表的な曲率としては以下の３つを覚えておけば良いと思います。ちなみにCS立体図で使われているのは曲率は基本的にはGeneral Curvatureです。

- General Curvature：３次元空間での斜面の凹凸形状を表します。Profile CurvatureとPlan Curvatureを合成することによっても求められます。
- Profile Curvature：斜面の縦断面（図の縦軸）に沿った凹凸を表します。
- Plan Curvature：斜面の横断面（図の横軸）に沿った凹凸を表します。（つまり等高線の凹凸）

- （補足）通常、CS立体図ではGeneral Curvatureのみを使っています。Plan Curvatureで作製すると水の流れを強調した図になります。私は、「水みち強調図」と呼んでいます。CS立体図を作製するときには曲率のレイヤを2回重ねますが、下のレイヤ（青白）のみをPlan Curvatureにすると湧水などを識別しやすい図になります。崩壊危険地抽出にはお勧めの図法です。（戸田）

これらの曲率を使うことによって上記のモデル図のような斜面形状が分類できることになっていますが、実際にはDEMを2回微分するため曲率の値は元のDEMの格子間隔や精度によって値がぶれやすく扱いずらい地形量です。なんの処理もしていないDEMから曲率を計算しCS立体図を作成すると下の例１のようになり、元のDEMデータの微細な形状や処理前のTIN構造などを拾ってしまいなんだかよくわからない絵となってしまいます。

CS立体図の作成にあたっては元のDEMに平滑化処理をかけることによってこの問題を回避しています。例２はガウシアンフィルタ（半径12、標準偏差3）を用いて平滑化処理を行ったDEMを用いて計算された曲率ですが、10～20mのスケールの地形の凹凸が旨く表現されています。

つまり美しいCS立体図を作成するためには、使用するDEMの格子間隔とターゲットとなる地形のスケールに合わせた平滑化範囲の設定、そしてDEMの精度に適したガウシアンフィルターの重みづけの設定が重要になると思います。

＜参考＞移動平均フィルタとガウシアンフィルタ
標準偏差の値が小さいほど平滑化の効果は小さくなり、大きいほど効果が大きくなる。



例１　未処理のDEMから計算されたGeneral Curvature（赤が凸、青が凹）とCS立体図



例２　平滑化されたDEMから計算されたGeneral Curvature（赤が凸、青が凹）とCS立体図

（補足）単純な移動平均ではなく、ガウシアンフィルタを使う理由
移動平均の場合、道路等の人工改変があると、そこから等距離の場所にゴーストの様な線が入ります。平均計算する際の検索円内に急激な地形変化が含まれるか？含まれないか？で値が大きく変わってしまうためです。ガウシアンフィルタを使うことで、距離が離れるにつれて影響が小さくなるようにウェイトをつけて、この様な急激な地形変化が図に反映されないようにしています。地理院の１０ｍメッシュなどでは、初めから人工改変の様な微地形は反映されていないため、移動平均でも大丈夫です。（戸田）

- 上の方に考え方を整理してみました。戸田さんなんか間違ってたら補足・修正お願いします！（勝部）
- 勝部様　有難うございます。いつか、ちゃんと論文書きます。．．．たぶん。（戸田）
- 共同で論文書いちゃいなよ！（いわさき）
- 是非、お願いします。ただし、私は筆が遅いのでご了承ください。（戸田）

### ８．Plugin
朝日さんがQGIS用のプラグインを作ってくれたよ！
https://github.com/waigania13/csmap_plugin

#### ８．１　動作確認

- QGIS 2.18.3 64btit版 Standalone（いわさき）
- ただし、通常起動では、SAGAプラグインが動かない。
- そのため、OSGeoコンソールから
qgis --confighpath d:\qgis
- として、空の設定をつくって起動するようにした。