我正在努力使用障碍来计算 Calenage 的adehabitatHR包中的 kernelUD 估计值。
我能够轻松地估计没有障碍物的家庭范围,但是由于鱼不能在陆地上游泳,我收到以下错误消息:
3 * h 中的错误:二元运算符的非数字参数
或者
.boundaryk(SpatialPoints(x, proj4string = CRS(as.character(pfs1))) 中的错误:不方便的边界:转动角度 > pi/2
取决于我是只使用湖的东侧(第一条错误消息),还是使用整个湖边界作为边界(第二条错误消息)。
下面是我正在使用的数据的片段。由于数据的庞大性,我只添加了一个简短的片段,但创建了一个指向所有其他文件的 github 链接。我的代码使用我的数据引导您完成直接在adehabitatHR示例中找到的步骤,您可以直观地看到湖泊屏障及其与估计的一条鱼的栖息地大小的关系。
我对多条鱼的家庭范围估计最终大于湖的实际面积。我想在包中解决这个问题,而不是仅仅将它们导出到 ArcGIS 并将它们剪切到湖边界。这不会那么准确。家庭范围估计的最大值应该是 3.25 平方公里,现在我可以得到几乎 3 倍大的估计值。本质上,鱼正在使用整个湖,但这些值扰乱了我的数据分析。
https://github.com/TRobin82/Home-Range.git
> head(burb)
coordinates TRANSMITTER DATETIME HPE TEMP DEPTH
1 (319278, 5223340) 366mm_M_a_t 2019-04-04 23:01:08 12.4 NA 10.6782
2 (319274, 5223340) 366mm_M_a_t 2019-04-04 23:17:39 13.0 4.1002 NA
3 (319275, 5223340) 366mm_M_a_t 2019-04-04 23:32:40 14.3 NA 10.6782
4 (319281, 5223340) 366mm_M_a_t 2019-04-05 00:03:33 10.3 NA 10.6782
5 (319282, 5223340) 366mm_M_a_t 2019-04-05 00:20:13 10.4 3.9433 NA
6 (319281, 5223340) 366mm_M_a_t 2019-04-05 00:35:25 10.3 NA 10.6782
Coordinate Reference System (CRS) arguments: +proj=utm +zone=15 +datum=WGS84 +units=m +no_defs
> HR_Burbs
TRANSMITTER SEX LENGTH WEIGHT BASIN HR_50 HR_75 HR_80 HR_85 HR_90 HR_95
X366mm_M_a_t 366mm_M_a_t M 366 556 1 0.088425 0.216225 0.255825 0.304650 0.370800 0.489825
X410mm_F_a_t 410mm_F_a_t F 410 710 1 1.086975 2.163825 2.502000 2.931300 3.525975 4.511700
X471mm_M_a_t 471mm_M_a_t M 471 980 2 0.141075 0.282600 0.325125 0.378225 0.447750 0.563850
X524mm_U_a_t 524mm_U_a_t U 524 1191 2 0.926100 2.023650 2.382300 2.841750 3.447675 4.352625
很抱歉没有包含太多示例,但代码和数据的性质不允许在堆栈溢出时进行有效和干净的帖子。
################# CODE FOR BARRIER JUST ON EAST SIDE OF LAKE ###################
################################################################################
library(rgdal)
library(maptools)
# Load KML coordinates
coords = getKMLcoordinates('E:\\BSUGradBurbot\\DATA\\BML Data\\BadMedEastBorderNEW.kml')
coords = SpatialPoints(coords, CRS('+proj=longlat'))
coords.proj = spTransform(coords, CRS=CRS('+init=epsg:32615 +units=m'))
a <- coords.proj@coords
b <- as.data.frame(a)
colnames(b) <- c("e","n","elev")
setwd("E:/BSUGradBurbot/DATA/BML Data")
#write.csv(b, "BML_east.csv")
Lake = read_csv("E:/BSUGradBurbot/DATA/BML Data/BML_east.csv")
str(Lake)
#convert read.table to list for bound layer
bound = structure(list(x = Lake$e, y = Lake$n), .Names = c("x","y"))
lines(bound, lwd=1)
str(bound)
## We convert bound to SpatialLines:
bound <- do.call("cbind",bound)
Slo1 <- Line(bound)
Sli1 <- Lines(list(Slo1), ID="frontier1")
barrier <- SpatialLines(list(Sli1), CRS("+init=epsg:32615"))
str(barrier)
## estimation of the UD with barrier
kud2 <- kernelUD(burb[,1], h="href", boundary = barrier)
> str(barrier)
Formal class 'SpatialLines' [package "sp"] with 3 slots
..@ lines :List of 1
.. ..$ :Formal class 'Lines' [package "sp"] with 2 slots
.. .. .. ..@ Lines:List of 1
.. .. .. .. ..$ :Formal class 'Line' [package "sp"] with 1 slot
.. .. .. .. .. .. ..@ coords: num [1:415, 1:2] 316441 316396 316397 316367 316377 ...
.. .. .. .. .. .. .. ..- attr(*, "dimnames")=List of 2
.. .. .. .. .. .. .. .. ..$ : NULL
.. .. .. .. .. .. .. .. ..$ : chr [1:2] "x" "y"
.. .. .. ..@ ID : chr "frontier1"
..@ bbox : num [1:2, 1:2] 316265 5219215 320108 5224210
.. ..- attr(*, "dimnames")=List of 2
.. .. ..$ : chr [1:2] "x" "y"
.. .. ..$ : chr [1:2] "min" "max"
..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
.. .. ..@ projargs: chr "+proj=utm +zone=15 +datum=WGS84 +units=m +no_defs"