Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

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

Open
dfifield opened this issue Feb 12, 2022 · 0 comments
Open

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

dfifield opened this issue Feb 12, 2022 · 0 comments

Comments

@dfifield
Copy link

Hi Phil,

I discovered that proj_rasters() fails if either ll$univariate or ll$analogue are NULL. I have a dataset where there is no univariate extrapolation so ll$univariate will be NULL. 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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

1 participant