野生動物の行動圏を推定するパッケージadehabitatHRとGISとの連携

# 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の概要
```r
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を指定する。
```r
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()でプロットが簡単にできる。
```r
ggplot(d_sf) + geom_sf()+theme_bw()
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/161102/e2bf2fd5-41b2-30ab-a747-72afafc1a007.png)

CRSを平面直角座標系であるWGS84/UTMに変換する。ここでは南アフリカの36Sゾーンを指定する。
```r
d.utm_sf = st_transform(d_sf, crs = 32736)
```
UTM変換もGISソフト上で実施すると、意外と手間がかかるものだがRでは一発である。
　そのほか、ここでは割愛するが、このsfパッケージはQ-GISなどで行っていたバッファーゾーンの発生等いろいろなポリゴン処理が簡単にできるので、非常に便利。

# 3．adehabitatHRを使った行動圏の解析
では、本題のadehabitatHRを使ってカーネル密度推定法を用いた行動圏推定を行う。
```math
\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オブジェクトから変換する必要がある。
```r
d.utm_sp = as(d.utm_sf, "Spatial")
```
カーネル密度推定（KDE)を使う
```r
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%行動圏を推定する。
```r
HR50 = getverticeshr(UD,50)
HR95 = getverticeshr(UD,95)
#再びsfへ変換
HR50sf = st_as_sf(HR50)
HR95sf = st_as_sf(HR95)
```
プロットしてみる
```r
ggplot(HR95sf) +geom_sf() + geom_sf(data=HR50sf) +theme_bw()
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/161102/44c2a75e-570c-6786-9612-f1f122ab17bc.png)

動物のコア行動圏と95%行動圏の推定が完了。
Q-GISで表示可能なESRI shapeファイルに変換し、保存する。append=TRUEにしておく。
```r
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で表示させるため）
```r
#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.
```
読み込んで見る
```r
d1 = read.csv("Utm.csv")
head(d1)
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/161102/22d5cbfa-b592-804a-144f-392d7e0c16b7.png)
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等で利用できない。


# 5.Enjoy