使用“移动”包计算突发利用率分布的体积面积

如何解决使用“移动”包计算突发利用率分布的体积面积

我将 brownian.bridge.dyn() 函数应用于“MoveBurst”对象。然后我重新标准化了这个对象，以便突发段总和为 1 并避免舍入问题。我现在想使用 getVolumeUD() 函数计算 UD 的体积，但这似乎不起作用，因为重新标准化的 UD 属于“list”类。

代码如下：

library(move)
    
move.ind <- move(x=ind$Lon,y=ind$Lat,time=ind$Datetime,proj=CRS("+proj=longlat +datum=WGS84"),data=ind,animal=ind$Tag,sensor="VR2W")

move.ind$LocationError = 1 

r.ind <- spTransform(move.ind,center=TRUE)

TimeDiff <- timeLag(r.ind,units="hours")
r.ind$TimeDiff <- append(0,TimeDiff)

e <- 5000
xUTM.ind <- raster(xmn=min(ind$NewEastingUTM)-e,xmx=max(ind$NewEastingUTM)+e,ymn=min(ind$NewNorthingUTM)-e,ymx=max(ind$NewNorthingUTM)+e,crs =CRS("+proj=utm +zone=17 +datum=WGS84"),resolution=50)

x.ind <- raster(xmn=min(DET$NewEastingUTM),xmx=max(DET$NewEastingUTM),ymn=min(DET$NewNorthingUTM),ymx=max(DET$NewNorthingUTM),crs=CRS("+proj=utm +zone=17 +datum=WGS84"))
proj4string(x.ind)

newTemplate <- projectExtent(xUTM.ind,proj4string(x.ind))

ones = rep(1,(newTemplate@ncols*newTemplate@nrows))
xAEQD.ind <- setValues(newTemplate,ones)

rNew.ind <- spTransform(r.ind,proj4string(xAEQD.ind))

bursted <- burst(rNew.ind,c('normal','long')[1+(timeLag(rNew.ind,units = 'hours')>4)]) 

bursted_dbbmm <- brownian.bridge.dyn(bursted,burstType='normal',raster = xAEQD.ind,location.error = "LocationError",time.step = 5,ext = 3,window.size = 29)

tmp <- calc(bursted_dbbmm,sum)
dbbmmlist <- new(".UD",tmp/sum(values(tmp)))

cont.ind <- getVolumeUD(dbbmmlist)

伪代码如下所示：

ind <- structure(list(Datetime = structure(c(1343535240,1343539560,1343543880,1343548200,1343622480,1343625120,1343628540,1343630280,1343634600,1343640660,1343644980,1343655360,1343656200,1343661360,1343665680,1343671740,1343707200,1343709780,1343711520,1343714100,1343715840,1343717520,1343718420,1343721000,1343721840,1343722740,1343723580,1343725320,1343727060),class = c("POSIXct","POSIXt"),tzone = "US/Eastern"),Tag = c(437L,437L,437L),Lat = c(25.73785,25.73779,25.73822,25.73803,25.7387,25.73772,25.73842,25.73759,25.73693,25.73765,25.73781,25.73819,25.73975,25.73766,25.73758,25.73799,25.74314,25.73943,25.74317,25.7407,25.74324,25.73739,25.74175,25.71287,25.73699,25.73749,25.74292,25.73809,25.743,25.74249
),Lon = c(-79.24489,-79.24543,-79.24039,-79.24479,-79.24606,-79.24551,-79.23763,-79.24604,-79.24535,-79.24622,-79.24678,-79.24583,-79.24725,-79.2454,-79.24459,-79.24632,-79.24553,-79.24629,-79.24731,-79.24647,-79.24724,-79.24652,-79.24577,-79.24749,-79.24636,-79.24537,-79.24547,-79.24742,-79.24582,-79.2464,-79.24771),NewEastingUTM = c(676052.593850702,675998.505006306,676503.524589617,676062.361342852,675933.957732871,675990.581885337,676780.134256345,675937.599356069,676007.79854611,675919.451865544,675863.032391998,675957.784408548,675813.022339476,676001.706417219,676082.426966281,675909.522121439,675988.177486275,675904.342388974,675807.473799857,675886.239923252,675812.627071404,675881.120633141,675964.982682709,675786.000042484,675941.884041032,676005.703553397,675994.93392188,675791.30042537,675958.935010342,675893.512994838,675762.839399625),NewNorthingUTM = c(2847824.10140759,2847816.73427449,2847871.10249762,2847844.17392274,2847916.6963574,2847808.8734514,2847896.95405131,2847793.76586291,2847721.57691275,2847800.1720647,2847817.14869788,2847860.50938921,2848031.42006168,2847802.37392041,2847844.44095791,2847792.28461383,2847838.75526749,2848408.21836108,2847995.89296004,2848411.30139552,2848136.66700018,2848418.98875582,2847771.97165951,2848252.64452269,2845055.05409167,2847728.1965547,2847783.44921703,2848382.34136664,2847849.44550836,2848392.56348938,2848334.32268777)),row.names = c(NA,31L),class = "data.frame")

解决方法

我已尝试重现您的问题。但是，您的代码似乎无法直接运行（DET 不存在，您的数据包含缺失数据），因此我使用一些不同的设置对其进行了修复。但是，这并没有显示您的问题（使用 move 的开发版本）：

library(move)
#> Loading required package: geosphere
#> Loading required package: sp
#> Loading required package: raster
#> Loading required package: rgdal
#> rgdal: version: 1.5-23,(SVN revision 1121)
#> Geospatial Data Abstraction Library extensions to R successfully loaded
#> Loaded GDAL runtime: GDAL 3.0.4,released 2020/01/28
#> Path to GDAL shared files: /usr/share/gdal
#> GDAL binary built with GEOS: TRUE 
#> Loaded PROJ runtime: Rel. 6.3.1,February 10th,2020,[PJ_VERSION: 631]
#> Path to PROJ shared files: /usr/share/proj
#> Linking to sp version:1.4-5
#> To mute warnings of possible GDAL/OSR exportToProj4() degradation,#> use options("rgdal_show_exportToProj4_warnings"="none") before loading rgdal.


ind <- structure(list(Datetime = structure(c(
  1343535240,1343539560,1343543880,1343548200,1343622480,1343625120,1343628540,1343630280,1343634600,1343640660,1343644980,1343655360,1343656200,1343661360,1343665680,1343671740,1343707200,1343709780,1343711520,1343714100,1343715840,1343717520,1343718420,1343721000,1343721840,1343722740,1343723580,1343725320,1343727060
),class = c(
  "POSIXct","POSIXt"
),tzone = "US/Eastern"),Tag = c(
  437L,437L,437L
),Lat = c(
  25.73785,25.73779,25.73822,25.73803,25.7387,25.73772,25.73842,25.73759,25.73693,25.73765,25.73781,25.73819,25.73975,25.73766,25.73758,25.73799,25.74314,25.73943,25.74317,25.7407,25.74324,25.73739,25.74175,25.71287,25.73699,25.73749,25.74292,25.73809,25.743,25.74249
),Lon = c(
  -79.24489,-79.24543,-79.24039,-79.24479,-79.24606,-79.24551,-79.23763,-79.24604,-79.24535,-79.24622,-79.24678,-79.24583,-79.24725,-79.2454,-79.24459,-79.24632,-79.24553,-79.24629,-79.24731,-79.24647,-79.24724,-79.24652,-79.24577,-79.24749,-79.24636,-79.24537,-79.24547,-79.24742,-79.24582,-79.2464,-79.24771
),NewEastingUTM = c(
  676052.593850702,675998.505006306,676503.524589617,676062.361342852,675933.957732871,675990.581885337,676780.134256345,675937.599356069,676007.79854611,675919.451865544,675863.032391998,675957.784408548,675813.022339476,676001.706417219,676082.426966281,675909.522121439,675988.177486275,675904.342388974,675807.473799857,675886.239923252,675812.627071404,675881.120633141,675964.982682709,675786.000042484,675941.884041032,676005.703553397,675994.93392188,675791.30042537,675958.935010342,675893.512994838,675762.839399625
),NewNorthingUTM = c(
  2847824.10140759,2847816.73427449,2847871.10249762,2847844.17392274,2847916.6963574,2847808.8734514,2847896.95405131,2847793.76586291,2847721.57691275,2847800.1720647,2847817.14869788,2847860.50938921,2848031.42006168,2847802.37392041,2847844.44095791,2847792.28461383,2847838.75526749,2848408.21836108,2847995.89296004,2848411.30139552,2848136.66700018,2848418.98875582,2847771.97165951,2848252.64452269,2845055.05409167,2847728.1965547,2847783.44921703,2848382.34136664,2847849.44550836,2848392.56348938,2848334.32268777
)),row.names = c(NA,31L),class = "data.frame")


ind<-ind[-28,][-10,]


move.ind <- move(
  x = ind$Lon,y = ind$Lat,time = ind$Datetime,proj = CRS("+proj=longlat +datum=WGS84"),data = ind,animal = ind$Tag,sensor = "VR2W"
)

move.ind$LocationError <- 1

r.ind <- spTransform(move.ind,center = TRUE)

rNew.ind <- spTransform(r.ind,CRS("+proj=utm +zone=17 +datum=WGS84"))

bursted <- burst(rNew.ind,c("normal","long")[1 + (timeLag(rNew.ind,units = "hours") > 4)])

bursted_dbbmm <- brownian.bridge.dyn(bursted,raster = 50,burstType = "normal",# raster = xAEQD.ind,location.error = "LocationError",time.step = 5,ext = 3,window.size = 15,margin = 7
)
#> Warning in proj4string(object): CRS object has comment,which is lost in output
#> Warning in .local(object,raster,location.error = location.error,ext = ext,:
#> Some burst are omitted for variance calculation since there are not segements of
#> interest
#> Warning in FUN(X[[i]],...): CRS object has comment,which is lost in output
#> Computational size: 9.3e+06
#> Warning in FUN(X[[i]],which is lost in output
#> Computational size: 2.3e+06

tmp <- calc(bursted_dbbmm,sum)
dbbmmlist <- new(".UD",tmp / sum(values(tmp)))

cont.ind <- getVolumeUD(dbbmmlist)

^{由 reprex package (v2.0.0) 于 2021 年 7 月 13 日创建}