とらりもんHOME  Index  Search  Changes  Login

とらりもん - 傾斜と方位・日射分布 Diff

  • Added parts are displayed like this.
  • Deleted parts are displayed like this.

筑波大学農林工学系 西田顕郎奈佐原顕郎
!傾斜と方位・日射分布

ひきつづき、数値地形モデルを処理しながら、GRASSの重要な概念を学ぶ。

以下、GRASS:~>はGRASSのコマンドラインのプロンプトのしるしであり、ユーザーは打ち込む必要は無い。(その後を打ち込むこと)

[[前回|数値地形モデル(GTOPO30)のダウンロードと読み込み]]のようにGRASSを起動せよ。


!傾斜と方位
GRASSには、数値地形モデル(DEM)から斜面傾斜・方位を求めるコマンドがある。日本の中部山岳地帯をd.zoomで選んで、それらのコマンドをためしてみよう。

    GRASS:~> d.rast GTOPO_E100N40
    GRASS:~> d.zoom

(日本の中部山岳地帯をマウスで選ぶ)

    GRASS:~> r.slope.aspect elevation=GTOPO_E100N40 slope=GTOPO_E100N40_slope aspect=GTOPO_E100N40_aspect
    GRASS:~> d.rast GTOPO_E100N40_slope
    GRASS:~> d.rast GTOPO_E100N40_aspect

{{attach_view(JP_GTOPO_slope.png)}}
{{attach_view(JP_GTOPO_aspect.png)}}

    左:傾斜(GTOPO_E100N40_slope)、右:アスペクト(GTOPO_E100N40_aspect)

傾斜は水平を0度、垂直を90度として斜面の傾きを表す。斜面方位(aspect)は、東を0度、北を90度、西を180度、南を270度とする反時計まわりの度数表示で斜面の向きを表す。

!斜面の受け取る日射量
上で行った、傾斜と方位の計算結果を、r.mapcalcで処理すれば、各地の斜面が太陽から受ける日射量を求めることができる。日射量の分布は、農林業や生態系管理・水文水資源管理の上で重要な情報であるが、GISを使えば、日射量に与える地形の影響を簡単に計算できる。

まず、原理を説明しよう。いま、日射が太陽方向からのみやってくるとし、周囲からの反射光の影響を考えないとすると、斜面が単位面積当たりにうけとる日射の量は、斜面の単位法線ベクトルnと太陽方向の単位ベクトルSの内積に比例する。

傾斜をθ、斜面方位をφとし、太陽の天頂角をθs、太陽の方位をφsとすると、

n = [ sinθ cosφ, sinθ sinφ, cosθ ]
S = [ sinθs cosφs, sinθs sinφs, cosθs ]

となる。従って、地形データから各地点のnを計算し、一方、所与の太陽位置(高度・方位)からSが計算し、 nとSの内積を計算すれば、各地表面の受け取る日射量が相対的に推定できる。この計算(上記の内積計算)は、GRASSの上で行なうことができる。ここでは、太陽の位置(θsとφs)をあとから任意に与えることができるように、シェルスクリプトを組んでみよう。viなどのテキストエディタで、rad.shという名前で以下のようなファイルを作る:

    GRASS:~> vi rad.sh

    #!/bin/sh
    r.mapcalc "nx = sin(GTOPO_E100N40_slope)*cos(GTOPO_E100N40_aspect)"
    r.mapcalc "ny = sin(GTOPO_E100N40_slope)*sin(GTOPO_E100N40_aspect)"
    r.mapcalc "nz = cos(GTOPO_E100N40_slope)"

    r.mapcalc "Sx = sin($1)*cos($2)"
    r.mapcalc "Sy = sin($1)*sin($2)"
    r.mapcalc "Sz = cos($1)"

    r.mapcalc "rad = nx*Sx + ny*Sy + nz*Sz"
    r.colors map=rad color=grey
    d.rast rad

これを保存し、いつものように実行権限を与える:

    GRASS:~> chmod +x rad.sh

そして、以下のように実行してみよう:

    GRASS:~> ./rad.sh 60 0    (太陽天頂角60度、つまり太陽高度30度で、方位角0度、つまり真東から日があたったとき: 下図左)
    GRASS:~> ./rad.sh 45 270  (太陽天頂角45度、つまり太陽高度45度で、方位角270度、つまり真南から日があたったとき: 下図右)

{{attach_view(JP_GTOPO_sunE.png)}}
{{attach_view(JP_GTOPO_sunS.png)}}


このように、r.mapcalcは、ラスターデータに関するかなり複雑な処理や演算もやってくれる。GRASSでラスターデータを解析する場合、このr.mapcalcコマンドを駆使することが重要である。

注: 試しに、d.soomコマンドを使って、表示範囲を拡大してみよう:

{{attach_view(JP_GTOPO_sunS_zo.png)}}



このように、中部山岳地帯以外の場所は空白になっていることがわかる。一般に、r.mapcalcを使った解析は、現在のregionの中でしか行われない。したがって、解析範囲を拡大したいときは、g.regionコマンドやd.zoomコマンドで、regionを適切に設定してからやり直す必要がある。(GRASSが現在のregionだけで解析を行うのは、処理時間を節約するためである)

注: 上のスクリプトでは、最終的なラスターマップ"rad"を作る前に、nx, ny, ...という6つのラスターマップを中間的に作っている。これは、みなさんにとって理解しやすいように、あえて段階をわけて処理したかったからである。本質的には、上の処理は、以下のようなスクリプトに集約できる。むしろこちらのほうが、中間データ(nx, ny, ...) を作らないぶん、処理は速いだろう。

    #!/bin/sh
    r.mapcalc "rad = sin(GTOPO_E100N40_slope)*cos(GTOPO_E100N40_aspect)*sin($1)*cos($2) + sin(GTOPO_E100N40_slope)*sin(GTOPO_E100N40_aspect)*sin($1)*sin($2) + cos(GTOPO_E100N40_slope)*cos($1)"
    r.colors map=rad color=grey
    d.rast rad

!実習課題
1. GTOPO30のページ(わからない場合は検索エンジンで検索して見つける!)から、北海道を含む領域をダウンロードせよ。

2. それを用いて、北海道の地形を適当なカラーテーブルによるカラーグラデーションで表示せよ。

3. 北海道に、南西から太陽高度30度で日射が当たったときの様子(日射量分布)を表示せよ。

!!!!!(ブラウザの「→ [[GIS入門]]に戻る」ボタンで戻ってください)

!!更新履歴
2011/11/08 筑波大学 田中健太郎