はじめに

仕事で野生動物に装着したGPS首輪のデータ解析をして行動圏の推定をしているのですが、MCPやKDEなどとは異なる新しい解析がしたいと思ったので最近の研究結果を漁ったところ、Dynamic Brownian Bridge Movement Model (dBBMM) という方法があることを知ったので備忘録で残します。

手法そのものの概念は省略し、以下にプログラムを記載します。
とりあえず解析したいのは以下のようなデータ。

言語はRを使って解析をしていきます。

dBBMM解析プログラム

qiita.rb

#必要なパッケージ
library(tidyverse)
library(sf)
library(move)

moveパッケージはdBBMMを実行するために必要なパッケージです。
KDEの解析などではadehabitatHRなどが有名ですが、ここではmoveパッケージを使用して解析を進めます。

データ読み込み

qiita.rb

setwd("")  # 作業ディレクトリを変更する
df <- read.csv("testdata.csv") # データ読み込み
df$date <- as.POSIXct(df$date, format = "%Y/%m/%d %H:%M") #date列を日付型として指定

dBBMMではGPSの日付データを解析に使用するので、日付の列を日付型として指定をしてあげます。

解析前処理

qiita.rb

#データをmove関数を使ってmove_objectへ渡す
move_object <- move(x = df$longitude, y = df$latitude, time = df$date, proj = CRS("+proj=longlat +datum=WGS84"), data = df, animal = df$collar_id)
# 座標変換 (UTM Zoneに変換)
move_object_utm <- spTransform(move_object, CRS("+proj=utm +zone=53 +datum=WGS84"))
# 異常値を除去
move_object_utm <- move_object_utm[coordinates(move_object_utm)[, 1] > 0 & coordinates(move_object_utm)[, 2] > 0]
# 座標範囲確認
summary(coordinates(move_object_utm))

解析の前処理として、使用するデータはmove関数を使ってmove_objectへ渡してあげます。
元データがWGS84の緯度経度データになっているので、設定するCRSは "+proj=longlat +datum=WGS84" を指定します。
その後、解析に座標変換が必要になるので、データが該当する地域のUTM系に変換をしてあげます。
基本的には、変換のみでよいのですが、この時点でデータに外れ値がある場合、解析に失敗することがあるので異常値の除去や解析に使うデータの座標範囲を確認します。

dBBMM実行

qiita.rb

bbmm_result <- brownian.bridge.dyn(
  object = move_object_utm, 
  location.error = 25, 
  dimSize = 300,       #ラスタグリッドのセル数。目安100~500で変える。値が大きいほど細かいが処理が重くなる。
  margin = 3, 
  window.size = 13, 
  ext = 1.2 #c(min_x - buffer, max_x + buffer, min_y - buffer, max_y + buffer) #地理的範囲指定。なくてもよい。適当なら0.5。正確に書くなら右の例のような書き方→ ext = extent(-118.5, -117.5, 33.5, 34.5) 
)
#動かなかったらlocation.error, dimsize, ext のパラメータを変更すると動く。

各種パラメータを設定してdBBMMを実行します。
パラメータは解析対象の動物によって変わってくるので、試行錯誤が必要になりますが、先行研究や類似の研究の論文などから引用してもよいと思います。
特に、ここでミソになってくるのが、location.errorとdimsize、extのパラメータです。
このあたりのパラメータが適切でない場合、解析が失敗してうまく出力がされませんので、微調整を繰り返す必要がでてきます。

extの調整があまく、失敗すると以下のエラーを吐かれます。

qiita.rb

Computational size: 4.4e+11
Error in .local(object, raster, location.error = location.error, ext = ext,  : 
  Higher y grid not large enough, consider extending the raster in that direction or enlarging the ext argument

Computational sizeが大きすぎて解析が失敗することもあります。そういった場合は、解析範囲やデータの中に極端な外れ値がないか確認をします。

extを細かく指定する場合には、以下のような記述をして、解析範囲を設定してあげるとうまく動きます。

qiita.rb

# 緯度経度の最大・最小値を取得
min_x <- min(coordinates(move_object_utm)[, 1]) # 東西方向（Easting）
max_x <- max(coordinates(move_object_utm)[, 1]) # 東西方向（Easting）
min_y <- min(coordinates(move_object_utm)[, 2]) # 南北方向（Northing）
max_y <- max(coordinates(move_object_utm)[, 2]) # 南北方向（Northing）
# バッファ値（適宜調整、例として100mを指定）
buffer <- 100

bbmm_result <- brownian.bridge.dyn( ,
   ext = c(min_x - buffer, max_x + buffer, min_y - buffer, max_y + buffer)

解析結果のプロット

qiita.rb

bbmm_df <- as.data.frame(bbmm_result, xy = TRUE)

p <- ggplot(data = bbmm_df, aes(x = x, y = y)) +
  geom_raster(aes(fill = layer)) +
  scale_fill_viridis_c() +
  labs(title = "BBMM Home Range", x = "Easting", y = "Northing", fill = "Utilization") +
  theme_minimal()

print(p)

解析したdBBMMの結果をプロットして確認します。
testdataの解析結果をプロットを実行すると以下の解析結果が出力されます。

上手く解析できているようです。

データエクスポート

qiita.rb

# dBBMM解析結果をrasterオブジェクトに変換
bbmm_raster <- rasterFromXYZ(bbmm_df[, c("x", "y", "layer")])

# 投影情報を設定
crs(bbmm_raster) <- CRS("+proj=utm +zone=53 +datum=WGS84")

# GeoTIFF形式でエクスポート
writeRaster(bbmm_raster, "bbmm_result.tif", format = "GTiff", overwrite = TRUE)

解析したdBBMMの結果をGISで使えるように、GeoTIFF形式でエクスポートします。
また、以下のようにすればベクタ形式で出力することもできます。
95%や50%などが必要な場合は以下のように出力するとよいでしょう。
ただし、変換や出力はデータが重いほど時間がかかります。

qiita.rb

###ベクタに変換###
###95%と50%に分けて出力###
# 利用分布 (Utilization Distribution, UD) の計算
bbmm_ud <- getVolumeUD(bbmm_result)

# 95%および50%行動圏の計算の描画
contour_95 <- raster::contour(bbmm_ud, levels=0.95, labels="", drawlabels=FALSE, lwd=2, col="blue")
contour_50 <- raster::contour(bbmm_ud, levels=0.50, labels="", drawlabels=FALSE, lwd=2, col="red")

# 投影情報を設定
crs(bbmm_ud) <- CRS("+proj=utm +zone=53 +datum=WGS84")

# 95%および50%行動圏のポリゴンを抽出
ud_95 <- rasterToPolygons(bbmm_ud, fun=function(x) x <= 0.95, dissolve=TRUE)
ud_50 <- rasterToPolygons(bbmm_ud, fun=function(x) x <= 0.50, dissolve=TRUE)

# sfオブジェクトに変換
ud_95_sf <- st_as_sf(ud_95)
ud_50_sf <- st_as_sf(ud_50)

# 投影情報を設定
st_crs(ud_95_sf) <- CRS("+proj=utm +zone=53 +datum=WGS84")
st_crs(ud_50_sf) <- CRS("+proj=utm +zone=53 +datum=WGS84")

# ジオパッケージとしてエクスポート
write_sf(obj = ud_95_sf, dsn = "bbmm_polygons_95.gpkg")
write_sf(obj = ud_50_sf, dsn = "bbmm_polygons_50.gpkg")

###95%50%結合版###
# 利用分布のポリゴン化（dimsizeが大きいと変換に時間がかかる）
bbmm_polygons <- rasterToPolygons(bbmm_raster, dissolve = TRUE)

#sf形式に変換
sf_bbmm_polygons <- st_as_sf(bbmm_polygons)

# ジオパッケージとしてエクスポート
write_sf(obj = sf_bbmm_polygons, dsn = "bbmm_polygons.gpkg")

結果の確認

最後に、解析したdBBMMの行動圏がGIS上で見られるか確認をします。
QGISでデータを読み込んであげ、スタイルを設定してあげると以下のようになりました。

おわりに

今回概念的なものは省略しましたが、ひとまず解析そのものは上記のプログラムで実行は可能なようです。
今後は野生動物それぞれでどのようなパラメータが必要か検討したいと思います。

日本では最近はなかなか野生動物の行動圏解析は流行っていませんが、海外ではかなり盛んなようです。
備忘録的に残すこの記事が行動圏解析をしたい学生さんなどに役立つと幸いです