r-spatial / rgee

Google Earth Engine for R
https://r-spatial.github.io/rgee/
Other
677 stars 146 forks source link

Guidance regarding writing functions #197

Closed bmaitner closed 2 years ago

bmaitner commented 2 years ago
library(reticulate)
python:         /home/rstudio/.virtualenvs/rgee/bin/python
libpython:      /home/rstudio/.local/share/r-miniconda/envs/r-reticulate/lib/libpython3.6m.so
pythonhome:     /home/rstudio/.virtualenvs/rgee:/home/rstudio/.virtualenvs/rgee
version:        3.6.13 | packaged by conda-forge | (default, Sep 23 2021, 07:56:31)  [GCC 9.4.0]
numpy:          /home/rstudio/.virtualenvs/rgee/lib/python3.6/site-packages/numpy
numpy_version:  1.19.5
ee:             /home/rstudio/.virtualenvs/rgee/lib/python3.6/site-packages/ee

Description

I'm trying to write a function to alter the MODIS NDVI band "DayOfYear" to give the day relative to the start of the UNIX era (1-1-1970). While the code isn't too terribly difficult to get working in general, I'm having trouble figuring out how the function needs to be structured in order to work with $map(). In the NDVI example for rgee, there are a few custom functions made to filter out the bad data using the QA band and these rely on base R functions without any problems. The custom function I wrote relies on other packages beyond base R, so perhaps this is the issue? Any general guidance or resources you can point me to would be much appreciated!

What I Did

ee_Initialize(drive = TRUE)

#Example code from the NDVI tutorial:
  modis_ndvi <- ee$ImageCollection("MODIS/006/MOD13A2")

  getQABits <- function(image, qa) {
    # Convert binary (character) to decimal (little endian)
    qa <- sum(2^(which(rev(unlist(strsplit(as.character(qa), "")) == 1))-1))
    # Return a mask band image, giving the qa value.
    image$bitwiseAnd(qa)$lt(1)
  }

  mod13A2_clean <- function(img) {
    # Extract the NDVI band
    ndvi_values <- img$select("NDVI")

    # Extract the quality band
    ndvi_qa <- img$select("SummaryQA")

    # Select pixels to mask
    quality_mask <- getQABits(ndvi_qa, "11")

    # Mask pixels with value zero.
    ndvi_values$updateMask(quality_mask)
  }

#As expected, the function works fine for either an image or image collection
  first_ndvi_clean <- mod13A2_clean(modis_ndvi$first())
  ndvi_clean <- modis_ndvi$map(mod13A2_clean)

# However, in building my own function I'm running into problems:
  get_integer_date <-function(img) {

    # extract the DayOfYear band
      day_values <- img$select("DayOfYear")

     # Get the first day of the year and convert to an integer.
      start_date <- rgee::ee_get_date_img(img)$time_start #can't use R code with map
      focal_origin <- lubridate::floor_date(x = start_date,unit = "years") #this gets the first date of the year.  all modis dates are relative to this
      focal_origin <- as.numeric(as.Date(focal_origin))-1 #1 is subtracted so that a value of 1 on the raster corresponds to a Jan 1, not Jan 2

    #Mask to only values greater than zero.
      mask <- day_values$gt(0)
      day_values <- day_values$updateMask(mask)

    # #Now, I just need to add the origin to the map
      day_values$add(focal_origin)

  }  

#This works for a single image:

  first_day_integer <- get_integer_date(img = modis_ndvi$first())

#But when attempting to map, I run into errors:

  all_days_integer <- modis_ndvi$map(get_integer_date)

# The example functions getQABits and mod13A2_clean both rely on R functions, so it isn't clear to me why some R functions can be used with $map() and others can't.  Any guidance here would be much appreciated!
csaybar commented 2 years ago

Hi @bmaitner you can not combine local-server objects into a map function. More details here https://developers.google.com/earth-engine/guides/client_server

The function below replicates your client-server function.

library(rgee)

ee_Initialize(drive = TRUE)

# However, in building my own function I'm running into problems:
get_integer_date <-function(img) {
  # 1. Extract the DayOfYear band
  day_values <- img$select("DayOfYear")

  # 2. Get the first day of the year and convert to an integer.
  start_date <- ee$Date(day_values$get("system:time_start"))
  base_date <- ee$Date("1970-01-01")

  # 3. Get the day diff betwen start_date and UNIX date
  daydiff <- start_date$difference(base_date, "day")

  # 4. Mask to only values greater than zero.
  mask <- day_values$gt(0)
  day_values <- day_values$updateMask(mask)

  # #Now, I just need to add the origin to the map
  day_values$add(daydiff)
}

modis_ndvi <- ee$ImageCollection("MODIS/006/MOD13A2")
all_days_integer <- modis_ndvi$map(get_integer_date)

gviz <- list(min=10980, max= 11080, palette = c("red", "green"))
Map <- R6Map$new()
Map$addLayer(all_days_integer$first(), gviz)
Map$addLegend(gviz)
Map
bmaitner commented 2 years ago

Ahhh, I think I see. Thanks so much for pointing me to that resource and for the assistance with the code! Apologies for asking a question that was answered elsewhere