densitymodelling / dsmextra

Extrapolation assessments in density surface models
GNU Lesser General Public License v3.0
4 stars 3 forks source link

proj_rasters fails if univariate (or analogue) are NULL #13

Open dfifield opened 2 years ago

dfifield commented 2 years ago

Hi Phil,

I discovered that proj_rasters() fails if eitherll$univariate or ll$analogue are NULL. I have a dataset where there is no univariate extrapolation so ll$univariate will beNULL. On another occasion, I messed up the dataframes passed to extrapolation_analysis() and ended up with no analogue either, and proj_rasters() failed in this case too.

As a workaround in my local copy of the package, I have passed the "types" object from map_extrapolation() to proj_rasters ().

Then within proj_rasters() I wrap two parts of the code in a conditional to test the value of types. Not sure if this is sensible or not, but it seems to fix the problem.

proj_rasters <- function(ll, coordinate.system, types){

  suppressWarnings(crs.webmerc <- sp::CRS("+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext  +no_defs"))

  llr <- ll # Copy list

  if(all(c("univariate", "analogue") %in% types)) {
      # Univariate extrapolation is negative by definition
      # When only a small number of cells are subject to UE, the resampling
      # may result in the loss of some of them.
      # By recording the indices of UE cells, we can perform a simplistic
      # correction to make sure they show up on the map.

      analogue.xy <- raster::as.data.frame(llr$analogue, xy = TRUE) %>% stats::na.omit(.)
      analogue.xy <- sp::SpatialPointsDataFrame(coords = analogue.xy[, c("x", "y")],
                                                data = analogue.xy,
                                                proj4string = coordinate.system)
      analogue.xy <- sp::spTransform(analogue.xy, CRSobj = crs.webmerc)

      univariate.ind <- raster::Which(llr$univariate < 0, cells = TRUE)
      univariate.values <- llr$univariate[univariate.ind]

      univariate.xy <- raster::as.data.frame(llr$univariate, xy = TRUE) %>% stats::na.omit(.)
      univariate.xy <- sp::SpatialPointsDataFrame(coords = univariate.xy[, c("x", "y")],
                                                  data = univariate.xy,
                                                  proj4string = coordinate.system)
      univariate.xy <- sp::spTransform(univariate.xy, CRSobj = crs.webmerc)
    }

    llr$all <- NULL
    llr <- purrr::discard(.x = llr, is.null)

    suppressWarnings(
    llr <- purrr::map(.x = llr, # Same extent as the full raster, allows correct alignment
                      .f = ~raster::projectRaster(from = .x,
                                                  to = ll$all,
                                                  method = 'ngb')) %>%
      purrr::map(.x = ., # CRS used by leaflet
                 .f = ~raster::projectRaster(from = .,
                                             crs = crs.webmerc,
                                             method = 'ngb')))

  if(all(c("univariate", "analogue") %in% types)) {
    llr.univariate.ind <- raster::cellFromXY(object = llr$univariate,
                                             xy = sp::coordinates(univariate.xy))

    llr$univariate[llr.univariate.ind[which(is.na(llr$univariate[llr.univariate.ind]))]] <-    univariate.values[which(is.na(llr$univariate[llr.univariate.ind]))]

    r1 <- raster::as.data.frame(llr$univariate, xy = TRUE)
    r2 <- raster::as.data.frame(llr$analogue, xy = TRUE)
    names(r1) <- names(r2) <- c("x", "y", "ExDet")

    duplicate.cells <- rbind(r1, r2) %>%
      stats::na.omit(.) %>%
      dplyr::select(., x, y) %>%
      .[duplicated(.),]

    llr.analogue.ind <- raster::cellFromXY(object = llr$analogue,
                                           xy = duplicate.cells)

    llr$analogue[llr.analogue.ind] <- NA
  }

  return(llr)}

Best! Dave