【QGISでのデータ整理】土地利用図から鳥類の個体数をマッピングする

#鳥類の個体数や種数をマッピングする前に

鳥類の個体数や種数をマッピングする際に用いられるモデルとして、一般化線形モデルや階層群集モデルが挙げられます。どちらのモデルでも、鳥類の種iの調査区jにおける個体数λijは,

Log⁡(λij)=hab ixj1 + βi1*xj2

のようにあらわされます。
xj1は調査区jにおける土地利用カテゴリを表し、xj2は調査区jの周囲の景観構造のように表されるでしょう。そして、hab i xj1はxj1についての、また、β_i1はそれぞれの説明変数についての種iの偏回帰係数です。

ここからは各グリッドにおける個体数を推定し、解析範囲全体で鳥類の個体数/種数のマッピングを試みていきます。実際に調査した区画、縦200m、幅100ｍのグリッドを解析範囲全体に発生させ、その周囲400mの開放地率を景観構造として扱います。

なお、このページでは解析については扱わず、あくまでもQGIS上での説明変数の算出・データ整理と、解析結果のマッピングに焦点を当てます。

#各グリッドの土地利用割合の算出

1. 土地利用ファイルは2km×2kmのベクタファイルで提供されるものとして、まず解析範囲全体に200m×100mのグリッドを発生させておきます。「ベクタ」の「調査ツール」からグリッドを作成できます。その後、各グリッドの面積を計算しておきます。

2. グリッドを発生させる際には注意事項があります。1. では、解析範囲よりも少し余分な範囲にグリッドを発生させましょう。解析範囲の最外縁グリッドにおいても、景観構造を計算できるようにするためです。この際に必ず、各グリッドに個別IDを格納しておきましょう。このIDがのちに非常に大事になります。

3. そして、2.で作成したグリッドから、解析範囲のグリッドを切り出します。この際に、再度各グリッドに個別のIDが振られているかを確認しましょう。

4. 1.で述べた2km×2kmの土地利用ファイルから、の土対象地利用カテゴリだけを抜き出したファイルを新たに作成します。これは、属性テーブルで対象土地利用を選択し、「選択した地物だけ保存」のチェックを入れてエクスポートすることで作成できます。

5. 1で作った200m×100mのファイルと4.のファイルを交差させます。交差させたら各グリッドの面積を計算させてあげて、属性に追加します。面積は属性テーブルから計算できます。

6. 1.のグリッドファイルに5.の属性を結合させます。結合は、1.のファイルの「レイヤプロパティ」→「結合」からファイルを選択。「＋」ボタンで結合フィールドとターゲットフィールドのどちらも共通IDを指定してやれば完成です。

7. 結合後に、グリッドの面積で各土地利用のグリッド内の面積を割り、属性に追加することで、各グリッド内の土地利用割合を計算できます。

#景観要素の算出
###ベクタファイルのラスタ化

景観要素の算出の際にはGRASSGISを使います。
ちなみにArcGISではフォーカル統計機能を用いれば算出可能なようです。

まずは、景観要素(ここでは周囲の開放地率)が格納されているベクタファイルをラスタファイルに変換します。ベクタファイルには複数の属性が含まれていますが、ここでは念のため用いる属性(周囲の開放地率)以外の属性を削除しておきます。

1. 「ラスタ→変換→ラスタ化」を選択
2. 「バーンイン値に使用するフィールド」には、周囲の開放地率を指定
3. 「燃焼する固定値」では、白く塗りつぶしたい値を設定(開放地率が高い場所を白く塗りつぶしたければ1を指定)
4. 「出力ラスターサイズの単位」では、ここでは200m×100mのセルにマッピングしたいので、「地理参照された単位」を選択
5. 幅を100、高さを200に設定
6. 出力領域は、適当に既存ファイルの範囲を選んでおく

この手順を用いることでベクタファイルをラスタファイルに変換することができます。

いい感じにラスタ化できました。
ここからは作成したラスタファイルをGRASSGISで作業するために準備していきます。

GRASSGISで開放地率を算出

1. QGIS desktop with GRASS版のQGISを起動する

2. QGISのプラグインでGRASSを有効にしておく

3. 「プラグイン→GRASS→新規マップセットを開く」

4. フォルダを指定したり、データの範囲を指定したりする

5. 「プラグイン→GRASS→マップセットを開く」でさきほど作成したファイルを開く

6. モジュールから「r.in.gdal.qgis」を選択、読み込みたいラスタファイルを選択し、「出力ラスターマップ名」を入力して実行、「出力を見る」

7. レイヤーにラスタファイルが加わる
   

8. 領域ではなるべくきれいに解析範囲を指定してあげた方が良いと思います。私は小数点まできっちり解析範囲に合うように指定しました。

9. 解像度では、同じく東西に100、南北に200を指定します。

10. 「適用」　ボタンを押す。
    

11. GRASSツールの「r.mapcalc」を指定する

12. ここでは隣接するセルを座標で指定することができ、モデルビルダーのように計算モデルを組むことができます。例えば、ラスタファイル名を「0609last」とすると、あるセルは「0609last[0,0]」とあらわすことができ、そのセルの左斜め上のセルは「0609last[-1,-1]」と表示されます。

13. この座標標記に従って、周囲400mのセルの開放地率を求めていきます。具体的に言うと、両脇に4個ずつのセルと上下に2個ずつのセルの数値を合計し、セル数で割ります。

14. こんなめちゃくちゃ汚いモデルができてしまいました...しんどい...。やっぱコードを自分で書いた方が良いかも。。点部分が赤→灰色に変化すると、そのノードが接続されたことを示します。
    

15. 最後に、「出力」にファイル名を入力します

16. GRASSGIS経由でも行うことができて(r.mapcalc)、その場合は以下のコードが上のモデルに相当します。これでもまだコード汚い。。線つなぎゲームをするよりましですが。。領域の設定をQGISでしっかりしておくことが大事で、ファイル名を数字だけで構成したりすると、Syntax errorがでてしまいます。

17. 解析を行う前には、「g.region」などでセルサイズを指定しておきましょう(「ラスターマップと揃うように領域セルを調整します」でファイルを指定)。 また、海域などの背景を考慮するために、NULL値の設定を「r.null」で設定しておきます(NULL値を0に置き換え)。

18. 境界部分はNULL値を含んでしまうためか、計算してくれないので境界部分に予め余裕を持たせておくといいと思います。ただし、解析範囲の全てのセルに値を含める必要はありません。

19. 最後に作成したラスタマップをベクタに変換し、土地利用割合を属性に含んだベクタファイルと交差させれば完了です。この際に微小ポリゴンが発生してしまうので、ある一定面積以下のポリゴンは削除しましょう。お疲れ様でした。

("ファイル名"@"マップ名" [-2,-4]+ #セルの指定
"ファイル名"@"マップ名" [-1,-4]+
......
"ファイル名"@"マップ名" [1,4]+
"ファイル名"@"マップ名" [2,4])
/44*100 
#セル数で割って平均値を新しいラスタに格納、
#ベクタに変換しなおすときに便利なように値を100倍しておく