使用 'move' 包计算突发利用率分布的体积面积
calculate volume area of a bursted utilization distribution using the 'move' package
我将 brownian.bridge.dyn() 函数应用于“MoveBurst”对象。然后我重新标准化了这个对象，以便突发段总和为 1 并避免舍入问题。我现在想使用 getVolumeUD() 函数计算 UD 的体积，但这似乎不起作用，因为重新标准化的 UD 属于 class ‘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, 1343634600, 1343640660, 1343644980, 1343655360, 1343656200,
1343661360, 1343665680, 1343671740, 1343707200, 1343709780, 1343711520,
1343714100, 1343715840, 1343717520, 1343718420, 1343721000, 1343721840,
1343722740, 1343722740, 1343723580, 1343725320, 1343727060), class = c("POSIXct",
"POSIXt"), tzone = "US/Eastern"), Tag = c(437L, 437L, 437L, 437L,
437L, 437L, 437L, 437L, 437L, 437L, 437L, 437L, 437L, 437L, 437L,
437L, 437L, 437L, 437L, 437L, 437L, 437L, 437L, 437L, 437L, 437L,
437L, 437L, 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.73803, 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 不存在，您的数据包含丢失的数据）所以我用一些不同的设置修复了它。然而，这并没有显示你的问题（移动的开发版本）：
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, 1343634600, 1343640660, 1343644980, 1343655360, 1343656200,
1343661360, 1343665680, 1343671740, 1343707200, 1343709780, 1343711520,
1343714100, 1343715840, 1343717520, 1343718420, 1343721000, 1343721840,
1343722740, 1343722740, 1343723580, 1343725320, 1343727060
), class = c(
"POSIXct",
"POSIXt"
), tzone = "US/Eastern"), Tag = c(
437L, 437L, 437L, 437L,
437L, 437L, 437L, 437L, 437L, 437L, 437L, 437L, 437L, 437L, 437L,
437L, 437L, 437L, 437L, 437L, 437L, 437L, 437L, 437L, 437L, 437L,
437L, 437L, 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.73803, 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]], ...): CRS object has comment, 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-07-13 创建
我将 brownian.bridge.dyn() 函数应用于“MoveBurst”对象。然后我重新标准化了这个对象，以便突发段总和为 1 并避免舍入问题。我现在想使用 getVolumeUD() 函数计算 UD 的体积，但这似乎不起作用，因为重新标准化的 UD 属于 class ‘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, 1343634600, 1343640660, 1343644980, 1343655360, 1343656200,
1343661360, 1343665680, 1343671740, 1343707200, 1343709780, 1343711520,
1343714100, 1343715840, 1343717520, 1343718420, 1343721000, 1343721840,
1343722740, 1343722740, 1343723580, 1343725320, 1343727060), class = c("POSIXct",
"POSIXt"), tzone = "US/Eastern"), Tag = c(437L, 437L, 437L, 437L,
437L, 437L, 437L, 437L, 437L, 437L, 437L, 437L, 437L, 437L, 437L,
437L, 437L, 437L, 437L, 437L, 437L, 437L, 437L, 437L, 437L, 437L,
437L, 437L, 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.73803, 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 不存在，您的数据包含丢失的数据）所以我用一些不同的设置修复了它。然而，这并没有显示你的问题（移动的开发版本）：
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, 1343634600, 1343640660, 1343644980, 1343655360, 1343656200,
1343661360, 1343665680, 1343671740, 1343707200, 1343709780, 1343711520,
1343714100, 1343715840, 1343717520, 1343718420, 1343721000, 1343721840,
1343722740, 1343722740, 1343723580, 1343725320, 1343727060
), class = c(
"POSIXct",
"POSIXt"
), tzone = "US/Eastern"), Tag = c(
437L, 437L, 437L, 437L,
437L, 437L, 437L, 437L, 437L, 437L, 437L, 437L, 437L, 437L, 437L,
437L, 437L, 437L, 437L, 437L, 437L, 437L, 437L, 437L, 437L, 437L,
437L, 437L, 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.73803, 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]], ...): CRS object has comment, 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-07-13 创建