Last active
January 23, 2021 22:03
-
-
Save csaybar/4f79c1dd63ea243d0d086a4bbd3829f3 to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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.
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Seems to be working but I cannot tell for sure until I get the for loops at the end sorted out.
As for the final chunk of code I think I might have to use a while loop or just ignore the '(var j = 0; j < nFires; j++)' statement and just loop through all fires or nBands..?? Not really sure yet.