系統解析備忘録: 3月 2022

前回までの記事で琵琶湖南湖の湖沼図から地形に関するラスターデータを作成しました。今回はこれを利用してMaxEnt解析を試みます。MaxEntはオープンソースのソフトでここからダウンロードできます。なお、ここではバージョン3.4.4を使用しています。サンプルとして、ブルーギルの分布調査¹から抜き出した2018年の分布データを使ってみます。

MaxEnt解析に必要なものはターゲットである種の分布データと地形・植生などその種の分布に関係があると思われる地理空間情報です。

1. 分布データの準備

ダウンロードした分布調査結果の図からブルーギルの分布情報（2018年6月）を抽出し、csv形式で保存しました（ファイルはこちら）。正確なものではありませんが、今回はとりあえずMaxEntを動かすことを目的としているので良しとします。

2. 地理空間情報データの準備

今回は地理空間情報データとしてこれまでの記事で作成した深度・斜度・底質のラスターデータおよびWorldClimからダウンロードした6月（分布調査が6月だったため）の平均気温・風速データを用います（空間解像度: 30 secondsを使用）。

地理空間情報データをMaxEntに読み込ませるためには、それぞれのラスターデータの範囲とセルサイズを一致させた上で、Esri ASCII ラスター形式（.asc）またはDIVA-GISグリッド形式（.grd）にする必要があります。

2.1. データの範囲

ラスターデータから特定の範囲を抜き出すための外枠を用意する必要があります。今回は琵琶湖が舞台なので、湖沼の水域データが活用できそうです。国土交通省の国土数値情報ダウンロードサイトからダウンロードします。

ダウンロードしたファイルを解凍し、QGISで「ベクタレイヤを追加」からシェープファイル（.shp）を読み込みます。QGISに日本の湖沼の水域が表示されます。他の湖沼は邪魔なので、まずは琵琶湖だけを別のレイヤに移します。「シングルクリックによる地物選択」をクリックし、琵琶湖の上でクリックします。すると図1のように琵琶湖全域が選択されます。

図1

続いてメニューの「編集」>「地物をコピー」から「編集」>「新規レイヤへの地物貼り付け」>「一時スクラッチレイヤ」をクリックし、レイヤ名を入力して「OK」を押します（ここでは「LakeBiwa」とした）。これにより琵琶湖だけを別レイヤにコピーできました。最初に読み込んだベクタレイヤを非表示にするか削除し、QGIS上に琵琶湖だけを表示します。

次に琵琶湖のうち、解析に使う南湖のみを切り出します。今はまだ琵琶湖全体が一つの地物となっているので、これを複数に分割します。レイヤを選択した上で、「編集モード切替」をクリックし、「地物を分割」をクリックします（図2）。

図2

この状態で、切りたい部分の開始点をクリックし、終点をクリックすると切断箇所が直線で表示されます。最後に右クリックすることで切断されます（図2の直線は切断後）。切断ができたら「部分の削除」をクリックし、消したい方の地物をクリックすることで地物が削除されます。これにより必要な領域を切り出します（図3）。

図3

これで完成に見えるのですが、このままだとラスターデータを切り抜く際にエラーとなります。これを回避するため、メニューの「ベクタ」>「ジオメトリツール」>「マルチパートをシングルパートに変換」をクリックし、「入力レイヤ」にさきほど切り出した琵琶湖南湖のレイヤを選択して実行します。「出力レイヤ」という名前のレイヤができるので、これを右クリックし、「エクスポート」>「新規ファイルに地物を保存」でシェープファイル形式（.shp）で適当な名前を付けて保存します（ここではSouthLakeBiwa_boundary.shpとした）。またこの際、「座標参照系」はラスターデータと一致させたほうが良いかと思います。

2.2. 地理情報空間データの編集

上で作成した琵琶湖南湖の範囲データを使用してラスターデータを切り抜きます。まずダウンロードした平均気温データ解凍し、「ラスタレイヤを追加」から読み込みます（データは.tif形式）。そして、メニューの「ラスタ」>「抽出」>「マスクレイヤで切り抜く」をクリックします。入力画面が表示されるので、「入力レイヤ」に平均気温データのレイヤ、「マスクレイヤ」に琵琶湖南湖の範囲データ（ここではSouthLakeBiwa_boundary）を選択し、「実行」をクリックします。完了すると「出力ファイル」という名前のレイヤに切り抜いたラスターデータが出力されます（図4）。このレイヤを右クリックし、「エクスポート」>「名前を付けて保存」でファイル形式をGeo TIFF形式（.tif）で保存します（ここではclipped_tavg06.tifとした）。この段階ではまだ（.asc）形式にはしないでください。

図4

繰り返しになりますが、MaxEnt解析を行うには全てのラスターデータの範囲とセルサイズを一致させなければなりません。そこで、先ほど作成したclipped_tavg06.tifを基準として他のラスターデータを変換していきます。QGISには「ラスタを揃える」という機能があり、これが使えるかと思ったのですが、なぜかエラーで動きませんでした。少々面倒ですが、Rを使って変換をしていきます。方法は以下の通りです。

>library(raster)　#使用するパッケージです。なければinstall.packages("raster")でインストールしてください。
>library(rgdal)　#同上。

>s = raster("clipped_tavg06.tif")　#基準とするラスターファイルの読み込み。
>ex = extent(s)　#基準ラスタの範囲を抽出。
>writeRaster(s, "clipped_tavg06.asc", "ascii")　#ここで基準としたclipped_tavg06.tifをEsri ASCII ラスター形式（.asc）で保存。

>r = raster("wc2.1_30s_wind_06.tif")　#変換するラスターファイルの読み込み。
>r = crop(r, ex)　#変換するラスターを基準ラスタの範囲で切り抜く。
>resample(r, s, "ngb", filename="clipped_wind06.asc", format="ascii")　#基準ラスタのセルサイズに調整。"ngb"は調整方法で最近傍法を意味する。他に"bilinear"も選択できる。また、filenameとformatを指定することで、ASCII ラスター形式で結果をファイルに保存している。

上記のうち、下の3行を繰り返すことで、他のラスターファイルも同じようにして範囲とセルサイズを調整します。完了したら作成した.asc形式のファイルをメモ帳などで開き、先頭の5行（NCOLS, NROWS, XLLCONER, YLLCORNER, CELLSIZE）が完全に一致しているか確認してください。できていれば新しいディレクトリを作り、これらのファイルをそこに移動させます。

3. MaxEnt解析の実行

MaxEntを起動し、サンプルに分布データ（.csv）、環境レイヤに.ascファイルを入れたディレクトリを指定します。下の枠にインプットしたデータのリストが表示されます（図5）。

図5

このうち、環境レイヤの底質データ（ここではclipped_sediment）はカテゴリカルデータなので、「Categorical」に設定します。他は「Continuous」のままで大丈夫です。最後に出力ディレクトリを指定し、「RUN」をクリックします。問題なければこのディレクトリに結果のファイルが出力されるはずです。このうちの.htmlファイルを開くと、モデルの評価、環境データごとの寄与度、予測された環境適地の地図などが確認できます。また、出力された.ascファイルをQGISに読み込むことで生息適地の地図を表示することもできます（図6）。

図6

参考文献

1. 田口貴史・石崎大介・根本守仁、令和元年(2019 年)6 月の琵琶湖南湖におけるブルーギルの分布、外来魚駆除対策研究（https://www.pref.shiga.lg.jp/file/attachment/5238608.pdf）

前回の記事では湖沼図の等深線を基に内挿を行うことで水深ラスターデータを作成しました。今回は同じ湖沼図から底質ラスターデータを作成します。

底質データの作成は以下の4つのステップからなります。

湖沼図の入手
ジオリファレンス
底質データの書き込み
内挿

このうち1と2については、前々回の記事を参照してください。ここでは3.底質データの書き込みと4.内挿について説明します。

3. 底質データの書き込み

QGISにジオリファリンス済みの地図を表示してください。等深線を書き込んだときは線を用いましたが、底質データは湖沼図上ではある地点の底質をアルファベットで示しているので点で書き込みます。なお、それぞれのアルファベットの意味はここにまとめました。

やることは単純で、各アルファベット上に点を打ち、底質の情報を入力するだけです。

まずは点を打つレイヤを準備します。画面左上にある「新規シェープファイルレイヤ」をクリックします（フロッピーディスクのようなアイコンの真下にあるやつ）。入力画面が開くので、ファイル名を入力し（ここでは、「sediment.shp」とした）、ジオメトリタイプで点（ポイント）を選択します。そして、新規フィールド欄の名前を「Sediment」とし（任意ですが、ここではSedimentとしました）、データ型を「テキストデータ」にしてから「フィールドリストに追加」をクリックします。すると下のフィールドリストに「Sediment」というフィールドが追加されます。下の「OK」をクリックすると「sediment」というレイヤが追加されます。

追加されたレイヤを選択し、「編集モード切替」をクリックし、続いて「点地物を追加する」をクリックします（図1）。

図1

その状態で地図上の底質を示すアルファベットにカーソルを合わせてクリックします（例えば図1の赤い円の中の「M」）。入力画面が現れるので、「Sediment」にその地点のアルファベットを入力してOKをクリックします。すると地図上に点が表示されます。これを繰り返して地図上の全ての底質データをレイヤに入力します。

4. 内挿

前回は水深という数値データを基に内挿を行いましたが、今回は文字データを基に内挿を行います。そのため、今回は「Multilevel b-spline interpolation for categories」を用います。「プロセシングツールボックス」の検索窓にこれを入力して選択します。

パラメータ入力画面が現れるので、以下のように設定して「実行」を押します（図2）。

Points: sediment（レイヤ名）
Attribute: Sediment （内挿を行う属性）
Output extent: 「レイヤから計算」>地図のレイヤを選択（内挿の範囲指定）
Cell size: 0.001（数字が小さいほど解像度が高い）
Fit: 0（よくわからないが、試した範囲では「1」にしても結果はほぼ変わらなかった）
Categories: sediment_cell001（アウトプットファイル名。任意）
Probability: sediment_cell001_prob（確率のアウトプットファイル名。任意）
図2

完了すると図3のような内装された底質レイヤが得られます。以上で底質ラスターデータの作製は完了です。

図3

＊エラーが出て内挿ができない場合

QGISは内部でSAGAというプログラムを呼び出して「Multilevel b-spline interpolation for categories」を実行しています。しかし、PCにインストールされているSAGAのバージョンがQGIS推奨のバージョンと一致しないとエラーになることがあります。使っているOSがLinuxの場合、以下のコマンドで直接SAGAによる内挿を実行できます（ここでは上記のパラメータ設定と同じ設定をしています）。コマンドを実行してできたファイルのうち拡張子が.sdatとなっているものをラスターとしてQGISに読み込ませれば結果を見ることができます。

saga_cmd grid_spline 7 -POINTS Biwa_Sedimentpoint_shape.shp -FIELD Sediment -TARGET_USER_SIZE 0.001 -TARGET_USER_XMIN 135.853200 -TARGET_USER_XMAX 136.006200 -TARGET_USER_YMIN 34.974000 -TARGET_USER_YMAX 35.173800 -TARGET_USER_FITS 0 -CATEGORIES sediment_cell001 -PROPABILITY sediment_cell001_prob

コマンド先頭の「saga_cmd grid_spline 7」はこの解析が「Multilevel b-spline interpolation for categories」であることを指定しているので変更する必要はありません。「-POINTS」は内挿を行うシェープファイル、「-FIELD」はQGISの入力画面で「Attribute」とされているもの、「-TARGET_USER_SIZE」は「Cell size」、「-TARGET_USER_XMIN」から「-TARGET_USER_YMAX」はQGISの「Output extent」で指定したものを別々に指定、「-TARGET_USER_FITS」は「Fit」、「-CATEGORIES」は「Categories」、「-PROBABILITY」は「Probability」をそれぞれ意味しています。一度QGISで内挿を走らせるとログが出力されるのでそれを参考にすると設定しやすいと思います。

系統解析備忘録

2022年3月22日火曜日

MaxEnt解析手順メモ