GRASS GIS による日照計算（初級編）

Last updated at 2025-09-07Posted at 2025-09-06

GRASS GIS を使用して日照計算を行い、対象日を特定して、地点ごとの日照時間や日射量を可視化した地図を作成します。

日照計算の機能は QGIS Desktop 本体には含まれておらず、プロセシング・ツール・ボックスを使って GRASS GIS の日照計算機能を利用するか、または、GRASS GIS そのものを起動して日照計算をおこなうかのどちらかになります。ここでは後者の GRASS GIS を直接利用する方法を解説します。中級編で解説するような広範囲かつ高精細な日照地図を作成するためには、GRASS GIS の仕組みについて理解しておく必要がありますので、最初から直接に GRASS GIS を扱って慣れておく方が良いでしょう。

1. 実際の作業の前に

1-1. 環境

現在の筆者の作業環境は以下の通りです。

Windows 11
GRASS GIS 8.4.1
- QGIS Desktop 3.40.9 とともにインストールされるバージョン
- 実行パスは C:\Program Files\QGIS 3.40.9\bin\grass84.bat

ここでは、この作業環境を前提として解説を進めますが、Windows 10 でも、また、linux でも Mac でも、GRASS GIS さえ起動できれば似たような手順で作業が出来るはずです。

1-2. GRASS GIS の起動

プログラム・グループ QGIS 3.40.9 の中に GRASS GIS 8.4.1 のショートカット・アイコンがありますので、それで GRASS GIS を起動します。よく使うので「スタートにピン留め」しておくのが良いでしょう。

1-3. GRASS GIS のデータベース

GRASS GIS は、必要な全てのデータを内蔵のデータベースに置いて作業をします。

QGIS では、GeoTIFF や shp ファイル、PostgreSQL 上の PostGIS テーブルなど、外部のデータを直接に読んで作業をすることが出来ますが、GRASS GIS ではそれらを直接に読むのではなく、自身のデータベースにインポートされたデータを読んで計算作業をする仕組みになっています。計算作業の途中で必要になる作業データも、最終的な成果物も、この内部データベースに出力されます。

そして、GRASS GIS の成果物を QGIS Desktop や Web GIS で利用する必要が生じたときには、それを GeoTIFF ファイルや shp ファイルなどの外部ファイルとしてエクスポートすることになります。

1-3-1. データベース・ディレクトリ

初期状態では、このデータベースを保管するディレクトリが C:\Users\xxxx\AppData\Roaming\GRASS8 の下に作られています。

これをそのまま使っても構いませんが、十分な空き容量があるドライブ上に新しいデータベース・ディレクトリを作ることをお薦めします。また、GRASS GIS ではディスク I/O の速度が計算速度に影響する場合が多々ありますので、可能であれば SSD を利用してください。

上の追加ボタンを押すと、フォルダ選択ダイアログが表示されますので、あらかじめ作成した空フォルダを選択してください。

このときに、ほかで作成した既存のデータベース・ディレクトリを選択することも可能です。

GRASS GIS のデータベースの実体は、ディレクトリ構造 + 管理用 SQLite DB ファイルなので、ディレクトリごとバックアップしたり、別の場所に移したりすることが簡単に出来ます。

1-3-2. プロジェクト

データベース・フォルダを新規に作成すると、"Do you want to create a new project (also known as location)?" と尋ねてきます。「新規のプロジェクト（ロケーションとも呼ばれておるぞ）を作成したいかね？」と言うことですね。「はい」を選んで次に進みましょう。

プロジェクトの名称と説明を入力して次に進みます。

プロジェクトの名称と同じ名前のサブ・ディレクトリがデータベース・ディレクトリの配下に作成されることになります
プロジェクト名は、後々、スクリプト内に記述することを考えて、英数字だけのシンプルなものにするのが無難です
説明は省略しても構いません

次に CRS を選択するように促されます。

一番上の "Select CRS from a list of EPSG or description" を選択した状態で次に進みます。

どのような CRS を選択するのが良いかは、扱う地図と目的に依存しますので、一概には言えません。

日照計算に使用する r.sun モジュールは、本来は地理座標系（経緯度）の利用を前提としています。ただし集落域や町村域の日照計算であれば、JGD2011/2000 の平面直角座標系で問題なく実用可能ですので、他の地図と合わせる形で JGD2011(2000) / Japan Plane Rectangular CS のタイプの CRS を選択するのが良いでしょう。県域を越える広範な地図を扱う場合は、地理座標系の選択を検討して下さい。

CRS を選択して次へ進むと、プロジェクトの作成が完了します。

プロジェクトの CRS 設定は決定的に重要です。

というのは、プロジェクトに属する全てのデータが、例外なく、この CRS を共有するからです。異なる CRS をもつ外部データは、この CRS によって再投影されたデータがインポートされます。計算によって新たに作成されるデータもこの CRS を持つものになります。

不適切な CRS 設定でプロジェクトを作成してしまった場合は、無理をせずに、適切な CRS 設定で新しく作成してください。不要なプロジェクトの削除方法については後述します。

1-3-3. マップセット

新規にプロジェクトを作成すると、自動的に PERMANENT という名前のデフォルトのマップセットがその下に作成されます。

地図のラスタ・データやベクタ・データは、すべて、いずれかのマップセットの下に置かなければなりません。

マップセットは自由に作成して追加することが出来ます。PERMANENT は read-only のマスター・データ、WORK は作業用、DIST は最終成果物などと、用途に応じてマップセットを使い分けることが可能です。

規模の小さいプロジェクトでは、PERMANENT のマップセットだけを使っても構いません。PERMANENT は read-only、と言うのは推奨される運用上の慣行に過ぎません。

当面はこの PERMANENT マップセットだけを使って作業を進めることにします。

GRASS GIS を操作する上では、カレント・マップセット という概念、および、マップセットの 計算領域 (region) と マスク について理解することがキーポイントになります。作業を進めながら、少しずつ解説します。

2. 標高データをインポートする

準備が出来たので日照計算に不可欠な標高データを GRASS GIS にインポートします。

2-1. カレント・マップセットを選択する

最初に、作業をする カレント・マップセット を選択します。

上の図では、"D:\data\i-gis\grass-test" データベースの "isg-sun" プロジェクトの "PERMANENT" マップセットがカレント（current）として選択されています。

別のマップセットを選びたい場合は、データ・ツリーを展開して、目的のマップセット上で右クリック・メニューを出して「マップセット切り替え」を実行して下さい。

2-2. ラスタ・データをインポートする

上の図の import raster data [r.import] ボタンを押します。あるいは、メニューから、ファイル > ラスターデータのインポート > Simplified raster import with reprojection [r.import] を選びます。

"Source input" として標高データのラスタ・データ・ファイルを選択して、「インポート」ボタンを押します。

インポートされたラスタ・データは、カレント・マップセットの下に配置され、レイヤ・ツリーに追加されて、メイン・ウィンドウ上に描画されます。

このとき、

データはプロジェクトの CRS に合わせて再投影されます
データの地理的領域と解像度は、特に指定しなければ、ソース・データのものが維持されます

ラスタ・データのメタデータの確認

インポートしたラスタ・データのメタデータを確認しておきましょう。

データ・ツリー上で選択 > 右クリック > Show metadata
または、レイヤ・ツリー上で選択 > 右クリック > メタデータ

どちらでも構いません。右側のコンソール・ペインに以下のような情報が表示されます。

(Sat Sep  6 13:00:08 2025)                                                      
r.info map=DEM_50CM_ISG@PERMANENT                                               
 +----------------------------------------------------------------------------+
 | Map:      DEM_50CM_ISG@PERMANENT         Date: Sat Sep 06 11:51:12 2025    |
 | Mapset:   PERMANENT                      Login of Creator: in-admin        |
 | Project:  isg-sun                                                          |
 | DataBase: D:\data\i-gis\grass-test                                         |
 | Title:    DEM_50CM_ISG                                                     |
 | Timestamp: none                                                            |
 |----------------------------------------------------------------------------|
 |                                                                            |
 |   Type of Map:  raster               Number of Categories: 0               |
 |   Data Type:    FCELL                Semantic label: (none)                |
 |   Rows:         12000                                                      |
 |   Columns:      8000                                                       |
 |   Total Cells:  96000000                                                   |
 |        Projection: JGD2011 / Japan Plane Rectangular CS V                  |
 |            N:     -93000    S:     -99000   Res:   0.5                     |
 |            E:      52000    W:      48000   Res:   0.5                     |
 |   Range of data:    min = 148.49  max = 1005.17                            |
 |                                                                            |
 |   Data Description:                                                        |
 |    generated by r.in.gdal                                                  |
 |                                                                            |
 |   Comments:                                                                |
 |    r.in.gdal -a -k input="D:\data\i-gis\hyogo50cm\DEM\TIFF-ISG\DEM_50CM\   |
 |    _ISG.tif" output="DEM_50CM_ISG" memory=300 offset=0 num_digits=0        |
 |                                                                            |
 +----------------------------------------------------------------------------+
(Sat Sep  6 13:00:08 2025) コマンド終了 (0 秒)

主な項目は、

南北 12,000 行 x 東西 8,000 桁 = 96,000,000 セルのデータ
CRS : JGD2011 / Japan Plane Rectangular CS V
北端 : -93,000 / 南端 : -99,000 / 南北の解像度 0.5m
東端 : 52,000 / 西端 : 48,000 / 東西の解像度 0.5m
標高データ最小値 : 148.49m / 標高データ最大値 : 1005.17m

という内容です。

2-3. 計算領域を設定する

日照計算に必要なデータのインポートが完了しましたので、次に 計算領域（region） を設定します。

ここで計算領域とは、データを読み書きする 範囲（東西南北の広がり） と 解像度 を意味します。

領域外のデータは無いものとして扱われます
入力データは、計算領域の解像度で再サンプリングされた値が計算に使用されます
出力データは計算領域の範囲と解像度を持ったものになります

計算領域を設定する方法として最も簡単なものは、ソースのラスタ・データに合わせる、という方法です。

メニューで、設定 > Computational region > 領域設定 [g.region] を選びます。

既存領域タブ
- Set region to match raster map(s) に先ほどインポートしたラスタ・データを選択

「実行」ボタンを押すと、ソースのラスタ・データに合致する計算領域が設定されます。

計算領域の確認

メニューで、設定 > Computational region > Show current region [g.region -p] を選びます。

すると、右側のコンソール・ペインに以下のような情報が表示されます。

(Sat Sep  6 13:45:02 2025)                                                      
g.region -p                                                                     
projection: 99 (JGD2011 / Japan Plane Rectangular CS V)
zone:       0
datum:      towgs84=0,0,0,0,0,0,0
ellipsoid:  grs80
north:      -93000
south:      -99000
west:       48000
east:       52000
nsres:      0.5
ewres:      0.5
rows:       12000
cols:       8000
cells:      96000000
(Sat Sep  6 13:45:03 2025) コマンド終了 (0 秒)

ソースのラスタ・データのメタデータと合致していることが分かりますね。

また、メインの地図表示ウィンドウでも、計算領域の範囲が赤い矩形で表示されます。

3. 傾斜と方位のマップを作成する

日照計算の準備として、標高データのマップから土地の傾斜と方位のマップを作成します。

これらのマップは、日照時間の計算には関係しませんが、日射量に影響します。つまり、南向きの斜面は北向きの斜面に比べて同じ日照時間でも陽光から受け取るエネルギー量が大きい、というようなことがあり、土地の傾斜と方位は日射量の計算で重要な要因となるわけです。

日照計算を一度しか行わないのであれば、傾斜と方位の計算を日照計算の中で実施しても構わないのですが、通常は、対象日を変えて何回か日照計算を繰り返しますので、事前に計算を済ませておいて、その結果を日照計算に渡すのが一般的な手法です。

この作業には r.slope.aspect を用います。

3-1. r.slope.aspect

メニューから、ラスター > 地形解析 > 傾斜と方位 [r.slope.aspect] を選びます。

必須タブ
- Name of input elevation raster map:（入力する標高ラスタ・マップの名前）... インポートした標高ラスタ・マップを選択
出力タブ
- Name for output slope raster map:（出力する傾斜ラスタ・マップの名前） ... 適当に
- Name for output aspect raster map:（出力する方位ラスタ・マップの名前） ... 適当に
- 残りの項目は入力不要
設定タブ
- このままでよい
オプション・タブ
- 既存のファイルに上書きする ... 必要に応じて
- Verbose モジュール出力 ... 好みで
- Number of threads for parallel computing ... PC の CPU 性能に合わせて（デフォルトは 1）
- Maximum memory to be used (in MB) ... PC の実装メモリ容量に合わせて（デフォルトは 300）

そして「実行」ボタンを押すと、現在の計算領域の設定に従って、傾斜マップと方位マップが計算されて作成され、カレント・マップセットである PERMANENT の下に追加されます。

4. 日照時間と日射量のデータを作成する

4-1. 対象日を決める

作成するマップの範囲は上記からおのずと決定されますが、計算対象とする日を決める必要があります。季節によって太陽の軌跡が変り、それに伴って日照時間と日射量も変化するためです。

筆者は次の5つの日を対象日として選びました。

節気	月日	通算日
冬至	12月21日	355
立春	2月4日	35
春分	3月20日	79
立夏	5月5日	125
夏至	6月21日	172

夏至から冬至に向う立秋・秋分・立冬は、それぞれ、立夏・春分・立春と太陽の軌跡が同じになりますので対象日から外しています。

通算日は、1月1日を1とする年間の通算日です。次に述べる r.sun の計算では対象日の指定に通算日を使用します。

では、いよいよ、r.sun を使って日照計算を行います。対象日は春分（通算日で言うと 79）とします。

4-2. r.sun

メニューから、ラスター > 太陽照度と影 > 太陽照度と照射 [r.sun] を選びます。

入力タブ
- Name of the input elevation raster map（入力する標高ラスタ・マップの名前） ... インポートした標高ラスタ・マップを選択
- Name of the input aspect map（入力する方位マップの名前） ... 作成済みの方位マップを選択
- Name of the input slope raster map（入力する傾斜ラスタ・マップの名前） ... 作成済みの傾斜マップを選択
- 他の項目はそのままにしておく
出力タブ
- Output of global(total) ... irradiance/irradiation raster map [Wh.m-2.day-1] (mode-2)（出力する全放射照度ラスタ・マップの名前） ... 計算対象日をファイル名に含める（e.g. GLOB_RAD_79）
- Output insolation time raster map [h] (mode-2 only)（出力する日照時間ラスタ・マップの名前） ... 計算対象日をファイル名に含める（e.g. INSOL_TIME_79）
- 他の項目は空白のままにしておく
テーマ・タブ
- No. of day of the year (1-365) ... 計算対象日を通算日で入力
- Time step when computing all-day radiation sums（一日の日照を積算するときの時間の刻み） ... 好みでより細かく（デフォルトは 0.5 時間刻み）
オプション・タブ
- 既存のファイルに上書きする ... 必要に応じて
- Verbose モジュール出力 ... 好みで
- Number of threads which will be used for parallel computing ... PC の CPU 性能に合わせて（デフォルトは 1）
- レイヤーツリーに作成したマップを追加 ... 好みで

最後に「実行」ボタンを押すと、コマンド出力タブに切り替って、計算の進捗状況が表示されます。

4-3. 時間がかかります

r.sun は鬼のように CPU 時間とメモリを食います。

コンソール出力の実例を示します。

(Sat Sep  6 15:55:36 2025)                                                      
r.sun --verbose elevation=DEM_50CM_ISG@PERMANENT aspect=DEM_ASPECT@PERMANENT slope=DEM_SLOPE@PERMANENT glob_rad=GLOB_RAD_79 insol_time=INSOL_TIME_79 day=79 nprocs=52
Number of threads <52>
Mode 2: integrated daily irradiation for a given day of the year
Using Linke constant: 3.000000
Using albedo constant: 0.200000
Using slope map <DEM_SLOPE@PERMANENT>
Using aspect map <DEM_ASPECT@PERMANENT>
(Sat Sep  6 16:40:56 2025) コマンド終了 (45 分 19 秒)

計算時間は、45 分 19 秒。条件は以下の通りです。

計算領域 : 12,000 x 8,000 = 96,000,000 セル
- 0.5m の解像度なので、6km x 4km の範囲
時間ステップ : 0.5 時間
並列スレッド : 52
PC のスペック :
- CPU : Xeon 8170M (26 cores/52 threads) x 2
- Meomory : 192 GB

これは、PC のスペックが高いから 45 分そこそこで完了している、と考える方が良いでしょう。一般の PC だと、もっともっと長い時間がかかるはずです。いや、時間がかかっても完走すれば良い方で、場合によってはメモリ不足のために途中で倒れる可能性も大いにあります。

4-4. 高速化の方策

r.sun の高速化のために取れる方策はそれほど多くありません。

PC にコア数の多い高性能な CPU を積む
PC にメモリを沢山積む
計算領域を狭くする
計算領域の解像度を落とす

最初のものが最も効果があります。r.sun は並列処理が可能な計算処理なので、コア数（スレッド数）の多い CPU の使用は非常に有効です。ただし、今すぐには対応できないこと、そして、お金がかかることが難点です。

メモリの追加は r.sun の安定動作に大いに寄与します。しかし、物理メモリの追加はお金もかかりますし、PC の仕様上の上限もありますので、思うようにできないことも多いでしょう。せめて十分な仮想メモリを確保できるようにするぐらいが、ただちに実施できる対策でしょうか。

後ろの二つは、対象となるデータの量を減らして計算量を抑えようという対策です。

ただし「計算領域を狭くする」については注意が必要です。太陽をさえぎる可能性のある地物、たとえば、北以外の方位にある高い山の稜線を計算領域外に追い出すと、計算の誤差が大きくなります。従って、対象となる地形によっては、むやみに計算領域を狭くすることは出来ません。

結局のところ、今すぐに出来ることと言えば、時間がかかっても辛抱強く待つか、計算領域の解像度を落とすか、ということになりそうです。

5. 成果物のエクスポート

下の図は、日照計算で出力された日射量マップ GLOB_RAD_79 です。

これをそのまま GeoTIFF としてエクスポートしても良いのですが、集落域が分かりやすいように、集落境界で切り出したものをエクスポートしたいと思います。

5-1. ベクタ・データのインポート

まず、集落域を示すベクタ・データをインポートします。

上の図の import vector data [v.import] ボタンを押します。あるいは、メニューから、ファイル > ベクトルデータのインポート > Simplified vector import with reprojection [v.import] を選びます。

"Source input" として集落域を示すベクタ・データ・ファイルを選択して、「インポート」ボタンを押します。

インポートされたベクタ・データは、カレント・マップセットの下に配置され、レイヤ・ツリーに追加されて、メイン・ウィンドウ上に描画されます。

5-2. マスクの作成

メニューから ラスター > マスク [r.mask] を選びます。

ベクトル・タブ
- Name of vector map to use as mask: ... インポートした集落域ベクタ・マップを指定

「実行」ボタンを押すとマスクが作成されて適用されます。

上の図はマスクを効かせた状態で日射量マップ GLOB_RAD_79 を表示しているところです。集落域に含まれるデータだけが表示されているのが確認できます。

なお、マスクが有効かどうかは、次の2つのことでも確認できます。

データ・ツリーに MASK というラスタ・データが有るか無いか
ウィンドウ右下に赤文字で MASK と表示されいるか否か（MASK インジケータ）

マスクを削除したいときは、MASK インジケータをクリックすると、MASK を削除するかどうかを尋ねてきますので、「はい」を選択して下さい。

面白い方法としては、データ・ツリー上の MASK というラスタ・データの名前を MASK 以外に変更するという手があります。そうするとマスクが消えます。そして、名前を MASK に戻すとマスクが復活します。

5-3. 計算領域の縮小

現在の計算領域はマスクの領域よりかなり広いものです。

このままエクスポートすると随分と空白部分（NoData）の多い、無駄に大きなファイルになりますので、計算領域をマスクの範囲まで縮小します。

ふたたび g.region を使用します。メニューで、設定 > Computational region > 領域設定 [g.region] を選びます。

境界タブ
- このラスターマップからの非NULLデータに会うまで領域を縮小します: ... ラスタ・データ MASK を選択

実行すると、下の図のように、計算領域がマスクの外接矩形まで縮小されます。

5-4. ラスタ・データのエクスポート

メニューから ファイル > ラスターマップのエクスポート > 一般的なエクスポートフォーマット [r.out.gdal] を選びます。

必須タブ
- Name of raster map (or group) to export ... エクスポートするラスタ・マップを指定
- Name for output raster file ... 保存するディレクトリとファイル名を指定
- Raster data format to write ... COG を推奨
  COG は、Cloud Optimized GeoTiff。通常の GeoTiff よりも部分的な読み出しに強く、Web GIS 等で使用するのに向いている
  小さいものなら GeoTIFF（ドロップダウン上の表記は GTiff）でもよい
作成タブ
- Do not write GDAL standard colortable ... チェックする
- Do not write non-standard metadata ... チェックする
- データタイプ ... Float32
  普通は Float64 の精度は不要
- Creation option(s) to pass to the output format driver: ... BIGTIFF=YES を推奨
- Number of overviews to create for the output dataset: ... 好みで
  COG が抱え持つ縮小画像の数
  この程度の大きさのラスタ・データなら 2 ぐらいか？ちょっと適当
オプション・タブ
- Force raster export despite any warnings of data loss ... 必要に応じて
  GRASS が Float64 で持っているデータを Float32 で書き出したい時などに有効
- 既存のファイルに上書きする ... 必要に応じて
- Verbose モジュールの出力 ... 好みで

以上のように設定して「実行」ボタンを押すと、ラスタ・データが COG 形式の GeoTIFF としてエクスポートされます。

6. 計算領域（region）とマスク（mask）

GRASS GIS 操作の要諦となる計算領域（region）とマスク（mask）についてまとめておきます。

6-1. 計算領域（region）

何か：GRASS が計算・入出力に使う範囲（N/S/E/W）と解像度の設定
単位：各マップセットごとに保持（$MAPSET/WIND）
影響：
- 領域外は NULL として扱われる
- 入力ラスタは region の解像度にリサンプルされて使われる
- 出力ラスタは region の範囲・解像度で書き出される
基本操作：
- GUI：設定 > Computational region > 領域設定 [g.region]
- CLI：
  - 既存ラスタに合わせる：g.region raster=DEM
  - 解像度だけ指定：g.region res=1（1m）
  - グリッド整合：g.region align=DEM（境界を DEM のセル格子に揃える）
  - 内容確認：g.region -p
  - 保存/復元：g.region save=full / g.region region=full
実務のコツ：
- 計算前に必ず領域を確認
- ラスタを重ねるなら align= を使って格子ずれを防ぐ
- 解像度を途中で変えると結果が変わる。再計算の意図がない限り触らない

6-1. マスク（mask） ... 計算対象の限定

何か：ラスタ計算で有効セル/無効セルを切り分ける仕組み。MASK という特別名のラスタで表現
影響：マスク外は NULL として扱われる（計算・出力とも除外）
基本操作：
- GUI：ラスター > マスク [r.mask]
- CLI：
  - ラスタで作成：r.mask raster=landuse
  - ベクタで作成：r.mask vector=area where="name='集落域'"
  - 反転：r.mask -i ...（外を残し内を除外）
  - 解除：r.mask -r（MASK を削除）
確認：データツリーに MASK がある／ウィンドウ右下に赤字 MASK 表示
実務のコツ：
- マスクはラスタ演算にのみ厳密に効く。ベクタ処理は原則として無関係
- MASK はマップセット単位で有効。別マップセットには影響しない

6-3. よくある落とし穴

region を確認せずに計算 ... 想定外の解像度・範囲で出力される
align= を使わない ... セル格子ずれでピクセル単位の誤差・ジャギー
マスク解除忘れ ... 次の処理・エクスポートが欠落する

7. 終りに

初級編は以上で終りです。GRASS GIS による日照計算の基本的な方法は解説出来たと思っています。

この先は、どのような地形と範囲を対象にするか、どの程度の精度の出力が必要か、どのようなスペックの PC を使うことが出来るかによって、いろいろと変ってきます。

中級編では、より広範囲な領域で日照計算を行う場合について、いくつかノウハウを共有したいと考えています。

ここで作成した日照地図を以下のサイトで掲載していますので、参考までにご覧ください。

i-GIS 岩座神ネット地理情報システム
- 岩座神農林業地図

You get articles that match your needs
You can efficiently read back useful information
You can use dark theme

What you can do with signing up