当我尝试按 ID 计算区域重叠时,group_by 函数不起作用
The group_by function is not working when I try to calculate area overlap by ID
我有 GPS 数据,可用于确定家庭范围。我想测量范围与受保护区域形状文件 (available here) 的重叠。当我有一个人时这很好用,但当我尝试使用多个 ID 时却不行。
问题似乎出在这段代码中的 group_by
函数上:
trk %>% group_by(id) %>%
hr_kde(., levels = c(.95)) %>%
hr_isopleths(.) %>%
st_intersection(., merged_Africa_tranform) %>%
st_area(.) / 1e6
它只会产生一个值 (29.77905),而我期望两个,每个 ID 一个(A 为 34.68964,B 为 38.99062)。
完整代码如下:
# load packages
library(tidyverse)
library(amt)
library(sf)
library(adehabitatHR)
#' load in the protected areas (working with Albers Equal Area)
merged_Africa = read_sf("shape//eswatini//WDPA_Apr2020_SWZ-shapefile-polygons.shp")
st_crs(merged_Africa) <- 4326
merged_Africa_tranform <- st_transform(merged_Africa, "ESRI:102022")
#' try a data frame with multiple IDs
x_ <- c(707692, 707589, 707998, 708407, 708916, 709415, 707743, 707429, 707971, 708143, 708981, 709156)
y_ <- c(-3030991,-3031423,-3031640,-3031750,-3032508,-3037158,-3030995,-3031723,-3031680,-3031755,-3032408,-3037758)
id <- c(rep("A",6), rep("B", 6))
mydata <- data.frame(x_,y_,id)
# transform to trk object
trk <-
mk_track(mydata,
.x = x_,
.y = y_,
id = id,
crs = 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"))
#' this should group by id, extract the home range, convert to polygon and measure the overlap
trk %>% group_by(id) %>%
hr_kde(., levels = c(.95)) %>%
hr_isopleths(.) %>%
st_intersection(., merged_Africa_tranform) %>%
st_area(.) / 1e6 #' only produces one value which is the same if id is ignored
我来回答你的问题..
# load packages
library(tidyverse)
library(amt)
library(sf)
library(adehabitatHR)
#loading and transforming shapefile can be simplified
merged_Africa <- read_sf("./temp/WDPA_Jun2020_SWZ-shapefile-polygons.shp") %>%
st_set_crs( 4326 ) %>%
st_transform( 102022 )
#mydata comes from sample data provided
#loop over the id's, crate track by ID, and calulate (and sum) overlapping area's
L <- lapply( unique( mydata$id ), function(x) {
track_id <- mk_track( mydata[ id == x, ],
.x = x_,
.y = y_,
id = id,
crs = 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"))
track_object <- track_id %>% hr_kde(., levels = c(.95)) %>% hr_isopleths(.)
area_totals <- track_object %>% st_intersection( merged_Africa ) %>% st_area()
sum( area_totals )
})
#set names
names(L) <- unique( mydata$id )
L
# $A
# 34689637 [m^2]
#
# $B
# 38990619 [m^2]
我有 GPS 数据,可用于确定家庭范围。我想测量范围与受保护区域形状文件 (available here) 的重叠。当我有一个人时这很好用,但当我尝试使用多个 ID 时却不行。
问题似乎出在这段代码中的 group_by
函数上:
trk %>% group_by(id) %>%
hr_kde(., levels = c(.95)) %>%
hr_isopleths(.) %>%
st_intersection(., merged_Africa_tranform) %>%
st_area(.) / 1e6
它只会产生一个值 (29.77905),而我期望两个,每个 ID 一个(A 为 34.68964,B 为 38.99062)。
完整代码如下:
# load packages
library(tidyverse)
library(amt)
library(sf)
library(adehabitatHR)
#' load in the protected areas (working with Albers Equal Area)
merged_Africa = read_sf("shape//eswatini//WDPA_Apr2020_SWZ-shapefile-polygons.shp")
st_crs(merged_Africa) <- 4326
merged_Africa_tranform <- st_transform(merged_Africa, "ESRI:102022")
#' try a data frame with multiple IDs
x_ <- c(707692, 707589, 707998, 708407, 708916, 709415, 707743, 707429, 707971, 708143, 708981, 709156)
y_ <- c(-3030991,-3031423,-3031640,-3031750,-3032508,-3037158,-3030995,-3031723,-3031680,-3031755,-3032408,-3037758)
id <- c(rep("A",6), rep("B", 6))
mydata <- data.frame(x_,y_,id)
# transform to trk object
trk <-
mk_track(mydata,
.x = x_,
.y = y_,
id = id,
crs = 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"))
#' this should group by id, extract the home range, convert to polygon and measure the overlap
trk %>% group_by(id) %>%
hr_kde(., levels = c(.95)) %>%
hr_isopleths(.) %>%
st_intersection(., merged_Africa_tranform) %>%
st_area(.) / 1e6 #' only produces one value which is the same if id is ignored
我来回答你的问题..
# load packages
library(tidyverse)
library(amt)
library(sf)
library(adehabitatHR)
#loading and transforming shapefile can be simplified
merged_Africa <- read_sf("./temp/WDPA_Jun2020_SWZ-shapefile-polygons.shp") %>%
st_set_crs( 4326 ) %>%
st_transform( 102022 )
#mydata comes from sample data provided
#loop over the id's, crate track by ID, and calulate (and sum) overlapping area's
L <- lapply( unique( mydata$id ), function(x) {
track_id <- mk_track( mydata[ id == x, ],
.x = x_,
.y = y_,
id = id,
crs = 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"))
track_object <- track_id %>% hr_kde(., levels = c(.95)) %>% hr_isopleths(.)
area_totals <- track_object %>% st_intersection( merged_Africa ) %>% st_area()
sum( area_totals )
})
#set names
names(L) <- unique( mydata$id )
L
# $A
# 34689637 [m^2]
#
# $B
# 38990619 [m^2]