最外郭法(Minimum Convex Polygon: MCP)で個体ごとの行動圏を推定し、
それぞれの個体ペアでどれくらい行動圏が重複するかを計算する。
カーネル推定の場合は adehabitatHR の kerneloverlaphr を使えばよいが、
MCP の場合は便利な関数がないような気がする。
というわけで試してみた。
とにかく推定&計算
いつものように UTM に変換した動物の位置情報を sp 形式に変換した後、
adehabitatHR の mcp、rgeos の gArea と gIntersection で簡単に計算できる。
library(adehabitatHR)
library(rgdal)
library(maptools)
library(rgeos)
# 1-2 列目に xy 座標、3 列目に ID を入れたファイルを読み込む
d<-read.csv(適切なファイル)
# sp 形式に変換
d_sp<-SpatialPoints(d[,c(1:2)], proj4string=CRS("+init=epsg:32652"))
d_sp<-data.frame(d_sp)
idd_sp<-data.frame(d[3])
coordinates(idd_sp)<-d_sp
# 早速 MCP
# デフォルトでは 95% MCP が推定される
hr<-mcp(idd_sp[,1],percent=100,unin="m",unout="ha")
# shp ファイルを出力
writePolyShape(hr,"homerange_mcp")
# 行動圏面積を計算
# byid = TRUE で ID ごとの行動圏面積を計算
hr_area<-gArea(hr,byid=TRUE)
# 重複面積を計算
# byid = TRUE で全 ID ペアの重複面積を計算
overlap<-gIntersection(hr,hr,byid=TRUE)
overlap_area<-gArea(overlap,byid=TRUE)
推定された行動圏
##データをひとつのデータフレームにまとめたい
これが意外とめんどうだった。
以下のようなデータフレームを得ることを目標とした。
ratio1 列は id1 の行動圏にしめる id1 と id2 の重複域の面積割合を示している。
# まずは overlap_area
# "id1 id2" というラベルがついているので、スペース前後で切り分けて split_name というリストを得る
split_name<-strsplit(names(overlap_area)," ")
# リスト内の各ベクトルの 1 番目と 2 番目の要素をそれぞれベクトルとしてまとめる
id1<-sapply(split_name,'[',1)
id2<-sapply(split_name,'[',2)
# データフレームでまとめる
d_overlap<-data.frame(overlap_area,id1,id2)
# 次に hr_area
# ラベルを取り出した後、データフレームでまとめる
d_hr<-data.frame(hr_area,names(hr_area))
# d_overlap と d_hr を id1 で merge したいので列の名前を変更
colnames(d_hr)[2]<-"id1"
# id1 で merge する
d_all<-merge(d_overlap,d_hr,by="id1")
# 行動圏面積にしめる重複面積の割合を計算
d_all$ratio1<-d_all$overlap_area/d_all$hr_area
(もっと丁寧な説明をする気力が出てきたら追記するかも)
(よりスマートにできる方法があったらおしえてください)
##参考にしたページ
adehabitatHR
https://cran.r-project.org/web/packages/adehabitatHR/adehabitatHR.pdf
Home Range Estimation in R: the adehabitatHR Package
https://cran.r-project.org/web/packages/adehabitatHR/vignettes/adehabitatHR.pdf
rgeos
https://cran.r-project.org/web/packages/rgeos/rgeos.pdf
apply系関数の使い方
http://takenaka-akio.org/doc/r_auto/list.html