Uso de Google Earth Engine para elaboración de mapas

En este ejemplo se creará un mapa utilizando directamente los datos en la nube de Google Earth Engine.

Cargar libreías

Antes de ejecutar el siguiente código, se debe instalar el paquete rgee siguiendo las instrucciones de su repositorio https://github.com/r-spatial/rgee. Los demás paquetes se pueden instalar con la función pacman::p_load() si es que no lo tiene.

## Carga de librerías
require(pacman)

pacman::p_load(rgee, mapview, sf, cptcity, raster, rgdal, tidyverse, rgeos,
               gtools, ggpubr, stringr, cowplot, foreign, ggmap, stars,
               extrafont, ggspatial)

## Consideraciones de la elaboración de este ejercicio
## Sisteama operativo: Windows 10
## R: 4.0.2
## rgee: 1.00
## mapview: 2.8.8
## sf: 0.9.5
## cptcity: 1.0.4
## raster: 3.3.13
## rgdal: 1.5.12
## tidyverse: 1.3.0
## rgeos: 0.5.3
## gtools: 3.8.2
## ggpubr: 0.4.0
## stringr: 1.4.0
## cowplot: 1.0.0
## foreign: 0.8.80
## ggmap: 3.0.0.902
## stars: 0.4.3
## extrafont: 0.17
## ggspatial: 1.1.4

Uso de rgee

Inicializar la sesión en Google Earth Engine en cada sesión.

## Iniciar la sesión de GEE
ee_Initialize()

## -- rgee 1.0.0 ------------------------------------------------- earthengine-api 0.1.226 -- 
##  v email: not_defined
##  v Initializing Google Earth Engine:
 v Initializing Google Earth Engine:  DONE!
## 
 v Earth Engine user: users/erikseras 
## ------------------------------------------------------------------------------------------

Crear las variables necesarias para el presente ejemplo.

## Área de interés
box <- ee$Geometry$Rectangle(coords = c(-76.97,-12.17,-76.90,-12.10),
                             ## WGS 84
                             proj = "EPSG:4326",
                             geodesic = FALSE)

## Escoger una imagen
image <- ee$Image("COPERNICUS/S2_SR/20191028T151711_20191028T152253_T18LTM")

## NDVI
ndvi <- image$normalizedDifference(c("B8", "B4"))

Preparación de las archivos a trabajar

Como se piensa trabajar sobre el área de Lomas de Villa María (Lima), se necesitará un archivo shape que contenga su contorno. En el siguiente link se podrá encontrar los archivos de las Áreas de Conservación Regional (ACR) https://drive.google.com/drive/folders/1yDB9oZBS6ZSZ-U1Rf7WsNlLdaHdHl-3z?usp=sharing.

## Shaoe del área de interés
loma_shape <- 
        ## Cambiar la ubicación  en donde usted guardó su archivo shapefile
        raster::shapefile("../data/ACR/ACR.shp") %>% 
        sf::st_as_sf(loma_shape) %>% 
        ## Obtener el shape de Lomas de Villa María
        dplyr::filter(acr_codi == "ACR23" & objectid == 2924) %>% 
        sf::as_Spatial()

## Raster de NDVI del área de interés
raster_ndvi <- ee_as_raster(
        image = ndvi,
        region = box,
        scale = 10
) %>%  
        raster::projectRaster(crs = '+proj=longlat +datum=WGS84 +no_defs') %>% 
        raster::crop(loma_shape) %>% 
        raster::mask(loma_shape)

## - region parameters
##  WKT      : POLYGON ((285644.2 8653871, 293263.5 8653925, 293209.5 8661670, 285588.2 8661616, 285644.2 8653871)) 
##  CRS      : 32718 
##  geodesic : FALSE 
##  evenOdd  : TRUE 
## region is too large ... creating  4  patches.
## Getting data from the patch: 1/4 
## Getting data from the patch: 2/4 
## Getting data from the patch: 3/4 
## Getting data from the patch: 4/4

Este gráfico muestra de manera sencilla los valores de NDVI de Lomas de Villa María.

## NDVI de Lomas de Villa María
plot(raster_ndvi)

Elaboración del mapa

Con el siguiente código se arreglará el gráfico anterior y se le dará un formato aceptable para su presentación.

## Elaboración del mapa final
mapa_final <- 
        ## Ploteo del raster
        rasterVis::gplot(raster_ndvi) +
        geom_tile(aes(fill = value))  +
        ## Establecer colores de los pixel
        scale_fill_gradientn(
                colours = RColorBrewer::brewer.pal(n = 8, name = "BuGn"),
                na.value = 'white'
        ) +
        ## Añadir el shape
        geom_polygon(
                data = loma_shape, aes(x=long, y = lat, group = group),
                     color = 'grey', fill='NA'
        ) +
        ## Arreglar el sistema de coordenadas del gráfico
        coord_sf(crs = crs(loma_shape)) +
        ## Título y subtítulo
        ggtitle(
                "NDVI de Lomas de Villa María (28/10/2019)",
                subtitle = "Procesado por: Erik Seras"
        ) +
        ## Añadir etiquetas de ejes
        labs(
                x="", y="", fill = "NDVI",
             caption ="Fuente: Sentinel-2 MSI: MultiSpectral Instrument, Level-2A"
        ) +
        ## Tema del gráfico
        theme_bw() +
        ## Personalizar el tema del gráfico
        theme(
                ## Establecer los colores del fonodo del gráfico y su contorno
                panel.background = element_rect(fill = 'white', colour = 'black'),
                ## No mostrar las grillas mayores
                panel.grid.major = element_blank(),
                ## No mostrar las grillas menores
                panel.grid.minor = element_blank(),
                ## Ubicar a la derecha la legenda
                legend.position = 'right',
                ## Ubicar al medio el título
                plot.title = element_text(hjust = 0.5),
                ## ubicar en el medio el título
                plot.subtitle = element_text(hjust = 0.5),
                ## Ajustar el centrado de las coordenadas de latitud del eje y
                axis.text.y = element_text(vjust = 0.5),
                ## Ajustar el centrado de las coordenadas de latitud del eje x
                axis.text.x = element_text(vjust = 0.5),
                ## Ancho de la barra de la leyenda
                legend.key.width = unit(1.5, 'line')
        )  +
        ## Añadir barra de escala
        annotation_scale(
                location = "br"
        ) +
        ## Añadir el norte
        annotation_north_arrow(
                location = "tr", which_north = "true",
                               style = north_arrow_fancy_orienteering, 
                               height = unit(1.2, "cm"),
                               width = unit(1.2, "cm")
        ) +
        ## Añadir texto extra
        annotate(
                geom = 'text', x = -76.935, y = -12.155,
                label = 'COPERNICUS/S2_SR/20191028T151711_20191028T152253_T18LTM',
                color = 'grey22', size = 3
        )

## Regions defined for each Polygons

Y por último, guardar el gráfico que se elaboró en el formato que desee, en este caso se escogió el formato png.

## Guardar el mapa
ggplot2::ggsave(
        ## Cambiar la direeción a donde usted quiera guardar el gráfico
        filename =  '../images/mapa_ndvi.png',
        plot = mapa_final,
        units = "cm", width = 20, 
        height = 20, dpi = 900
)

Este sería el mapa generado luego de haber corrido todo el código. Usted puede modificarlo y adecuarlo a sus necesidades; y con su imaginación, representar una infinidad de cosas. Espero que le haya sido útil este breve código. Muchas gracias por leerlo.

Ejemplo de elaboración de un mapa usando GEE

Uso de Google Earth Engine para elaboración de mapas

Erik Seras (er.seras.perez@gmail.com)

26/7/2020

Cargar libreías

Uso de rgee

Preparación de las archivos a trabajar

Elaboración del mapa