Google Earth Engine con R

Índices diferenciales normalizados

Usando GEE en R

Inicializando rgee

Google earth engine (GEE) es una plataforma computacional gratis en la nube que está diseñada para el análisis de datos espaciales a una escala planetaria.

Podemos explorar, descargar, manipular, y procesar cientos de imágenes satelitales (Sentinel, Landsat, Modis, etc.) de cualquier parte del mundo en cuestión de minutos mediante la integración del ecosistema de r spatial con Google Earth Engine y viceversa. Actualmente, Google ofrece soporte sólo para Python y JavaScript, es aquí donde rgee llenará el vacío comenzando a brindar soporte a R!, con una sintaxis muy similar a las otras dos bibliotecas cliente compatibles con Google (Aybar et al., 2020).

Es necesario tener instalado “RTolls”, así poder instalar tanto los paquete rgee y googledrive. La inicialización se dá desde Google drive, por tanto, es importante una cuenta Gmail.

# devtools::install_github("tidyverse/googledrive")
# install.packages('rgee')

## Importar librerias
library(googledrive)
library(rgee)

## Paquetes de rgee
# ee_install()

# Inicializando rgee
ee_Initialize('ronaldoag1007@gmail.com',drive=T)

## -- rgee 1.1.0 --------------------------------------- earthengine-api 0.1.283 -- 
##  v user: ronaldoag1007@gmail.com 
##  v Google Drive credentials:
 v Google Drive credentials:  FOUND
##  v Initializing Google Earth Engine:
 v Initializing Google Earth Engine:  DONE!
## 
 v Earth Engine account: users/RonaldoAnticona 
## --------------------------------------------------------------------------------

Selección del directorio de trabajo y asignamos geométria al shp cargado:

library(sf)
library(rgdal)
setwd("C:/Users/RONALDO/Desktop/rgee2021/ndwi")
basin_ichu<- st_read("shp/ichu.shp")%>%
  sf_as_ee()

## Reading layer `ichu' from data source 
##   `C:\Users\RONALDO\Desktop\rgee2021\ndwi\shp\ichu.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 1 field
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -75.23861 ymin: -13.08958 xmax: -74.76798 ymax: -12.53413
## Geodetic CRS:  WGS 84

geom_ichu<- basin_ichu$geometry()

Manejo de imágen satelital

Empleo de la imágen de Landsat 8, colección 1, nivel 1 de reflectancia calibrada en la parte superior de la atmósfera (TOA)(Chander et al., 2009), consisten en Números Digitales (DN) escalados, cuantificados y calibrados que representan los datos de imágenes multiespectrales. Los datos de los productos Landsat 8 adquiridos tanto por Operational Land Imager (OLI) como por el sensor térmico de infrarrojos (TIRS) se entregan en formato de entero sin signo de 16 bits. Las escenas Landsat con la mayor calidad de datos disponible se colocan en el Nivel 1 y se consideran adecuadas para el análisis de series de tiempo. Ver más:https://developers.google.com/earth-engine/datasets/catalog/LANDSAT_LC08_C01_T1_TOA#description

El snippet correspondiente es: ee.ImageCollection(“LANDSAT/LC08/C01/T1_TOA”), y contiene las bandas:

Name	Pixel Size	Wavelength	Description
B1	30 meters	0.43 - 0.45 µm	Coastal aerosol
B2	30 meters	0.45 - 0.51 µm	Blue
B3	30 meters	0.53 - 0.59 µm	Green
B4	30 meters	0.64 - 0.67 µm	Red
B5	30 meters	0.85 - 0.88 µm	Near infrared
B6	30 meters	1.57 - 1.65 µm	Shortwave infrared 1
B7	30 meters	2.11 - 2.29 µm	Shortwave infrared 2
B8	15 meters	0.52 - 0.90 µm	Band 8 Panchromatic
B9	15 meters	1.36 - 1.38 µm	Cirrus
B10	30 meters	10.60 - 11.19 µm	Thermal infrared 1, resampled from 100m to 30m
B11	30 meters	11.50 - 12.51 µm	Thermal infrared 2, resampled from 100m to 30m
BQA	30 meters		Landsat Collection 1 QA Bitmask

Cargar ImageCollection de GEE y verificar las imágenes menores al 10% de nubosidad, en el rango de fecha establecido:

coll<-ee$ImageCollection('LANDSAT/LC08/C01/T1_TOA')$
  filterDate('2019-01-01','2020-06-30')$
  filterBounds(ee$Geometry$Point(-74.97,-12.78))$
  filterMetadata('CLOUD_COVER','less_than',10)

ee_get_date_ic(coll)

##                                             id          time_start
## 1 LANDSAT/LC08/C01/T1_TOA/LC08_006069_20190626 2019-06-26 15:05:03
## 2 LANDSAT/LC08/C01/T1_TOA/LC08_006069_20190712 2019-07-12 15:05:06
## 3 LANDSAT/LC08/C01/T1_TOA/LC08_006069_20200511 2020-05-11 15:04:34
## 4 LANDSAT/LC08/C01/T1_TOA/LC08_006069_20200527 2020-05-27 15:04:37
## 5 LANDSAT/LC08/C01/T1_TOA/LC08_006069_20200612 2020-06-12 15:04:47
## 6 LANDSAT/LC08/C01/T1_TOA/LC08_006069_20200628 2020-06-28 15:04:55

Selección de imágen de interés:

L8<-'LANDSAT/LC08/C01/T1_TOA/LC08_006069_20190712' %>%
  ee$Image() %>%
  ee$Image$select(c('B2','B3','B4','B5','B6','B7'),
                  c('BLUE','GREEN','RED','NIR','SWIR 1','SWIR 2'))

Visualización en combinación de bandas (B6,B5,B4), el cuál muestra una mayor diferenciación entre el suelo y el agua. La vegetación se muestra en diversas tonalidades de verde y rosa, que varían en función del tipo y de las condiciones de ubicación; las áreas urbanas y el suelo expuesto se presentan en tonos rosados; el agua, independiente de la cantidad de sedimentos en suspensión, aparece en negro. Más sobre combinación de bandas: https://www.hidraulicafacil.com/2016/03/Com.Landsat7.html

viz_654<- list(bands=c('SWIR 1','NIR','RED'),
               min=0.1,
               max=0.8,
               gamma=1.6)

Map$centerObject(L8, zoom = 8)
Map$addLayer(eeObject = L8, visParams =viz_654)

Recorte en zona de interés y visualización:

ichu_clip<- L8$clip(basin_ichu)

Map$centerObject(ee$Geometry$Point(list(-75.012367,-12.771628)), zoom =11)
Map$addLayer(ichu_clip, visParams =viz_654)

Cálculo de NDWI

El Índice Diferencial de Agua Normalizado (NDWI) derivados de imágenes de satélite se utiliza con frecuencia y con éxito en la detección y mapeo de masas de agua superficial (JRC, 2011; Özelkan, 2020).

Más sobre índices diferenciales:https://acolita.com/lista-de-indices-espectrales-en-sentinel-2-y-landsat/

Función para cálculo de NDWI

# Método de McFeeters (1996)
F_NDWI<- function (image){
  ndwi<-image$expression('float(GREEN-NIR)/(GREEN+NIR)',
                         opt_map=list("GREEN"=image$select("GREEN"),
                                      "NIR"=image$select("NIR")))$rename("NDWI")
  return(ndwi)
}

ndwi_ichu_clip <- F_NDWI(ichu_clip)

Visualización de NDWI:

viz_ndwi<- list(min=-1,
                max=1,
                palette=c('FF2A00','FF6619','FFCC65','FFEE99',
                          'CCFFFF','99EDFF','3299FF','075AFF'))
Map$centerObject(ee$Geometry$Point(list(-75.012367,-12.771628)), zoom=11)
Map$addLayer(ndwi_ichu_clip,visParams =viz_ndwi)

Selección de NDWI>=0.2 y enmascaramiento:

ndwi_agua<- ndwi_ichu_clip$gte(0.2)
ndwi_mask <- ndwi_ichu_clip$updateMask(ndwi_agua)
# Reproyectar mascara
ndwi_mask_utm <- ndwi_mask$reproject(crs="EPSG:32718", scale=30)

Visualización de máscara NDWI>=0.2:

Map$centerObject(ee$Geometry$Point(list(-75.0546,-12.69039)), zoom = 13)
Map$addLayer(ndwi_mask_utm,visParams =viz_ndwi)

Serie de tiempo diária de NDVI

El producto MOD13A2 Versión 6 proporciona valores de índice de vegetación (VI) por píxel a una resolución espacial de 1 kilómetro (km). Hay dos capas de vegetación primarias. El primero es el Índice de Vegetación de Diferencia Normalizada (NDVI), que se conoce como el índice de continuidad del NDVI derivado del Radiómetro Avanzado de Muy Alta Resolución de la Administración Nacional Oceánica y Atmosférica existente (NOAA-AVHRR). La segunda capa de vegetación es el Índice de Vegetación Mejorado (EVI), que ha mejorado la sensibilidad en las regiones de alta biomasa. El algoritmo de este producto elige el mejor valor de píxel disponible de todas las adquisiciones del período de 16 días. Los criterios utilizados son nubes bajas, ángulo de visión bajo y el valor más alto de NDVI / EVI (USGS, 2021).

Las principales capas presentan las siguientes características:

SDS Name	Description	Units	Data Type	Fill Value	No Data Value	Valid Range	Scale Factor
1 km 16 days NDVI	1 km 16 days NDVI	NDVI	16-bit signed integer	-3000	N/A	-2000 to 10000	0.0001
1 km 16 days EVI	1 km 16 days EVI	EVI	16-bit signed integer	-3000	N/A	-2000 to 10000	0.0001
1 km 16 days VI Quality	VI Quality indicators	Bit Field	16-bit unsigned integer	65535	N/A	0 to 65534	N/A

El fragmento de EE viene a ser: ee.ImageCollection(“MODIS/006/MOD13A2”).

Cargar banda NDVI desde el productor MODIS MOD13A2.006, desde 01/01/16 hasta el 30/06/21 en el punto de interés:

NDVI_d<-ee$ImageCollection('MODIS/006/MOD13A2')$
  select('NDVI')$
  filterDate(rdate_to_eedate("2016-01-01"),rdate_to_eedate("2021-06-30"))$
  filterBounds(ee$Geometry$Point(-74.97,-74.97))

Reescalar valores de NDVI en el rango de 0-1 en la colección de imágenes a partir de una función.

rename_scale<-function(image){
  return(image$multiply(.0001))
}

# Aplicamos la función en la ImageCollection
NDVI_scale<- NDVI_d$map(rename_scale)

Extracción de valores de píxel en un punto de interés:

point<- ee$Geometry$Point(c(-74.97,-12.78))
ndvi_ts_p<- ee_extract(NDVI_scale, ee_as_sf(point),sf=FALSE)

Más sobre pivot_longer: https://tidyr.tidyverse.org/reference/pivot_longer.html Más sobre manejo de fechas: https://estadistica-dma.ulpgc.es/cursoR4ULPGC/6h-Fechas.html

library(tidyverse)
library(lubridate)
ndvi_ts_p<- ndvi_ts_p %>% 
  pivot_longer(everything(),names_to="date",values_to="NDVI")%>%
  mutate(date=as_date(date, format = "X%Y_%m_%d_NDVI"))
ndvi_ts_p

## # A tibble: 127 x 2
##    date        NDVI
##    <date>     <dbl>
##  1 2016-01-01 0.416
##  2 2016-01-17 0.245
##  3 2016-02-02 0.376
##  4 2016-02-18 0.286
##  5 2016-03-05 0.596
##  6 2016-03-21 0.392
##  7 2016-04-06 0.257
##  8 2016-04-22 0.358
##  9 2016-05-08 0.251
## 10 2016-05-24 0.313
## # ... with 117 more rows

Serie de tiempo mensual de NDVI en el punto de interés:

ndvi_ts_p%>% 
  ggplot() +
  geom_line(mapping = aes(x=date, y=NDVI ,group=1 ,colour = "NDVI"), size = 0.5) +
  scale_colour_manual("", 
                      breaks = c("NDVI" ),
                      values = c("#5FEF16")) +
  labs(title='NDVI: MODIS/006/MOD13A2')+
  labs(x = 'FECHA', y = expression('NDVI'))+
  geom_point(mapping = aes(x=date, y=NDVI),size=0.7,shape = 1)

Creación de NDVI en Gift

La creación de animaciones nos permite dar una idea de los cambios en los fenómenos sucitados en la zona de interes; la secuencia en la realización de este Gift, tiene bastante apoyo del siguiente url: https://github.com/r-spatial/rgee

library(ps) 
library(sf)
library(rgdal)
setwd("C:/Users/RONALDO/Documents/R/win-library/4.0/rgee/shp/")
mask <- ("ichu.shp") %>% 
  st_read(quiet = TRUE) %>% 
  sf_as_ee()
region <- mask$geometry()$bounds()
class(mask)

## [1] "ee.featurecollection.FeatureCollection"
## [2] "ee.collection.Collection"              
## [3] "ee.element.Element"                    
## [4] "ee.computedobject.ComputedObject"      
## [5] "ee.encodable.Encodable"                
## [6] "python.builtin.object"

Conjunto de datos MODIS Terra Vegetation Indices 16-Day Global 1km como un ee.ImageCollection y seleccione la banda NDVI:

col <- ee$ImageCollection('MODIS/006/MOD13A2')$select('NDVI')

Agrupar imágenes por fecha compuesta:

col <- col$map(function(img) {
  doy <- ee$Date(img$get('system:time_start'))$getRelative('day', 'year')
  img$set('doy', doy)
})
distinctDOY <- col$filterDate('2019-01-01', '2021-01-01')

Definir un filtro que identifique las imágenes de colección completa:

filter <- ee$Filter$equals(leftField = 'doy', rightField = 'doy');

Definir una combinación; convertir el FeatureCollection resultante en un ImageCollection.

join <- ee$Join$saveAll('doy_matches')
joinCol <- ee$ImageCollection(join$apply(distinctDOY, col, filter))

Aplicar la reducción media entre las colecciones DOY coincidentes:

comp <- joinCol$map(function(img) {
  doyCol = ee$ImageCollection$fromImages(
    img$get('doy_matches')
  )
  doyCol$reduce(ee$Reducer$median())
})

Definir los parámetros de visualización RGB:

visParams = list(
  min = 0.0,
  max = 9000.0,
  bands = "NDVI_median",
  palette = c(
    'FFFFFF', 'CE7E45', 'DF923D', 'F1B555', 'FCD163', '99B718', '74A901',
    '66A000', '529400', '3E8601', '207401', '056201', '004C00', '023B01',
    '012E01', '011D01', '011301'
  )
)

Crear imágenes de visualización RGB para usarlas como cuadros de animación:

rgbVis <- comp$map(function(img) {
  do.call(img$visualize, visParams) %>% 
    ee$Image$clip(mask)
})

Definir los parámetros de visualización GIF:

gifParams <- list(
  region = region,
  dimensions = 600,
  crs = 'EPSG:3857',
  framesPerSecond = 5
)

Renderiza la animación GIF en la consola:

print(rgbVis$getVideoThumbURL(gifParams))

## [1] "https://earthengine.googleapis.com/v1alpha/projects/earthengine-legacy/videoThumbnails/b8e4ddeb15f0182c5f8555865fc74715-bce5db0616dfa51887ecb42c040619d4:getPixels"

browseURL(rgbVis$getVideoThumbURL(gifParams))

Obtener nombres de meses:

dates_modis_mabbr <- distinctDOY %>%
  ee_get_date_ic %>% # Get Image Collection dates
  '[['("time_start") %>% # Select time_start column
  lubridate::month() %>% # Get the month component of the datetime
  '['(month.abb, .) # subset around month abbreviations

ee_utils_gif_*

Use las funciones ee_utils_gif_ * para renderizar la animación GIF y agregar algunos textos. https://rdrr.io/github/csaybar/rgee/src/R/ee_utils.R

Gif_creator

ee_utils_gif_creator <- function(ic, parameters, quiet = FALSE, ...) {
  # check packages
  #ee_check_packages("ee_utils_gif_creator", "magick")

  if (!quiet) {
    message("1. Creating gif ... please wait ....")
  }
  animation_url <- ee$ImageCollection$getVideoThumbURL(ic, parameters)
  temp_gif <- tempfile()
  if (!quiet) {
    message("1. Downloading GIF from: ", animation_url)
  }
  download.file(
    url = animation_url,
    destfile = temp_gif,
    quiet = quiet,
    ...)
  magick::image_read(path = temp_gif)
}

Gif_annotate

ee_utils_gif_annotate <- function(image,
                                  text,
                                  gravity = "northwest",
                                  location = "+0+0",
                                  degrees = 0,
                                  size = 20,
                                  font = "sans",
                                  style = "normal",
                                  weight = 400,
                                  kerning = 0,
                                  decoration = NULL,
                                  color = NULL,
                                  strokecolor = NULL,
                                  boxcolor = NULL) {
  # check packages
  #ee_check_packages("ee_utils_gif_annotate", "magick")

  if (length(text) == 1) {
    image <- magick::image_annotate(image, text, gravity = gravity,
                                    location = location, degrees = degrees, size = size,
                                    font = font, style = style, weight = weight,
                                    kerning = kerning, decoration = decoration,
                                    color = color, strokecolor = strokecolor,
                                    boxcolor = boxcolor)
  } else if(length(text) == length(image)) {
    image <- magick::image_annotate(image, text, gravity = gravity,
                                    location = location, degrees = degrees, size = size,
                                    font = font, style = style, weight = weight,
                                    kerning = kerning, decoration = decoration,
                                    color = color, strokecolor = strokecolor,
                                    boxcolor = boxcolor)
  } else {
    stop(
      "The text argument has not the same length as the magick-image object",
      "\nActual:",length(text),
      "\nExpected:", length(image)
    )
  }
  image
}

Gif_save

ee_utils_gif_save <- function(image,
                              path = NULL,
                              format = NULL,
                              quality = NULL,
                              depth = NULL,
                              density = NULL,
                              comment = NULL,
                              flatten = FALSE) {

  # check packages
  #ee_check_packages("ee_utils_gif_save", "magick")
  magick::image_write(image = image, path = path, format = format,
                      quality = quality, depth = depth, density = density,
                      comment = comment, flatten = flatten)
}

Realizando la animación

Animación realizada desde el ‘2019-01-01’ hasta ‘2021-01-01’

animation <- ee_utils_gif_creator(rgbVis, gifParams, mode = "wb")
animation %>%
  ee_utils_gif_annotate(
    text = "NDVI: MODIS/006/MOD13A2",
    size = 15, color = "white",
    location = "+10+10"
  ) %>%
  ee_utils_gif_annotate(
    text = dates_modis_mabbr,
    size = 30,
    location = "+240+510",
    color = "white",
    font = "arial",
    boxcolor = "#000000"
  ) 

ee_utils_gif_save(animation, path = "raster_ndvi_ichu.gif")

Referencias

Aybar, C., Wu, Q., Bautista, L., Yali, R., & Barja, A. (2020). Rgee: An r package for interacting with google earth engine. Journal of Open Source Software. https://github.com/r-spatial/rgee/

Chander, G., Markham, B. L., & Helder, D. L. (2009). Summary of current radiometric calibration coefficients for landsat mss, tm, etm+, and eo-1 ali sensors. Remote Sensing of Environment, 113(5), 893–903. https://doi.org/https://doi.org/10.1016/j.rse.2009.01.007

JRC. (2011). In Edo.jrc.ec.europa.eu. https://edo.jrc.ec.europa.eu/documents/factsheets/factsheet_ndwi.pdf

Özelkan, E. (2020). Water body detection analysis using ndwi indices derived from landsat-8 oli. Polish Journal of Environmental Studies, 29(2), 1759–1769. https://doi.org/10.15244/pjoes/110447

USGS. (2021). MOD13A2 v006. https://lpdaac.usgs.gov/products/mod13a2v006/