0

我有 5 只动物的 GPS 位置。我正在尝试计算 50% 和 95% 的内核分布，以确定以 km2 为单位的核心和外围区域的大小（作为表格），并且我还想绘制这些数据的等高线（50% 和 95%）。有人可以帮忙写代码吗？

加载所需的包后，我已经能够使用下面的代码生成光栅图像，但不能生成这些区域的轮廓或大小。

dat <- read.csv("data/GW1.csv", stringsAsFactors = FALSE)
dat <- dat[, 1:6]

dat1 <- dat %>% mutate(date = ymd_hm(Time)) %>% 
  make_track(Longitude, Latitude, date, id = ID, crs = CRS("+init=epsg:4326")) %>% 
  filter(y_ < 10) %>% 
  group_by(id) %>% nest() %>% 
  mutate(hr_kde_ud = map(data, hr_kde), 
         hr_kde_area = map_dbl(hr_kde_ud, hr_area), 
         hr_mcp = map_dbl(data, ~ hr_mcp(.) %>% hr_area %>% pull(area)))

# get the uds for first animal
raster::plot(dat1$hr_kde_ud[[1]]$ud)

我是一个非常基础的 R 用户，之前没有使用过 dput 函数。所以在这里，我提供了我的列标题和前 23 行：

Time                Latitude    Longitude
2019-06-12 15:00    4.8708      97.3669
2019-06-12 14:00    4.87185     97.37221667
2019-06-12 13:00    4.86825     97.37365
2019-06-12 12:00    4.874516667 97.3715
2019-06-12 11:00    4.875483333 97.369
2019-06-12 10:00    4.8749      97.3695
2019-06-12 9:00     4.873416667 97.3693
2019-06-12 8:01     4.868933333 97.37481667
2019-06-12 7:00     4.872683333 97.37523333
2019-06-12 6:00     4.8725      97.37618333
2019-06-12 5:00     4.872466667 97.37611667
2019-06-12 4:00     4.872316667 97.37696667
2019-06-12 3:00     4.875716667 97.37618333
2019-06-12 2:00     4.876083333 97.37525
2019-06-12 1:00     4.877633333 97.37145
2019-06-12 0:00     4.8786      97.37205
2019-06-11 23:00    4.879683333 97.37051667
2019-06-11 22:00    4.8795      97.37171667
2019-06-11 21:00    4.877583333 97.37091667
2019-06-11 20:00    4.878233333 97.36876667
2019-06-11 19:01    4.872666667 97.36978333
2019-06-11 18:00    4.87035     97.37046667
2019-06-11 17:00    4.86995     97.37108333

我希望得到一张表，上面写着：

50% 100 km2
95% 210 km2

以及显示轮廓的光栅图像（请注意，我的图像中没有显示轮廓） 在此处输入图像描述

4

1 回答 1

0

我使用 adehabitatHR 包对您的 MCP 数据进行了家庭范围分析（我添加了一个 id 变量作为第一列并排除了 Time 列）。您可能需要与您的数据集更相关的不同 CRS。KDE 应该遵循相同的逻辑。

library(adehabitatHR)

setwd("C:\\Users\\User\\desktop")

mydata <- read.csv("test.csv", header = T)
head(mydata)


mydata.spdf <-
  SpatialPointsDataFrame(
    coords = as.data.frame(cbind(mydata$Longitude,
                                 mydata$Latitude)),
    data = mydata,
    proj4string = CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
  )

#'  The data are in a degree-based coordinate system. transform them to a meter-based coordinate system, for purposes of area calculations and such with home ranges
mydata.spdf <- spTransform(
  mydata.spdf,
  CRS(
    "+proj=aea +lat_1=20 +lat_2=-23 +lat_0=0 +lon_0=25 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"
  )
)

# First: 95% MCP
mydata.mcp.95 <- mcp(mydata.spdf[, 1],
                     percent = 95,
                     unin = "m",
                     unout = "km2")
mydata.mcp.95

# Now: 50% MCP
mydata.mcp.50 <- mcp(mydata.spdf[, 1],
                     percent = 50,
                     unin = "m",
                     unout = "km2")
mydata.mcp.50
plot(mydata.mcp.95)
plot(mydata.mcp.50, add = T)

#' Kernel Density Estimate
mydata.Khref <- kernelUD(mydata.spdf[, 1], h = "LSCV")

#' measure the size of the home range using this kernel estimator
mydata.Khref.poly <- getverticeshr(mydata.Khref, percent = 95, unin = "m", unout = "km2")
print(mydata.Khref.poly)  # returns the area of each polygon
image(mydata.Khref)
plot(mydata.Khref.poly, add = T)
points(mydata.spdf, pch = 16)
于 2019-09-03T13:24:14.277 回答