In RETURN-project/makeDataCube: Prepare raster datasets for upscaling recovery indicators

knitr::opts_chunk$set(echo = TRUE, 
                      eval = FALSE)
library(makeDataCube)
library(zoo)
library(lubridate)
library(raster)
library(tidyverse)
library(rgdal)
library(terra)

General

The goal of this vignette is to generate a mask file and add it to an existing data cube with optical data. This cube should be generated with the FORCE software and follow the folder structure as defined in make_Landsat_cube.Rmd. If no data cube is available, please generate one first.

The mask equals 1 for pixels that meet the following requirements:

• Have a tree cover percentage higher than Ttree in the year 2000
• Being burned with a confidence level higher than Tconf
• Not being classified as a non-natural land cover type over the study period
• Being classified as forest in the user-defined year Tyr

Other pixels are set to 0.

Inputs

Rtmpfolder <- normalizePath(file.path('..', 'data'))# folder where all temporary R files should be stored Rtmpfolder <- '/project/return/Data/tempR'
forcefolder <- normalizePath(file.path('..', 'data'))#'/home/wanda/Documents/data/force_small'#file.path('..', 'data')#forcefolder <- '/home/wanda/Documents/data/force' forcefolder <- '/project/return/Data/force'
ext <- c(-43.38238361637443,-43.27938679020256,-4.555765244985907,-4.451717415449725)#extent, vector with xmin, xmax, ymin, ymax in degrees
startyear <- 1998 #start year
endyear <- 2018 #end year

# settings for the mask
Ttree <- 60# threshold for the tree cover percentage (in year 2000)
Tconf <- 95 # threshold on confidence level to detect a fire
Tyr <- 2000 # year for which land cover should be forest; typically beginning of disturbance monitoring period
tempRes <- 'monthly'# desired temporal resolution of the fire data
resMask <- 30 # desired resolution of the raster in the cube (nearest neighbor resampling)
resFire <- 30 # desired resolution of the raster in the cube (nearest neighbor resampling)

starttime <- c(startyear,01,01)#start date: year, month, day
endtime <- c(endyear,12,31)#end date: year, month, day

Generate the directory structure

A standard structure of the data folders is assumed.

# set the folder to store temporary data
unixtools::set.tempdir(Rtmpfolder)
# set the upper memory limit for raster data. Raster will move data to disk if it exceeds the upper limit (defaults to 100MB). 
# options(rasterMaxMemory = 1e10)
# generate folder structure
fldrs <- setFolders(forcefolder)
tmpfolder <- fldrs['tmpfolder']  # temporary files
l2folder <- fldrs['l2folder']# level 2 data
queuefolder <- fldrs['queuefolder']# queue file
queuefile <- fldrs['queuefile']# name queue file
lcfolder <- fldrs['lcfolder']# raw land cover data
lclogfile <- fldrs['lclogfile']# land cover log file
tcfolder <- fldrs['tcfolder']# raw tree cover data
tclogfile <- fldrs['tclogfile']
firefolder <- fldrs['firefolder']# raw fire data
firelogfile <- fldrs['firelogfile']
paramfolder <- fldrs['paramfolder'] # parameter file
paramfile <- fldrs['paramfile']# name parameter file

generate a shapefile of the grid

iterate over the tiles and prepare mask and fire data

# for which continents do we have data available?
continents <- file.path(fldrs[['l2folder']],c('n-america', 's-america', 'europe', 'asia', 'africa', 'oceania', 'antartica'))
continents <- continents[continents %in% list.dirs(fldrs[['l2folder']], recursive = F)]

for(gi in 1:length(continents)){
  # generate a list of tiles to process
  elist <- getGrid(continents[gi], ext)
  if (length(elist)>0){
    for(ti in 1:length(elist)){
      exti <- elist[[ti]]# extent of the tile of interest
      tilename <- names(elist)[ti]# name of the tile of interest
      # folder with data of interest
      tilefolder <- file.path(continents[gi], tilename)
      # check if the tile folder of interest exists
      if(dir.exists(tilefolder)){
        # generate folder for temporary data of study area of interest
        extfolder <- file.path(forcefolder, 'temp', paste0('Area', exti[1],'_', exti[2],'_', exti[3],'_', exti[4]))
        if(!dir.exists(extfolder)){dir.create(extfolder)}

        # ----------------------
        ## Prepare land cover data
        # First MapBiomas land cover data are downloaded if they are not available yet. Then the land cover data are cropped to the extent of the study area, and the study period of interest. Currently, MapBiomas data are only available for the years 1985 - 2019.

        # Download MapBiomas land cover data
        dllLandcover(lcfolder, lclogfile)
        # Crop the land cover data to extent and time period of interest 
        lc_rst <- list()
        lcregions <- c('AMAZONIA', 'PANTANAL', 'CAATINGA', 'MATAATLANTICA', 'CERRADO', 'PAMPA')
        for(i in 1:length(lcregions)){
          lcfiles <- paste0(lcregions[i], '-', 1985:2019, '.tif')
          lc_rst[[i]] <- terra::rast(file.path(ifolder, lcfiles)) # open land cover data file
        }
        prepLandcover(lc_rst, extfolder, exti, fname = 'landcover.tif', as.Date(paste0(starttime[1],'-',starttime[2],'-',starttime[3])), as.Date(paste0(endtime[1],'-',endtime[2],'-',endtime[3])))

        # ----------------------
        ## Prepare tree cover data
        # Download tree cover data
        tcfls <- dllTreecover(tcfolder, exti, tclogfile)
        # Crop the tree cover data to exent of interest and mask nodata values
        hanfiles <- list()
        hanmaskfiles <- list()
        for(i in 1:length(tcfls$hanfiles)){
          hanfiles[[i]] <- terra::rast(file.path(ifolder, tcfls$hanfiles[i]))
          hanmaskfiles[[i]] <- terra::rast(file.path(ifolder, tcfls$hanmaskfiles[i]))
        }
        prepTreecover(tcfolder, extfolder, exti, fname = 'treecover.tif', hanfiles, hanmaskfiles)

        # ----------------------
        ## Generate mask without fire information
        han <- terra::rast(file.path(extfolder,'treecover.tif'))
        lc <- terra::rast(file.path(extfolder,'landcover.tif'))
        lcDates <- as.Date(loadRData(file.path(extfolder,'lcDates')))

        msk <- makeMaskNoFire(lc, lcDates, han, extfolder)
        rm(lc)

        # ----------------------
        ## Prepare fire data
        # Download fire data
        firefls <- dllFire(firefolder, firelogfile)
        # Crop the fire data to exent and time period of interest, resample and generate a binary time series (1 = fire, 0 = no fire) at the temporal resolution of interest
        fcl <- terra:rast(file.path(firefolder,firefls$fireclfiles))
        fjd <- terra::rast(file.path(firefolder,firefls$firejdfiles))
        fdts <- as.Date(firefls$fireclfiles, format = "%Y%m%d-ESACCI-L3S_FIRE-BA-MODIS-AREA_2-fv5.1-CL.tif")
        tsFire2 <- prepFire(fcl, fjd, fdts, han, msk, tempRes, Tconf, starttime, endtime, extfolder)
        # Export the fire data to the temporary folder
        nms <-str_replace_all(names(tsFire2),'-','_')
        terra::writeRaster(tsFire2, file.path(extfolder, paste0('fire_',nms,'.tif')), bylayer=TRUE, overwrite=TRUE)

        # ----------------------
        ## update mask with locations that experienced fire
        msk <- terra::app(tsFire2, sum, na.rm=T, filename = file.path(extfolder, 'fresum.tif'), overwrite=T)
        msk <- (msk > 0) # areas that have experienced a fire get the value 1, other areas the value 0
        # Export the mask layer tot the temporary folder
        terra::writeRaster(msk, file.path(extfolder, 'mask.tif'), overwrite=TRUE)

        # ----------------------
        ## Add mask and fire data to data cube
        system(paste0("force-cube ",file.path(extfolder, 'mask.tif')," ", file.path(continents[gi])," near ", resMask), intern = TRUE, ignore.stderr = TRUE)

        for(i in 1:dim(tsFire2)[3]){
          system(paste0("force-cube ",file.path(extfolder,  paste0('fire_',nms[i],'.tif'))," ", file.path(continents[gi])," near ", resFire), intern = TRUE, ignore.stderr = TRUE)
        }

        # remove temporary folder and its contents
        unlink(extfolder, recursive=TRUE)
      }
    }
  }
}

RETURN-project/makeDataCube documentation built on Feb. 11, 2022, 3:04 p.m.