gistfile1.R

library(rgee)

# -------------------------------------------------------------------------
# 1. Load Earth Engine ----------------------------------------------------
# -------------------------------------------------------------------------

ee_Initialize(drive = TRUE)

# -------------------------------------------------------------------------
# 2. Auxiliary functions --------------------------------------------------
# -------------------------------------------------------------------------

cloudMaskL573 <- function(image) {
  # Select qa band
  qa <- image$select("pixel_qa")
  
  # Create a mask
  cloud <- qa$bitwiseAnd(bitwShiftL(1, 5))$
    And(qa$bitwiseAnd(bitwShiftL(1, 7)))$
    Or(qa$bitwiseAnd(bitwShiftL(1, 3)))
  
  # Estimate the cloud percentage according to the study area
  pixel_qav <- cloud %>% 
    ee$Image$reduceRegion(
      reducer = ee$Reducer$sum(),
      geometry = region) %>% 
    ee$Dictionary$get("pixel_qa") %>% 
    ee$Number() %>% 
    ee$Number$divide(426.83)
  
  # Apply the mask, estimate the NDVI and transform to Int16 to make the download faster. 
  mask2 <- image$mask()$reduce(ee$Reducer$min())
  image$updateMask(cloud$Not())$updateMask(mask2)$
    normalizedDifference(c("B4", "B3"))$
    multiply(1000)$toInt16()$
    copyProperties(image, list("system:time_start", "system:id"))$ #save some properties from the original image
    set(list(cloud_percentage = pixel_qav))
}

# Function to cloud mask from the pixel_qa band of Landsat 8 SR data.
cloudMaskL53 <- function(image) {
  # Select qa band
  qa <- image$select("pixel_qa")
  
  # Create a mask
  cloudShadowBitMask <- bitwShiftL(1, 3)
  cloudsBitMask <- bitwShiftL(1, 5)
  cloud <- qa$bitwiseAnd(cloudShadowBitMask)$eq(0)$
    And(qa$bitwiseAnd(cloudsBitMask)$eq(0))
  
  # Estimate the cloud percentage according to the study area
  pixel_qav <- cloud$Not() %>% 
    ee$Image$reduceRegion(
      reducer = ee$Reducer$sum(),
      geometry = region) %>% 
    ee$Dictionary$get("pixel_qa") %>% 
    ee$Number() %>% 
    ee$Number$divide(426.83)
  
  # Apply the mask, estimate the NDVI and transform to Int16 to make the download faster. 
  image$updateMask(cloud)$
    normalizedDifference(c("B5", "B4"))$
    multiply(1000)$toInt16()$
    copyProperties(image, list("system:time_start", "system:id"))$
    set(list(cloud_percentage = pixel_qav))
}

# -------------------------------------------------------------------------
# 3. Define a region ------------------------------------------------------
# -------------------------------------------------------------------------

region <- ee$Geometry$Rectangle(coords = c(-72.03,-13.41,-71.95,-13.37)) 
# Map$centerObject(region)
# Map$addLayer(region)



# -------------------------------------------------------------------------
# 4. Download Landsat5 - NDVI ---------------------------------------------
# -------------------------------------------------------------------------

# Map the function over the collection.
l5_composite_ndvi <- ee$ImageCollection("LANDSAT/LT05/C01/T1_SR") %>%
  ee$ImageCollection$filterBounds(region) %>% 
  ee$ImageCollection$map(cloudMaskL573) %>% 
  ee$ImageCollection$filterMetadata("cloud_percentage", "less_than", 10)

# Download!
setwd("/home/csaybar/Desktop/carlos/l5")
dates_ic <- basename(ee_get_date_ic(l5_composite_ndvi)$id)
ee_imagecollection_to_local(
  ic = l5_composite_ndvi, 
  region = region, 
  dsn =  dates_ic
)


# -------------------------------------------------------------------------
# 5. Download Landsat8 - NDVI ---------------------------------------------
# -------------------------------------------------------------------------

# Map the function over the collection.
l8_composite_ndvi <- ee$ImageCollection("LANDSAT/LC08/C01/T1_SR") %>% 
  ee$ImageCollection$filterBounds(region) %>% 
  ee$ImageCollection$map(cloudMaskL53) %>% 
  ee$ImageCollection$filterMetadata("cloud_percentage", "less_than", 10)

# Download!
setwd("/home/csaybar/Desktop/carlos/l8")
dates_ic <- basename(ee_get_date_ic(l8_composite_ndvi)$id)
ee_imagecollection_to_local(
  ic = l8_composite_ndvi, 
  region = region, 
  dsn =  dates_ic
)