1. はじめに

QGISで可視領域の分析をしたかったので、その実施の経緯と結果をまとめます。
ここでは、北海道百年記念塔を例として取り上げます。
この塔の頂上から見える範囲を可視領域として、塔の周囲からこの頂上を確認できる範囲を分析します。

2. 前提条件

以下に使用したソフト、ハード、データを記します。

ソフトウェア
- QGIS ver. 3.28.1-Firenze
- GRASS GIS ver. 8.3
ハードウェア
- ThinkPad X1 Carbon 4th
  - OS: Windows 10 Pro
  - CPU: Intel(R) Core(TM) i5-6200U CPU @ 2.30GHz 2.40 GHz
  - メモリ: 8GB
データ
- 国土地理院基盤地図情報
- 国土交通省 PLATEAU

3. 方法

可視領域の分析手法として、以下の2つの方法を使用します。

QGISのプロセシングツールボックスの見通し領域(Viewshed)ツール
GRASS GISのr.viewshed

見通し領域を知る上で、いずれか1つで良いのですが、両方試してみて違いを見てみます。

3.1 前準備

前準備として、対象範囲の数値標高モデル(DEM: Digital Elevation Model)をダウンロードし、ラスタデータ（GeoTiff）を作成しておきます。

国土地理院基盤地図情報からGeoTiffの作成方法については、福山大学工学部金子邦彦先生のWebサイトに分かりやすい解説記事があるので、そちらをご参照ください。

ここでは国土地理院の基盤地図情報ダウンロードサービスから、北海道百年記念塔のある付近の数値標高モデルをダウンロードします。
以下の範囲の数値標高モデルをダウンロードしました。

地域メッシュ	範囲（市町村）
644132～644135	札幌市南区～長沼町
644142～644145	札幌市中央区～南幌町、長沼町
644152～644155	札幌市北区～江別市、南幌町、岩見沢市
644162～644165	石狩市～新篠津村、岩見沢市

作成したGeoTiffをQGIS上に表示したものを以下に示します。
GeoTiffレイヤーの下に国土地理院標準地図を重ねています。

座標参照系(CRS)には、日本測地系2011(JGD2011)平面直角座標系12系(EPSG 6680)を使用しました。
平面直角座標系を使用した理由は、塔からの距離を正確に求めたかったからです。

3.2 方法1. QGISの見通し領域（Viewshed）ツールを使う方法

QGISの見通し領域ツールを使って分析を行います。
可視領域の解析には、QGISのプロセシングツールボックスを開き、"viewshed"と検索します。
いくつかの候補が出て来ますが、この中の「見通し領域」をダブルクリックします。

ダイアログの項目の内、入力レイヤには数値標高モデルから作成したGeoTiffレイヤーを指定します。
以下に使用したパラメタを示します。

視点の場所：塔のある地点を地図上で指定
視点の高さ：塔の高さである100mに設定
ターゲットの高さ：そのまま（1.0m）にしました。
視点からの最大距離：5,000m としました。

「実行」をクリックすると処理が始まります。そして処理が成功すると、以下のメッセージが表示されます。

ダイアログを閉じて出力結果を見ると、以下のように表示されます。
白い部分が塔から見える範囲、黒い部分は見えない範囲です。何となくそれっぽいですね。

このままだと見づらいので、レイヤーの透明度を変更します。レイヤーの透明度は、レイヤスタイルの透明度で変更できます。

この結果を見ると、塔から西側に向かっては見える場所が多いですが、反対に東側～南東側は見えない場所が多いようです。
東側～南東側は丘陵地になっており、この丘陵地によって塔からの視界が遮られてしまうようです。
反対に視点を塔周辺に置き換えた場合、塔の西側からは塔は見え、東側～南東側から見ることのできる範囲は限定されるようです。

3.3 方法2. GRASS GISのr.viewshedを使う方法

次に、GRASS GISのr.viewshed を使用して可視領域分析を行います。
何故、QGISではなくGRASS GISを使用するかですが、QGISのプロセシングツールボックスにもGRASS GISのr.viewshed があるのですが、私の環境ではこれを実行してもエラーになってしまったため、別途GRASS GISを実行して、その上で直接r.viewshedを使ってみます。

なお、ここではGRASS GIS ver.8.3をダウンロードして使用しましたが、QGISに付属するGRASS GIS ver.7.8.7でも使えるのではないかと思います。（確かめていないので、間違っていたらすみません。）

GRASS GISを起動します。

3.3.1 GRASS GISデータベースの作成

まず、GRASS GISを利用するために、まずデータベースを作成します。
GRASS GISを起動すると、レイヤーマネジャウィンドウとマップウィンドウが表示されます。
この内、レイヤーマネジャウィンドウ上で、ツールバーの左から4つ目にある"Add existing or create new database"ボタンをクリックします。

すると、"Create location?"と尋ねられるので、「はい」をクリックします。

次に、ロケーションの名前を入力します。Description欄はオプションです。Nextをクリックします。

CRS（座標参照系）を聞かれるので、一番上の"Select CRS from a list by EPSG or description"にチェックを入れて、Nextをクリックします。

CRSのリストが表示されます。ここではテキスト入力ボックスに6680（JGD2011 平面直角座標系12系のEPSG番号）を入力します。すると、リストの中に"JGD2011 / Japan Plane Rectangular CS XII"が表示されるので、これを選択してNextをクリックします。

最後に、ここまでの入力内容を確認するダイアログが表示されるので、問題なければFinishをクリックします。

メッセージダイアログが表示され、GRASS GISのデータベース、ロケーション、マップセットが作成されたことが表示されます。

3.3.2 数値標高モデルの取り込み

r.viewshedを実行するために、数値標高モデルをマップセットに取り込みます。ここでは先に作成したGeoTiffファイルをインポートします。

レイヤーマネジャのツールバーにある"Import raster data [r.import]"をクリックします。

ラスタデータをインポートするためのダイアログが表示されます。ここで、"Source input"に予め作成しておいたGeoTiffファイルを指定し、Importをクリックします。

レイヤーマネジャのConsoleタブにr.importの処理結果が表示されます。

マップ表示ウィンドウにはインポートされたGeoTiffファイルが表示されます。綺麗ですね。

ちなみにマップ表示モードは"2D view"になっていますが、これを"3D view"に変えると、3D表示ができます。

上図は石狩平野を北東側から眺めた3D表示ですが、山や谷、丘、川の流れが見えて、これも綺麗です。

3.3.3 r.viewshedの実行

さて、これで下準備が整いました。いよいよr.viewshedを実行します。

まず、viewshedの対象範囲を指定します。レイヤーマネジャのConsoleタブを開き、以下のコマンドを実行します。

g.region rast=（GeoTiffのレイヤー名）

対象範囲の指定が上手く行ったかを確認するには、以下のコマンドを実行します。

g.region -p

日本測地系平面直角座標系12系の座標系で、対象範囲の東西南北範囲が設定されています。

次に、Toolsタブを開き、テキストボックスに"viewshed"と入力します。すると、r.viewshedが表示されるので、これをダブルクリックします。

r.viewshedのダイアログが表示されます。各項目には以下を設定しました。

"Name of input elevation raster maps"：数値標高モデルのGeoTiffレイヤーを指定します。
"Name for output raster maps"：出力するレイヤー名を入力します。
"Coordinates of viewing position"：塔の平面直角座標の位置を入力します。位置の指定には、予め塔の緯度経度を調べておき、平面直角座標系12系でのXY位置を調べておくと良いです。

なお、位置の調べ方は、QGISのCoordinate Captureプラグインを利用すると、塔の位置の平面直角座標の位置をキャプチャできます。簡単に手順を記します。

Coordinate Captureプラグインをインストールし、QGISのベクタメニューから、"Coordinate Capture"を選択します。
塔の位置をクリックすると、Coordinate Captureのパネルに緯度経度と、平面直角座標系12系での位置が表示されます。この位置の右側にあるアイコンをクリックすると、平面直角座標での塔の位置がクリップボードにコピーされます。

別の方法としては、塔の緯度経度を求め、それを平面直角座標に変換する方法があります。
その場合、まず、塔の緯度経度を求めます。
緯度経度から平面直角座標系への変換にはK'z Labの平面直角座標→緯度・経度変換を利用させて頂きました。

塔の緯度経度は以下としました。

座標	値
緯度	43.05650
経度	141.49635

上の緯度経度は平面直角座標系12系では以下の値になりました。

座標軸	値
Y	-61391.268
X	-104539.960

塔の位置が分かり、値を入力したので、次に進みます。

Settingsタブを開きます。各項目に値を入力します。

視点の高さである"Viewing elevation above the ground"には100(m)を入力
ターゲットの高さである"Offset for target elevation above the ground"には1.0(m)を入力
視界の到達距離である"Maximum visibility radius"には5000(m)を入力

最後にRunをクリックします。諸々メッセージが表示されて、最後に"Command finished"が表示されたら、Closeをクリックします。

マップ表示ウィンドウに可視領域が表示されます。

なお、この処理結果は見える見えないの2値ではなく、連続量になっています。このことは、@ishiijunpei 氏の「GRASS GISの「r.viewshed」コマンドで出力されるラスタ画像」に解説されています。
こちらの記事を参考にさせて頂き、r.viewshedの処理結果を見える見えないの2値に変換します。

Toolsタブを開き、"reclass"と入力。候補中の"Reclassify [r.reclass]"を選択します。

r.reclassのダイアログが表示されます。以下の値を入力します。

"Name of raster maps to be reclassified"：先ほど作成したr.viewshedの結果レイヤー名を指定
"Name for output raster map"：出力先のレイヤー名を入力
"File containing reclass rules or enter value directory"：ishiijunpei氏の解説では別ファイルに記載しておくことが紹介されていますが、ここでは直接"0 thru 180 = 1", "180 thru 1000 = NULL"を入力

最後にRunをクリックします。

処理が終わると、マップ表示ウィンドウにr.reclassの結果が表示されます。

作成したr.reclassの結果をQGIS上で表示するため、この結果をGeoTiffファイルにエクスポートします。

Layersタブを開き、r.reclassの結果レイヤー上で右クリックし、"Export"を選択します。

r.out.gdalのダイアログが表示されます。以下の値を入力します。

"Name of raster map (or group) to export"：r.reclassで作成した結果レイヤー名を指定
"Name for output raster file"：出力先のGeoTiffファイル名を入力

最後にRunをクリックします。処理が上手く行けば、Command outputタブに、"Command finished"と表示されます。

作成したGeoTiffファイルをQGISに取り込みします。
QGISのレイヤメニューから[レイヤを追加]→[ラスタレイヤを追加]をクリックします。

先ほど作成したGeoTiffファイルを指定して、レイヤに追加します。

QGISにGRASS GISで作成したGeoTiffファイルが取り込まれました。
このページで表示可能な画像サイズが一杯になったため、QGISの画面は掲載できませんが、QGIS上でも表示できています。

この後、作成した可視領域の検証を行います。
続きは別記事に記載します。