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.288 -- 
##  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 
## --------------------------------------------------------------------------------

Análisis de MDE

Trabajaremos con los datos de Shuttle Radar Topography Mission (SRTM), cuyos modelos digitales de elevación (DEM) son a una escala casi global. Este producto SRTM V3 (SRTM Plus) es proporcionado por NASA JPL con una resolución de 1 segundo de arco (aproximadamente 30m.). Este conjunto de datos se ha sometido a un proceso de llenado de vacíos utilizando datos de código abierto (ASTER GDEM2, GMTED2010 y NED), a diferencia de otras versiones que contienen vacíos o se han llenado de vacíos con fuentes comerciales (Farr et al., 2007 ) . Para obtener más información: https://developers.google.com/earth-engine/datasets/catalog/USGS_SRTMGL1_003#description

Para más conjuntos de datos DEM en Earth Engine : https://developers.google.com/earth-engine/datasets/tags/dem o https://github.com/csaybar/rgee/tree/examples/Datasets/Terrain

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

## Reading layer `basinIchu' from data source 
##   `C:\Users\RONALDO\Desktop\rgee2021\shp\basinIchu.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 0 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 474113.3 ymin: 8552958 xmax: 525192.8 ymax: 8614381
## Projected CRS: WGS 84 / UTM zone 18S

geom_ichu<- basin_ichu$geometry()

Modelo digital de elevación a escala mundial

elev <- ee$Image("USGS/SRTMGL1_003")$
  select('elevation')

# Características de la ee$Image
cat("Band names: ", paste(elev$bandNames()$getInfo(), collapse = " "))

## Band names:  elevation

cat("type: ", elev$projection()$getInfo()$type)

## type:  Projection

cat("crs: ", elev$projection()$getInfo()$crs)

## crs:  EPSG:4326

cat("geotransform: ", paste0(elev$projection()$getInfo()$transform, " "))

## geotransform:  0.000277777777777778  0  -180.000138888889  0  -0.000277777777777778  60.0001388888889

# Resolución de la ee$Image
cat("SRTM resolution: ", elev$projection()$nominalScale()$getInfo())

## SRTM resolution:  30.92208

Visualización

Map$centerObject(geom_ichu, zoom = 3)
vizPrm<- list( min = 0,
               max = 6000,
               palette =c("006600", "002200", "fff700",
                          "ab7634", "c4d0ff", "ffffff"))
Map$addLayer(elev,visParams =vizPrm)

# Recortamos en zona de interés
dem_ichu_clip<- elev$clip(basin_ichu)

Map$centerObject(geom_ichu, zoom = 10)
Map$addLayer(dem_ichu_clip,visParams =vizPrm)

Descarga de imágen

# Descarga desde URL de la imágen.

geom_params <- list(
  scale = 30,
  crs = 'EPSG:32718',
  region = geom_ichu
)

path <- elev$getDownloadUrl(geom_params)
print(path)

## [1] "https://earthengine.googleapis.com/v1alpha/projects/earthengine-legacy/thumbnails/6f5b53842106d99b3f07a483b2509755-8fd40860cb3fd9e439ca5cb8713776c9:getPixels"

browseURL(path)

Mapa de pendiente, aspecto y sombra

Establacemos las siguientes funciones (Visite : https://github.com/csaybar/rgee/blob/examples//Visualization/hillshade.R)

Radians

Radians <- function(img) {
  img$toFloat()$multiply(base::pi)$divide(180)
}

Hillshade

# Compute hillshade for the given illumination az, el.
Hillshade <- function(az, ze, slope, aspect) {
  azimuth <- Radians(ee$Image(az))
  zenith <- Radians(ee$Image(ze))
  azimuth$subtract(aspect)$cos()$
    multiply(slope$sin())$
    multiply(zenith$sin())$
    add(zenith$cos()$multiply(slope$cos()))
}

Mapa de pendientes en grados

Para visualizar, gestionamos una sencilla paleta de colores y establecemos valores mínimos y máximos de 0 y 90, respectivamente.

terrain <- ee$Algorithms$Terrain(dem_ichu_clip)
slope_ichu_d<- terrain$select("slope")

Map$centerObject(geom_ichu, zoom = 10)
viz_slope<- list( min = 0,
                  max = 90,
                  palette=c('3ae237','b5e22e','d6e21f','fff705',
                      'ffd611','ffb613','ff8b13','ff6e08','ff500d',
                      'ff0000','de0101','c21301','0602ff','235cb1',
                      '307ef3','269db1','30c8e2','32d3ef','3be285',
                      '3ff38f','86e26f'))
                        
Map$addLayer(slope_ichu_d,visParams =viz_slope)

Mapa de pendientes en radianes

slope_ichu_r <- Radians(terrain$select("slope"))

Map$centerObject(geom_ichu, zoom = 10)
viz_slope<- list( min = 0,
                  max = 3.1416/2,
                  palette=c('3ae237','b5e22e','d6e21f','fff705',
                      'ffd611','ffb613','ff8b13','ff6e08','ff500d',
                      'ff0000','de0101','c21301','0602ff','235cb1',
                      '307ef3','269db1','30c8e2','32d3ef','3be285',
                      '3ff38f','86e26f'))
                        
Map$addLayer(slope_ichu_r,visParams =viz_slope)

Mapa de orientación de laderas (Aspect)

aspect_ichu_d<- terrain$select("aspect")

Map$centerObject(geom_ichu, zoom = 10)
viz_asp<- list( min = 0,
                max = 360,
                palette=c('3ae237','b5e22e','d6e21f','fff705',
                          'ffd611','ffb613','ff8b13','ff6e08',
                          'ff500d','ff0000','de0101','c21301', 
                          '0602ff','235cb1','307ef3','269db1',
                          '30c8e2','32d3ef','3be285','3ff38f',
                          '86e26f'))

Map$addLayer(aspect_ichu_d,visParams =viz_asp)

Mapa de sombras (Hillshade)

More Colors in R : https://r-charts.com/

#Definiendo paleta de colores para hillshade
col = colorRampPalette(c(rgb(0, 50, 100, maxColorValue = 255), rgb(255, 240, 170, maxColorValue = 255)))(25)

aspect_ichu_r<- Radians(terrain$select("aspect"))

hillShade_ichu<- Hillshade(0, 60,slope_ichu_r, aspect_ichu_r)

viz_hls<- list( palette=col)
Map$centerObject(geom_ichu, zoom = 10)
Map$addLayer(hillShade_ichu, visParams =viz_hls)

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/

Farr, T. G., Rosen, P. A., Caro, E., Crippen, R., Duren, R., Hensley, S., Kobrick, M., Paller, M., Rodriguez, E., Roth, L., Seal, D., Shaffer, S., Shimada, J., Umland, J., Werner, M., Oskin, M., Burbank, D., & Alsdorf, D. (2007). The shuttle radar topography mission. Reviews of Geophysics, 45(2). https://doi.org/https://doi.org/10.1029/2005RG000183

Google Earth Engine con R

Manejo de Modelo Digital de Elevación (DEM)

Autor: Ronaldo Anticona Guillen