stack_q.R

library(dplyr)     
library(rgee)
library(sf) 

library(googledrive)
library(googleCloudStorageR)
library(lwgeom)

library(reticulate)
library(jsonlite)
library(tidyverse)

# use for all other accounts
ee_Initialize()


df <- structure(list(FIRENAME = c("Ash", "Bitumul"), Year = c(1985L, 1985L), 
                     Fire_ID = c("1985-AZCNF-000071", "1985-AZASF-000148"), 
                     geometry = structure(list(structure(list(structure(c(-212345.986249991, -212789.545625001, -213239.380625002, -212691.334124997, -212345.986249991, 3816853.32925, 3816529.052375, 3816894.243125, 3817495.70187501, 3816853.32925), .Dim = c(5L, 2L))), class = c("XY", "POLYGON", "sfg")), 
                                               structure(list(structure(c(-245447.64837499, -245857.643499993, -245697.833750002, -245154.533999994, -245447.64837499, 4037099.307625, 4037754.51575, 4038240.66825, 4037862.15787501, 4037099.307625), .Dim = c(5L, 2L))), class = c("XY", "POLYGON", "sfg"))), 
                                          n_empty = 0L, crs = structure(list(input = "North_America_Lambert_Conformal_Conic", wkt = "PROJCRS[\"North_America_Lambert_Conformal_Conic\",\n    BASEGEOGCRS[\"NAD83\",\n        DATUM[\"North American Datum 1983\",\n            ELLIPSOID[\"GRS 1980\",6378137,298.257222101,\n                LENGTHUNIT[\"metre\",1]],\n            ID[\"EPSG\",6269]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"Degree\",0.0174532925199433]]],\n    CONVERSION[\"unnamed\",\n        METHOD[\"Lambert Conic Conformal (2SP)\",\n            ID[\"EPSG\",9802]],\n        PARAMETER[\"Latitude of false origin\",0,\n            ANGLEUNIT[\"Degree\",0.0174532925199433],\n            ID[\"EPSG\",8821]],\n        PARAMETER[\"Longitude of false origin\",-108,\n            ANGLEUNIT[\"Degree\",0.0174532925199433],\n            ID[\"EPSG\",8822]],\n        PARAMETER[\"Latitude of 1st standard parallel\",32,\n            ANGLEUNIT[\"Degree\",0.0174532925199433],\n            ID[\"EPSG\",8823]],\n        PARAMETER[\"Latitude of 2nd standard parallel\",36,\n            ANGLEUNIT[\"Degree\",0.0174532925199433],\n            ID[\"EPSG\",8824]],\n        PARAMETER[\"Easting at false origin\",0,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8826]],\n        PARAMETER[\"Northing at false origin\",0,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8827]]],\n    CS[Cartesian,2],\n        AXIS[\"(E)\",east,\n            ORDER[1],\n            LENGTHUNIT[\"metre\",1,\n                ID[\"EPSG\",9001]]],\n        AXIS[\"(N)\",north,\n            ORDER[2],\n            LENGTHUNIT[\"metre\",1,\n                ID[\"EPSG\",9001]]]]"), class = "crs"), class = c("sfc_POLYGON", "sfc"), precision = 0, bbox = structure(c(xmin = -245857.643499993, ymin = 3816529.052375, xmax = -212345.986249991, ymax = 4038240.66825), class = "bbox"))), row.names = c(NA, -2L), class = c("sf", "data.frame"), sf_column = "geometry", agr = structure(c(FIRENAME = NA_integer_, Year = NA_integer_, Fire_ID = NA_integer_), class = "factor", .Label = c("constant", "aggregate", "identity")))





## Script to develop and export spectral index metrics (see list below) for fires in North America. 
##    Script by Morgan Voss and Than Robinson - University of Montana
##    with revisions by Lisa Holsinger - U.S. Forest Service Rocky Mountain Research Station
##    and Ellen Whitman - Canadian Forest Service, Northern Forestry Centre
##    Oct 2019
##    For updates or questions on this code, please contact:  
##        Ellen Whitman, ellen.whitman@canada.ca
##        Lisa Holsinger, lisa.holsinger@usda.gov
##        Sean Parks, sean.parks@usda.gov



## This script backfills 'No Data' areas, which occur due to clouds, cloud-shadows, snow with imagery from
##    up to five years prior to fire for 'pre' fire imagery, and up to two years after fire for 'post' fire
##    imagery.  Note, data are mosaiced such that if 'clear' imagery exists for one year prior and one year after fire,
##    only that imagery is used, and additional imagery is used to fill 'no data' areas, in a successive manner.

##--------------------       INPUTS       ------------------------------##
## import shapefile with fire polygons as a feature collection - these must have the following attributes:
##    Fire_ID, Year
fires <-  ee$FeatureCollection(sf_as_ee(df)); ##Add your google earth engine account name between the two slashes. 
##Upload and name appropriately your 'Fires' shapefile as an asset in your directory  


## specify fire severity metrics to create
bandList <-  ('rdnbr_w_offset'); 

## Enter beginning and end days for imagery season as julian dates, adapt for your local fire season
startday <-  91;
endday   <-  273;

##  visualize fire perimeters
Map$centerObject(fires);
Map$addLayer(fires);

##--------------------    END OF INPUTS   ----------------------------##


##--------------------     PROCESSING     ----------------------------##
##-------- Initialize variables for fire perimeters  -----------------##
## create list with fire IDs  
fireID    <- ee$List(fires$aggregate_array('Fire_ID'))$getInfo();
fireName    <- ee$List(fires$aggregate_array('FIRENAME'))$getInfo(); 

nFires <- length(fireID);
nNames <- length(fireName);

##------------------- Image Processing for Fires Begins Here -------------##
## Landsat 5, 7, and 8 TOA Reflectance Tier 1 collections
##     Only include SLC On Landat 7
ls8SR <- ee$ImageCollection('LANDSAT/LC08/C01/T1_SR')
ls7SR <- ee$ImageCollection('LANDSAT/LE07/C01/T1_SR')
ls5SR <- ee$ImageCollection('LANDSAT/LT05/C01/T1_SR')
ls4SR <- ee$ImageCollection('LANDSAT/LT04/C01/T1_SR')

##///////////////////////////////
## FUNCTIONS TO CREATE NBR
##///////////////////////////////


## Returns vegetation indices for LS8
ls8_Indices <- function(lsImage){
  nbr <- lsImage$normalizedDifference(c("B5", "B7"))$toFloat()
  qa <- lsImage$select("pixel_qa")
  lsImage$select(list(0,1), list("nbr", "pixel_qa"))$
    copyProperties(lsImage, list('system:time_start'))
}

## Returns vegetation indices for LS4, LS5 and LS7
ls4_7_Indices <- function(lsImage){
  nbr <- lsImage$normalizedDifference(c("B4", "B7"))$toFloat()
  qa <- lsImage$select("pixel_qa")
  nbr$addBands(qa)
  lsImage$select(list(0,1), list("nbr", 'pixel_qa'))$
    copyProperties(lsImage, list('system:time_start'))
}

## Mask Landsat surface reflectance images
## Creates a mask for clear pixels 
lsCfmask <- function(lsImg){
  quality <-lsImg$select('pixel_qa')
  clear <- quality$bitwiseAnd(8)$eq(0)$ ## cloud shadow
    And(quality$bitwiseAnd(32)$eq(0))$ ## cloud
    And(quality$bitwiseAnd(4)$eq(0))$ ## water
    And(quality$bitwiseAnd(16)$eq(0)) ## snow
  lsImg$updateMask(clear)$
    select(0)$                                  
    copyProperties(lsImg, list('system:time_start'))
}

## Map functions across Landsat Collections

ls8 <- ls8SR$map(ls8_Indices)$map(lsCfmask)
ls7 <- ls7SR$map(ls4_7_Indices)$map(lsCfmask)
ls5 <- ls5SR$map(ls4_7_Indices)$map(lsCfmask)
ls4 <- ls4SR$map(ls4_7_Indices)$map(lsCfmask)

## Merge Landsat Collections
lsCol <- ee$ImageCollection(ls8$merge(ls7)$merge(ls5)$merge(ls4))


                            
                            
## ------------------ Create and Export Fire Severity Imagery for each fire -----------------##
indices <- ee$ImageCollection(fires$map(function(ft){
  ## use 'Fire_ID' as unique identifier
  fName    <-  ft$get("Fire_ID")
  ## select fire
  fire <-  ft
  fireBounds <-  ft$geometry()$bounds()
  
  ## create pre- and post-fire imagery
  fireYear <-  ee$Date$parse('YYYY', ee$Number$format(fire$get('Year')))
  
  ## Pre-Imagery
  preFireYear <-  fireYear$advance(-1, 'year')
  preFireIndices <-  lsCol$filterBounds(fireBounds)$
    filterDate(preFireYear, fireYear)$
    filter(ee$Filter$dayOfYear(startday, endday))$
    mean()$
    rename('preNBR')
  
  ## if any pixels within fire have no data due to clouds, go back two years (and up to five) to fill in no data areas
  ## two years back
  preFireYear2    <-  fireYear$advance(-2, 'year')
  preFireIndices2 <-  lsCol$filterBounds(fireBounds)$
    filterDate(preFireYear2, fireYear)$
    filter(ee$Filter$dayOfYear(startday, endday))$
    mean()$
    rename('preNBR')
  
  pre_check <-  preFireIndices$clip(fire)
  pre_filled <- pre_check$unmask(preFireIndices2)
                              
  ## three years back
  preFireYear3    <-  fireYear$advance(-3, 'year')
  preFireIndices3 <-  lsCol$filterBounds(fireBounds)$
    filterDate(preFireYear3, fireYear)$
    filter(ee$Filter$dayOfYear(startday, endday))$
    mean()$
    rename('preNBR')
  pre_filled2<- pre_filled$unmask(preFireIndices3)
  
  ## four years back
  preFireYear4 <-  fireYear$advance(-4, 'year')
  preFireIndices4 <-  lsCol$filterBounds(fireBounds)$
    filterDate(preFireYear4, fireYear)$
    filter(ee$Filter$dayOfYear(startday, endday))$
    mean()$
    rename('preNBR')
  pre_filled3 <- pre_filled2$unmask(preFireIndices4)
                              
  ## five years back
  preFireYear5 <-  fireYear$advance(-5, 'year')
  preFireIndices5 <-  lsCol$filterBounds(fireBounds)$
    filterDate(preFireYear5, fireYear)$
    filter(ee$Filter$dayOfYear(startday, endday))$
    mean()$
    rename('preNBR')
  pre_filled4 <- pre_filled3$unmask(preFireIndices5) 
                              
  ## Post-Imagery
  postFireYear <-  fireYear$advance(1, 'year')
  postFireIndices <-  lsCol$filterBounds(fireBounds)$
    filterDate(postFireYear, fireYear$advance(2, 'year'))$
    filter(ee$Filter$dayOfYear(startday, endday))$
    mean()$
    rename('postNBR')
                              
  ## if any pixels within fire have only one 'scene' or less, add additional year for imagery window
  postFireIndices2 <-  lsCol$filterBounds(fireBounds)$
    filterDate(postFireYear, fireYear$advance(3, 'year'))$
    filter(ee$Filter$dayOfYear(startday, endday))$
    mean()$
    rename('postNBR')
  post_check <-  postFireIndices$clip(fire)
  post_filled <-  post_check$unmask(postFireIndices2)
  fireIndices <-  pre_filled4$addBands(post_filled)
                              
  ## create fire severity indices    
  ## calculate dNBR 
                              
  burnIndices <-  fireIndices$expression(
    "(b('preNBR') - b('postNBR')) * 1000")$
    rename('dnbr')$toFloat()$addBands(fireIndices)
  
  ## calculate dNBR with Offset developed from 180-m ring outside the fire perimeter
  ring   <-  fire$buffer(180)$difference(fire, maxError=0.1)
  burnIndices2 <-  ee$Image$constant(ee$Number(burnIndices$select('dnbr')$reduceRegion(
    reducer = ee$Reducer$mean(),
    geometry = ring$geometry(),
    scale = 30,
    maxPixels = 1e9
  )$get('dnbr')))$rename('offset')$toFloat()$addBands(burnIndices) 
                                              
  burnIndices3 <-  burnIndices2$expression(
    "b('dnbr') - b('offset')")$
    rename('dnbr_w_offset')$toFloat()$addBands(burnIndices2)
                              
  ## calculate RBR
  burnIndices4 <-  burnIndices3$expression(
    "b('dnbr') / (b('preNBR') + 1.001)")$
    rename('rbr')$toFloat()$addBands(burnIndices3)
                              
  ## calculate RBR with offset  
  burnIndices5 <-  burnIndices4$expression(
    "b('dnbr_w_offset') / (b('preNBR') + 1.001)")$
  rename('rbr_w_offset')$toFloat()$addBands(burnIndices4)
                              
                               
  burnIndices7 <-  burnIndices6$expression(
    "(b('dnbr') / sqrt(b('preNBR2')))")$
    rename('rdnbr')$toFloat()$addBands(burnIndices6)
  
  ## calculate RdNBR with offset
  burnIndices8 <-  burnIndices7$expression(
    "b('dnbr_w_offset') / sqrt(b('preNBR2'))")$
    rename('rdnbr_w_offset')$toFloat()$addBands(burnIndices7)
  
  burnIndices8 <-  burnIndices8$select(bandList)
  burnIndices8$set(
    list(
      Fire_ID=ft$get('Fire_ID'),
      FIRENAME=ft$get('FIRENAME'),
      Year=ft$get('Year')
    )
  )
}))

In trying to further debug the indices imagecollection I tried using an example you have at https://r-spatial.github.io/rgee/reference/ee_imagecollection_to_local.html to export it locally. I came up with the following error

Error in py_call_impl(callable, dots$args, dots$keywords) : 
  EEException: Image.load: Image asset 'users/mtd25/burnIndices7' not found.

which would seem to me that there may be a problem with burnIndices6 not importing correctly to #7.