1. モチベーション
野生動物の調査としてGPSを装着し位置情報を収集する場合があるが、得られたデータをできれば、点で表示するだけでなく、行動圏範囲の推定など高度な解析利用にも使いたい。
そのための神パッケージ”adehabitat"がRには以前から存在し、その界隈の人たちに利用されていた。
やがて時代の波はRにおいて、位置情報はspオブジェクトで取り扱うようになり”adehabitat"はCRANから消えさり、spオブジェクトを取り扱う”adehabitatHR”に進化した。しかし、時代の波は非情にもspオブジェクトも今度はsfオブジェクトに取り変わられようとしている。
また”adehabitatHR"を使うときには、”maptools"を組み合わせることで、adehabitatHRが吐き出した解析結果のspオブジェクトを他のGISアプリケーションで標準的に使うshpファイルに変換させ、例えばQ-GISなどで表示可能になるため、連携して利用されてきた。
（Q-GISに結果が表示できなければ使い勝手がかなり悪い）
しかし、いまや時代の波はspオブジェクトとGISの架け橋であった”maptools"もCRANから消え去った・・
そこで、今、Rの地理空間情報界の王者として君臨するsfパッケージと連携し、”adehabitatHR"をGISで利用できるように試みた。
環境
R version 4.4.1 (2024-06-14) -- "Race for Your Life"
OS: windows11 or ubuntu22.04LT どちらでもOK
2．DATAの概要
library(adehabitatHR)
library(tidyverse)
library(sf)
rm(list = ls())
toni = read.csv("https://raw.githubusercontent.com/kenkenvw/data/main/toni.csv")
toni
2005-06年に南アフリカのクルーガー国立公園で、1頭のバッファローのバッファローに装着されたGPS首輪からの位置情報である
Source
MoveBank http://www.movebank.org
Name: Kruger African Buffalo, GPS tracking, South Africa
まず、csvをsfオブジェクトに変換してGPSデータの概要を見てみる。
位置情報を抜き出して、CRSを指定する。WRS84のEPSGコードである4326を指定する。
d_sf = st_as_sf(toni, coords = c("long", "lat"), crs = 4326)
d_sf
Simple feature collection with 6371 features and 3 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 31.71222 ymin: -24.35659 xmax: 31.93429 ymax: -24.13321
Geodetic CRS: WGS 84
First 10 features:
X id timestamp.utc geometry
1 17930 toni 2005-08-23 06:35:00.000 POINT (31.75345 -24.1695)
2 17931 toni 2005-08-23 07:34:00.000 POINT (31.73884 -24.15402)
3 17932 toni 2005-08-23 08:34:00.000 POINT (31.73969 -24.15359)
4 17933 toni 2005-08-23 09:35:00.000 POINT (31.73874 -24.15329)
5 17934 toni 2005-08-23 10:34:00.000 POINT (31.73946 -24.15336)
6 17935 toni 2005-08-23 11:35:00.000 POINT (31.73898 -24.15363)
7 17936 toni 2005-08-23 12:35:00.000 POINT (31.73947 -24.15346)
8 17937 toni 2005-08-23 13:35:00.000 POINT (31.74052 -24.1545)
9 17938 toni 2005-08-23 14:36:00.000 POINT (31.74146 -24.16009)
10 17939 toni 2005-08-23 15:35:00.000 POINT (31.74665 -24.17546)
位置情報は、geometry列にまとめられている。
spオブジェクトではこれまで位置情報はX,Y座標に区別されていたが、sfではまとめて情報が保持される。
ではプロットしてみよう。ggplotではgeom_sf()でプロットが簡単にできる。
ggplot(d_sf) + geom_sf()+theme_bw()
CRSを平面直角座標系であるWGS84/UTMに変換する。ここでは南アフリカの36Sゾーンを指定する。
d.utm_sf = st_transform(d_sf, crs = 32736)
UTM変換もGISソフト上で実施すると、意外と手間がかかるものだがRでは一発である。
そのほか、ここでは割愛するが、このsfパッケージはQ-GISなどで行っていたバッファーゾーンの発生等いろいろなポリゴン処理が簡単にできるので、非常に便利。
3．adehabitatHRを使った行動圏の解析
では、本題のadehabitatHRを使ってカーネル密度推定法を用いた行動圏推定を行う。
\hat{f}_h(x)=\frac{1}{Nh}\sum_{i=1}^NK(\frac{x-x_i}{h})
Kはカーネル関数（通常、ガウス関数）
ｈはsmoothing factor とよばれるパラメータ。
値を大きくすればするほど、なめらかな図形になる。
ここでは、行動圏推定にカーネル密度関数を適用することについては、これ以上ふれません。
次の関数により、動物のカーネル行動圏が推定できる。
kernelUD(xy, h = "href")
h = "href" とする以下の式から h が計算される。
$h = σ × n^{1/6}$
$σ = 0.5 × (σ_x + σ_y)$
n は GPS ポイント数、σx と σy はそれぞれ x 座標と y 座標の標準偏差。
（上記は2変量正規分布を仮定する。hの決め方には、もう一つ"lscv"がある）
adehabitatHRは、spオブジェクトを用いるため、sfオブジェクトから変換する必要がある。
d.utm_sp = as(d.utm_sf, "Spatial")
カーネル密度推定（KDE)を使う
UD = kernelUD(d.utm_sp, h="href")
kernelUD(d.utm_sp, h = "href") で警告がありました:
xy should contain only one column (the id of the animals)
id ignored
警告が出ましたが、これは無視してOK。
動物の50%行動圏と95%行動圏を推定する。
HR50 = getverticeshr(UD,50)
HR95 = getverticeshr(UD,95)
#再びsfへ変換
HR50sf = st_as_sf(HR50)
HR95sf = st_as_sf(HR95)
プロットしてみる
ggplot(HR95sf) +geom_sf() + geom_sf(data=HR50sf) +theme_bw()
動物のコア行動圏と95%行動圏の推定が完了。
Q-GISで表示可能なESRI shapeファイルに変換し、保存する。append=TRUEにしておく。
st_write(HR50sf, "toni50.shp", append = TRUE)
st_write(HR95sf, "toni95.shp", append = TRUE)
Updating layer `toni50' to data source `toni50.shp' using driver `ESRI Shapefile'
Writing 1 features with 2 fields and geometry type Polygon.
Updating layer `toni95' to data source `toni95.shp' using driver `ESRI Shapefile'
Writing 1 features with 2 fields and geometry type Polygon.```
さらに位置座標をUTM座標に変換したsfオブジェクトをX-Y座標別に分けて、csvにアウトプットする。（Q-GISで表示させるため）
#X列の名前が重なるため削除
st_write(d.utm_sf[,-1], "Utm.csv", layer_options = "GEOMETRY=AS_XY",
append=TRUE)
Updating layer `Utm' to data source `Utm.csv' using driver `CSV'
options: GEOMETRY=AS_XY
Writing 6371 features with 2 fields and geometry type Point.
読み込んで見る
d1 = read.csv("Utm.csv")
head(d1)
WGS84/UTMに変換したCSVが得られた。
4.まとめ
GPSデータから変換する場合のsf パッケージのポイント
○ st_as_sf( x , coords = c("long", "lat"), crs = 4326)
XY座標を含むデータフレームからsfオブジェクトに変換。
またspオブジェクトもsfへ変換できる。
○ st_transform(x , crs = 32736)
sfオブジェクトのcrsを変換する。
○ as(x, "Spatial")
sfオブジェクトをspオブジェクトに変換する。
○ st_write(x, "yyy.shp", append = TRUE)
sfオブジェクトをファイル名後の拡張子で自動的に変換されてセーブされる。
csvでセーブする場合は、layer_options = "GEOMETRY=AS_XY", をつけないと位置情報が１セル内に格納されてしまうので、その後のQ-GIS等で利用できない。