error in scales::colour_ramp

[eb-bruns]

Hello! I'm trying to run the SummaryHTML function by following the Cucurbita example but using occurrence points and raster files for my own target species. Everything prior in the Cucurbita example has worked using my data, but I get this error within the "Calculating ERSex gap map" step:

Quitting from lines 439-458 (summaryHTML.Rmd) 
Error in scales::colour_ramp(colors, alpha = alpha) : 
  Must provide at least one colour to create a colour ramp

I'm using my own occurrence data and your raster files for J. californica and J. hindsii from the 2020 PNAS paper, found here: https://github.com/dcarver1/cwr_pnas_results/tree/main/speciesLevelData20210706/Juglans

Thanks and let me know what other information you may need!

[dcarver1]

@esbeckman Hi, thanks for the details. I would recommend grabbing the .rmd file that the summaryHTML() function is calling directly. run system.file(package = "GapAnalysis") to get the folder path for the package. The rmd will be in the data/preloaded_data folder

This error is connected to the Insitu ERS maps generation base on the line reference in your error. I can't tell you the exact cause but my guess is that no ecoregions are being captured in the creation of the gEco variable. What's a bit odd is that the code to create the map is wrapped in a try() function so it shouldn't be quitting the creation of the html.

Please run the rmd outside of the summaryHTML function and report back with addational context on the exact location in which the process is failing. Should be better able to get at a work around from there.

[eb-bruns]

Hi @dcarver1 - thanks for the instructions. I've tried my example directly in the summaryHTML.Rmd file, and the error appears in this section:

if(length(gEco)>30){
  try(tm_shape(r1, name = "Native area and distribution model where ecoregions conserved in situ")+
  tm_raster(style = "pretty", palette = ersBPal, title = "Map of ERSin",
            labels = c("0: Native area", "1: Potential distribution where ecoregions conserved in situ"))+
tm_shape(g1, name = "Ecoregions not conserved in situ")+
  tm_raster(style = "cat", palette = ecoPal, title = "Ecoregions not conserved in situ"))
}else{
  try(tm_shape(r1, name = "Native area and distribution model where ecoregions conserved in situ")+
  tm_raster(style = "cat", palette = ersBPal, title = "Map of ERSin",
            labels = c("0: Native area", "1: Potential distribution where ecoregions conserved in situ"))+
tm_shape(g1, name = "Ecoregions not conserved in situ")+
  tm_raster(style = "cat", palette = ecoPal, title = "Ecoregions not conserved in situ"))
}

My length(gEco) is 1, so it runs the 'else' section and throws the 'Error in scales::colour_ramp' I quoted above. I checked gEco and it is indeed NA.

[dcarver1]

@esbeckman

We need to trace this error back a bit. the gEco object is define on line 267 and reference object g1

g1 <- exERS$gap_maps[[1]]
gEco <- unique(raster::values(g1))

g1 is define from the exERS function. You can view the details of this function on github or access it in R folder of the system.file(package = "GapAnalysis") file.

Try running this function in isolation and pay special attention to the lines 197 to 219

  if(isTRUE(Gap_Map)){
          message(paste0("Calculating ERSex gap map for ",as.character(Species_list[i])),"\n")

          # ERSex Gap Map
          # select all ecoregions present in ecoVal(all points) but absent in ecoValG(g buffers)
          ecoGap <- ecoVals[!ecoVals %in% ecoValsG]
          if(length(ecoGap) == 0){
            GapMapEx_list[[i]] <- paste0("All ecoregions within the model are within ", Buffer_distance,
                                         "km of G occurrence. There are no gaps")

          }else{
            # pull selected ecoregions and mask to presence area of the model
            eco2 <- Ecoregions_shp[Ecoregions_shp$ECO_ID_U %in% ecoGap,]
            #convert to sf object for conversion using fasterize
            eco2a <- sf::st_as_sf(eco2, SdmMask)
            # generate a ecoregion raster keeping the unique id.
            eco3 <- fasterize::fasterize(eco2a, SdmMask, field = "ECO_ID_U")
            # mask so only locations within the predicted presence area is included.
            gap_map <- eco3 * SdmMask
            GapMapEx_list[[i]] <- gap_map
            names(GapMapEx_list[[i]] ) <- Species_list[[i]]
          }
        }

If a NA value is being returned I expect something within the else statement if failing. Get to the bottom of that and we can take the next steps.

[eb-bruns]

Thanks, @dcarver1! I'm still trying to figure out exactly what's going wrong, but I did discover that the error is in the ERSin function (in situ rather than ex situ). It looks like all ecoregions have a protected area, so I believe there are no gaps and therefore the NA value for g1 is correct? The try() function is meant to catch this exception? I went through the ERSin function line by line and no errors are thrown. Here is the output of the ERSin function when I run it in isolation:

$ERSin
              species ERSin
1 Juglans_californica   100
2     Juglans_hindsii   100

$gap_maps
$gap_maps[[1]]
class      : RasterLayer 
dimensions : 180, 360, 64800  (nrow, ncol, ncell)
resolution : 0.04189815, 0.08819444  (x, y)
extent     : -124.4167, -109.3333, 27.83333, 43.70833  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs 
source     : memory
names      : layer 
values     : NA, NA  (min, max)

$gap_maps[[2]]
class      : RasterLayer 
dimensions : 180, 360, 64800  (nrow, ncol, ncell)
resolution : 0.03576389, 0.08703704  (x, y)
extent     : -124.4167, -111.5417, 28.04167, 43.70833  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs 
source     : memory
names      : layer 
values     : NA, NA  (min, max)

[dcarver1]

@esbeckman ,

Thanks for the troubleshooting on your end. Here are a few thoughts.

1. 100% ERSin species is great news! -- it also makes for a very boring map. You can see some examples here; . No gaps to show. Given that you might consider dropping the creation of the map for those species. You can do that by setting Gap_Map = FALSE, in the generation of the inERS dataset.

   inERS <- GapAnalysis::ERSin(Species_list = Sl, Occurrence_data = Od, 
                           Raster_list = Rl, Pro_areas = Pro_areas, Ecoregions_shp = Ecoregions_shp,
                           Gap_Map = TRUE)

   This should allow the rest of the content to generate and there will be no map for that feature.

2. We really need an additional clause to capture this condition. You should try replacing the code in your summaryHTML for ERSin function with the following. I added a ifelse statement to render a map without any ecoregion elements if the only value present is NA. I did not test this locally, so please reach back out and let me know if it worked for you.

g1 <- inERS$gap_maps[[1]]
gEco <- unique(raster::values(g1))
if(length(gEco)==1 & is.na(gEco)){
  gEco
}
g1[g1 ==0] <- NA

if(length(gEco)>30){
  try(tm_shape(r1, name = "Native area and distribution model where ecoregions conserved in situ")+
  tm_raster(style = "pretty", palette = ersBPal, title = "Map of ERSin",
            labels = c("0: Native area", "1: Potential distribution where ecoregions conserved in situ"))+
tm_shape(g1, name = "Ecoregions not conserved in situ")+
  tm_raster(style = "cat", palette = ecoPal, title = "Ecoregions not conserved in situ"))
  }else{
    ### adding new conditions 
    ifelse(test = length(gEco)== 1 & is.na(gEco[1]),
           yes = try(tm_shape(r1, name = "Native area and distribution model where ecoregions conserved in situ")+
  tm_raster(style = "cat", palette = ersBPal, title = "Map of ERSin",
            labels = c("0: Native area", "1: Potential distribution where ecoregions conserved in situ"))),
          no = try(tm_shape(r1, name = "Native area and distribution model where ecoregions conserved in situ")+
  tm_raster(style = "cat", palette = ersBPal, title = "Map of ERSin",
            labels = c("0: Native area", "1: Potential distribution where ecoregions conserved in situ"))+
tm_shape(g1, name = "Ecoregions not conserved in situ")+
  tm_raster(style = "cat", palette = ecoPal, title = "Ecoregions not conserved in situ"))
           )
  }

[dcarver1]

Need additional condition when ERSin and ERSex to account for when all ecoregions are effectively conserved. Should be able to use the ERS score of 100 as the testing condition.

[eb-bruns]

Hi @dcarver1 - sorry for the delay! I’ve tried a few things and wanted to make sure I mostly understand what’s happening before asking more questions!

First, a separate but related issue: I noticed that the raster resolution on the maps wasn’t looking right, and realized the GetDatasets() function had downloaded the 10 arc minutes version of the protected areas instead of the 2.5 arc mins version which matches the resolution of the SDMs I’m using. So I downloaded the wdpa_reclass.tif from here and used that instead – map resolution looked good after that.

After realizing the resolution-difference issue, I decided to try the SummaryHTML function again with 2.5 arc mins protected areas layer. This time, there was not the scales::colour_ramp error! It may indeed still be a logical loophole, but not for my test species anymore :). Before I fixed the resolution issue, I did try running your fix to the ERSin function. It did let the SummaryHTML run without an error, but outputted a block of script instead of the map and still did that when I added a species that I don’t think fit the exception (i.e. should have produced a map). See image below for the output.

[eb-bruns]

So, I have another question now @dcarver1. I'll put it in this thread but you can move to a separate issue as desired!

With the SummaryHTML function running, I've looked more closely at the output and realized the map labels are not quite right. After a deeper look, I've determined that the rasters I’m using are missing the “Native area (ecoregions within countries where species has been observed in the field)” designation – they only have the modeled potential distribution from Khoury et al 2020. So, my maps are missing one element in the legend, which shifts the remaining parts up so they are incorrect. For example, the purple designation is missing in the GRSex map – below is my Juglans major map and then one of the Cucurbita examples from the README for comparison.

I (wrongly) assumed that the "native area" designation was added somewhere within the GapAnalysis package workflow, since it wasn't mentioned explicitly in the README. Wondering what you suggest as the easiest workaround? I'm guessing I should just add the "native area" myself? I'm not as proficient with rasters, so if you have relevant code to share, that would be amazing. Otherwise, I'm sure I can figure it out.

Thanks, as always!

[dcarver1]

@eb-bruns, thanks for provided a detailed example.

The native area designation was a product of the modeling processing associated with the original modeling effort and truthful should have have been brought into this project.

If you want to do some direct edits on your end I can provide an example of what changes would be need to be made to make the corrections. There is a lot of repeated code between the map elements so it'll be time consuming but not conceptually difficult to apply those changes to other features.

We're working on trying to find some time to address some of these larger issues that you are running into but it's unlike to happen in 2022.