写给于欣的笔记

于欣 动物迁徙

---

球面几何学若干

假如地球上一点$A$的经度是$\alpha_A$，纬度是$\delta_A$，另一点$B(\alpha_B,\delta_B)$
两点的经度差是$\Delta\alpha$，纬度差是$\Delta\delta$
两点在地球上的角距离（以弧度为单位）

$$\Delta\theta=\sqrt{(\Delta\alpha\cos\delta_A)^2+(\Delta\delta)^2}$$
如果要求距离长度，乘以地球半径$R$即可，即是$R\Delta\theta$。
注意：此公式仅在两点经纬度分别相差不大时（$\Delta\alpha$，$\Delta\delta$较小）使用，当经纬度有一个相差较大时，需要使用更严格的公式，根据球面三角形边的余弦定理，可得如下公式：
$$\cos\Delta\theta=\sin\delta_A\sin\delta_B+\cos\delta_A\cos\delta_B\cos(\alpha_A-\alpha_B)$$
然后利用反三角函数解出$\Delta\theta=\arccos(\cos\Delta\theta)$即可。

经纬度的“面积”与球面上真实的面积换算

假设$A$点附近有一个区域的经纬度面积是$dS_{l}=d\alpha d\delta$，其在球面上的真实面积是$dS_{r}=R\cos\delta_{A}d\alpha \cdot Rd\delta$，它们之间的换算关系如下：

$$dS_{r} = dS_l\cdot R^2\cos\delta_{A}.$$
如果区域的面积很大，真实面积只能由积分得到
$$S_{r}=\int_{S} dS_{r}=\int_{S}dS_{l}R^2\cos\delta_{A}$$
关于纬度的说明：
北纬：$\delta>0$
南纬：$\delta<0$
注意以上计算是粗估距离，实际地球半径在不同处不一样，大地基准面不是球形，严格的计算建议用Mathematica软件。

1.23445e-06=$1.23445\times10^{-6}$

地图投影

关于地图投影，小距离的范围内用墨卡托投影没有问题。
球形墨卡托投影的公式
经纬度换算为横纵坐标
$(\alpha,\delta)\rightarrow(x,y)$

$$x=R\alpha$$
$$y=R\ln\left(\tan\left(\frac{\pi}{4}+\frac{\delta}{2}\right)\right)$$
严格地说来，要保证面积不变，最好用高斯投影

adehabitatHR软件包的说明

首先要计算Utilization Distribution（UD），再根据UD估算Home Range
BRD是基于Biased Random Bridges对进行UD估计
kernelbb is the estimation of Kernel Brownian Bridge Home-Range
kernelUD是基于kernel estimation方法对UD进行估计

UD

A utilization distribution is a probability distribution giving the probability density that an animal is found at a given point in space. It is estimated from data sampling the location of an individual or individuals in space over a period of time using, for example, telemetry or GPS based methods.
来自Wiki

Kernel estimation方法

kud <- kernelUD(puechabonsp$relocs[,1], h="href")
UD

$$\hat{f}(\boldsymbol{x})=\frac{1}{nh}\sum_{i=1}^{n}K\left(\frac{\boldsymbol{x}-\boldsymbol{X}_i}{h}\right)$$
$h$是光滑参数，reference bandwidth $h=\sigma n^{-1/6}$， $\sigma=\frac{\sigma_x+\sigma_y}{2}$
用likelihood cross-validation来找
grid控制分辨率，grid越高分辨率越高，extent控制显示的范围

Brownian bridge kernel方法

tata <- kernelbb(x, sig1 = 6.23, sig2 = 58, grid = 50) #UD
image(tata)


Biased random bridge kernel方法(MBKE)

原理

$$f(\boldsymbol{r},t_i|\boldsymbol{r}_1,\boldsymbol{r}_2)=\frac{t_2-t_1}{4\pi D(t_2-t_i)}\exp\left[\frac{\boldsymbol{r}_m\boldsymbol{D}\boldsymbol{r}_m}{4p_i(t_2-t_i)}\right]$$
$\boldsymbol{r}_m=(x_1+p_i(x_2-x_1),y_1+p_i(y_2-y_1))$, $p_i=(t_i-t_1)/(t_2-t_1)$
协方差矩阵$\boldsymbol{D}$可以简化为扩散系数$D$
光滑参数
$$h_i^2=h_{\min}^2+4p_i(1-p_i)(h_{\max}^2-h_{\min}^2)T_i/T_{m}ax$$
这里$h_{\max}^2=h_{\min}^2+D\times T_{\max}/2$或者$h_{\max}^2=D\times T_{\max}/2$
$h_{\min}$至少要达到定位的误差
Code
vv <- BRB.D(tr, Tmax=180*60, Lmin=50, habitat=buffalo$habitat, activity="act")
ud <- BRB(tr, D = vv, Tmax = 180*60, tau = 300, Lmin = 50, hmin=100,
habitat = buffalo$habitat, activity = "act", grid = 100, b=0,
same4all=FALSE, extent=0.01)

比较
kernel estimation仅考虑了采样点附近的UD, Brownian bridge加上了时间，考虑任意两采样点之间的每一点的UD。

Home range

找home range有多种方法

vector mode
homerange <- getverticeshr(UD)
plot(homerange)
getverticeshr用来找home range的范围
plot(getverticeshr(UD, 95), add=TRUE, lwd=2)
画出home range
raster mode
vud <- getvolumeUD(UD)
image(vud)
contour(vud, add=TRUE, nlevels=1, levels=95, lwd=3, drawlabels=FALSE)

kernel.area用于求home range的面积
ii <- kernel.area(UD, percent=seq(50, 95, by=5))

Overlap

相关系数$r_s$
VI

$$\iint\min[{\rm UD_1}(x,y),{\rm UD_2}(x,y)]dxdy$$
BA
$$\iint\sqrt{{\rm UD_1}(x,y){\rm UD_2}(x,y)}dxdy$$
UDOI
$$\iint{\rm UD_1}(x,y){\rm UD_2}(x,y)dxdy$$
HD

数值积分一览

$$\iint f(x,y)dxdy\approx\sum_{i,j}f(x_{ij},y_{ij})\Delta x\Delta y$$
假如$x_{ij},y_{ij}$是$N\times N$的矩阵，$\Delta x=\frac{x_{\max}-x_{\min}}{N}$, $\Delta y=\frac{y_{\max}-y_{\min}}{N}$

潮汐水位与时间关系的计算

可参考文章《潮汐调和分析与应用研究》中的模型，以及书《潮汐原理与计算》
有现成的python程序像UTide, ttide, pytides，可参见 http://pyoceans.github.io/sea-py/

R笔记

R语言报错的一些处理
变量不是矩阵类型的，可以转化为矩阵
就保留它的所有数据，
[[1]]代表单个个体

报错记录

出现报错
non unique datas for a given burst
的可能原因：
在排除数据表格里的数据本身重复之外，出现该报错的常见原因
1. 日期（或时间）输入格式不对, 年Y写成y, -写成/之类的，时间的名称不对time或t等等。
2. 没有告知时间，只有日期，同时告知时间信息。
3. 读入数据时，读出大量空白行，检查行数是否与实际数据一致，否则删去空白行。
4. 多个个体时，没有告知ID。

蛋坨子
经度122°4'41''=122.078°
纬度39°10'47''=39.1797°