Issue in function prepTC

[Ohm-Np]

While running the code below:

library(sf)
library(terra)
remotes::install_github("mapme-initiative/mapme.forest", ref = "terra")
library(mapme.forest)

# ---- test download function and prepare data -----
# load raster files
raster_files <- rast(paste0("../../datalake/mapme.protectedareas/input/global_forest_watch/",
                            list.files("../../datalake/mapme.protectedareas/input/global_forest_watch/")))
# load region of interest
aoi <-
  read_sf("../../datalake/mapme.protectedareas/input/wdpa_kfw/wdpa_kfw_spatial_latinamerica_2021-02-01_supportedPAs_unique.gpkg")
# we need to reproject the aoi object
aoi <-
  st_transform(aoi, crs = st_crs(raster_files))

# create function
area_comp <- function(aoi) {

  # crop the rasters to the extent of the aoi
  rasters <- terra::crop(raster_files, aoi)
  # assigning proper names to the layers and simply plot the results
  names(rasters) <- c("co2_emission", "lossyear", "treecover")

  #----- test prepTC -----
  prep_treeCover <- prepTC(
    inputForestMap = rasters$treecover,
    thresholdCover = 10,
    thresholdClump = 6)

  # ----- test Yearly forest mask -----
  yearlyForestMask <- getTM(
    inputForestMap = prep_treeCover,
    inputLossMap = rasters$lossyear,
    years = 2001:2020)

  # ---- test Area Calc -----
  result_latlon <- AreaCalc(yearlyForestMask,
                            studysite = aoi,
                            latlon = TRUE,
                            polyName = "WDPA_PID",
                            ncores = 1,
                            years = 2001:2020)
  # ---- remove other columns ----
  result <- result_latlon %>%
    select(-c(2:50))
  # ---- return results ----
  return(result)
}

### run function for first polygon **WORKS FINE**
area_comp(aoi[1, ])

### run function for the fourth polygon **THROWS ERROR**
area_comp(aoi[4, ]) 

# Error: [%in%] no matches supplied 

# 8. stop("[", f, "] ", emsg, ..., call. = FALSE) 
# 7. error(f, x@ptr$getError()) 
# 6. messages(x, "%in%") 
# 5. clummy %in% clumpTmp 
# 4. clummy %in% clumpTmp 
# 3. `[<-`(`*tmp*`, clummy %in% clumpTmp, value = 0) 
# 2. prepTC(inputForestMap = rasters$treecover, thresholdCover = 10, 
#        thresholdClump = 6) 
# 1. area_comp(aoi[4, ])

Seems like the error comes from the prepTC function.
@goergen95 , could you please help me here in this part?

[goergen95]

Hi @Ohm-Np , difficult to reproduce the issue because I don't have access to your aoi object. I suspect something weird about the 4th geometry but the error message is quite cryptic. Could you try an aoi <- st_make_valid(aoi) and maybe inspect the the cropped rasters for the 4th element manually, e.g. a plot of the rasters? Maybe the geometry is empty/invalid or there are no values in the resulting rasters.

[Ohm-Np]

Thank you @goergen95.

The geometry is valid and also the cropped rasters have some valid values (min 0 - max 100).

However, I found one weird thing while processing that fourth polygon though. The script works fine when I set the values of thresholdCover to 75 and thresholdClump to 20 but doesn't work when cover is 10 & clump is 6.

[Ohm-Np]

When we provide thresholdCover value - 10, almost every cell of clummy raster for this 4th polygon was shown as forest due to which the dataframe clumpTmp shows only one row. Thus, the clump removal could not be applied in this case and throws this error: # Error: [%in%] no matches supplied

To avoid such ideal condition of having 100% forest cover, I made the following change (add if-else statement) in the prepTC function and it seems to work as of now. But I am not 100% sure if the logic I am following here is correct or not!

    clummy = patches(inputForestMap, directions = 8, zeroAsNA=TRUE)
    clumpTmp = data.frame(freq(clummy))

    # added if statement here in this line - if there are multiple rows, only then it would process further lines of code
    if (nrow(clumpTmp) > 1) {

      clumpTmp = clumpTmp[which(clumpTmp$count < thresholdClump), ]
      clumpTmp = as.vector(clumpTmp$value)
      clummy[clummy %in% clumpTmp] = 0
      clummy[clummy != 0] = 1
      clummy[is.na(clummy)] = 0
      inputForestMap = clummy

    } 
   # if there is only one row, then whole clummy raster would be returned as inputForestMap i.e. in an ideal case of 100% forest cover
   else {  
      inputForestMap = clummy
    }

[goergen95]

Thank you for figuring this out! Very interesting edge case!! Your solution seems fine, though I guess it would also return the raster as is if it only has values of 0. Which in any case is also fine, just that it is the opposite edge case. Let's figure out soon when and where to put your changes (but you could open a PR anytime).

[Ohm-Np]

I guess it would also return the raster as is if it only has values of 0

Thank you @goergen95 for pointing this out. I will try to adapt this case too.