【QGIS×Python】DEMから起伏量図を作成しよう

るみ

2024年8月16日 21:59

QGISで起伏量図を作成してみよう。
起伏量図とは、ある区域において設定したグリッドごとの標高差（起伏量）の分布図であり、地形の起伏の激しさを表したものです。

この記事の最後に起伏量図を作成するpythonスクリプトを置いておきます。

手順

１：DEMのラスタを用意
２：グリッドの作成
３：グリッド内の起伏量を計算（ゾーン統計量）
４：表示方法の設定（シンポロジ、ラベル）

手書きで起伏量図を作るときは、グリッドごとに等高線と睨めっこする必要があるので非常に時間がかかりますが、QGISを用いることでより速く簡単に起伏量図を作成できます。

1.DEMの用意

調査地域のDEMをダウンロードして、ラスタに変換しましょう。
ここではDEMの解説を省略します。DEMについては以下のサイトの解説が分かりやすいです。

2.グリッドの作成

DEMラスタを用意したら、DEMに合わせたグリッドを作成しましょう。

グリッド（ベクタ）は「ベクタ」→「調査ツール」→「グリッドを作成」から作成できます。

続いて、グリッドの設定をします。

グリッドタイプ：方眼なら「長方形」を選択
グリッドの範囲：地図上で選択もしくはラスタレイヤを選択してグリッド範　　　　　囲を設定します。
水平/垂直方向の間隔：グリッドの大きさを設定します。正方形のグリッドを作る場合は、同じ長さを入力。
出力グリッドのCRS：グリッドの座標系と投影法を設定。DEMレイヤおよび地図のレイヤに合わせる。

＊「重なり」は何なのかよく分からない。

上記のように入力したら保存先を指定し「実行」を押してください。

こんな感じになります。500m*500mのグリッドを作成しました。

3.グリッド内の起伏量を計算

ゾーン統計量（ベクタ）を用いてグリッド内の起伏量を計算します。

プロセシングツール内の「ゾーン統計量(ベクタ)」を開きます。

入力レイヤ：グリッドレイヤ
ラスタレイヤ：DEMレイヤ
計算する統計量：範囲(range)を選択する

「範囲」はラスタレイヤの最大値-最小値を計算した値です。DEMレイヤには標高データが入っているので、「範囲」を計算することで起伏量が求められます。

起伏量を求めたレイヤが出来上がりました。

4.表示方法の設定

最後に、起伏量に応じてグリッドの色を塗り分けてみましょう。
左のレイヤボックスにあるグリッドのレイヤをダブルクリックして、「シンポロジ」を開きます。
シンポロジで「連続値による定義」を選択。

値：「range」を選択。ここでは表示に用いる値を選択します。
カラーランプ：表示するグラデーションの種類を選択。
モード、間隔サイズ：数値の間隔を設定。ここでは30mの固定間隔に設定しています。
レイヤレンダリング：不透明度などの設定

上記の設定を済ませたら「分類」を押して判例が表示されるのを確認してください。

完成！！！

Pythonスクリプト

以下のスクリプトをQGISのPythonコンソールにコピペすれば使えます。なお、パラメータ部分を各自で置き換えてください。

from qgis.core import QgsProcessingFeatureSourceDefinition,QgsCoordinateReferenceSystem,QgsProject, QgsVectorLayer
from qgis import processing



#パラメータ設定
#グリッドに関するパラメータ設定
input_layer_name = "DEMレイヤ名"  #グリッドを作成する範囲のレイヤ名
grid_length_EW = 500 #長方形の東西方向の辺の長さ(メートル)
grid_length_NW = 500 #長方形の南北方向の辺の長さ(メートル)
output_path_grid = "パス"  #グリッドを保存する出力ファイルのパス
EPSG_code = 32654 #作成するグリッドのEPSGコード。例として32654を置いてます

#起伏量レイヤに関するパラメータ設定
DEM_path = "パス" #ゾーン統計量で用いるDEMレイヤのパス
output_path_KIHUKU = "パス" #起伏量を含むベクタレイヤの出力パス
layer_name_KIHKU = "起伏量を含むベクタレイヤのレイヤ名"



#グリッド作成
#入力レイヤを取得
input_layer = QgsProject.instance().mapLayersByName(input_layer_name)[0]
if not input_layer:
    raise Exception(f"レイヤ '{input_layer_name}' が見つかりません。入力したレイヤ名が正しいか確認してください")

#「native:creategrid」の実行
processing.run("native:creategrid", {
    'TYPE': 2,  #グリッドタイプの指定。2=長方形
    'EXTENT': QgsProcessingFeatureSourceDefinition(input_layer.source(), selectedFeaturesOnly=False),
    'HSPACING': grid_length_EW,
    'VSPACING': grid_length_NW,
    'CRS': QgsCoordinateReferenceSystem(EPSG_code),  #グリッドのCRSの設定
    'OUTPUT': output_path_grid  #出力ファイルのパス
})

print(f"グリッドが作成され、{output_path_grid} に保存されました。")

#作成したグリッドレイヤをプロジェクトに追加
grid_layer = QgsVectorLayer(output_path_grid, "grid500500", "ogr") 
if grid_layer.isValid():
    QgsProject.instance().addMapLayer(grid_layer)
    print("グリッドレイヤを作成しました。")
else:
    print("グリッドレイヤが無効です。パスや設定を確認してください。")



#ゾーン統計量（ベクタ）で起伏量を計算する。
processing.run("native:zonalstatisticsfb", {
    'INPUT': output_path_grid,
    'INPUT_RASTER': DEM_path,
    'RASTER_BAND':1, #DEMラスタの標高データを持つバンド値
    'COLUMN_PREFIX':'_range',
    'STATISTICS':[7], #計算する統計量。７は範囲(=最大値-最小値)。
    'OUTPUT': output_path_KIHUKU
})

#作成した起伏量図レイヤ(ラスタ)をプロジェクトに追加
KIHUKU_layer = QgsVectorLayer(output_path_KIHUKU, layer_name_KIHKU, "ogr") 
if grid_layer.isValid():
    QgsProject.instance().addMapLayer(KIHUKU_layer)
    print("起伏量レイヤを作成しました。")
else:
    print("起伏量レイヤが無効です。パスや設定を確認してください。")

試しに鳥海山の起伏量図を作成してみました。山頂付近の起伏量が大きいことが分かりますね。

それでは～