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

Last updated at 2025-09-11Posted at 2025-09-11

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

1. 日照計算の対象領域を広げたい

基本的には同じ事をやろうとしていますが、初級編では比較的狭い集落域を対象としたのに対して、中級編ではかなり広い町域を対象とします。

対象範囲	南北	東西	解像度	行	桁	セル
集落	6 km	4 km	50 cm	12,000	8,000	96,000,000
町	30 km	16 km	50 cm	60,000	32,000	1,920,000,000

町域は集落域のちょうど 20 倍の広さです。

1-1. 時間がかかりすぎる

最初は簡単に考えていました。集落域が 1 日分 1 時間未満で計算が完了したから、町域も 1 日分 20 時間あれば出来るだろう、それぐらいは辛抱して待てばよい、と。

ところが、実際にやってみると、一昼夜 r.sun を回しても戻ってこない。二昼夜たってもまだ回っている。三日目になると、どこかで永久ループに陥っているんじゃないかと心配になってくる。結局、計算が完了するのに 85 時間を要しました。

え？何で想定の 4 倍以上も時間がかかるの？という話は後でまたしますが、r.sun の計算量は対象領域の面積だけでは決まらないのですね。

ともかく 85 時間も待つのは我慢できない（計算対象日がまだ 4 日残っていた）ので対策を考えたのですが、すぐに思いついたのは次の二つです。

スパコンを買う
対象領域の解像度を下げる

1-2. マシン性能はほぼ限界

前者はもちろん即座に却下ですが、現在、筆者が使用している PC は、ヤフオクで買った激安の中古品とは言え、CPU は Xeon 8170M (26 cores/52 threads) が 2 基、メモリは 192 GB という贅沢なマシンなので、ハードウェアにテコ入れをしようとすると、後はもうスパコンぐらいしか考えられない訳です。

時間貸しのスパコンを借ることも考えましたが、それはそれでハードルが高そうです。

1-3. 解像度を落とすのが現実的

その結果、解像度を下げることを考えて作り直した表がこれです。

対象範囲	南北	東西	解像度	行	桁	セル
集落	6 km	4 km	50 cm	12,000	8,000	96,000,000
町	30 km	16 km	50 cm	60,000	32,000	1,920,000,000
			1 m	30,000	16,000	480,000,000
			2 m	15,000	8,000	120,000,000
			5 m	6,000	3,200	19,200,000

解像度を 2 m まで落とせば、集落の場合とあまり変らない時間で計算が出来るだろう、ということです。実用上もそんなに遜色はないでしょう。

解像度を落として対応するということであれば、話はここで終了です。

以下は、解像度を落とさずに広域の日照計算をしようと試みた作業の記録です。テーマも作業環境も非常に特殊なものですから、具体的な技術的覚え書きとしては、筆者以外にはほとんど無意味かもしれません。

しかし、GRASS GIS で日照計算を行う場合に役立つ一般的な考え方や、有用性の高いテクニックを解説しますので、興味のある方は是非とも読み進めてください。

2. コマンド・ラインで GRASS GIS を使う

初級編では日照計算に関わる作業を GRASS の GUI で行いました。そのときに気付いたことは、GRASS GIS は慣れてくるとコマンド・ライン・インタフェイスで使いたくなる ということです。

たとえば、メニューには必ず [g.region] や [r.mask] のように対応するコマンドが記されています。また、機能別のダイアログにおいても、設定項目に対応するコマンド・オプションが分かるように記されていますし、「コピー」ボタンを押せば、これから実行する（または既に実行した）コマンドの詳細な内容が文字列としてクリップ・ボードにコピーされるなど、あたかも、コマンド・ラインで使ってくれと言わんばかりの配慮がなされています。

2-1. PowerShell から GRASS を使う

いろいろな方法があると思いますが、筆者は Windows の PowerShell からスクリプトを通じてGRASS を呼び出して使っています。

例えば、下記は、集落の日照計算に使用した isg-sun.ps1 というスクリプトです。（実は、初級編の最後の方から、コマンド・ラインを使うようになりました）

isg-sun.ps1

# 作業環境設定
$GISDBASE = "D:\data\i-gis\grassdata"    # GRASS データベース
$LOCATION = "isg2447"                    # プロジェクト(Location)
$BASEMAPSET = "PERMANENT"                # ベースとなるマップセット
$GRASS = "C:\Program Files\QGIS 3.40.9\bin\grass84.bat"  # 環境に依存するので要注意

# 対象日リスト
$days = 35, 79, 125, 172, 355

$start = Get-Date
Write-Host "#### $start #### Started processing ..."

# region を masks400 に設定
Write-Host "Setting the region to mapset ($BASEMAPSET) ..."
& $GRASS "$GISDBASE\$LOCATION\$BASEMAPSET" --exec g.region region=masks400m@PERMANENT

# --- 逐次処理
foreach ($day in $days) {

    $t1 = Get-Date
    Write-Host "== $t1 == Started processing $day ..."

    & $GRASS "$GISDBASE\$LOCATION\$BASEMAPSET" --exec r.sun --overwrite --verbose `
        elevation=DSM@PERMANENT `
        aspect=DSM_ASPECT@PERMANENT `
        slope=DSM_SLOPE@PERMANENT `
        glob_rad=GLOB_$day `
        insol_time=INSOL_$day `
        day=$day step=0.2 `
        nprocs=52

    $t2 = Get-Date
    $tt = New-TimeSpan -Start $t1 -End $t2
    Write-Host "== $t2 == Completed processing $day"
    Write-Host "== Time taken : $tt =="
}

$end = Get-Date
$time = New-TimeSpan -Start $start -End $end
Write-Host "#### $end #### Completed processing."
Write-Host "#### Time Taken : $time ####"

このスクリプトを PowerShell から呼び出します。

PS D:\data\i-gis\grassdata> ./isg-sun.ps1

すると、対象日を 35, 79, 125, 172, 355 と自動的に入れ替えて、日照計算が 5 回繰り返して実行されます。

スクリプトを書く手間はかかりますが、毎回 GUI のダイアログで設定を行うことに比べれば、確実かつ手早く作業が出来ます。また、作業の正確な記録にもなり、後日の再現性も高いので、何かと便利です。

& $GRASS "$GISDBASE\$LOCATION\$BASEMAPSET" --exec は、後に続くコマンドを実行するための決まり文句です。

$GRASS は、GRASS を起動するバッチ・ファイル（GRASS のインストールされた環境に依存）
"データベース\プロジェクト\マップセット" は、カレント・マップセットの指定

r.sun など、後に続くコマンドは、GUI の機能別ダイアログにある「コピー」ボタンを使用して、実行する（実行した）コマンドとパラメータを取得して貼り付けると、楽に編集できます。筆者は、GUI で試して、行けそうだと思ったら、コピー & ペーストでスクリプトに移すようにしています。

コマンドの詳細なオプションについて確認したいときは、ウェブ・ブラウザで grass [command] 例えば grass r.sun や grass r.horizon を検索して、最新のマニュアル・ページを開くのがお薦めです。

3. 日照計算と地平角度について

日照計算で最も CPU 時間を食うのは、地平角度（horizon angle）の計算です。

3-1. 地平角度（horizon angle）

地平角度は、ある地点からある方位に向いたときに見える、地物と空（そら）の境界の仰俯角です。大洋や広大な砂漠の真ん中であれば、おそらくどの方位の地平角度も 0 度でしょう。逆に筆者が住んでいる集落のような中山間地域では、いずれの方位の地平角度であっても、軽く 10 度や 20 度に達します。そして、平野部であっても、解像度を高くして細かく見ると、背の高い樹木や建造物の近くでは、方位によっては地平角度が 90 度近くになります。

3-2. 日照計算における地平角度の役割

日照計算においては、この地平角度が、ある日時においてその地点に太陽光が到達するか否かを決定します。

以下に、日照計算の計算処理フローを記します。

1 計算領域を特定
2 対象日を指定
3 開始時刻（水平線日出時刻）と終了時刻（水平線日没時刻）を対象日から計算して決定
4 計算領域内の対象地点を変えながらループ (※1)
    A 日照時間を 0 で初期化
    B 日射量を 0 で初期化
    C 対象地点の傾斜と方位を取得 (※2)
    D 開始時刻から終了時刻まで step 刻みで時刻を変えながらループ
        a 対象地点と時刻における太陽の方位と高度を計算して取得
        b 対象地点における太陽方位の地平角度を取得 (※3)
        c 太陽高度 > 地平角度なら陽光が当たると判断
            - 日照時間 ... step 分の時間を追加
            - 日射量 ... 太陽の方位・高度と対象地点の方位・傾斜から計算し、step 分の時間を掛けて追加
        d 次の時刻へ
    E 次の対象地点へ
5 日照時間のマップを書き出し
6 日射量のマップを書き出し

(※1) 対象地点のループはそれぞれ独立して計算できるので並列化が可能
- nprocs オプション
- コア数の多い CPU の使用が有効
(※2) 対象地点の傾斜と方位は計算済みのマップを与えることが可能
- slope, aspect のオプション
(※3) 地平角度は計算済みのマップを与えることが可能
- horizon_basename, horizon_step のオプション
- 多数のマップを準備する必要がある

対象地点の傾斜・方位と地平角度はループの内側で何度も繰り返して参照されるため、その場で計算するのでなく、計算済みの値を持つマップから読み出すことが出来るようになっています。

ただ、地点の傾斜・方位については、それぞれ一つずつマップがあれば済むのに対して、地平角度のマップは、季節と時刻によって変化する太陽の方位に対応するために、非常に沢山のものを用意しておく必要があります。全方位を満遍なくカバーするのが理想ですが、通常は代表的な方位のマップだけを準備しておいて、その間に位置する方位については補間値で対応するようにします。例えば、10 度刻みなら 36 枚、5 度刻みなら 72 枚のマップを準備する訳です。

3-3. 地平角度の計算

地平角度の計算は、与えられた標高データ・マップを使って、対象地点から計算領域の外縁まで対象方位のラインに沿って標高データをスキャンして仰角を計算し、最も角度の高いものを選択するという方法で行います。

従って、計算領域が広くなれば広くなるほど、スキャンするラインが長くなって、計算に要する時間が長くなっていきます。

地平角度の計算にかかる時間は、

r.sun 内部で計算する場合 ... [計算領域のセル数] x [時間刻み step 総数] x [標高データ・スキャン長] に比例
r.horizon で事前計算する場合 ... [計算領域のセル数] x [方位刻み step 総数] x [標高データ・スキャン長] に比例

となります。

筆者の環境で観測した結果も、この関係式で説明できます。最初の試みにおいて、面積比 20 倍の領域での計算時間が 80 倍以上になったのは、標高データ・スキャン長がほぼ 4 倍になったためだと推測されます。

3-4. 地平角度の事前計算の是非

そして、ここが問題なのですが、この地平角度のマップを準備するために時間をかけても報われるとは限らない、ということがあります。

事前に計算済みのマップを作成するためには、実際には参照されない方位の分も含めて、かなり大量の計算をしなければなりません。そのため、トータルで考えると、日照計算の内部で地平角度の計算をする方が速かったりします。また、補間が入らない分、内部計算方式の方が精度は高くなります。

4. 地平角度計算の合理化

地平角度の計算については、次のような工夫をすることが出来ます。

4-1. 参照されない方位は計算しない

北緯 35 度あたりの地方では、真北を中心に東西 60 度ずつぐらいの範囲の地平角度は、日照計算に影響しませんので、計算をする必要がありません。太陽が最も北上する夏至でも、太陽の方位がその範囲まで入ることが無いからです。このことを考慮すると、計算量を約 2/3 に減らすことが出来ます。

この無視して良い方位の範囲は緯度によって異なりますので、対象となる土地の夏至の日出・日没の方位を確認してください。

なお、r.sun に計算済みの地平角度マップを渡す場合は、参照されない範囲のマップも渡す必要がありますので、0 で埋めたダミーのマップを作成することで対応します。

4-2. 高解像度版と低解像度版を合成する

これは、近くにある樹木や建造物は高解像度版で捉え、遠方の山は低解像度版で捉えて、二つを合成しようというアイデアです。

地平角度を計算する r.horizon には、maxdistance というオプションがあり、対象地点から標高データをスキャンする距離を限定することが出来ます。このオプションを使って、高解像度版では計算領域の周縁まで長々とスキャンすることをやめて、比較的近く、例えば 100 m 範囲内にある樹木や建造物だけを考慮に入れた計算を行い、高速化をはかります。

それだけでは地平角度が不当に低くなる箇所がたくさん出てきますので、別途作成した対象領域全体をカバーする低解像度版と突き合せ、値の高い方を取って合成して、近くの遮蔽物と遠方の山の両方を捉えた地平角度の合成マップを作成します。

このテクニックは、遠方に高い山があって、日の出や日没に影響を及ぼしているけれども、その山を計算領域に含めると、とんでもなく広い範囲になってしまう、という場合に問題を解決する手法として利用できます。高解像度版は最終的に日照計算をしたい狭い範囲で作成し、低解像度版はその山を含んだ広い範囲で作成して合成すればよいのです。

富士山が見える場所とか、阿蘇の外輪山の中とかで役に立つはずです。

4-3. PowerShell の並列処理を利用する

GRASS で地平角度を計算する r.horizon には、r.sun の nprocs のような並列化のためのオプションはありません。

かわりに PowerShell の並列処理機能を使うことが出来ます。（PowerShell 7.0 以降が必要です）

5. サンプル・コード

5-1. 地平角度マップ作成

上で述べた合理化策を採用した地平角度マップ作成のサンプル・コードを示します。

5-1-1. 高解像度版の作成

horz-high-res.ps1

# 作業環境設定
$GISDBASE = "D:\data\i-gis\grassdata"    # GRASS データベース
$LOCATION = "taka2447"                   # プロジェクト(Location)
$BASEMAPSET = "PERMANENT"                # ベースとなるマップセット
$GRASS = "C:\Program Files\QGIS 3.40.9\bin\grass84.bat"  # 環境に依存するので要注意

$HIGH_RES_MAPSET = "HIGH_RES"            # 高解像度マップセット (50cm)
$HIGH_RES_ELEV = "DSM_50CM"              # 高解像度用標高マップ
$HIGH_RES_HORZ_BASENAME = "HORZ_50CM"    # 高解像度地平角度マップのベース名

# パラメータ設定
$ex_start = 35    # 除外範囲開始
$ex_end   = 145   # 除外範囲終了
$step     = 5     # 角度刻み

# スレッド数（同時実行数）
$throttle = 8

# マップセット初期設定
function Init-Mapset {
    # 高解像度マップセットを作成
    Write-Host "Creating a high resolution mapset ..."
    & $GRASS "$GISDBASE\$LOCATION\$BASEMAPSET" --exec g.mapset -c mapset=$HIGH_RES_MAPSET

    # 高解像度マップセットの計算領域を設定 ... PERMANENT にある DSM_50CM に合わせる
    Write-Host "Setting the region to the high resolution mapset ..."
    & $GRASS "$GISDBASE\$LOCATION\$HIGH_RES_MAPSET" --exec g.region raster=$HIGH_RES_ELEV@$BASEMAPSET
}

# 処理開始
$start = Get-Date
Write-Host "#### $start #### Started creating high resolution horizon maps ..."

# マップセット初期設定
Init-Mapset

# 角度リストを生成
$angles = @()
for ($a = 0; $a -lt 360; $a += $step) { $angles += $a }

# 角度ごとの並列実行
$angles | ForEach-Object -Parallel {
    # 外側の変数は $using: で参照
    $angle = $_
    $tag = $angle.ToString("000")
    $GISDBASE = $using:GISDBASE
    $LOCATION = $using:LOCATION
    $BASEMAPSET = $using:BASEMAPSET
    $GRASS = $using:GRASS
    $HIGH_RES_MAPSET = $using:HIGH_RES_MAPSET
    $HIGH_RES_ELEV = $using:HIGH_RES_ELEV
    $HIGH_RES_HORZ_BASENAME = $using:HIGH_RES_HORZ_BASENAME

    if ($angle -ge $using:ex_start -and $angle -le $using:ex_end) {
        # 除外範囲 ... 何もしない
    }
    else {
        # 地平計算実行 高解像度
        $t1 = Get-Date
        Write-Host "== $t1 == Started creating the high resolution horizon map for angle [$angle] ..."
        
        & $GRASS "$GISDBASE\$LOCATION\$HIGH_RES_MAPSET" --exec r.horizon `
            elevation=$HIGH_RES_ELEV `
            direction=$angle `
            horizon="${HIGH_RES_HORZ_BASENAME}_$tag" `
            maxdistance=100 `
            --overwrite --verbose
    
        $t2 = Get-Date
        Write-Host "== $t2 == Completed creating the high resolution horizon map for angle [$angle]."
        $ts = New-TimeSpan -Start $t1 -End $t2
        Write-Host "== Time taken : $ts =="
    }

} -ThrottleLimit $throttle

# 処理終了
$end = Get-Date
Write-Host "#### $end #### Completed creating high resolution horizon maps"
$span = New-TimeSpan -Start $start -End $end
Write-Host "#### Time taken : $span ####"

元になる標高マップは $BASEMAPSET (='PERMANENT') にあると想定しています
- 高解像度で日照計算をしたいのは樹木や建造物の影響も可視化するためですので、DEM ではなく DSM を使っています
成果物は $HIGH_RES_MAPSET (='HIGH_RES') に HORZ_50CM_xxx という名前の一連のマップとして作成されます
- xxx は 000 から 5 度刻みで 355 まで、全 72 枚
- ただし、除外範囲分は欠落
並列実行数の上限は、CPU コア数よりも、計算領域の広さと実装メモリ容量に影響を受けます
- 小さめに設定する必要があります
- 筆者の PC (192 GB) では、12 ぐらいが上限でした
かなりの時間を食います

5-1-2. 低解像度版の作成

horz-low-res.ps1

# 作業環境設定
$GISDBASE = "D:\data\i-gis\grassdata"    # GRASS データベース
$LOCATION = "taka2447"                   # プロジェクト(Location)
$BASEMAPSET = "PERMANENT"                # ベースとなるマップセット
$GRASS = "C:\Program Files\QGIS 3.40.9\bin\grass84.bat"  # 環境に依存するので要注意

$HIGH_RES_MAPSET = "HIGH_RES"            # 高解像度マップセット (50cm)
$LOW_RES_MAPSET = "LOW_RES"              # 低解像度マップセット (2m)
$HIGH_RES_ELEV = "DSM_50CM"              # 高解像度用標高マップ
$LOW_RES_ELEV = "DSM_2M"                 # 低解像度用標高マップ
$HIGH_RES_HORZ_BASENAME = "HORZ_50CM"    # 高解像度地平角度マップのベース名
$LOW_RES_HORZ_BASENAME = "HORZ_2M"       # 低解像度地平角度マップのベース名
$COMP_HORZ_BASENAME = "HORZ_COMP"        # 合成された地平角度マップのベース名

# パラメータ設定
$ex_start = 35    # 除外範囲開始
$ex_end   = 145   # 除外範囲終了
$step     = 5     # 角度刻み

# スレッド数（同時実行数）
$throttle = 32

# マップセット初期設定
function Init-Mapset {
    # 低解像度マップセットを作成
    Write-Host "Creating a low resolution mapset ..."
    & $GRASS "$GISDBASE\$LOCATION\$BASEMAPSET" --exec g.mapset -c mapset=$LOW_RES_MAPSET

    # 低解像度マップセットの計算領域を設定 ... PERMANENT にある DSM_50CM に合わせる
    Write-Host "Setting the region to the low resolution mapset ..."
    & $GRASS "$GISDBASE\$LOCATION\$LOW_RES_MAPSET" --exec g.region raster=$HIGH_RES_ELEV@$BASEMAPSET
    # 低解像度マップセットの計算領域の解像度を 2m に変更
    & $GRASS "$GISDBASE\$LOCATION\$LOW_RES_MAPSET" --exec g.region -a res=2 align=$HIGH_RES_ELEV@$BASEMAPSET
}

# 低解像度用標高データを作成
function Make-Low-Res-Elev {
    $t1 = Get-Date
    Write-Host "== $t1 == Started creating a low resolution elevation map ..."

    & $GRASS "$GISDBASE\$LOCATION\$LOW_RES_MAPSET" --exec r.resamp.stats method=average input=$HIGH_RES_ELEV output=$LOW_RES_ELEV

    $t2 = Get-Date
    Write-Host "== $t2 == Completed creating a low resolution elevation map."
    $ts = New-TimeSpan -Start $t1 -End $t2
    Write-Host "== Time taken : $ts =="
}

# 処理開始
$start = Get-Date
Write-Host "#### $start #### Started creating low resolution horizon maps ..."

# マップセット初期設定
Init-Mapset

# 低解像度用標高データを作成
Make-Low-Res-Elev

# 角度リストを生成
$angles = @()
for ($a = 0; $a -lt 360; $a += $step) { $angles += $a }

# 角度ごとの並列実行
$angles | ForEach-Object -Parallel {
    # 外側の変数は $using: で参照
    $angle = $_
    $tag = $angle.ToString("000")
    $GISDBASE = $using:GISDBASE
    $LOCATION = $using:LOCATION
    $BASEMAPSET = $using:BASEMAPSET
    $GRASS = $using:GRASS
    $LOW_RES_MAPSET = $using:LOW_RES_MAPSET
    $LOW_RES_ELEV = $using:LOW_RES_ELEV
    $LOW_RES_HORZ_BASENAME = $using:LOW_RES_HORZ_BASENAME

    if ($angle -ge $using:ex_start -and $angle -le $using:ex_end) {
        # 除外範囲 ... 何もしない
    }
    else {
        # 地平計算実行 低解像度
        $t1 = Get-Date
        Write-Host "== $t1 == Started creating the low resolution horizon map for angle [$angle] ..."
        
        & $GRASS "$GISDBASE\$LOCATION\$LOW_RES_MAPSET" --exec r.horizon `
            elevation=$LOW_RES_ELEV `
            direction=$angle `
            horizon="${LOW_RES_HORZ_BASENAME}_$tag" `
            --overwrite --verbose
    
        $t2 = Get-Date
        Write-Host "== $t2 == Completed creating the low resolution horizon map for angle [$angle]."
        $ts = New-TimeSpan -Start $t1 -End $t2
        Write-Host "== Time taken : $ts =="
    }

} -ThrottleLimit $throttle

# 処理終了
$end = Get-Date
Write-Host "#### $end #### Completed creating low resolution horizon maps"
$span = New-TimeSpan -Start $start -End $end
Write-Host "#### Time taken : $span ####"

元になる標高マップを $LOW_RES_MAPSET (='LOW_RES') に作成します
成果物は $LOW_RES_MAPSET に HORZ_2M_xxx という名前の一連のマップとして作成されます
- xxx は 000 から 5 度刻みで 355 まで、全 72 枚
- ただし、除外範囲分は欠落
並列実行数の上限は、高解像度版に比べて、かなり大きく取ることが出来ます
比較的短い時間で完了します

5-1-3. 合成版の作成

horz-merge.ps1

# 作業環境設定
$GISDBASE = "D:\data\i-gis\grassdata"    # GRASS データベース
$LOCATION = "taka2447"                   # プロジェクト(Location)
$BASEMAPSET = "PERMANENT"                # ベースとなるマップセット
$GRASS = "C:\Program Files\QGIS 3.40.9\bin\grass84.bat"  # 環境に依存するので要注意

$HIGH_RES_MAPSET = "HIGH_RES"            # 高解像度マップセット (50cm)
$LOW_RES_MAPSET = "LOW_RES"              # 低解像度マップセット (2m)
$HIGH_RES_ELEV = "DSM_50CM"              # 高解像度用標高マップ
$HIGH_RES_HORZ_BASENAME = "HORZ_50CM"    # 高解像度地平角度マップのベース名
$LOW_RES_HORZ_BASENAME = "HORZ_2M"       # 低解像度地平角度マップのベース名
$COMP_HORZ_BASENAME = "HORZ_COMP"        # 合成された地平角度マップのベース名

# パラメータ設定
$ex_start = 35    # 除外範囲開始
$ex_end   = 145   # 除外範囲終了
$step     = 5     # 角度刻み

# スレッド数（同時実行数）
$throttle = 32

# マップセット初期設定
function Init-Mapset {
    # 高解像度マップセットの計算領域を設定 ... PERMANENT にある DSM_50CM に合わせる
    Write-Host "Setting the region to the high resolution mapset ..."
    & $GRASS "$GISDBASE\$LOCATION\$HIGH_RES_MAPSET" --exec g.region raster=$HIGH_RES_ELEV@$BASEMAPSET
}

# 処理開始
$start = Get-Date
Write-Host "#### $start #### Started creating composite horizon maps ..."

# マップセット初期設定
Init-Mapset

# 角度リストを生成
$angles = @()
for ($a = 0; $a -lt 360; $a += $step) { $angles += $a }

# 角度ごとの並列実行
$angles | ForEach-Object -Parallel {
    # 外側の変数は $using: で参照
    $angle = $_
    $tag = $angle.ToString("000")
    $GISDBASE = $using:GISDBASE
    $LOCATION = $using:LOCATION
    $BASEMAPSET = $using:BASEMAPSET
    $GRASS = $using:GRASS
    $HIGH_RES_MAPSET = $using:HIGH_RES_MAPSET
    $HIGH_RES_HORZ_BASENAME = $using:HIGH_RES_HORZ_BASENAME
    $LOW_RES_MAPSET = $using:LOW_RES_MAPSET
    $LOW_RES_HORZ_BASENAME = $using:LOW_RES_HORZ_BASENAME
    $COMP_HORZ_BASENAME = $using:COMP_HORZ_BASENAME

    if ($angle -ge $using:ex_start -and $angle -le $using:ex_end) {
        # 除外範囲 ... ダミーを作成
        $t1 = Get-Date
        Write-Host "== $t1 == Started creating dummy horizon map for angle [$angle] ..."
        
        $calc = "${COMP_HORZ_BASENAME}_$tag = 0"
        & $GRASS "$GISDBASE\$LOCATION\$HIGH_RES_MAPSET" --exec r.mapcalc "$calc"
    
        $t2 = Get-Date
        Write-Host "== $t2 == Completed creating dummy horizon map for angle [$angle]."
        $ts = New-TimeSpan -Start $t1 -End $t2
        Write-Host "== Time taken : $ts =="
    }
    else {
        # 高解像度版と低解像度版を合成する
        $1 = Get-Date
        Write-Host "== $t1 == Started creating composite horizon map for angle [$angle] ..."
        
        $calc = "${COMP_HORZ_BASENAME}_$tag = max(${HIGH_RES_HORZ_BASENAME}_$tag, ${LOW_RES_HORZ_BASENAME}_$tag@$LOW_RES_MAPSET)"
        & $GRASS "$GISDBASE\$LOCATION\$HIGH_RES_MAPSET" --exec r.mapcalc "$calc"
    
        $t2 = Get-Date
        Write-Host "== $t2 == Completed creating composite horizon map for angle [$angle]."
        $ts = New-TimeSpan -Start $t1 -End $t2
        Write-Host "== Time taken : $ts =="
    }

} -ThrottleLimit $throttle

# 処理終了
$end = Get-Date
Write-Host "#### $end #### Completed creating composite horizon maps"
$span = New-TimeSpan -Start $start -End $end
Write-Host "#### Time taken : $span ####"

最終成果物は $HIGH_RES_MAPSET (='HIGH_RES') に HORZ_COMP_xxx という名前の一連のマップとして作成されます
- xxx は 000 から 5 度刻みで 355 まで、全 72 枚

5-2. 地平角度マップを使う日照計算

上で作成した地平角度マップを使用する日照計算のサンプル・コードを示します。

sun.ps1

# 作業環境設定
$GISDBASE = "D:\data\i-gis\grassdata"    # GRASS データベース
$LOCATION = "taka2447"                   # プロジェクト(Location)
$BASEMAPSET = "PERMANENT"                # ベースとなるマップセット
$GRASS = "C:\Program Files\QGIS 3.40.9\bin\grass84.bat"  # 環境に依存するので要注意

$HIGH_RES_MAPSET = "HIGH_RES"            # 高解像度マップセット (50cm)
$ELEV = "DSM_50CM@$BASEMAPSET"           # 標高マップ
$ASPECT = "ASPECT@$BASEMAPSET"           # 方位マップ
$SLOPE = "SLOPE@$BASEMAPSET"             # 傾斜マップ
$HORZ_BASENAME = "HORZ_COMP"             # 地平角度マップのベース名
$HORZ_STEP = 5                           # 地平角度マップの方位ステップ

# 角度リスト
$days = 35, 79, 125, 172, 355

$start = Get-Date
Write-Host "#### $start #### Started processing ..."

# $HIGH_RES_MAPSET に region を設定
Write-Host "Setting the region to mapset ($HIGH_RES_MAPSET) ..."
& $GRASS "$GISDBASE\$LOCATION\$HIGH_RES_MAPSET" --exec g.region raster=$ELEV

# --- 逐次処理
foreach ($day in $days) {

    $t1 = Get-Date
    Write-Host "== $t1 == Started processing day [$day] ..."
    & $GRASS "$GISDBASE\$LOCATION\$HIGH_RES_MAPSET" --exec r.sun --overwrite --verbose `
        npartitions=128 `
        elevation=$ELEV `
        aspect=$ASPECT `
        slope=$SLOPE `
        horizon_basename=$HORZ_BASENAME `
        horizon_step=$HORZ_STEP `
        glob_rad=GLOB_RAD_$day `
        insol_time=INSOL_TIME_$day `
        day=$day step=0.2 `
        nprocs=52

    $t2 = Get-Date
    Write-Host "== $t2 == Completed processing day [$day]."
    $tt = New-TimeSpan -Start $t1 -End $t2
    Write-Host "== Time taken : $tt =="
}

$end = Get-Date
Write-Host "#### $end #### Completed."
$time = New-TimeSpan -Start $start -End $end
Write-Host "#### Time taken : $time ####"

標高マップに対応する方位マップと傾斜マップは作成済みであると仮定しています
作業は地平角度マップが格納されている $HIGH_RES_MAPSET で行います
- horizon_basename オプションにはマップセットの指定を含めることが出来ませんので、$HIGH_RES_MAPSET をカレントにする必要があります
- 最終成果物の GLOB_RAD_xxx および INSOL_TIME_xxx も $HIGH_RES_MAPSET 下に作成されます
地平角度マップのベース名、方位刻みは、事前に作成したものと合致させる必要があります
非常にメモリを食います
- これは、地平角度マップを使わない場合に比べて、格段に入力マップ数が増えるためです
  - 地平角度マップを使わない場合
    - 入力マップ = 標高マップ、方位マップ、傾斜マップ
    - 出力マップ = 日照時間マップ、日射量マップ
  - 地平角度マップを使う場合
    - 入力マップ = 標高マップ、方位マップ、傾斜マップ、ステップ数分の地平角度マップ
    - 出力マップ = 日照時間マップ、日射量マップ
- 通常は npartitions オプションによって入力マップを分割読み込みして、メモリ負担を軽減する必要があります
計算所要時間はそれなりに短縮されます
- 地平角度マップを使わずに 85 時間かかった計算が、9 時間程度で終了しました

結果から言うと、地平角度を事前計算する方式が圧倒的に優れているとは言えない状況です。事前準備の手間も考えると、地平角度の事前計算をせずに、単純に日照計算を行う方が、作成する対象日が少ない場合は、トータルの所要時間が短いかも知れません。

地平角度を事前計算する方式のメリットは

遠方にある遮蔽物の効果を考慮に入れることが出来る（低解像度版と高解像度版を合成するテクニックが使える）
npartitions オプションを使える（地平角度を内部計算する方法では npartitions オプションは使えない）

という二つの点でしょう。

6. 終りに

中級編は以上で終りです。主として、地平角度を事前計算する方法について解説しました。

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

i-GIS 岩座神ネット地理情報システム
- 岩座神農林業地図 ... 初級編レベル
- 多可町の地図 ... 中級編レベル

実際には、どちらの地図にも、まだ解説していないテクニックが適用されています。

DSM のノイズ除去
- 兵庫県が提供しているメッシュ・データの DSM には、送電線の架線の高さが含まれています。これは厳密に言えば外れ値（ノイズ）ではないのですが、そのままにしておくと、最終的な日照地図に巨大な万里の長城の影が表れることになります
山林エリアの解像度をわざと粗くする
- 山林エリアの日照地図は、50cm 解像度の DSM で作成すると細かすぎる影が出来て、却って見づらいものになりますので、住宅・農地エリアとは別に低解像度で日照地図を作って、二つを合成しています
- 山林エリアの日照地図は、本来なら DSM でなく DEM で作成するほうが有用ですが、一枚の地図にまとめるために DSM の低解像度版を採用しました（DEM を使うと、繋ぎ目で不自然な段差が出来る）
- 山林エリアと宅地・農地エリアを分けるマスクの作成が作業のキモになります

これらについても、上級編ならぬマニアック編を書いて解説するかもしれません。

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