I am attempting to run a script written for Google Earth Engine within Rstudio using the rgee package. I don't really know anything about Javascript and converting it to work within r has basically been a bash my head against a wall until something works experience.
Anyways in working through the below code I get to the indices/ft function section (4th code chunk) and receive the following message.
Error in .Call(_reticulate_py_str_impl, x) : reached elapsed time limit
I have used the setTimeLimit and setSessionTimeLimit functions with everything on inf (removed from the code examples below) but it did not seem to help and I receive the same message. I also used gc() to see if I could clear up some space but that had no effect. What can I do to make it stop freezing up at this point?
Also, I would appreciate any help in translating all of this script to work within Rstudio using the rgee package. So if you see anywhere else that I made an error that I am not catching or know how to approach the last bonus chunk at the end please let me know.
Here's the original script in case anyone wants to look https://code.earthengine.google.com/?scriptPath=users%2Fmtd25%2FFire_severity%3AFire_atlas
library(rgee)
earthengine_python <- Sys.getenv("EARTHENGINE_PYTHON", unset = NA)
print(earthengine_python)
Sys.setenv(RETICULATE_PYTHON = earthengine_python)
ee_current_version <- system.file("python/ee_utils.py", package = "rgee")
ee_utils <- rgee:::ee_source_python(ee_current_version)
print(ee_utils$ee$'__version__')
library(dplyr)
library(sf)
library(googledrive)
library(googleCloudStorageR)
library(lwgeom)
library(reticulate)
library(jsonlite)
library(tidyverse)
# use for all other accounts
ee_Initialize()
Dataframe
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")))
Inputs and setup
## 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, [email protected]
## Lisa Holsinger, [email protected]
## Sean Parks, [email protected]
## 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")
nbr$addBands(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)
Now here is the part I have problems with. When debugging line by line I get the "Error in .Call(_reticulate_py_str_impl, x) : reached elapsed time limit" message once I get to burnIndices3 and after each subsequent bunindices 4 - 8.
## ------------------ 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', 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)
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)
## calculate RdNBR
burnIndices6 <- burnIndices5$expression(
"(b('preNBR') < 0.001) ? 0.001 +
: b('preNBR')")$
sqrt()$rename('preNBR2')$toFloat()$addBands(burnIndices5)
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({
ft$get('Fire_ID')
ft$get('FIRENAME')
ft$get('Year')}
)
}))
Bonus chunk of code I haven't gotten to yet.
## ## Export fire indices for each fire
nBands <- bandList$length;
for (j = 0; j < nFires; j++){
id <- fireID[j];
Name <- id;
fireExport <- ee$Image(indices.filterMetadata('fireID', 'equals', id).first());
fireBounds <- ee$Feature(fires.filterMetadata('Fire_ID', 'equals', id).first()).geometry().bounds();
firesname <- fireName[j];
Nameoffire <- firesname;
fireExport <- ee$Image(indices.filterMetadata('fireName', 'equals', firesname).first());
fireBounds <- ee$Feature(fires.filterMetadata('FIRENAME', 'equals', firesname).first()).geometry().bounds();
for (i == 0; i < nBands; i++) {
bandExport <- bandList[i];
exportImg <- fireExport.select(bandExport);
Export.image.toDrive({
image: exportImg,
##description: Name + '_' + bandExport,
##fileNamePrefix: Name + '_' + bandExport,
description: Name + '_' + Nameoffire + '_' + bandExport,
fileNamePrefix: Name + '_' + Nameoffire + '_' + bandExport,
maxPixels: 1e13,
scale: 30,
crs: "EPSG:4326",
folder: 'fires',
region: fireBounds
})
}}
I see a couple of problems here
## create pre- and post-fire imagery
fireYear <- ee$Date$parse('YYYY', ee$Number$format(fire$get('Year')))
## calculate RdNBR
burnIndices6 <- burnIndices5$expression(
"(b('preNBR') < 0.001) ? 0.001: b('preNBR')")$
# if("b('preNBR') < 0.001"){
# "b('preNBR') + 0.001"
# } else {
# "b('preNBR')"})$
sqrt()$rename('preNBR2')$toFloat()$addBands(burnIndices5)
{} to List:burnIndices8$set(
list(
Fire_ID=ft$get('Fire_ID'),
FIRENAME=ft$get('FIRENAME'),
Year=ft$get('Year')
)
)
To avoid problems like reached elapsed time limit I highly recommend taking a look at https://r-earthengine.github.io/intro_03/ Section 18. See a reproducible example here: https://gist.github.com/csaybar/4f79c1dd63ea243d0d086a4bbd3829f3