AdehabitatHR 的 MCP 家庭范围太小

Posted 2023-03-17

【中文标题】AdehabitatHR 的 MCP 家庭范围太小

【英文标题】：MCP Home range with AdehabitatHR too small

【发布时间】：2021-11-09 19:50:47

【问题描述】：

我是使用 R 进行空间分析的新手。我正在尝试主要使用 R 中的 adehabitatHR 包计算每日家庭范围（MCP 100、后来的 95 和 50%）。我使用了我的同事提供的脚本和它一直为她工作。现在，当我运行它时，它可以工作，但是我收到的 e-xx 值太小了。 我真的不知道哪里出了问题，我认为它与投影有关，但我太缺乏经验，不知道哪个应该更好。

您可以看到下面的代码。也许有人知道如何解决这个问题？

library (adehabitatHR)
library (sp)
library (scales)  # helps make Polygons partly transparent using the alpha argument below 
library(dplyr)

# Read data frame into R 
tab1<-read.csv(file="animals.csv",sep=",",header=TRUE)

# Remove two rows with NA´s 
tab1 <- tab1 [!is.na(tab1$x_) & !is.na(tab1$y_),]

# Only include three columns (AnimalID + Exp, x and y coordinates) for making MCP´s 
# Creating new ID to group days and AnimalID
tab2 <- tab1 [, c("id", "x_", "y_", "day", "month", "year")] # week, month, day + month
tab2$Time <- paste(tab2$day, tab2$month, sep= " ") #  day + month
tab2$Time2 <- paste(tab2$Time, tab2$year)
tab2$Mix <- paste(tab2$AnimalID, tab2$Time2, sep= " ") # id + time 
tab3 <- tab2 [, c("Mix", "x_", "y_")]
tibble(tab3)

# Check in tab2 which Animal has less than 5 GPS location
tab4 <- tab3 %>%
  group_by(Mix) %>%
  summarise(n = n())

# Filter data with less than 5 
tab5<-filter(tab4, n > 5)
tab6<-filter(tab4, n < 5)
  
# Remove rows from data set which have less than 5 GPS location
tab7 <- tab3 %>% 
  filter(!grepl(" 31 May 2019", Mix))
tab8 <- tab7 %>% 
  filter(!grepl(" 29 Jun 2019", Mix))

# Create a SpatialPointsDataFrame by defining the coordinates
coordinates(tab8) <- c("x_", "y_")
str(tab8)

# Set the coordinate reference system, (CRS) 
# The sample data are UTM points in WGS84 from zone 33N 
proj4string(tab8) <- CRS( "+proj=utm +zone=33N +datum=WGS84 +units=m +no_defs" )

animals.mcp <- mcp(tab8[, "Mix"], percent = 100, unin = "m", unout = "m2") 
as.data.frame(animals.mcp)

输出数据如下所示：

                      id         area
 1 Dec 2019    1 Dec 2019 5.180087e-03
 1 Feb 2020    1 Feb 2020 7.171711e-04
 1 Jan 2020    1 Jan 2020 1.692209e-03
 1 Jul 2019    1 Jul 2019 3.384402e-07
 1 Jun 2019    1 Jun 2019 1.829374e-01
 1 Nov 2019    1 Nov 2019 7.794313e-04
...

提前致谢！

【问题讨论】：

嗨，贾斯汀——欢迎您！如果您可以提供一些示例数据，如本常见问题解答中所述，关于如何提供获得最佳答案所需的所有信息，将会有所帮助：***.com/questions/5963269/…。没有它，就很难确定。

看来你是正确的，投影是可能的罪魁祸首。幸运的是，您可以将面积单位从测量单位转换为所需单位（例如，通过除以 1e6 从 m^2 到 km^2）。

谢谢！我能够解决它，不确定它是否是一种优雅的方式，但它现在可以工作：coordinates(tab7) <- c("x_", "y_") proj4string(tab7) <- CRS("+init=epsg:4326") tab8 <- spTransform(tab7, CRS("+proj=utm +zone=33N +datum=WGS84 +units=m +no_defs")) 所以基本上不得不告诉他这是 WGS84，然后将它投影到我的区域。在 GIS 中导出，现在看起来不错！

太好了！我很高兴=）你应该把上面写的作为答案（这样做是完全好的形式！）以便其他人可以找到他们需要的答案。

【参考方案1】：

所以我可以通过添加来解决它

coordinates(tab7) <- c("x_", "y_") 
proj4string(tab7) <- CRS("+init=epsg:4326") 
tab8 <- spTransform(tab7, CRS("+proj=utm +zone=33N +datum=WGS84 +units=m +no_defs"))

所以首先告诉 R 使用什么坐标系，然后将其投影到需要的区域。

希望这可以帮助其他有同样问题的人！