GRASS GISで日照時間を計算してみる！！

はじめに

この記事は、FOSS4G Advent Calendar 2020 の　12/24の記事です。　・・・えええ、
すでに26日で、遅れておりますが（汗

事の始まりは、QGIS 初心者質問グループの任意の月日,時刻の日射分布図という投稿でした。こうした機能の多くは、QGISそのものよりも、GRASS GISやSAGAの機能を利用することが多いです。
今回は、GRASS GISを使って、任意の日時の日照時間、および、それを累積した1日での累積日照時間の算出方法について、説明します。

・・・すみません、QGISでスマートにやる方法がわからなかttので、GRASSでやりたいと思います（汗

なお、手順については、こちらのword ファイルにもまとめてあります。併せてごらんください。

処理の流れ

GISにより、何かの結果を得たいと考える場合、多くの場合は単一の機能して実装されておらず、複数の機能を組み合わせる必要があります。今回の場合は、ある地域について、一日のうちに何時間太陽が当たっているかを計算することとしました。GRASS GISには多彩な機能があるものの、さすがにこの値を直接出してくれる機能はありませんでした。
関連する機能としては、r.summaskがあります。これは、

太陽の位置と高度のラスターマップから影の領域を計算します。
正確な太陽の位置(A)を指定するか、太陽の位置(B)をr.sunmask自身で計算する日付/時刻を指定します。みんなの自動翻訳で翻訳

という機能を持つものになります。
この機能を使って分析する場合の流れを、実施の手順と逆（バックキャスト、という奴）に説明すると、

ある日の日照時間を計算する ← 総時間から、r.sunmaskで計算された日陰の時間を引くことで、総日照時間とする ← r.maskである時間の影か否かを計算し、それを積算することで総日陰時間データを作成する

という手順を取ることにしました。

今回は、サンプル地域として以下の範囲について処理を行いました。

なお、DEMのデータは、国土地理院の基盤地図情報からダウンロードしたものを使っていますが、5mや10mの解像度だと処理に時間がかかるため、100mに解像度を落として行っています。解像度を落とす方法は、r.sun.incidoutの処理高速化は可能でしょうか？を参考にして下さい。

動作環境

動作確認をした環境は、

Windows 10 Pro (64bit)
GRASS GIS 7.8 (OSGeo4W 64bit版)

です。OSGeo4Wのコマンドラインを使うので、QGIS付属のGRASSでは動きませんでした ¹

データ処理

GRASSへのデータのインポート

さて、実際の処理を行う前に、まず初めにデータをGRASS GISのデータベースにインポートする必要があります。これは、GRASS GISがGeoTiffやSpatiaLiteといったファイルを直接扱うのではなく、独自のデータベースにデータを格納し、処理するものだからです。GRASSを使おうとして初めに躓くのがここなのですが、ここを乗り越えれば、全行程の半分は突破したものと思って下さい！
重要な概念は、以下の2枚のスライドの通りになります。（FOSS4G Tokyo 2011 GRASS GISハンズオンセッションよりこぴぺ）

具体的な手順については、word ファイルにステップbyステップで書いてありますので、そちらをご覧下さい。

r.sunmaskの実行

つぎに、実際にr.sunmaskを実行してみます。さてここで、Mapsetの変更を行います。これまでの手順でDEMをインポートしましたが、その際のMapsetは、PERMANENTマップセットになっています。ですがこの先、複数のラスタデータ作成する予定です。そのため、これをPERMANENTマップセットとは別のマップセットに作成することにします。ここでは、sunny_shade と言うマップセットを作成し、そこで作業することにします。

実行には、Layer Manager ウィンドから、

設定→GRASS作業環境→マップセット作成　[g.mapsets]

を選択します。

表示されたg.mapsetウィンドの「Name of mpaset」に「sunny_shade」と入力します。

次に「作成」タブをクリックし、「マップセットがない場合はマップセットを作成」にチェックを入れます。

その上で実行をクリックすると、下の画面のようにマップセットの作成が完了します。

作成が完了したら、閉じるをクリックしてください。これで、これから作業するデータは、「sunny_shade」に作成されます。

次に、r.summaskで日向、日陰を計算します。実行には、Layer Manager ウィンドから、

ラスター→太陽照度と影→陰影図　[r.sunmask]

をクリックします。

以下の画面が表示されるので、まず、計算に用いる標高マップを選択します。ここでは、PERMANENTマップセットに含まれるmerge_50を選択します²。

次に、テーマタブをクリックします。ここで、計算する日時を指定します。例えば、日本における2020年5月1日午前6時を指定する場合、以下のようになります。この際、一度、Minutesの所に「0」を入力して下さい（何か数字を入力しないと、上手くいかないようです）。また、一番最後のタイムゾーンに「9」を入力するのを忘れないでください

最後に、オプションタブを選択し、出力する地図を選択します。ここでは、「SS_20200501_06」としました。

計算が上手くいけば、以下のように表示されると思います。

このうち、黄色いところが影になっているところです。また、Layer Managerのチェックボックスを外すと、

のように表示されます。
次に、12時のマップを作成してみます。テーマの時間を12、出力ファイル名をSS_20200501_12とします。上手く動くと、以下のようになります（SS_20200501_06はチェックを外して表示しません）。
一見すると、何もないように見えますが、これは日陰がないという事なので、問題ありません。

18時について計算すると以下の通りです。

改めて3枚の地図を並べると

となり、太陽の位置に合わせて影の領域が変化していることが分かります。

複数の日陰マップを一括して作る

さて、これである時間の日陰の地図を作ることができました。次にこの地図を毎分ごとに作成する必要があります。2020年5月1日の高知の日の出は5:17、日の入りは18:49（国立天文台 > 暦計算室 > 各地のこよみ > 高知(高知県)参照）となります。そうすると、5時から19時まで、毎分ごとのデータを作ろうとすると、13時間×60分＝780回も上記の入力をする必要があります。それは、なかなかキツいです。そこでこの計算を自動的に何度も行ってもらうことが必要になります。この際に利用できるのが、黒い画面、コマンドウインドです。

そこに行く前に、まず、r.sunmuskウィンドにあるコピーボタンをクリックし、Windowsであれば、メモ帳等にペーストしてみて下さい。すると、以下の画面のような文字列が貼り付けられると思います。

テキストにすると、以下の通りです。

r.sunmask elevation=merge_50@PERMANENT output=SS_202005010601 year=2020 month=5 day=1 hour=06 minute=01 timezone=9

これが、GRASSがr.sunmuskを実行する際に、実際に行っている命令になります。この命令をそのままコマンドウィンドにペースとしてリターンを押すと、以下の通り実行が可能です（赤枠をご覧下さい）。

ですので、上記のコマンドのhour=12 minute=0が自動的に変わり、さらに出力するマップ名SS_202005010601の時・分を示す0601の数字が自動的に変わるようにプログラムを書けばいいことになります。今回はLinuxのshellscriptの機能を用いて、この機能を持つスクリプトを作成しました。³

sum_mask_hou_min.sh

# 5時から9時までの計算
for hou in {5..9}
 do
  # 0分から9分までの計算
  for min in {0..9}
   do
    r.sunmask elevation=merge_50@PERMANENT output=SS_202005010${hou}0${min} year=2020 month=5 day=1 hour=${hou} minute=${min} timezone=9
   done
  # 10分から59分までの計算
  for min in {10..59}
   do
    r.sunmask elevation=merge_50@PERMANENT output=SS_202005010${hou}${min} year=2020 month=5 day=1 hour=${hou} minute=${min} timezone=9
   done
done

# 10時から18時までの計算
for hou in {10..18}
 do
  # 0分から9分までの計算
  for min in {0..9}
   do
    r.sunmask elevation=merge_50@PERMANENT output=SS_20200501${hou}0${min} year=2020 month=5 day=1 hour=${hou} minute=${min} timezone=9
   done
  # 10分から59分までの計算
  for min in {10..59}
   do
    r.sunmask elevation=merge_50@PERMANENT output=SS_20200501${hou}${min} year=2020 month=5 day=1 hour=${hou} minute=${min} timezone=9
   done
  done

実際のファイルは、こちらです(sun_mask_hou_min.shへのリンク)
この内容を、sum_mask_hou_min.shとして、c:\GRASSフォルダに保存して下さい。このコマンドを実行するには、

C:\c:
C:\cd c:\grass
C:\g.region res=100
C:\sh sum_mask_hou_min.sh

と入力します（入力する際は、c:\は除いて下さい）。成功すると、以下のように毎分ごとの日陰マップが生成されます。

日照時間マップの作成

残りは、わずかです！　ここまでに作成した分ごとの日陰地図から、1日の間の日照時間マップを作成します。考え方としては、毎分日陰マップを作成した5時から18時、つまり780分のうち、何分間、ひなたになっているかを計算することにます。こうした、地図間の計算を行う場合、GRASS GISでは、r.mapcalcというコマンドを使います（r.mapcalcのマニュアル）具体的な計算の内容は、

全時間分のデータを持つ地図を作る → 毎分後のと日陰地図を引く

という流れです。つまり、日陰地図には、日陰の所に1という数字が入っています。これを全時間（780分）から引くことで、残りの時間が、その日の日照時間という事になります。この計算についても、日陰地図を作ったときと同様に、スクリプトを作成して行いました。以下、スクリプトの内容です。

sunny_min.sh

# 日陰時間を計算する元ファイルを作成
r.mapcalc expression="shade_min=0"

# 総日陰時間を計算
# 5時から9時までの計算
for hou in {5..9}
  do
   for min in {0..9}
    do
     # null値を0で置き換え
     r.null map=SS_202005010${hou}0${min}@sunny_shade null=0

     # tmp_mapに日陰時間を加える
     r.mapcalc expression="tmp_map=shade_min+SS_202005010${hou}0${min}@sunny_shade"

     # 古い日陰時間マップを削除
     g.remove -f type=raster name=shade_min

     # tmp_mapを新しい日陰時間マップに変更
     g.rename raster=tmp_map,shade_min
    done
   for min in {10..59}
    do
     r.null map=SS_202005010${hou}${min}@sunny_shade null=0
     r.mapcalc expression="tmp_map=shade_min+SS_202005010${hou}${min}@sunny_shade"
     g.remove -f type=raster name=shade_min
     g.rename raster=tmp_map,shade_min
    done
   done

# 10時から18時までの計算
for hou in {10..18}
 do
   for min in {0..9}
    do
     r.null map=SS_20200501${hou}0${min}@sunny_shade null=0
     r.mapcalc expression="tmp_map=shade_min+SS_20200501${hou}0${min}@sunny_shade"
     g.remove -f type=raster name=shade_min
     g.rename raster=tmp_map,shade_min
    done
   for min in {10..59}
    do
     r.null map=SS_20200501${hou}${min}@sunny_shade null=0
     r.mapcalc expression="tmp_map=shade_min+SS_20200501${hou}${min}@sunny_shade"
     g.remove -f type=raster name=shade_min
     g.rename raster=tmp_map,shade_min
    done
 done

# 総時間から総日陰時間を引いて、日向時間マップを作成
r.mapcalc expression="sunny_min=780-shade_min"

この内容を、sunny_min.shとして、c:\GRASSフォルダに保存して下さい。
実際のファイルは、こちらです(sunny_min.shへのリンク)

このコマンドを実行するには、コマンド画面に

C:\sh sum_mask_hou_min.sh

と入力すればOKです。計算した結果は、以下の通りとなります。

この画像をQGIS等で見たい場合には、GeoTiffとしてエクスポートするといでしょうか。Layer Manage ウィンドから、

ファイル→ラスターマップのエクスポート→一般的なアクスポートフォーマット　[r.out.gdal]

をクリックして下さい。

ここで、

Name of raster map → sunnmy_min
Name for output raster file → C:\GRASS\sunny_min.tif
Raster data format to write → GTiff

とします。

これで実行を押すと、GRASSで作られたデータがGeoTiffファイルとして出力されます。実際のファイルはこちらです（sunny_min.tifへのリンク）

上記のファイルをQGISに読み込み、彩色等を少し修正すると、以下のように表示できます。レイヤを選択し手情報を表示すると、日照時間が分で表示されることが確認できます（右上）。

以上で、手順の説明は終わりです。

おわりに

さて、今回の記事はこれで終わりですが実はやろうとして出来なかったことがあります。

生成した日陰マップを連続してエクスポートし、アニメーションGIF化する
- エクスポート自体は出来るのですが、何故かカラーテーブルの設定がおかしく、あんまり美しくなくなってしまう
生成したマップではなく、GRASSのマップウィンドに表示されたデータをエクスポートする
- これは、エクスポートできなかった。Windowsだからダメなのかしら？

このへんは、誰かがっつりやってくれるとありがたいですw

また、今回は自分が使い慣れている環境と言うことで、GRASS GIS + Shell Scriptでの処理について説明しました。FOSS4Gのいいところは、色々な選択肢がある事です。QGISでダメならGRASS GISで。GRASSでもダメならPostGISで等など、目的の応じてツールの使い分けることで、様々なことが実現できると思います。なので、これをお読みになって頂いたかにも、出来れば様々なFOSS4Gにさわってみて、新しい使い方を発見して頂ければと思います。

$\huge{それではみなさん、よいお年を！}$ ⁴

私がよく知らんだけかも... ↩
すみません、100mにリサンプルする前です ↩
桁数をそろえるために、かなり美しくない実装になっています。WindwosのOSGeo4Wでは、私ではここまでが限界でした。誰か、美しい実装を・・・。 ↩
アドベンドカレンダーだったのではないか・・・・？ ↩