转换 GEE 脚本以与 rgee

编程问答 2022-05-29

问题描述

我正在尝试使用 rgee 包在 Rstudio 中运行为 Google Earth Engine 编写的脚本。我对 Javascript 一无所知，将它转换为在 r 中工作基本上是我的头撞墙，直到有一些工作经验为止。

无论如何，通过下面的代码，我进入了索引/ft 函数部分（第 4 个代码块）并收到以下消息。

.Call(_reticulate_py_str_impl,x) 中的错误：已达到已用时间限制

我对 inf 上的所有内容都使用了 setTimeLimit 和 setSessionTimeLimit 函数（从下面的代码示例中删除），但它似乎没有帮助，我收到了相同的消息。我还使用了 gc() 来查看是否可以清理一些空间，但这没有任何效果。我该怎么做才能让它在此时停止冻结？

此外，如果您能使用 rgee 包将所有脚本翻译为在 Rstudio 中工作，我将不胜感激。因此，如果您发现我在其他任何地方犯了错误而我没有发现或知道如何处理最后的奖励块，请告诉我。

这是原始脚本，以防有人想看 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()

数据框

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),"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\",\n            ANGLEUNIT[\"Degree\",0.0174532925199433]]],\n    CONVERSION[\"unnamed\",\n        METHOD[\"LAmbert Conic Conformal (2SP)\",9802]],\n        ParaMETER[\"Latitude of false origin\",0.0174532925199433],8821]],\n        ParaMETER[\"Longitude of false origin\",-108,8822]],\n        ParaMETER[\"Latitude of 1st standard parallel\",32,8823]],\n        ParaMETER[\"Latitude of 2nd standard parallel\",36,8824]],\n        ParaMETER[\"Easting at false origin\",\n            LENGTHUNIT[\"metre\",1],8826]],\n        ParaMETER[\"northing at false origin\",8827]]],\n    CS[Cartesian,2],\n        AXIS[\"(E)\",east,\n            ORDER[1],1,\n                ID[\"epsg\",9001]]],\n        AXIS[\"(N)\",north,\n            ORDER[2],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")
  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",'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)

现在这是我遇到问题的部分。逐行调试时，我得到 “.Call(_reticulate_py_str_impl,x) 中的错误：已达到已用时间限制” 一旦我到达burnIndices3 以及随后的每个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,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,endday))$
    mean()$
    rename('preNBR')
  pre_filled2<- pre_filled$unmask(preFireIndices3)
  
  ## four years back
  preFireYear4 <-  fireYear$advance(-4,'year')
  preFireIndices4 <-  lsCol$filterBounds(fireBounds)$
    filterDate(preFireYear4,endday))$
    mean()$
    rename('preNBR')
  pre_filled3 <- pre_filled2$unmask(preFireIndices4)
  
  ## five years back
  preFireYear5 <-  fireYear$advance(-5,'year')
  preFireIndices5 <-  lsCol$filterBounds(fireBounds)$
    filterDate(preFireYear5,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,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')}
  ) 

}))

我还没有使用的额外代码块。

## ## 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',id).first()).geometry().bounds();
  firesname <-  fireName[j];
  Nameoffire <-  firesname;
  fireExport <-  ee$Image(indices.filterMetadata('fireName',firesname).first());
  fireBounds <-  ee$Feature(fires.filterMetadata('FIRENAME',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
    })
  }}

解决方法

我发现这里有几个问题

  1. 在解析日期之前将 ee.Number 对象转换为 ee.String。
## create pre- and post-fire imagery
fireYear <-  ee$Date$parse('YYYY',ee$Number$format(fire$get('Year')))
  1. 删除 RdNBR 中的加号运算符。
  ## 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)
  1. 将 {} 更改为列表：
burnIndices8$set(
  list(
    Fire_ID=ft$get('Fire_ID'),FIRENAME=ft$get('FIRENAME'),Year=ft$get('Year')
  )
)

为避免出现达到已用时间限制之类的问题，我强烈建议您查看https://r-earthengine.github.io/intro_03/第 18 节。请在此处查看可重现的示例：https://gist.github.com/csaybar/4f79c1dd63ea243d0d086a4bbd3829f3

google-earth-enginejavascriptrrrgee