jmsigner / amt

37 stars 13 forks source link

Problem with `hr_isopleths()` when aKDE has >1 level #56

Closed bsmity13 closed 2 years ago

bsmity13 commented 2 years ago

Demo of problem:

library(amt)

dat <- amt_fisher %>% 
  filter(id == "F1")

# Works fine:
akde <- hr_akde(dat)
hr_isopleths(akde)

# Doesn't work (levels argument ignored):
hr_isopleths(akde, levels = c(0.95, 0.5))

# Doesn't work (bug in hr_isopleth.akde)
akde2 <- hr_akde(dat, levels = c(0.95, 0.5))
hr_isopleths(akde2) # Error

# amt:::hr_isopleths.akde
x <- akde2
conf.level <- 0.95
# function (x, conf.level = 0.95, ...) 
# {
  # This is fine
  checkmate::assert_number(conf.level, lower = 0, upper = 1)
  res <- ctmm::SpatialPolygonsDataFrame.UD(x$akde, level.UD = x$levels, 
                                           conf.level = conf.level)
  res1 <- sf::st_as_sf(res)
  res1 <- sf::st_transform(res1, akde$crs)
  res1 <- res1[, setdiff(names(res1), "name")]

  # Error here
  # Should each always be 3? (lci, estimate, uci)
  # Alternatively, should we keep res1$name to extract "level" and "what"?
  res1$level <- rep(x$level, each = nrow(res1))
  res1$what <- rep(c(paste0("lci (", conf.level, ")"), 
                     "estimate", paste0("uci (", conf.level, ")")), 
                   length(x$levels))

  # Rest of code is fine
  res1$area = sf::st_area(res1)
  res1[, c("level", "what", "area", "geometry")]
# }

Potential solution (might not work if length(conf.level) > 1)

# function (x, conf.level = 0.95, ...) 
# {
  checkmate::assert_number(conf.level, lower = 0, upper = 1)
  res <- ctmm::SpatialPolygonsDataFrame.UD(x$akde, level.UD = x$levels, 
                                           conf.level = conf.level)
  res1 <- sf::st_as_sf(res)
  res1 <- sf::st_transform(res1, akde$crs)

  ## Proposed fix
  # Extract level and what from raw names
  split_names <- strsplit(res1$name, " ", fixed = TRUE)

  level_perc <- sapply(split_names, getElement, 2)
  est <- sapply(split_names, getElement, 3)

  res1$level <- as.numeric(gsub(pattern = "%", replacement = "", 
                     x = level_perc, fixed = TRUE))/100
  res1$what <- ifelse(est == "low", paste0("lci (", conf.level, ")"),
                      ifelse(est == "est", "estimate",
                      paste0("uci (", conf.level, ")")))
  res1$name <- NULL
  row.names(res1) <- NULL
  ##

  res1$area = sf::st_area(res1)
  res1[, c("level", "what", "area", "geometry")]
# }
jmsigner commented 2 years ago

Thanks for reporting, this now fixed 38177f7