En este documento se ejemplifica brevemente el uso del paquete rgee para calcular el Indice de Vegetación mejorado (EVI) para la ciudad de Asunción (Paraguay) durante el periodo comprendido entre el año 2015 y 2019 como parte de un proyecto que busca replicar la propuesta de Tuanmu y Jetz (2015) quienes caracterizaron la heterogeneidad del hábitat para la biodiversidad a escala global por medio de 14 métricas de paisaje a varias resoluciones espaciales y las cuales estan basadas en características texturales de imágenes del EVI, las cuales forman parte de los productos derivados del Moderate Resolution Imaging Spectroradiometer (MODIS).
El objetivo de esto es calcular las mismas métricas pero a partir de imágenes del satélite Landsat 8, las cuales cuentan con una mejor resolución espacial que MODIS. El primer paso para esto es el cálculo del EVI y luego utilizar los resultados para obtener las 14 métricas propuestas. Este script fue elaborado primero en Google earth Engine para probarlo y usarlo como guía al transcribir los parámetros y funciones al entorno de R.
En primer lugar, se deben instalar todas las librerías necesarias. Tidyverse es una de las más utilizadas y cuenta con la mayoría de los paquetes que pueden ser de utilidad para correr este script o trabajar con los datos. También las librerías “sf” y “raster” pueden ser de utilidad. Por último, se debe instalar el paquete “rgee” ya que servirá de nexo entre la API de GEE (en Python), GEE (js) y R.
Para comenzar primero cargamos las librerías necesarias para los cálculos.
Cabe destacar que es necesario contar con una cuenta en GEE para que el paquete rgee funcione. Si contamos con rgee, Python (también Anaconda3 o miniconda) instalados correctamente y una cuenta en GEE, el proceso es muy sencillo. rgee cuenta con instrucciones bastante sencillas de seguir en el repositorio de github pero algunas personas pueden tener problemas conectando python con R.
El paquete cuenta con la función “ee_install()” para instalar y configurar todo lo necesario. En el caso que no encuentre a Python, sugerirá la instalación de miniconda.
Dicho esto, iniciamos rgee y nos conectamos a GEE:
Definimos las variables: limites del área de interés (AOI), periodo de estudio, porcentaje máximo de nubosidad permitida en las imágenes y los parámetros de visualización para los resultados:
#cargar AOI
mask <- st_read("D:/DATA/R/Rprojects/EVI_calculation/shape/ASUNCION.SHP",
quiet = TRUE) %>%
sf_as_ee()
#Convertimos el shapefile en una geometría
region <- mask$geometry()
#Parámetros de visualización
visP = list(
min = 0,
max = 1,
palette = c("FFFFFF", "CE7E45", "DF923D", "F1B555", "FCD163", "99B718",
"74A901", "66A000", "529400", "3E8601", "207401", "056201",
"004C00", "023B01", "012E01", "011D01", "011301")
)
#Fechas para las imágenes de base
startDate <- "2015-01-01"
endDate <- "2019-12-31"
#Cobertura de nubes
cc <- 60L
Cargamos la colección de imágenes Landsat 8 a utilizar junto con los datos de JRC. La colección de L8 utilizada cuenta con la reflectancia superficial de los sensores OLI/TIRS corregida atmosféricamente. Por otro lado, los datos del JRC Global Surface Water Mapping Layers (v1.2) contiene mapas de la ubicación y distribución temporal de las aguas superficiales desde 1984 hasta 2019, además proporciona estadísticas sobre la extensión y el cambio de esas superficies de agua; estos datos usaremos como máscara para los cuerpos de agua permanentes.
Definimos las funciones a utilizar: máscara de nubes, máscara para zonas de agua permanente, máscara para cortar las imágenes por el AOI, cálculo de EVI2 para la colección completa, máscara para valores de EVI2 > 1 < 0 (outliers)
#Mascara de nubes
maskClouds <- function(img) {
clear <- img$select('pixel_qa')$bitwiseAnd(2)$neq(0)
clear <- img$updateMask(clear)
return(clear)
}
#Máscara de agua
water_mask <- gsw$select('seasonality')$lt(11)$unmask(1)$clip(region)
#l8 evi
AddEVI <- function (img) {
evi <- img$expression(
expression = '2.5 * ((nir - red) / (nir + 2.4 * red + 1))', opt_map = list(
nir = img$select("B5")$multiply(0.0001),
red = img$select("B4")$multiply(0.0001)))
return (img$addBands(evi))
}
#Clip por AOI
cliplim <- function (img) {return (img$clip(region))}
#Función que enmascara valores de EVI definidos como > 1 and < 0
maskExcess <- function(img) {
hi <- img$lte(1)
lo <- img$gte(0)
masked <- img$mask(hi$and(lo))
return (img$mask(masked))
}
Aplicamos filtros espaciales, temporales y de porcentaje de nubosidad máximo a la colección de imágenes L8, luego aplicamos las funciones para la máscara de nubes y además la función para cortar todas las escenas de la colección por los límites del AOI.
Aplicamos la función para el cálculo del EVI a todas las imágenes de la colección. Una vez terminado reducimos la imágen al percentil 90 y finalmente aplicamos una máscara para eliminar todos los pixeles que correspondan a aguas permanentes.
La extracción del percentil 90 de EVI se realiza para capturar el “verdor” de la superficie terrestre en el pico de una temporada de crecimiento, evitando al mismo tiempo la inclusión de valores de EVI demasiado altos.
La selección de la serie temporal de 5 años se realizó para equilibrar los potenciales cambios de cobertura que puedan darse y un número suficiente de valores de EVI libre de nubes durante el pico de crecimiento.
También el uso de una máscara para las aguas permanentes se realiza debido a los valores artificialmente altos de EVI que suelen tener (Solano et al. 2010).
#Aplicar función para EVI
col1 <- coleccion$map(AddEVI)
#Aplicar función para percentil y función de máscara de agua
col1 <- col1$reduce(ee$Reducer$percentile(list(90)))$mask(water_mask)
Por último, corroboramos los nombres asignados a las bandas de la imagen resultante y anotamos el que corresponde al EVI ya que posteriormente lo utili- zaremos (también podemos renombrar la banda modificando la función para el cálculo de EVI antes de aplicarla a la colección)
#Inspeccionar el nombre de las bandas
bandNames <- col1$bandNames()
#ee_help(bandNames)
#cat("Bandas disponibles: ", paste(bandNames$getInfo(), collapse = " "))
Separamos la banda que corresponde al percentil 90 del EVI calculado y guardamos en una nueva variable:
Por último, centramos la visualización sobre los límites del AOI, elegimos el porcentaje de transparencia de la capa, le damos el nombre que queremos y visualizamos el resultado para controlar que todo este en orden:
El resultado interactivo puede observarse haciendo clic acá.
Consultando la documentación de “rgee” es posible escribir una función para exportar los resultados a Google Drive o a los Assets de GEE. Este paquete resulta en una poderosa herramienta para análisis espaciales utilizando las bondades de GEE acortando la brecha entre los diferentes lenguajes utilizados, y permite a muchos usuarios de R utilizar GEE sin la necesitad de aprender todo un lenguaje nuevo al momento de analizar sus datos.
Aybar, Cesar. 2020. Rgee: R Bindings for Calling the ’Earth Engine’ Api.
Aybar, Cesar, Quisheng Wu, Lesly Bautista, Roy Yali, and Antony Barja. 2020. “Rgee: An R Package for Interacting with Google Earth Engine.” Journal of Open Source Software.
Bivand, Roger S., Edzer Pebesma, and Virgilio Gomez-Rubio. 2013. Applied Spatial Data Analysis with R, Second Edition. Springer, NY. https://asdar-book.org/.
Didan, Kamel, Armando Barreto, Ramón Solano, and Alfredo Huete. 2015. “MODIS Vegetation Index User’s Guide (MOD13 Series).” Tucson, Arizona: Vegetation Index; Phenology Lab, University of Arizona. https://vip.arizona.edu/documents/MODIS/MODIS_VI_UsersGuide_June_2015_C6.pdf.
Henry, Lionel, and Hadley Wickham. 2020. Purrr: Functional Programming Tools. https://CRAN.R-project.org/package=purrr.
Hijmans, Robert J. 2020. Raster: Geographic Data Analysis and Modeling. https://CRAN.R-project.org/package=raster.
Müller, Kirill, and Hadley Wickham. 2020. Tibble: Simple Data Frames. https://CRAN.R-project.org/package=tibble.
Pebesma, Edzer. 2018. “Simple Features for R: Standardized Support for Spatial Vector Data.” The R Journal 10 (1): 439–46. https://doi.org/10.32614/RJ-2018-009.
———. 2020. Sf: Simple Features for R.
Pebesma, Edzer, and Roger Bivand. 2020. Sp: Classes and Methods for Spatial Data. https://CRAN.R-project.org/package=sp.
Pebesma, Edzer J., and Roger S. Bivand. 2005. “Classes and Methods for Spatial Data in R.” R News 5 (2): 9–13. https://CRAN.R-project.org/doc/Rnews/.
R Core Team. 2020. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. https://www.R-project.org/.
Tuanmu, Mao-Ning, and Walter Jetz. 2015. “A Global, Remote Sensing-Based Characterization of Terrestrial Habitat Heterogeneity for Biodiversity and Ecosystem Modelling: Global Habitat Heterogeneity.” Global Ecology and Biogeography 24 (11): 1329–39. https://doi.org/10.1111/geb.12365.
Wickham, Hadley. 2016. Ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York. https://ggplot2.tidyverse.org.
———. 2019a. Stringr: Simple, Consistent Wrappers for Common String Operations. https://CRAN.R-project.org/package=stringr.
———. 2019b. Tidyverse: Easily Install and Load the ’Tidyverse’. https://CRAN.R-project.org/package=tidyverse.
———. 2020. Forcats: Tools for Working with Categorical Variables (Factors). https://CRAN.R-project.org/package=forcats.
Wickham, Hadley, Mara Averick, Jennifer Bryan, Winston Chang, Lucy D’Agostino McGowan, Romain François, Garrett Grolemund, et al. 2019. “Welcome to the tidyverse.” Journal of Open Source Software 4 (43): 1686. https://doi.org/10.21105/joss.01686.
Wickham, Hadley, Winston Chang, Lionel Henry, Thomas Lin Pedersen, Kohske Takahashi, Claus Wilke, Kara Woo, Hiroaki Yutani, and Dewey Dunnington. 2020. Ggplot2: Create Elegant Data Visualisations Using the Grammar of Graphics. https://CRAN.R-project.org/package=ggplot2.
Wickham, Hadley, Romain François, Lionel Henry, and Kirill Müller. 2020. Dplyr: A Grammar of Data Manipulation. https://CRAN.R-project.org/package=dplyr.
Wickham, Hadley, and Lionel Henry. 2020. Tidyr: Tidy Messy Data. https://CRAN.R-project.org/package=tidyr.
Wickham, Hadley, Jim Hester, and Romain Francois. 2018. Readr: Read Rectangular Text Data. https://CRAN.R-project.org/package=readr.