GRASS GIS r.sunで日射量を計算したい

[sho_he...@yahoo.co.jp]

GRASS GIS r.sunというツールで日射量の計算が可能らしく、試していますがうまくいきません。

QGISは2.4.0です。古いです。3.4.0でも試しましたが同様の結果となりました。
必要なデータは標高、傾斜方位、傾斜です。標高データは10mメッシュのものを基盤数値情報から入手しGeoTIFF形式に変換しました。傾斜方位、傾斜ラスタも作成済みで、CRSはEPSG:2445です。プロジェクトも同じCRSにしています。

GRASS GIS r.sunから、標高、傾斜方位、傾斜に入力するラスタを選択して、Output global(total) irradiance/irradiation layerに出力先を入力しました。
RUNすると、「r.out.gdal.exeは動作を停止しました。」と表示されます。出力先にTIFF形式のファイル(8バイト)が作成されていますが、ワールドファイルはありません。

サイズが大きすぎるのかとも思い、小さなファイルでも試しましたがダメでした。
また、2バイト文字は出入力先ともに使っていません。

どなたかわかる人がいらっしゃったら教えてください。

[sho_he...@yahoo.co.jp]

Google groupでスマホからファイルを添付する方法がわからず、ファイルは添付できませんでした。
基盤地図情報から入手したDEMは、基盤地図情報ビューアで、シェープファイルにエクスポート→QGISでベクタのラスタ化→ラスタ解析ツールのワープ(再投影)でCRSの変更→ラスタ解析ツールのDEM(地形モデル)で傾斜と傾斜方位に変換しています。

[福岡]

こんばんは

DEMの作成方法なのですが、ちょっと違和感があります。

ダウンロードしたxmlファイルからgeotiffのDEMを作成するのにはいくつか方法があると思いますが、提示なさった方法はあまり聞いたことがありません。

エコリスさんで公開している「基盤地図情報 標高DEMデータ変換ツール」で作成したDEMで試してみてはどうでしょうか？

基盤地図情報 標高DEMデータ変換ツール

http://www.ecoris.co.jp/contents/demtool.html

[sho_he...@yahoo.co.jp]

たしかに今回のやり方だと、工数も多くて面倒でした。明日職場で試してみます。
ありがとうございます。

[adachi]

正規のはずのビューアで変換しようとするとなぜか直接ラスタにはできませんよね。なかなか使い勝手が悪い…

2.Xをお使いであれば、「fgddemImporter」というプラグインも使えるかと思います。

zipファイルから直接取り込めるのでそこは便利ですが、以前のスレッドでも議論されていたように

読み込み直後にデフォルトで緯度経度座標系になるので、平面直角座標系等に再投影する必要があります。

そして件のr.sunですが、3.4.4でやってみようとしたところ

・「ラスタ」メニューから作成した傾斜ラスタ、傾斜方位ラスタを使うとエラーが出る

・GRASSツールのr.slope.aspectで作成した傾斜ラスタ、傾斜方位ラスタを使うと（時間はかなりかかるが）ちゃんと作成される

という現象が起こりました。

同じGRASSツール内で作成したデータのほうが整合がとれる？のでしょうか。

同じ目的物を作るツールのはずなのに、できたものが違う挙動をするのは何か気持ちわるいですが…

[sho_he...@yahoo.co.jp]

ありがとうございます。

エコリスさんのツールで作成し、試しましたが同様の結果となりました。

それから、GRASSツールのr.slope.aspectで、傾斜と傾斜方位を作成しようとすると、「r.out.gdal.exeは動作を停止しました。」と表示されてしまいます。

GRASSツール自体が使えないということでしょうか？

[福岡]

こんにちは

adachiさんの情報を含めて、自分の環境で検証してみました。

QGIS3.4.4では、r.sunはご相談と同じようにエラーが発生して、正しく動作しませんでした。

また、r.slope.aspectでエラーメッセージが出るものの傾斜と傾斜方位は出力されましたが、これを使ってr.sunで処理しようとするとエラーが発生しました。

どうも自分の環境ではQGIS3.4.4ではうまく処理することができませんでした。

そこでQGIS2.18.28もインストールされているので、こちらでも試してみました。

DEMとGGIS3.4.4で出力した傾斜図と傾斜方位図を使用してr.sunを実行しました。

範囲が広すぎたためと思われますが、2時間ほど掛かりましたが実行が完了しました。出力もOKのようです。

自分の感触としては、3.4.4に何らかの問題があるように思います。安定版を使った方が良いように思いました。

[Shohei Murata]

皆さん本当にありがとうございます。

3.2.2でも試しましたが、同様の結果でした。

ちなみに、この投稿でPCからファイルを添付する方法を教えていただけますか。

実行後のエラーのログを追うと、青字部分が

「Starting GRASS GIS...

警告: マップセットの平行ロックはウィンドウズではサポートされていません。」で始まります。

その後、コードが進んで

「エラー：insol_time and incidout are incompatible options」とあります。

その後は数行続いて「エラー：ラスタマップ＜・・・＞が見つかりません」というのが繰り返されます。

[福岡]

入力欄外の左上に「ファイルを添付」とありますので、こちらをクリックすると添付できると思います。

これから打ち合わせに出かけますので、すぐに返信はできないと思いますがよろしくお願いいたします。

[Shohei Murata]

添付できない理由がわかりました。

今朝メンバー登録申請をしたばかりなので、まだメンバーになっていないからと思われます。

メンバー登録が完了したら添付できるようになるかと。

[水谷貴行]

こんにちは。水谷です。

とあるQGIS講習で作成した以下の資料（抜粋）に日射量図の作り方がありますので、ご参考にしていただければと思います。

https://docs.google.com/presentation/d/1QLJTiV-QrPeaIkZFMEqfrPbc33KaJ-n63K8y-4TgfIQ/edit?usp=sharing 

 

ただ、 QGIS3.4.2では、r.sunは実行してもエラーとなるようですが、

なぜだかグラフィカルモデラーから呼び出すと大丈夫でした。

現バージョンでは、まだGRASSとの連携が上手くいかない場合があるようです。

[Shohei Murata]

水谷様

まさに水谷様が書かれました、「［オープンデータ＋QGIS］ 統計・防災・環境情報がひと目でわかる地図の作り方」を参考に作成を試みています。
内容は非常にわかりやすく、いつも手元に置いています。

メンバーとして承認されましたので、ファイルを添付します。
もしお時間がありましたら試してみてください。

[Shohei Murata]

傾斜と傾斜方位も添付します。

[adachi]

Murataさま

アップロードいただいたファイルでの、以下の環境での動作確認です。

Windows 10, 64 bit

QGIS 3.4.4

1. r.sunツールでElevation、Aspect、Slopeをそれぞれ設定

2-1. その他はすべてデフォルト値で実行

　→ 「レイヤが生成されませんでした」エラーが発生

2-2. アウトプットファイル設定のうち、一番上にある「incidence angle raster map」のみ、右側の「…」から「出力をスキップ」に設定して実行

　→ １～２分で完了

となりました。

この結果から「incidence angle raster map」というファイルを作成するプロセスがエラーを発生させていたように見えます。

他の方の環境でも再現性があるのか、わかりませんが…

[Shohei Murata]

adachi様

原因を特定して下さってありがとうございます。
明日は出張のため、明後日職場に着いたらすぐに試してみます。

[福岡]

皆様へ

adachiさんの投稿に、なるほどと思いました。。

グラフィカルモデラーではアウトプットに指定しているもの以外がスキップされるので、処理されたということなのですね？

私も「incidence angle raster map」をスキップにしたところ、r.sunで出力できました

ただ、ログを見ると1カ所エラーメッセージがあるようです。

ERROR 6: SetColorTable() only supported for Byte or UInt16 bands in TIFF format.

メッセージは出ているものの出力は正常そうなので、今のところ無視しても大丈夫かなと思っています。

ところで今回の検証のために地理座標系のDEMから書籍の手順に従ってr.sunの出力を行うグラフィカルモデラーを作って、会社のPCの保存しています。

添付されたファイルから位置がわかりますので、明日にでも基盤地図情報からデータをダウンロードして実行してみようと思います。

確か書籍では、再投影する際にセルサイズを変更するなどの処理も行っていたと思いますので、出力にどのような差がでるのかを確認してみたいと思います。

[福岡]

皆様へ

昨日お約束していたにもかかわらず、遅くなって申し訳ありません。m(__)m

グラフィカルモデラーを作成して処理する方法について、手順書を作りましたので添付いたします。

参考にして頂けると幸いです。

[Shohei Murata]

皆さまありがとうございました。
3.2.2でadachiさまが教えてくださった方法で実行すると計算できました。

福岡さま
このような便利な方法があるとは知りませんでした。
QGIS初心者でもすぐに使うことができ、業務に非常に役に立ちます。

[sakuya_izayoi]

izayoiです。

お世話になります。

QGIS3.4では確認できていたr.sunなのですが、現在の安定板であるQGISの3.10からは消えちゃっているようです。

代わりに、r.sun.incidoutというコマンドがありました。これが代替手段になるのかなと。ヘルプをクリックするとr.sunにつながるという不思議な状態になっておりますが・・・。

(incidoutってどういう意味なんだろう？)

[sakuya_izayoi]

前回、現在のQGISではr.sunはr.sun.incidoutで対応できそうだと書いたのですが、

現在のQGISから操れるのは、一日当たりの総日射量(Wh/m^2/day)ではなく、放射照度(Wh/m^2)になっているようです。

以前のように使うなら

https://grass.osgeo.org/grass78/manuals/r.sun.html

を読むと、モード2に切り替えなければならないみたいですが、モードの切り替え方法ご存じのかた、お見えになりますでしょうか？

grassを直接回すしかない？

[福岡]

izayoi 様

良く検証していない状況での投稿です。間違っていたらごめんなさい。m(__)m

QGISの言語を英語にして「r.sun.incidout」を立ち上げて、項目を比較してみました。

そうすると「通日(1-365)」とあるのは、「No. of day of the year(1-365)」となっています。

ここでは、1月1日をから何日目という数値を設定することになっているものだと思います。

総放射照度[W.m-2]は、Global(total)irradiance/irradiation[W.m-2]。

単位には「1day」が付いていませんが、そもそも計算する日を設定していることが引っかかっています。

もしかして計算結果を比較すると同じくなるのではないか、と疑っています。

[sakuya_izayoi]

福岡様

ありがとうございます。

izayoiです。

「No. of day of the year(1-365)」の部分は、太陽の高度や、日照時間をシミュレートしている部分です。そのため、これは放射照度及び日射量どちらにも必要な値となっております。

もともと、r.sunには、モード1とモード2がありまして。

モード1が総放射照度(W×m^-2)、モード2が総日射量(Wh×m^-2×day^-1)

となっております。そのため、両者は単位が違い、当然結果が変わります。

先ほど、ありたさまが、チャットの方で

r.sun.incidout と r.sun.insoltime

に分岐していると教えていただけました。

まだ確認できていないのですが、r.sun.insoltimeが総日射量の方になるそうです。

[ありた]

こんにちは

昼間、携帯から返信しようと思ったらうまくいかなかったのでチャットで

いざよいさんに伝えましたが、こちらでも補足しておきます。

GRASS manuals r.sun にあるように、モードが1と2があります。

https://grass.osgeo.org/grass78/manuals/r.sun.html

QGIS の ChangeLog には、プロセシング上で r.sun を2つに分割したとあります。

>  [processing] split GRASS r.sun algorithm into two (fix #21637)

https://github.com/qgis/QGIS/blob/master/ChangeLog

テストしていないので、ただしいかわかりませんが、おそらく r.sun は

r.sun.incidout と r.sun.insoltime に分割され、モード2は後者と思われます。

https://github.com/qgis/QGIS/commit/ab646c4d75c6eac3c37cfa6a284191040048ffe8

「通日」はある日を起点にしたときの経過日数で、 day of year といったときは、通常

1月0日からの経過日数を指します。ですので、50を指定すると2月19日を指します。

r.sun.incidout には、通日とともに地方時を指定します。

おそらくですが、ここでいう地方時は指定経度における時刻（太陽の位置）を指します。

東経140度における地方時12時は日本標準時の11時40分になります。こうして、設定した

通日と地方時における瞬間の、単位面積当りの日射の強さ（ W/m^2 ）を計算します。

r.sun.insoltime は通日で指定した日における、朝から晩までの単位面積当りの日射量の

総量（ Wh/m^2/day ）を計算するものと思われます。

[福岡]

こんにちは

先日は検証もしないで、投稿して申し訳ありませんでした。m(__)m

そこで、私なりの検証をした結果を報告いたします。

同じデータ（標高、傾斜方位、傾斜）で、「No. of day of the year(1-365)」を「172」にして、QGIS3.4で「r.sun」を、QGIS3.10で「r.sun.insoltime」を実行してみました。

「r.sun」から出力された「Global(total)irradiance/irradiation」と「r.sun.insoltime」から出力された「glob_rad」を比較したところ、値が一致しないようです。

値の範囲ですが「Global(total)irradiance/irradiation」は「1530.83～7375.58」、「glob_rad」は「602.919～9443.17」となりました。

全くでなくとも、ほぼ同じ値を得られることを期待したのですが、なぜかうまく行きませんでした。

どこかでポカをしているかもしれませんが、もしかすると必要な情報かもしれないと思いましたので、投稿いたします。

[sakuya_izayoi]

福岡さま

izayoiです。

お世話になります。

私の方でもやってみたのですが、同じく値が近くなりませんね。

まだ全部のパラメータを確認しきれていないのですが、なんだかちょっと気になりますね。

どちらかが垂直面で、どちらかが水平面の総日射量を示している可能性があります。もう少し見ないと。

結果として出力される値としてはそんなもんかなとも。（年間の平均日射量としては3000(Wh×m^-2×day^-1)程度のオーダーになるからどっちも範囲内。一応、岐阜で学生時代に計測していたからその程度かなと。)

取り急ぎ、ご報告まで。

[omid mohammadi]

大家好，我如何通過這個軟件計算我的盆地的寬度函數？
謝謝，'

[sakuya_izayoi]

關於盆寬函數，你說的是哪個參數？是輸出文件的分辨率嗎？

2022年3月14日月曜日 5:30:58 UTC+9 omid mohammadi: