QGIS で DTM から段彩図を作成する方法

# はじめに

QGIS（Quantum GIS）は、地理情報システム（GIS）のためのオープンソースソフトウェアです。QGIS は、多様な地理空間データの作成、表示、解析を可能にし、ユーザーに豊富なツールセットを提供します。

この記事では、デジタル標高モデル（DTM）を使用して段彩図を作成する手順について解説します。まずは QGIS を使用した段彩図の作成手順を説明し、その後、同様の処理をPythonを使用して PyQGIS で実装する方法を紹介します。

段彩図は、地形の標高を色分けして視覚的に表現する地図の一種です。段彩図を使用することで、地形の高低差や傾斜の変化を簡単に把握することができます。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2659159/551b3d2b-29d8-4add-c76f-5df204941659.png)
**段彩図に DTM で高さを適用した例**

# QGIS
## 段彩図の作成手順
QGIS でデジタル地形モデル（DTM）から陰影段彩図を作成する手順は以下の通りです。

### 1. DTM データの読み込み
1. QGIS を起動し、「レイヤ」→「レイヤを追加」→「ラスターレイヤを追加」を選択します。
    <img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2659159/03cfa5fd-6075-2f66-ca19-f68fc2fe2507.png" width="550px" />

1. DTM データ（例: GeoTIFFファイルなど）を選択し、「追加」をクリックします。
    <img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2659159/bdf33d5f-91ce-d033-09c1-29ca39a6d8dd.png" />

### 2. 段彩図の作成
1. プロセシングツールボックスで「ラスターレイヤの段彩図」を検索し、「ラスターレイヤの段彩図」ツールを選択します。
    <img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2659159/d8e4d271-a609-5698-aa25-12a3dd6be357.png" width="350px" />

1. 「入力ラスタ」には DTM データを選択し、「出力ファイル名」を指定します。
1.  「階層数」や「階層の色」を設定し、「OK」をクリックして段彩図を生成します。
    <img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2659159/9ae91674-6772-9a35-e24d-a709706dc93c.png" />

1. 出力された段彩図を確認します。
    <img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2659159/f8a875b7-abb8-148a-3713-6d72e89abfb6.png" width="550px" />

# PyQGIS
PyQGIS を使って GeoTIFF 形式の DTM を読み込み、段彩図を作成して、結果を GeoTIFF 形式で出力する手順を以下に示します。

## 事前セットアップ

1. **QGIS のインストール**：PyQGIS は QGIS に含まれている Python ライブラリです。まずはじめに QGIS をインストールしておく必要があります。
2. **Python スクリプトの準備**：PyQGIS スクリプトを作成するために Python 環境が必要となります。
3. **環境変数の設定**: QGIS の Python ライブラリを使用するには、環境変数を設定する必要があります。QGIS のインストールディレクトリを `C:\Program Files\QGIS 3.34.6` とした場合の設定例を以下に示します。
    
    | 環境変数 | 設定例 |
    |:--|:--|
    | Path | C:\Program Files\QGIS 3.34.6\bin<br/>C:\Program Files\QGIS 3.34.6\apps\qgis-ltr\bin<br/>C:\Program Files\QGIS 3.34.6\apps\grass\grass83\bin<br/>C:\Program Files\QGIS 3.34.6\apps\qt5\bin<br/>C:\Program Files\QGIS 3.34.6\apps\Python312\Scripts<br/>C:\Program Files\QGIS 3.34.6\apps\Python31 |
    | PYTHONHOME | C:\Program Files\QGIS 3.34.6\apps\Python312 |
    | PYTHONPATH | C:\Program Files\QGIS 3.34.6\apps\qgis-ltr\python |
    | QGIS_PLUGINPATH | C:\Program Files\QGIS 3.34.6\apps\qt5\plugins<br/>C:\Program Files\QGIS 3.34.6\apps\qgis-ltr\qtplugins |
    | QGIS_PREFIX_PATH | C:\Program Files\QGIS 3.34.6\apps\qgis-ltr |
    | QGIS_PYTHON_PATH | C:\Program Files\QGIS 3.34.6\apps\Python312\Lib\site-packages |
    
4. Python パッケージのインストール
    ```bash
    pip install pyqt5
    ```

## 実装例

### 1. 必要なモジュールのインポート

```python
from qgis.core import (QgsApplication,
                       QgsProject,
                       QgsRasterFileWriter,
                       QgsProcessingException)
from qgis.analysis import QgsRelief
import os
```

- `qgis.core.QgsApplication`: QGIS アプリケーションの初期化と終了に使用。
- `qgis.core.QgsProject`: QGIS プロジェクトの管理に使用。
- `qgis.core.QgsRasterFileWriter`: ラスタファイルの書き出しに使用
- `qgis.core.QgsProcessingException`: 処理例外を扱うために使用
- `qgis.analysis.QgsRelief`: DEM データから段彩図を生成するために使用

### 2. QGIS 環境の初期化
まず初めに `setPrefixPath` で QGIS のインストールパスを設定します。
その後、QGIS アプリケーションのインスタンスを作成し (`QgsApplication([], False)`)、QGIS アプリケーションの初期化を完了します (`initQgis()`)。

```python
QgsApplication.setPrefixPath(r"C:\Program Files\QGIS 3.34.6", True)  # QGIS のインストールパスを設定
qgs = QgsApplication([], False)  # QGIS アプリケーションを初期化
qgs.initQgis()  # QGIS アプリケーションの初期化を完了
```

### 3. プロジェクトの作成
`QgsProject.instance()` で、QGIS プロジェクトのインスタンスを取得します。

```python
project = QgsProject.instance()
```

### 4. 入力ファイルと出力ファイルのパス設定

入力ファイル (`dem_path`) と出力ファイル (`output_path`) の設定を行います。
また、後続の `QgsRelief` で引数で利用するため、`QgsRasterFileWriter.driverForExtension` を使い、出力ファイルの拡張子から GDAL ドライバ名 (`outputFormat`) を取得しておきます。

```python
dem_path = r"data\FG-GML-5338-04-DEM.tif"
output_path = r"data\DTM_relief.tif"
outputFormat = QgsRasterFileWriter.driverForExtension(os.path.splitext(output_path)[1])  # GTiff
```

:::note info
[GDAL ドライバ](https://gdal.org/drivers/raster/index.html) は、特定の地理空間データフォーマットを読み書きするためのモジュールであり、GDAL はこれらのドライバを提供するライブラリです。
QGIS は、GDAL を利用して、多種多様な地理空間データフォーマットの操作、変換、解析を行います。
:::

### 5. 段彩図の生成
DEM ファイルから段彩図を生成するためのインスタンスを作成します。

```python
relief = QgsRelief(dem_path, output_path, outputFormat)
```

### 6. 段彩の区分を自動的に生成
最適化された段彩区分を計算するよう、段彩区分を設定します。

```python
reliefColors = relief.calculateOptimizedReliefClasses()  # 最適化された段彩区分を計算
relief.setReliefColors(reliefColors)
```

### 7. z 係数の設定
z 係数（標高データの垂直倍率）を設定します。

```python
relief.setZFactor(1.000000)
```

### 8. 頻度分布の出力
頻度分布を CSV ファイルとして出力したい場合は、以下の通り記載します。

```python
relief.exportFrequencyDistributionToCsv(r"data\fd.csv")
```

### 9. ラスタ処理の実行
ラスタ処理を実行します。

```python
relief.processRaster()
```

### 10. QGIS 環境の終了
QGIS アプリケーションを終了します。

```python
qgs.exitQgis()  # QGIS アプリケーションを終了
```

## 実行
ソースコードを実行すると、QGIS を起動することなく、段彩図（`data\DTM_relief.tif`）を生成することができます。
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2659159/ece2dd01-e3aa-3bc0-d442-55e6cad95f1c.png" width="550px" />

<details><summary>参考：ソースコード全文</summary>

```python
from qgis.core import (QgsApplication,
                       QgsProject,
                       QgsRasterLayer,
                       QgsRasterFileWriter,
                       QgsCoordinateReferenceSystem,
                       QgsProcessingException)
from qgis.analysis import QgsRelief
import os

# QGIS 環境の初期化
QgsApplication.setPrefixPath(r"C:\Program Files\QGIS 3.34.6", True)  # QGIS のインストールパスを設定
qgs = QgsApplication([], False)  # QGIS アプリケーションを初期化
qgs.initQgis()  # QGIS アプリケーションの初期化を完了

# プロジェクトの作成
project = QgsProject.instance()

# DEMファイルのパス
dem_path = r"data\FG-GML-5338-04-DEM.tif"

# 出力ファイルのパス
output_path = r"data\DTM_relief.tif"
outputFormat = QgsRasterFileWriter.driverForExtension(os.path.splitext(output_path)[1])  # GTiff

 # DEMファイルから段彩図を生成するためのインスタンスを作成
relief = QgsRelief(dem_path, output_path, outputFormat)

# 段彩の区分を自動的に生成
reliefColors = relief.calculateOptimizedReliefClasses()  # 最適化された段彩区分を計算
relief.setReliefColors(reliefColors)

# z 係数
relief.setZFactor(1.000000)

## 頻度分布の出力
#relief.exportFrequencyDistributionToCsv(r"data\fd.csv")

# ラスタ処理を実行
res = relief.processRaster()

if res == 1:
    raise QgsProcessingException('入力ファイルを開けません。(' + dem_path + ')')
elif res == 2:
    raise QgsProcessingException('出力ファイルの GDAL ドライバを取得できません。')
elif res == 3:
    raise QgsProcessingException('出力ファイルを作成できません。(' + output_path + ')')
elif res == 4:
    raise QgsProcessingException('入力バンドを取得できません。')
elif res == 5:
    raise QgsProcessingException('出力バンドを作成できません。')
elif res == 6:
    raise QgsProcessingException('出力ラスタサイズが小さすぎます。')
elif res == 7:
    print('処理がキャンセルされました。')

# QGIS 環境の終了
print('処理が完了しました。(' + output_path + ')')
qgs.exitQgis()  # QGIS アプリケーションを終了
```
</details>