我正试图计算整个物种的活动范围,并根据性别将两组分开。我使用mcp(),它运行,但输出是有问题的。
这是我的数据
library(raster)
library(dismo)
library(rgdal)
library(dplyr)
library(tidyverse)
library(tidyselect)
library(sp)
library(adehabitatHR)
library(scales)
cat.shark.data <- read.csv2("C:/Users/jcdeb/OneDrive/Bureaublad/SASC/hopefully last file ever.csv", header=T)
head(cat.shark.data)
Common.name Sex..M.F. Long Lat
1 Pyjama Catshark F 19.34785 -34.75671
2 Pyjama Catshark F 19.29512 -34.61850
3 Pyjama Catshark M 19.29512 -34.61850
4 Pyjama Catshark M 19.29512 -34.61850
5 Pyjama Catshark M 19.29512 -34.61850
6 Pyjama Catshark F 19.34581 -34.55419
# converting Lat long spatialpoints to UTM zone 34 spatialPoints
cord.dec = SpatialPoints(cbind(cat.shark.data$Long, cat.shark.data$Lat), proj4string=CRS("+proj=longlat"))
cord.dec
cord.dec@coords
zone <- 34
# used example from
# https://stackoverflow.com/questions/38621660/unexpected-convertion-output-from-latlong-to-utm-in-r
cord.UTM <- spTransform(cord.dec, CRS(paste("+proj=utm +south +zone=",zone,"ellps=WGS84",sep='')))
cord.UTM
par(mfrow = c(1, 2))
plot(cord.dec, axes = TRUE, main = "Lat-Long Coordinates", cex.axis = 0.95)
plot(cord.UTM, axes = TRUE, main = "UTM Coordinates", col = "red", cex.axis = 0.95)
# replace column with Lat and Long for UTM coordinates
cord.dec.coords <- as.data.frame(cord.UTM@coords)
cat.shark.data$Lat <-cord.dec.coords$coords.x1
cat.shark.data$Long <- cord.dec.coords$coords.x2
plot(cat.shark.data$Long, cat.shark.data$Lat)
这部分应该是好的,因为我检查了longlat到UTM的转换,坐标确实与正确的UTM值相对应。因此,我认为问题出在下一部分:
par(mfrow = c(1, 1))
x<- cat.shark.data[,"Long"]
y<- cat.shark.data[,"Lat"]
plot(x,y)
shark.sp<- cat.shark.data[,c("Sex..M.F.", "Long", "Lat")]
coordinates(shark.sp)<- c("Long", "Lat")
class(shark.sp)
slot(shark.sp, "proj4string") <- CRS( "+proj=utm +zone=34 +datum=WGS84 +units=m +no_defs" )
#proj4string(shark.sp) <- CRS( "+proj=utm +zone=34H +datum=WGS84 +units=m +no_defs" ) # should be the same right ?
shark.MCP<- mcp(shark.sp, percent = 95, unout = c( "km2"))
shark.MCP
Object of class "SpatialPolygonsDataFrame" (package sp):
Number of SpatialPolygons: 2
Variables measured:
id area
F F 49.26988
M M 49.26988
plot(shark.sp, col = as.factor(shark.sp@data$Sex..M.F.), pch = 16)
plot(shark.MCP, col = alpha(1:5, 0.5), add = TRUE)
hrs <- mcp.area(shark.sp, percent = seq(50, 100, by = 5))
hrs
F M
50 50.91919 50.91919
55 50.91919 50.91919
60 50.91919 50.91919
65 50.91919 50.91919
70 50.91919 50.91919
75 233.67845 57.53127
80 299.54642 87.06809
85 301.38459 127.67519
90 633.39131 606.42969
95 4926.98764 4926.98764
100 34146.77787 20543.01517
当我绘制数据点时,它们看起来有实际的间距,对于物种来说,如果95%的范围在50平方公里左右,这不是不可能的。然而,女性和男性不太可能有完全相同的范围,当我看hrs的结果时,它们根本不一致。不幸的是,这是我第一次尝试,所以我不知道问题出在哪里。
在mcp()代码行之后我也得到了这个警告:
In proj4string(xy) :
CRS object has comment, which is lost in output; in tests, see
https://cran.r-project.org/web/packages/sp/vignettes/CRS_warnings.html
和11个警告,当我运行hrs()代码。然而,我把这个网站作为参考,并尝试了他们的脚本与他们的数据:https://jamesepaterson.github.io/jamespatersonblog/03_trackingworkshop_homeranges并且得到了相同的警告,同时仍然得到与示例相同的结果。
我查看了警告,甚至调整了一些代码来解决它,但我仍然得到警告和相同的结果。
编辑:在玩了更多的代码,我只得到警告,如果我运行这个代码行
slot(shark.sp, "proj4string") <- CRS( "+proj=utm +zone=34 +datum=WGS84 +units=m +no_defs" )
有没有人知道我哪里错了,因为我已经尝试了很多,似乎没有工作。
提前感谢您!
在运行mcp之前,我会尝试将鲨鱼点分成雄性和雌性,看看这是否会给出更真实的结果。