R/nightlights.R

In chrisvwn/Rnightlights: Satellite Nightlight Data Extraction

######################## processNLCountry ###################################

#' Processes nightlights for an individual country in a particular nlPeriod
#'
#' Given a \code{countryCode}, \code{yearMonth} and preferred processing
#'     methods \code{cropMaskMethod} and \code{extractMethod}, this function
#'     will first check if the data already exists in the cache. First it
#'     will check if the data file exists and if it does not it will create
#'     a dataframe of the country data containing only the administrative
#'     properties and move to processing. If the data file exists it will
#'     check to see if the particular year month already exists. If it
#'     exists, it will exit with a message. If it does not exist, it will
#'     load the country data file and move on to processing.
#'
#'     Processing consists of:
#'     \itemize{
#'        \item Reading in the country polygon in ESRI Shapefile format
#'        \item Reading in the tiles that the particular country intersects
#'            with and then clipping
#'        the tile(s) to the country boundary
#'        \item Extract the data from the clipped raster and compute various
#'            statistics at the lowest admin level in the country.
#'        \item Finally, these stats are appended to the data frame and
#'            written to the data file.
#'     }
#'
#'     NOTE: \code{processNLCountry()} assumes that all inputs are available
#'     and will not attempt to download them. It should ideally be called
#'     from the function \code{processNlData()} which does all the
#'     preparation for processing. \code{processNlData()} which can process
#'     multiple countries and time periods will download all the required
#'     tiles and polygons prior to calling \code{processNlCountry}.
#'     \code{getCtryNlData} can also be used with the option
#'     \code{ignoreMissing=FALSE} which will call \code{processNlData}
#'     in the background.
#'
#' @param ctryCode \code{character} The ctryCode of interest
#'
#' @param admLevel \code{character} The country admin level in the given
#'     ctryCode at which to calculate stats
#'
#' @param nlType \code{character} The nlType of interest
#'
#' @param configName character the config shortname of raster being processed
#' 
#' @param extension character the extension of raster being processed
#'
#' @param multiTileStrategy character How to handle multiple tiles per nlPeriod
#'
#' @param multiTileMergeFun character The function to use to merge tiles
#'
#' @param removeGasFlaresMethod character The method to use to perform gas flare removal
#'     or NULL to disable
#'     
#' @param nlPeriod \code{character} The nlPeriod of interest
#'
#' @param nlStats the statistics to calculate. If not provided will calculate
#'     the stats specified in \code{pkgOptions("nlStats")}
#'
#' @param downloadMethod The method used to download polygons and rasters
#'
#' @param cropMaskMethod \code{character} Whether to use rasterize or
#'     gdal-based functions to crop and mask the country rasters
#'
#' @param extractMethod ("rast" or "gdal") Whether to use rasterize or
#'     gdal-based functions to crop and mask the country rasters
#'
#' @param gadmVersion The GADM version to use
#'
#' @param gadmPolyType The format of polygons to download from GADM
#'
#' @param custPolyPath Alternative to GADM. A path to a custom shapefile zip
#'
#' @return None
#'
#' @examples
#'
#' #calculate only the sum of monthly VIIRS radiances for Dec 2014 using gdal
#' #for both cropMask and extraction for KEN
#' \dontrun{
#' Rnightlights:::processNLCountry("KEN", "KEN_adm2", "VIIRS.M", "201412", "gdal", "gdal", "sum")
#' }
#'
processNLCountry <- function(ctryCode,
                             admLevel,
                             nlType,
                             configName = pkgOptions(paste0("configName_", nlType)),
                             extension,
                             multiTileStrategy = pkgOptions("multiTileStrategy"),
                             multiTileMergeFun = pkgOptions("multiTileMergeFun"),
                             removeGasFlaresMethod = pkgOptions(paste0("removeGasFlaresMethod_", nlType)),
                             nlPeriod,
                             nlStats = pkgOptions("nlStats"),
                             downloadMethod = pkgOptions("downloadMethod"),
                             cropMaskMethod = pkgOptions("cropMaskMethod"),
                             extractMethod = pkgOptions("extractMethod"),
                             gadmVersion = pkgOptions("gadmVersion"),
                             gadmPolyType = pkgOptions("gadmPolyType"),
                             custPolyPath = NULL)
{
  if (missing(ctryCode))
    stop(Sys.time(), "Missing required parameter ctryCode")
  
  if (missing(nlPeriod))
    stop(Sys.time(), "Missing required parameter nlPeriod")
  
  if (missing(nlType))
    stop(Sys.time(), "Missing required parameter nlType")
  
  if (!validCtryCodes(ctryCodes = ctryCode))
    stop(Sys.time(), "Invalid ctryCode: ", ctryCode)
  
  if (!allValidNlPeriods(nlPeriods = nlPeriod, nlTypes = nlType))
    stop(Sys.time(),
         "Invalid nlPeriod: ",
         nlPeriod,
         " for nlType ",
         nlType)
  
  if (!validNlTypes(nlTypes = nlType))
    stop(Sys.time(), "Invalid nlType: ", nlType)
  
  if (missing(admLevel))
    admLevel <- getCtryShpLowestLyrNames(
      ctryCodes = ctryCode,
      gadmVersion = gadmVersion,
      gadmPolyType = gadmPolyType,
      custPolyPath = custPolyPath
    )
  
  if (is.list(nlStats) &&
      length(nlStats) > 1 &&
      all(sapply(2:length(nlStats),
                 function(i)
                   ! is.list(nlStats[[i]]) &&
                 (grepl("=", nlStats[i]) ||
                  length(names(nlStats[i])) > 0))))
    nlStats <- list(nlStats)
  
  funArgs <- match.call()[-1]
  
  message(
    Sys.time(),
    ": ProcessNlCountry: ",
    paste(
      names(funArgs),
      sapply(names(funArgs), function(x)
        eval(parse(text = x))),
      sep = "=",
      collapse = ", "
    ),
    "****"
  )
  
  wgs84 <- getCRS()
  
  message(Sys.time(), ": Check for existing data file")
  
  ctryNlDataFnamePath <- getCtryNlDataFnamePath(
    ctryCode = ctryCode,
    admLevel = admLevel,
    gadmVersion = gadmVersion,
    gadmPolyType = gadmPolyType,
    custPolyPath = custPolyPath
  )
  
  if (existsCtryNlDataFile(
    ctryCode = ctryCode,
    admLevel = admLevel,
    gadmVersion = gadmVersion,
    gadmPolyType = gadmPolyType,
    custPolyPath = custPolyPath
  ))
  {
    message(Sys.time(), ": Data file found: ", ctryNlDataFnamePath)
    
    existStats <- NULL
    
    for (i in seq_along(nlStats))
    {
      existStat <- existsCtryNlData(
        ctryCode = ctryCode,
        admLevel = admLevel,
        nlTypes = nlType,
        configNames = configName,
        extensions = extension,
        multiTileStrategy = multiTileStrategy,
        multiTileMergeFun = multiTileMergeFun,
        removeGasFlaresMethod = removeGasFlaresMethod,
        nlPeriods = nlPeriod,
        nlStats = nlStats[i],
        gadmVersion = gadmVersion,
        gadmPolyType = gadmPolyType,
        custPolyPath = custPolyPath
      )
      
      existStats <- c(existStats, existStat)
    }
    
    if (all(existStats))
    {
      message(
        Sys.time(),
        ": **All stats exist for ",
        paste(ctryCode, admLevel, nlPeriod, sep = " "),
        ". Skipping"
      )
      
      return(-1)
    }
    else
    {
      nlStats <- nlStats[!existStats]
      message(Sys.time(),
              ": Processing stats: ",
              nlStatSignature(nlStats))
    }
    
    message(Sys.time(), ": Load country data file")
    ctryNlDataDF <-
      utils::read.csv(
        ctryNlDataFnamePath,
        header = TRUE,
        check.names = FALSE,
        encoding = "UTF-8"
      )
    
    message(Sys.time(), ": Load country polygon admin level")
    
    # ctryExtent <- raster::extent(x = ctryPolyAdm0)
    #
    # raster::projection(ctryPolyAdm0) <- sp::CRS(projargs = wgs84)
    
  } else
  {
    message(Sys.time(), ": Data file not found. Creating ...")
    
    ctryNlDataDF <- createCtryNlDataDF(
      ctryCode = ctryCode,
      admLevel = admLevel,
      gadmVersion = gadmVersion,
      gadmPolyType = gadmPolyType,
      custPolyPath = custPolyPath
    )
    
    message(Sys.time(), ": Data file not found. Creating ... DONE")
  }
  
  message(Sys.time(), ": Begin processing ", nlPeriod)
  
  #get the path we will use to save the cropped raster
  ctryRasterOutputFnamePath <-
    getCtryRasterOutputFnamePath(
      ctryCode = ctryCode,
      nlType = nlType,
      configName = configName,
      extension = extension,
      multiTileStrategy = multiTileStrategy,
      multiTileMergeFun = multiTileMergeFun,
      removeGasFlaresMethod = removeGasFlaresMethod,
      nlPeriod = nlPeriod,
      gadmVersion = gadmVersion,
      gadmPolyType = gadmPolyType,
      custPolyPath = custPolyPath
    )
  
  if (!file.exists(ctryRasterOutputFnamePath))
  {
    message(Sys.time(), ": Country output raster not found. Creating")
    
    message(Sys.time(), ": Load polygon layer for crop")
    ctryPolyAdm0 <- readCtryPolyAdmLayer(
      ctryCode = ctryCode,
      admLevel = unlist(
        getCtryShpLyrNames(
          ctryCodes = ctryCode,
          lyrNums = 0,
          gadmVersion = gadmVersion,
          gadmPolyType = gadmPolyType,
          custPolyPath = custPolyPath
        )
      ),
      gadmVersion = gadmVersion,
      custPolyPath = custPolyPath
    )
    
    # if(removeGasFlaresMethod == "OGP")
    # {
    #   message(Sys.time(), ": Gas flare removal (OLS Gas-Flare Polygons) is ON")
    #   
    #   #read in the mosaiced gas flare polys
    #   #does this country require gas flare removal
    #   if(hasNlCtryGasFlares(ctryCode = ctryCode,
    #                         gadmVersion = gadmVersion,
    #                         gadmPolyType = gadmPolyType,
    #                         custPolyPath = custPolyPath))
    #   {
    #     message(Sys.time(), ": Reading in gas-flare polygons")
    #     
    #     #replace ctryPolyAdm0 with a gas-flare corrected
    #     #ctryPolyAdm0 (from DMSP-OLS) which will be used to crop the 
    #     #raster
    #     ctryPolyAdm0 <- getNlCtryGasFlaresPoly(
    #       ctryCode = ctryCode,
    #       gadmVersion = gadmVersion,
    #       gadmPolyType = gadmPolyType,
    #       custPolyPath = custPolyPath
    #     )
    #     
    #   } else
    #   {
    #     message(Sys.time(),
    #             ": Gas flare removal not required for: ",
    #             ctryCode)
    #   }
    # }
    
    if(!missing(removeGasFlaresMethod) && !is.null(removeGasFlaresMethod) && 
       removeGasFlaresMethod == "OGP")
    {
      gfRemoved <- removeGasFlares(removeGasFlaresMethod = removeGasFlaresMethod,
                                     ctryPolyAdm0 = ctryPolyAdm0,
                                     ctryCode = ctryCode, nlType = nlType,
                                     nlPeriod = nlPeriod,gadmVersion = gadmVersion,
                                     gadmPolyType = gadmPolyType, custPolyPath = custPolyPath)
      
      ctryPolyAdm0 <- gfRemoved$ctryGFPolyAdm0
    }

    message(Sys.time(), ": Reading in the raster tiles")
    
    tileList <-
      getCtryTileList(ctryCodes = ctryCode, nlType = nlType)
    
    ctryRastCropped <- NULL
    
    for (tile in tileList)
    {
      rastFilename <- getNlTileTifLclNamePath(
        nlType = nlType,
        configName = configName,
        extension = extension,
        nlPeriod = nlPeriod,
        tileNum = tileName2Idx(tileName = tile,
                               nlType = nlType)
      )
      
      rastTile <- raster::raster(x = rastFilename)
      
      raster::projection(rastTile) <- sp::CRS(projargs = wgs84)
      
      ctryPolyAdm0 <- sp::spTransform(ctryPolyAdm0, sp::CRS(SRS_string = wgs84))
      
      message(Sys.time(), ": Cropping the raster tiles ")
      
      #extTempCrop <- crop(rastTile, ctryExtent)
      
      message(Sys.time(), ": Cropping tile = ", tile)
      
      #we crop using raster for both rast and gdal
      #Note: gdalwarp is not used for cropping because the crop_to_cutline option
      #causes a shift in the cell locations which then affects the stats extracted.
      #A gdal-based crop to extent would be highly desirable for performance reasons
      #though so seeking other gdal-based workarounds
      tempCrop <-
        raster::crop(x = rastTile,
                     y = ctryPolyAdm0,
                     progress = 'text')
      
      #will only be non-null if there are multiple tiles to be mosaiced i.e. a country
      #that straddles multiple tiles. So in the 2nd+ loop ctryRastCropped will not be null
      if (is.null(ctryRastCropped))
      {
        ctryRastCropped <- tempCrop
      }
      else
      {
        #if we are here we are mosaicing multiple tiles
        ctryRastMerged <- ctryRastCropped
        
        ctryRastCropped <- NULL
        
        #mosaic the tiles and store back in ctryRastCropped
        ctryRastCropped <-
          raster::merge(x = ctryRastMerged, y = tempCrop)
        
        rm(ctryRastMerged)
      }
      
      rm(tempCrop)
    }
    
    #unload the tile from memory
    rm(rastTile)
    
    #release unused memory
    gc()
    
    message(Sys.time(), ": Masking the raster ")
    
    if (cropMaskMethod == "rast")
    {
      #RASTERIZE
      message(Sys.time(), ": Convert to raster object using fasterize ")
      
      
      # Use fasterize to mask
      # Convert sf object to raster
      
      ctryPolyAdm0_raster <-
        fasterize::fasterize(sf = sf::st_as_sf(ctryPolyAdm0),
                             raster = ctryRastCropped)
      message(Sys.time(), ": Masking ")
      
      ctryRastCropped <-
        raster::mask(x = ctryRastCropped, mask = ctryPolyAdm0_raster)
      
      message(Sys.time(), ": Writing the raster to disk ")
      
      raster::writeRaster(
        x = ctryRastCropped,
        filename = ctryRasterOutputFnamePath,
        overwrite = TRUE,
        progress = "text"
      )
      
      message(Sys.time(), ": Crop and mask using rasterize ... Done")
      
    }
    else if (cropMaskMethod == "gdal")
    {
      message(Sys.time(), ": Crop and mask using gdalwarp ... ")
      
      #GDALWARP
      rstTmp <-
        file.path(getNlDir(dirName = "dirNlTemp"),
                  paste0(basename(tempfile()), ".tif"))
      
      message(Sys.time(),
              ": Writing merged raster to disk for gdalwarp masking")
      
      raster::writeRaster(x = ctryRastCropped,
                          filename = rstTmp,
                          progress = "text")
      
      ctryRastCropped <- NULL
      
      gc()
      
      #the polygon will already correspond to having gas flares removed or not
      ctryPolyAdm0TmpDir <- tools::file_path_sans_ext(rstTmp)
      
      rgdal::writeOGR(
        obj = methods::as(ctryPolyAdm0, "SpatialPolygonsDataFrame"),
        dsn = ctryPolyAdm0TmpDir,
        driver = "ESRI Shapefile",
        layer = "GID_0_IDX"
      )
      
      outputFileVrt <-
        file.path(
          getNlDir(dirName = "dirNlTemp"),
          paste0(ctryCode, "_", nlType, "_", nlPeriod, ".vrt"),
          fsep =
        )
      
      if (file.exists(outputFileVrt))
        file.remove(outputFileVrt)
      
      message(Sys.time(), ": gdalwarp masking to VRT")
      
      gdalUtils::gdalwarp(
        srcfile = rstTmp,
        dstfile = outputFileVrt,
        s_srs = wgs84,
        t_srs = wgs84,
        cutline = ctryPolyAdm0TmpDir,
        cl = "GID_0_IDX",
        multi = TRUE,
        wm = pkgOptions("gdalCacheMax"),
        wo = paste0("NUM_THREADS=",
                    pkgOptions("numThreads")),
        q = FALSE
      )
      
      message(Sys.time(), ": gdal_translate converting VRT to TIFF ")
      gdalUtils::gdal_translate(co = "compress=LZW",
                                src_dataset = outputFileVrt,
                                dst_dataset = ctryRasterOutputFnamePath)
      
      message(Sys.time(), ": Deleting the component rasters ")
      
      file.remove(rstTmp)
      file.remove(outputFileVrt)
      unlink(ctryPolyAdm0TmpDir,
             recursive = T,
             force = T)
      
      ctryRastCropped <- raster::raster(ctryRasterOutputFnamePath)
      
      #GDALWARP
      message(Sys.time(), ": Crop and mask using gdalwarp ... DONE")
    }
  }
  else
  {
    rastFilename <- ctryRasterOutputFnamePath
    
    message(Sys.time(), ": ", rastFilename, " already exists")
    
    ctryRastCropped <- raster::raster(x = rastFilename)
    
    raster::projection(x = ctryRastCropped) <-
      sp::CRS(SRS_string = wgs84)
  }
  
  if(!missing(removeGasFlaresMethod) && !is.null(removeGasFlaresMethod) && 
     removeGasFlaresMethod %in% c("OTM","VTM"))
  {
    gfRemoved <- removeGasFlares(removeGasFlaresMethod = removeGasFlaresMethod,
                                                ctryRastCropped = ctryRastCropped,
                                                ctryCode = ctryCode, nlType = nlType,
                                                nlPeriod = nlPeriod,gadmVersion = gadmVersion,
                                                gadmPolyType = gadmPolyType, custPolyPath = custPolyPath)
    
    ctryRastCropped <- gfRemoved$ctryGFMaskCropped
  }
  
  #message(Sys.time(), ": Create web version of raster")
  
  #gdal_translate -co COMPRESS=JPEG -co PHOTOMETRIC=YCBCR -co TILED=YES 5255C.tif 5255C_JPEG_YCBCR.tif
  
  #gdalUtils::gdal_translate(src_dataset = rastFilename,
  #                          dst_dataset = rastWebFilename,
  #                          co = "COMPRESS=JPEG PHOTOMETRIC=YCBCR TILED=YES")
  
  #gdaladdo --config COMPRESS_OVERVIEW JPEG --config PHOTOMETRIC_OVERVIEW YCBCR
  #--config INTERLEAVE_OVERVIEW PIXEL -r average 5255C_JPEG_YCBCR.tif 2 4 8 16
  
  #rastWebFilename <- file.path(getNlDir("dirRasterWeb"), basename(rastFilename))
  
  #gdalUtils::gdaladdo(filename = rastWebFilename, r = "average", levels = c(2, 4, 8, 16))
  
  #message(Sys.time(), ": End create web raster ")
  #system(cmd)
  
  message(Sys.time(), ": Begin extracting data from the raster ")
  
  ctryPoly <- readCtryPolyAdmLayer(
    ctryCode = ctryCode,
    admLevel = admLevel,
    gadmVersion = gadmVersion,
    gadmPolyType = gadmPolyType,
    custPolyPath = custPolyPath
  )
  
  if (extractMethod == "rast")
    sumAvgRad <- fnAggRadRast(
      ctryPoly = ctryPoly,
      ctryRastCropped = ctryRastCropped,
      nlType = nlType,
      configName = configName,
      extension = extension,
      nlStats = nlStats,
      custPolyPath = custPolyPath
    )
  else if (extractMethod == "gdal")
    sumAvgRad <- fnAggRadGdal(
      ctryCode = ctryCode,
      admLevel = admLevel,
      ctryPoly = ctryPoly,
      nlType = nlType,
      configName = configName,
      extension = extension,
      multiTileStrategy = multiTileStrategy,
      multiTileMergeFun = multiTileMergeFun,
      removeGasFlaresMethod = removeGasFlaresMethod,
      nlPeriod = nlPeriod,
      nlStats = nlStats,
      gadmVersion = gadmVersion,
      gadmPolyType = gadmPolyType,
      custPolyPath = custPolyPath
    )
  
  for (i in seq_along(nlStats))
  {
    nlStatName <- nlStats[[i]][[1]]
    
    ctryNlDataDF <- insertNlDataCol(
      ctryNlDataDF = ctryNlDataDF,
      dataCol = sumAvgRad[, nlStatName],
      nlStat = nlStats[i],
      nlPeriod = nlPeriod,
      nlType = nlType,
      configName = configName,
      extension = extension,
      multiTileStrategy = multiTileStrategy,
      multiTileMergeFun = multiTileMergeFun,
      removeGasFlaresMethod = removeGasFlaresMethod
    )
    
    #finally we are sure the nlStat is good and working,
    #save it
    message(Sys.time(),
            ": Saving nlStat '",
            nlStatSignature(nlStats[[1]]),
            "' ...")
    
    statSaved <- saveNlStat(nlStat = nlStats[[i]])
    
    if (statSaved)
      message(Sys.time(), ": Saving nlStat '", nlStatName, "' ... DONE")
  }
  
  message(Sys.time(), ": DONE processing ", ctryCode, " ", nlPeriod)
  
  message(Sys.time(), ": **COMPLETE. Writing data to disk")
  
  #Write the country data dataframe to disk
  saveCtryNlData(
    ctryNlDataDF = ctryNlDataDF,
    ctryCode = ctryCode,
    admLevel = admLevel,
    gadmVersion = gadmVersion,
    gadmPolyType = gadmPolyType,
    custPolyPath = custPolyPath
  )
  
  #release the cropped raster
  rm (ctryRastCropped)
  
  #delete temporary raster files
  raster::removeTmpFiles(h = 0)
}

######################## getCtryRasterOutputFname ###################################

#' Constructs the name of the output raster
#'
#' Constructs the name of the output raster
#'
#' @param ctryCode the ctryCode of interest
#'
#' @param nlType the nlType of interest
#'
#' @param configName character the config shortname of raster being processed
#' 
#' @param extension character the extension of raster being processed
#'
#' @param multiTileStrategy character How to handle multiple tiles per nlPeriod
#'
#' @param multiTileMergeFun character The function to use to merge tiles
#'
#' @param removeGasFlaresMethod logical Whether to perform gas flare removal pre-processing
#'
#' @param nlPeriod the nlPeriod of interest
#'
#' @param gadmVersion The GADM version to use
#'
#' @param gadmPolyType The format of polygons to download from GADM
#'
#' @param custPolyPath The path to a custom polygon as an alternative
#'     to using GADM polygons
#'
#' @return Character the name of country raster for a country and a given
#'     nlType and nlPeriod
#'
#' @examples
#'
#' Rnightlights:::getCtryRasterOutputFname(ctryCode = "KEN", nlType = "VIIRS.M", nlPeriod = "201412")
#'
getCtryRasterOutputFname <- function(ctryCode,
                                     nlType,
                                     configName = pkgOptions(paste0("configName_", nlType)),
                                     extension,
                                     multiTileStrategy = pkgOptions("multiTileStrategy"),
                                     multiTileMergeFun = pkgOptions("multiTileMergeFun"),
                                     removeGasFlaresMethod = pkgOptions(paste0("removeGasFlaresMethod_", nlType)),
                                     nlPeriod,
                                     gadmVersion = pkgOptions("gadmVersion"),
                                     gadmPolyType = pkgOptions("gadmPolyType"),
                                     custPolyPath = NULL)
{
  if (missing(ctryCode))
    stop(Sys.time(), "Missing required parameter ctryCode")
  
  if (missing(nlType))
    stop(Sys.time(), "Missing required parameter nlType")
  
  if (missing(nlPeriod))
    stop(Sys.time(), "Missing required parameter nlPeriod")
  
  if (!validCtryCodes(ctryCodes = ctryCode))
    stop(Sys.time(), "Invalid ctryCode: ", ctryCode)
  
  if (!validNlTypes(nlTypes = nlType))
    stop(Sys.time(), "Invalid nlType: ", nlType)
  
  if (!allValidNlPeriods(nlPeriods = nlPeriod, nlTypes = nlType))
    stop(Sys.time(),
         "Invalid nlPeriod: ",
         nlPeriod,
         " for nlType: ",
         nlType)
  
  if (missing(custPolyPath))
    custPolyPath <- NULL
  
  #configName <- toupper(configName)
  
  extension <- toupper(extension)
  
  fname <- if (is.null(custPolyPath))
    paste0(
      "NL_",
      ctryCode,
      "_",
      nlType,
      "_",
      nlPeriod,
      "_",
      configName,
      "_",
      extension,
      "-MTS",
      toupper(multiTileStrategy),
      "-",
      toupper(multiTileMergeFun),
      "-RGF",
      removeGasFlaresMethod,
      "_GADM-",
      gadmVersion,
      "-",
      toupper(gadmPolyType),
      ".tif"
    )
  else
    paste0(
      "NL_",
      ctryCode,
      "_",
      nlType,
      "_",
      nlPeriod,
      "_",
      configName,
      "_",
      extension,
      "-MTS",
      toupper(multiTileStrategy),
      "-",
      toupper(multiTileMergeFun),
      "-RGF",
      removeGasFlaresMethod,
      "_CUST-",
      basename(custPolyPath),
      "-SHPZIP.tif"
    )
  return (fname)
}


######################## getCtryRasterOutputFnamePath ###################################

#' Get the full path to the file where the cropped VIIRS country raster is stored.
#'
#' Get the full path to the file where the cropped VIIRS country raster is stored. This file is created
#'     when processing the country before extracting the data. It can be used to re-process a country much faster
#'
#' @param ctryCode (character) the ctryCode of interest
#'
#' @param nlType (character) the nlType of interest
#'
#' @param configName character the config shortname of raster being processed
#' 
#' @param extension character the extension of raster being processed
#'
#' @param multiTileStrategy (character) How to handle multiple tiles per nlPeriod
#'
#' @param multiTileMergeFun (character) The function to use to merge tiles
#'
#' @param removeGasFlaresMethod (character) The method to use to perform gas flare removal
#'
#' @param nlPeriod (character) the nlPeriod of interest
#'
#' @param gadmVersion (character) The GADM version to use
#'
#' @param gadmPolyType (character) The format of polygons to download from GADM
#'
#' @param custPolyPath (character) The path to a custom polygon as an alternative
#'     to using GADM polygons
#'
#' @return (character) full path to the cropped VIIRS country raster for a country and a given year and month
#'
#' @examples
#' \dontrun{
#' getCtryRasterOutputFnamePath("KEN","VIIRS.M", "201412")
#' }
#'
#'#export for gui() shiny app
#'@export
getCtryRasterOutputFnamePath <- function(ctryCode,
                                         nlType,
                                         configName,
                                         extension,
                                         multiTileStrategy = pkgOptions("multiTileStrategy"),
                                         multiTileMergeFun = pkgOptions("multiTileMergeFun"),
                                         removeGasFlaresMethod = pkgOptions(paste0("removeGasFlaresMethod_", nlType)),
                                         nlPeriod,
                                         gadmVersion = pkgOptions("gadmVersion"),
                                         gadmPolyType = pkgOptions("gadmPolyType"),
                                         custPolyPath = NULL)
{
  if (missing(ctryCode))
    stop(Sys.time(), "Missing required parameter ctryCode")
  
  if (missing(nlType))
    stop(Sys.time(), "Missing required parameter nlType")
  
  if (missing(nlPeriod))
    stop(Sys.time(), "Missing required parameter nlPeriod")
  
  if (missing("configName") ||
      is.na(configName) || length(configName) == 0 || configName == "")
    configName <- pkgOptions(paste0("configName_", nlType))
  
  if (missing("extension") ||
      is.na(extension) || length(extension) == 0 || extension == "")
    extension <- pkgOptions(paste0("extension_", nlType))
  
  if (!validCtryCodes(ctryCodes = ctryCode))
    stop(Sys.time(), "Invalid ctryCode: ", ctryCode)
  
  if (!validNlTypes(nlTypes = nlType))
    stop(Sys.time(), "Invalid nlType: ", nlType)
  
  if (!allValidNlPeriods(nlPeriods = nlPeriod, nlTypes = nlType))
    stop(Sys.time(),
         "Invalid nlPeriod: ",
         nlPeriod,
         " for nlType: ",
         nlType)
  
  return (file.path(
    getNlDir("dirRasterOutput"),
    getCtryRasterOutputFname(
      ctryCode = ctryCode,
      nlType = nlType,
      configName = configName,
      extension = extension,
      multiTileStrategy = multiTileStrategy,
      multiTileMergeFun = multiTileMergeFun,
      removeGasFlaresMethod = removeGasFlaresMethod,
      nlPeriod = nlPeriod,
      gadmVersion = gadmVersion,
      gadmPolyType = gadmPolyType,
      custPolyPath = custPolyPath
    )
  ))
}

######################## getCtryMaskOutputFname ###################################

#' Constructs the name of the output raster
#'
#' Constructs the name of the output raster
#'
#' @param ctryCode the ctryCode of interest
#'
#' @param nlType the nlType of interest
#'
#' @param nlPeriod the nlPeriod of interest
#'
#' @param gadmVersion The GADM version to use
#'
#' @param gadmPolyType The format of polygons to download from GADM
#'
#' @param custPolyPath The path to a custom polygon as an alternative
#'     to using GADM polygons
#'
#' @return Character the name of country raster for a country and a given
#'     nlType and nlPeriod
#'
#' @examples
#'
#' Rnightlights:::getCtryGFMaskOutputFname(ctryCode = "KEN", nlType = "VIIRS.M", nlPeriod = "201412")
#'
getCtryGFMaskOutputFname <- function(ctryCode,
                                     nlType,
                                     nlPeriod,
                                     gadmVersion = pkgOptions("gadmVersion"),
                                     gadmPolyType = pkgOptions("gadmPolyType"),
                                     custPolyPath = NULL)
{
  if (missing(ctryCode))
    stop(Sys.time(), "Missing required parameter ctryCode")
  
  if (missing(nlType))
    stop(Sys.time(), "Missing required parameter nlType")
  
  if (missing(nlPeriod))
    stop(Sys.time(), "Missing required parameter nlPeriod")
  
  if (!validCtryCodes(ctryCodes = ctryCode))
    stop(Sys.time(), "Invalid ctryCode: ", ctryCode)
  
  if (!validNlTypes(nlTypes = nlType))
    stop(Sys.time(), "Invalid nlType: ", nlType)
  
  if (!allValidNlPeriods(nlPeriods = nlPeriod, nlTypes = nlType))
    stop(Sys.time(),
         "Invalid nlPeriod: ",
         nlPeriod,
         " for nlType: ",
         nlType)
  
  if (missing(custPolyPath))
    custPolyPath <- NULL
  
  fname <- if (is.null(custPolyPath))
    paste0(
      "NL_GFMask_",
      ctryCode,
      "_",
      nlType,
      "_",
      nlPeriod,
      "_",
      gadmVersion,
      "-",
      toupper(gadmPolyType),
      ".tif"
    )
  else
    paste0(
      "NL_GFMask_",
      ctryCode,
      "_",
      nlType,
      "_",
      nlPeriod,
      "_",
      "_CUST-",
      basename(custPolyPath),
      "-SHPZIP.tif"
    )
  return (fname)
}

######################## getCtryGFMaskOutputFnamePath ###################################

#' Get the full path to the file where the cropped VIIRS country mask is stored.
#'
#' Get the full path to the file where the cropped VIIRS country mask is stored. This file is created
#'     when processing the country before extracting the data. It can be used to re-process a country much faster
#'
#' @param ctryCode (character) the ctryCode of interest
#'
#' @param nlType (character) the nlType of interest
#'
#' @param nlPeriod (character) the nlPeriod of interest
#'
#' @param gadmVersion (character) The GADM version to use
#'
#' @param gadmPolyType (character) The format of polygons to download from GADM
#'
#' @param custPolyPath (character) The path to a custom polygon as an alternative
#'     to using GADM polygons
#'
#' @return (character) full path to the cropped VIIRS country raster for a country and a given year and month
#'
#' @examples
#' \dontrun{
#' getCtryRasterOutputFnamePath("KEN","VIIRS.M", "201412")
#' }
#'
getCtryGFMaskOutputFnamePath <- function(ctryCode,
                                         nlType,
                                         nlPeriod,
                                         gadmVersion = pkgOptions("gadmVersion"),
                                         gadmPolyType = pkgOptions("gadmPolyType"),
                                         custPolyPath = NULL)
{
  if (missing(ctryCode))
    stop(Sys.time(), "Missing required parameter ctryCode")
  
  if (missing(nlType))
    stop(Sys.time(), "Missing required parameter nlType")
  
  if (missing(nlPeriod))
    stop(Sys.time(), "Missing required parameter nlPeriod")

  if (!validCtryCodes(ctryCodes = ctryCode))
    stop(Sys.time(), "Invalid ctryCode: ", ctryCode)
  
  if (!validNlTypes(nlTypes = nlType))
    stop(Sys.time(), "Invalid nlType: ", nlType)
  
  if (!allValidNlPeriods(nlPeriods = nlPeriod, nlTypes = nlType))
    stop(Sys.time(),
         "Invalid nlPeriod: ",
         nlPeriod,
         " for nlType: ",
         nlType)
  
  return (file.path(
    getNlDir("dirRasterOutput"),
    getCtryGFMaskOutputFname(
      ctryCode = ctryCode,
      nlType = nlType,
      nlPeriod = nlPeriod,
      gadmVersion = gadmVersion,
      gadmPolyType = gadmPolyType,
      custPolyPath = custPolyPath
    )
  ))
}


######################## processNlData ###################################

#' Downloads nightlight tiles and country polygons and calls the function to process them
#'
#' Downloads nightlight tiles and country polygons in preparation for the appropriate functions
#'     to process them. Given a list of countries and nlPeriods and an nlType, processNlData
#'     will first determine which countries and periods do not actually have data i.e. have not
#'     been previously processed. From the list of countries and prediods that need processing,
#'     it determines which nightlight tiles require to be downloaded. At the same time, it
#'     downloads any country polygons which have not already been downloaded.
#'
#'     processNlData then calls processNlCountry with the nlType supplied as a parameter. The
#'     processing is essentially the same for all nlTypes.
#'
#'     This is the main entry-point to the package and the main user-facing function. However, it
#'     works in batch mode and caches all data without returning any data to the user. However, it
#'     will return TRUE/FALSE depending on whether it completed successfully. Since it saves state
#'     regularly it can be run multiply in case of failure until it finally returns TRUE. This is where
#'     the only constraints are downloading and processing errors due to bandwidth or other
#'     resource constraints. A good example is running a long-running batch on an AWS EC2
#'     spot-priced machine where the server may be decommissioned at any time. See more in the examples.
#'
#' @param ctryCodes (character) the list of countries to be processed
#'
#' @param admLevels (character) the list of admin levels in the given countries at
#'     which to calculate stats
#'
#' @param nlTypes (character) the types of nightlights to process
#'
#' @param configNames character the config shortnames of rasters being processed
#' 
#' @param extensions character the extensions of rasters being processed
#'
#' @param multiTileStrategy (character) How to handle multiple tiles per nlPeriod
#'
#' @param multiTileMergeFun (character) The function to use to merge tiles
#'
#' @param removeGasFlaresMethod character The method to use to perform gas flare removal
#'     or NULL to disable
#'     
#' @param nlPeriods (character) the nlPeriods of interest
#'
#' @param nlStats (character) the statistics to calculate. If not provided will calculate
#'     the stats specified in \code{pkgOptions("nlStats")}
#'
#' @param gadmVersion (character) The GADM version to use
#'
#' @param gadmPolyType (character) The format of polygons to download from GADM
#'
#' @param custPolyPath (character) Alternative to GADM. A path to a custom shapefile zip
#'
#' @param downloadMethod (character) The method used to download rasters and polygons
#'
#' @param cropMaskMethod (character) The method used to crop and mask satellite rasters
#'
#' @param extractMethod (character) The method used to extract and perform functions on raster data
#'
#' @return None
#'
#' @examples
#'
#' #long running examples which may require large downloads
#' \dontrun{
#' #Example 1: process monthly VIIRS nightlights for all countries at the
#'     lowest admin level and for all nlPeriods available e.g. to create a
#'     local cache or repo
#'
#'     processNlData() #process VIIRS nightlights for all countries and all periods
#'
#' #Example 2: process monthly VIIRS nightlights for all countries in 2014 only
#'
#'     nlPeriods <- getAllNlPeriods("VIIRS.M") #get a list of all nightlight periods to present-day
#'
#'     nlPeriods <- nlPeriods[grep("^2014", nlPeriods)] #filter only periods in 2014
#'
#'     processNlData(nlTypes="VIIRS.M", nlPeriods=nlPeriods)
#'
#' #Example 3: process OLS nightlights for countries KEN & RWA from 1992
#' #     to 2000
#'
#'     cCodes <- c("KEN", "RWA")
#'
#'     nlPeriods <- getAllNlPeriods("VIIRS.M")
#'
#'     nlPeriods <- nlRange("1992", "2000", "OLS.Y")
#'
#'     processNlData(ctryCodes=cCodes, nlPeriods=nlPeriods)
#'
#' #Example 4: process VIIRS nightlights for countries KEN & RWA in 2014 Oct to 2014 Dec only
#'
#'     processNlData(ctryCodes=c("KEN", "RWA"), nlTypes="VIIRS.M",
#'         nlPeriods=c("201410", "201411", "201412"))
#'
#' #Example 5: process all nightlights, all countries, all stats in one thread
#'
#'    processNlData()
#'
#' #Example 6: process all VIIRS monthly nightlights, all countries, all stats with each
#' #   year in a separate thread. Create a separate R script for each year as follows:
#'
#'     library(Rnightlights)
#'
#'     nlPeriods <- getAllNlPeriods("VIIRS.M")
#'
#'     nlPeriods_2014 <- nlPeriods[grep("^2014", nlPeriods)]
#'
#'     processNlData(nlPeriods=nlPeriods_2014)
#'
#'     #Run the script from the command line as:
#'
#'     #R CMD BATCH script_name_2014.R
#'     }
#' @export
processNlData <- function (ctryCodes,
                           admLevels,
                           nlTypes,
                           configNames,
                           extensions,
                           multiTileStrategy = pkgOptions("multiTileStrategy"),
                           multiTileMergeFun = pkgOptions("multiTileMergeFun"),
                           removeGasFlaresMethod = pkgOptions(paste0("removeGasFlaresMethod_", nlType)),
                           nlPeriods,
                           nlStats = pkgOptions("nlStats"),
                           custPolyPath = NULL,
                           gadmVersion = pkgOptions("gadmVersion"),
                           gadmPolyType = pkgOptions("gadmPolyType"),
                           downloadMethod = pkgOptions("downloadMethod"),
                           cropMaskMethod = pkgOptions("cropMaskMethod"),
                           extractMethod = pkgOptions("extractMethod"))
{
  if (missing(ctryCodes))
    ctryCodes <- getAllNlCtryCodes(omit = "error")
  
  if (missing(admLevels))
    admLevels <- "adm0"
  
  if (missing(nlTypes))
    nlTypes <- getAllNlTypes()
  
  if (missing(configNames))
    configNames <- pkgOptions(paste0("configName_", nlTypes))
  
  if(missing(extensions))
    extensions <- pkgOptions(paste0("extension_", nlTypes))
  
  #if the period is not given process all available periods
  if (missing("nlPeriods") ||
      is.na(nlPeriods) || length(nlPeriods) == 0 || nlPeriods == "")
  {
    nlPeriods <- getAllNlPeriods(nlTypes = nlTypes)
  }
  
  if (!allValid(testData = ctryCodes, testFun = validCtryCodes))
    stop(Sys.time(), ": Invalid ctryCode detected")
  
  #try to understand admLevels provided  and put into a format
  #that we can test later i.e. each ctryCode with corresponding admLevel entries
  if (length(admLevels) != length(ctryCodes))
  {
    #if admLevels are named we expect each country to have an entry
    #for admLevels. Simple check that lengths must be the same
    if (!is.null(names(admLevels)))
      stop(Sys.time(), ": AdmLevel names do not match ctryCodes")
    
    #otherwise assume the admLevels refer to each country
    #combine the admLevels and repeat for each ctryCode
    #any wrong matches should be caught later
    admLevels <-
      stats::setNames(rep(list(unlist(admLevels)), length(ctryCodes)), ctryCodes)
  }
  
  funArgs <- match.call()[-1]
  
  message(
    Sys.time(),
    ": **** START PROCESSING: ",
    paste(
      names(funArgs),
      sapply(names(funArgs), function(x)
        eval(parse(text = x))),
      sep = "=",
      collapse = ", "
    ),
    "****"
  )
  
  #uppercase vars
  ctryCodes <- toupper(x = ctryCodes)
  
  #Ensure we have all polygons before checking admLevels
  message(Sys.time(), ": Downloading country polygons ...")
  
  #cl <- snow::makeCluster(pkgOptions("numThreads"))
  
  #doSNOW::registerDoSNOW(cl = cl)
  
  pb <- utils::txtProgressBar(min = 0,
                              max = length(ctryCodes),
                              style = 3)
  
  progress <- function(n)
    utils::setTxtProgressBar(pb, n)
  
  #download all country polygons if they don't already exist
  #foreach::foreach (ctryCode = ctryCodes,
  #                  #.export = c("gadmVersion", "gadmPolyType", "custPolyPath", "downloadMethod"),
  #                  .options.snow = list(progress=progress)) %dopar%
  for (ctryCode in ctryCodes)
  {
    message(Sys.time(), ": Downloading polygon: ", ctryCode)
    
    dnldCtryPoly(
      ctryCode = ctryCode,
      gadmVersion = gadmVersion,
      gadmPolyType = gadmPolyType,
      custPolyPath = custPolyPath,
      downloadMethod = downloadMethod
    )
    
    progress(n = 1)
  }
  
  close(pb)
  
  #snow::stopCluster(cl)
  
  message(Sys.time(), ": Downloading country polygons ... DONE")
  
  admLevels <-
    sapply(seq_along(unlist(ctryCodes)), function(cCodeIdx)
    {
      ctryCode <- ctryCodes[[cCodeIdx]]
      
      #ctryAdmLevels <- admLevels[[cCodeIdx]]
      #if(!is.null(names(admLevels)[cCodeIdx]))
      #  return(admLevels[ctryCode])
      
      ctryAdmLevels <- if (is.list(admLevels))
        unlist(admLevels[cCodeIdx])
      else
        admLevels[cCodeIdx]
      
      unlist(sapply(ctryAdmLevels, function(admLevel)
      {
        #allow keywords/aliases instead of admLevels
        if (length(admLevel) == 1)
        {
          if (admLevel == "country")
            admLevel <- getCtryShpLyrNames(
              ctryCodes = ctryCode,
              lyrNums = 0,
              gadmVersion = gadmVersion,
              gadmPolyType = gadmPolyType,
              custPolyPath = custPolyPath
            )
          else if (admLevel %in% c("bottom", "lowest"))
            admLevel <- getCtryShpLowestLyrNames(
              ctryCodes = ctryCode,
              gadmVersion = gadmVersion,
              gadmPolyType = gadmPolyType,
              custPolyPath = custPolyPath
            )
          else if (admLevel %in% c("top", "highest"))
            admLevel <- getCtryShpLyrNames(
              ctryCodes = ctryCode,
              lyrNums = 1,
              gadmVersion = gadmVersion,
              gadmPolyType = gadmPolyType,
              custPolyPath = custPolyPath
            )
          else if (admLevel == "all")
            admLevel <- getCtryShpAllAdmLvls(
              ctryCodes = ctryCode,
              gadmVersion = gadmVersion,
              gadmPolyType = gadmPolyType,
              custPolyPath = custPolyPath
            )
          else
          {
            tmpAdmLevel <- searchAdmLevel(
              ctryCodes = ctryCode,
              admLevelNames = admLevel,
              downloadMethod = downloadMethod,
              gadmVersion = gadmVersion,
              gadmPolyType = gadmPolyType,
              custPolyPath = custPolyPath
            )
            
            admLevel <-
              ifelse(is.na(tmpAdmLevel), admLevel, tmpAdmLevel)
          }
          
          # #after processing admLevels if any are not in proper format e.g. KEN_adm0
          # #check if they might have been supplied as e.g. adm0 or e.g. 0
          # if(!length(grep(paste0(ctryCode,"_adm\\d+"), admLevel, ignore.case = T)) == length(admLevel))
          # {
          #   if(length(grep("^adm\\d+$", admLevel, ignore.case = T)) > 0)
          #     admLevel <- paste(ctryCode, admLevel, sep="_")
          #   else if(length(grep("^\\d+$", admLevels, ignore.case = T)) > 0)
          #     admLevel <- paste(ctryCode, paste0("adm", admLevel), sep="_")
          # }
        }
        
        admLevel
      }))
    })
  
  #if the admLevels are input as "adm0" convert to proper admLevel e.g. "KEN_adm0"
  # for(i in 1:length(ctryCodes))
  # if(!length(grep(paste0(ctryCodes[i],"_adm\\d+"), admLevels[i], ignore.case = T)) == length(admLevels[i]))
  # {
  #   if(length(grep("^adm\\d+$", admLevels[i], ignore.case = T)) == length(admLevels[i]))
  #     admLevels[i] <- paste(ctryCodes[i], admLevels[i], sep="_")
  #   else if(length(grep("^\\d+$", admLevels[i], ignore.case = T)) == length(admLevels[i]))
  #     admLevels[i] <- paste(ctryCode, paste0("adm", admLevels[i]), sep="_")
  # }
  
  # if(!allValidCtryAdmLvls(ctryCode = ctryCodes, admLevels = admLevels, gadmVersion = gadmVersion, custPolyPath = custPolyPath))
  #   stop(Sys.time(), "Invalid admLevels detected")
  
  #traverse and check all admLevel component lists/vectors
  if (!all(sapply(seq_along(ctryCodes),
                  function(i)
                  {
                    if (is.list(admLevels))
                    {
                      allValidCtryAdmLvls(
                        ctryCode = ctryCodes[i],
                        admLevels = admLevels[[i]],
                        gadmVersion = gadmVersion,
                        gadmPolyType = gadmPolyType,
                        custPolyPath = custPolyPath
                      )
                    } else
                    {
                      if (length(ctryCodes) == 1)
                        allValidCtryAdmLvls(
                          ctryCode = ctryCodes[i],
                          admLevels = admLevels,
                          gadmVersion = gadmVersion,
                          gadmPolyType = gadmPolyType,
                          custPolyPath = custPolyPath
                        )
                      else
                        allValidCtryAdmLvls(
                          ctryCode = ctryCodes[i],
                          admLevels = admLevels[i],
                          gadmVersion = gadmVersion,
                          gadmPolyType = gadmPolyType,
                          custPolyPath = custPolyPath
                        )
                    }
                  })))
  
  stop(Sys.time(), "Invalid admLevels detected")
  
  if (!allValidNlPeriods(nlTypes = nlTypes, nlPeriods = nlPeriods))
  {
    stop(Sys.time(), "Invalid nlPeriod(s) detected")
  }
  
  nlTiles <- NULL
  
  ##First step: Determine which tiles are required for processing. This is determined by the
  #list of ctryCodes. Since theoretically we need the polygons of all countries to determine which
  #tiles to download we take the opportunity to download all required shapefiles
  #Practically, we can stop as soon as all nlTiles for an nlPeriod are flagged for download.
  
  ##If any tiles cannot be found/downloaded then abort and try for next country
  #we probably need to flag failed downloads so we don't try to process them and report back to user
  
  #extract nlPeriods ensuring they match with nlPeriods
  for (idxNlType in seq_along(nlTypes))
  {
    nlType <- nlTypes[idxNlType]
    
    configName <- configNames[idxNlType]
    
    extension <- extensions[idxNlType]
      
    if (length(nlTypes) == 1)
    {
      if (is.list(nlPeriods))
      {
        nlTypePeriods <- unlist(nlPeriods)
      }
      else
      {
        nlTypePeriods <- nlPeriods
      }
    }
    else
    {
      if (is.list(nlPeriods))
      {
        nlTypePeriods <- unlist(nlPeriods[idxNlType])
      }
      else
      {
        nlTypePeriods <- nlPeriods[[idxNlType]]
      }
    }
    
    #if the tile mapping does not exist create it
    if (!exists("nlTiles") || is.null(nlTiles))
      nlTiles <- getNlTiles(nlType = nlType)
    
    #for all nlPeriods check if the tiles exist else download
    for (nlPeriod in nlTypePeriods)
    {
      message(
        Sys.time(),
        ": **** PROCESSING nlType:",
        nlType,
        " | configName: ",
        configName,
        " | extension: ",
        extension,
        " | nlPeriod:",
        nlPeriod,
        "****"
      )
      message(Sys.time(),
              ": Checking tiles required for ",
              paste(nlType, nlPeriod))
      
      #init the list of tiles to be downloaded
      tileList <- NULL
      
      #determine tiles to download
      ctryTiles <- NULL
      
      tileList <- NULL
      
      #For each country
      for (idxCtryCode in seq_along(ctryCodes))
      {
        ctryCode <- ctryCodes[idxCtryCode]
        
        if (is.list(admLevels))
          ctryAdmLevels <- unlist(admLevels[idxCtryCode])
        else
          ctryAdmLevels <- admLevels[idxCtryCode]
        
        existAdmLvlStats <- unlist(lapply(ctryAdmLevels,
                                          function(admLevel)
                                            allExistsCtryNlData(
                                              ctryCodes = ctryCode,
                                              admLevels = admLevel,
                                              nlTypes = nlType,
                                              configNames = configNames,
                                              extensions = extensions,
                                              multiTileStrategy = multiTileStrategy,
                                              multiTileMergeFun = multiTileMergeFun,
                                              removeGasFlaresMethod = removeGasFlaresMethod,
                                              nlPeriods = nlPeriod,
                                              nlStats = nlStats,
                                              gadmVersion = gadmVersion,
                                              gadmPolyType = gadmPolyType,
                                              custPolyPath = custPolyPath
                                            )))
        
        
        #Check if all stats exist for the ctryCode
        if (all(existAdmLvlStats))
        {
          message (Sys.time(), ": **", ctryCode, ": All stats exist")
          
          next
        }
        
        rasterOutputFnamePath <-
          getCtryRasterOutputFnamePath(
            ctryCode = ctryCode,
            nlType = nlType,
            configName = configName,
            extension = extension,
            multiTileStrategy = multiTileStrategy,
            multiTileMergeFun = multiTileMergeFun,
            removeGasFlaresMethod = removeGasFlaresMethod,
            nlPeriod = nlPeriod,
            gadmVersion = gadmVersion,
            gadmPolyType = gadmPolyType,
            custPolyPath = custPolyPath
          )
        
        #if the cropped raster is found do not list tile
        if (file.exists(rasterOutputFnamePath))
        {
          message(Sys.time(),
                  ": ",
                  ctryCode,
                  ": cropped raster exists. Tile not required")
          
          next
        }
        
        message(Sys.time(),
                ": ",
                ctryCode,
                ": Stats missing. Adding tiles")
        
        #get the list of tiles required for the ctryCode
        ctryTiles <-
          getCtryTileList(ctryCodes = ctryCode, nlType = nlType)
        
        #combine the list of unique tiles across all ctryCodes in tileList
        tileList <-
          c(tileList, setdiff(x = ctryTiles, y = tileList))
        
        #if all the unique tiles have been listed no need to proceed checking
        if (length(tileList) == nrow(nlTiles))
        {
          message (Sys.time(),
                   ": All tiles have been listed. No need to check other country tiles")
          
          break
        }
      }
      
      message(
        Sys.time(),
        ": numTiles: ",
        length(tileList),
        ", Required tiles: ",
        paste(tileList, collapse = ",")
      )
      
      if (length(tileList) == 0)
      {
        message(
          Sys.time(),
          ": No tiles needed for ",
          nlPeriod,
          ". Processing countries with missing stats and existing cropped rasters"
        )
        
        #next
      }
      else
      {
        #download all required tiles
        if (!downloadNlTiles(
          nlType = nlType,
          configName = configName,
          extension = extension,
          multiTileStrategy = multiTileStrategy,
          nlPeriod = nlPeriod,
          tileList = tileList
        ))
        {
          message(Sys.time(),
                  ": Something went wrong with the tile downloads. Aborting ...")
          
          break
        }
      }
      
      #processNlCountry for all countries
      #Run through all countries since not guaranteed to check all countries in the
      #loop above checking existing stats. Premature exit allowed once all tiles listed
      for (idxCtryCode in seq_along(ctryCodes))
      {
        message(Sys.time(), ": Processing ctryCode: ", ctryCodes[idxCtryCode])
        
        ctryCode <- ctryCodes[idxCtryCode]
        
        if (is.list(admLevels))
          ctryAdmLevels <- unlist(admLevels[idxCtryCode])
        else
          ctryAdmLevels <- admLevels[idxCtryCode]
        
        for (admLevel in ctryAdmLevels)
        {
          message(Sys.time(),
                  ": Processing ctryCode: ",
                  ctryCode,
                  " AdmLevel: ",
                  admLevel)
          
          processNLCountry(
            ctryCode = ctryCode,
            admLevel = admLevel,
            nlType = nlType,
            configName = configName,
            extension = extension,
            multiTileStrategy = multiTileStrategy,
            multiTileMergeFun = multiTileMergeFun,
            removeGasFlaresMethod = removeGasFlaresMethod,
            nlPeriod = nlPeriod,
            nlStats = nlStats,
            downloadMethod = downloadMethod,
            cropMaskMethod = cropMaskMethod,
            extractMethod = extractMethod,
            gadmVersion = gadmVersion,
            gadmPolyType = gadmPolyType,
            custPolyPath = custPolyPath
          )
        }
      }
      
      #post-processing. Delete the downloaded tiles to release disk space
      if (pkgOptions("deleteTiles"))
      {
        for (tile in tileList)
        {
          tileTif <- getNlTileTifLclNamePath(
            nlType = nlType,
            nlPeriod = nlPeriod,
            tileNum = tileName2Idx(tile, nlType)
          )
          
          message(Sys.time(), ": Deleting tile TIF: ", tileTif)
          
          #del the tif file
          if (file.exists(tileTif))
            unlink(file.path(tileTif), force = TRUE)
          
          tileZip <-
            file.path(
              getNlTileZipLclNamePath(
                nlType = nlType,
                nlPeriod = nlPeriod,
                tileNum = tileName2Idx(tileName = tile,
                                       nlType = nlType)
              )
            )
          
          message(Sys.time(), ": Deleting tile ZIP: ", tileZip)
          
          #del the zip file
          if (file.exists(tileZip))
            unlink(x = file.path(tileZip), force = TRUE)
        }
      }
    }
  }
  
  message(
    Sys.time(),
    ": **** COMPLETED PROCESSING :",
    paste(
      names(funArgs),
      sapply(names(funArgs), function(x)
        eval(parse(text = x))),
      sep = "=",
      collapse = ", "
    ),
    "****"
  )
}