R

English version

R

El paquete de R monitoreo llama, con reticulate, a las funciones del paquete de Python: calcula exactamente lo mismo que el plugin de QGIS y el servidor MCP. Devuelve data.frame, rutas de GeoTIFF para terra y, con rgee, imágenes de Earth Engine.

Las opciones (regiones, periodos, métricas, filtros) se explican en Conceptos. Los ejemplos de esta página son reales y usan el Parque Nacional de Sierra Nevada (ID de WDPA 389011).

Instalación

1. El Python de MonitorEO. Lo más sencillo es el entorno de conda de Primeros pasos (y earthengine authenticate, solo para usar tu cuenta de Earth Engine).

2. El paquete de R, desde GitHub:

   install.packages(c("remotes", "reticulate", "ggplot2", "sf", "terra", "leaflet"))
   remotes::install_github("monitoreo-obsnev/monitoreo", subdir = "r/monitoreo")

   Para una versión concreta (la que se cita con su DOI), añade ref = "v0.3.2"; desde un clon del repositorio, install.packages("r/monitoreo", repos = NULL, type = "source").

3. Decirle a R qué Python usar, antes de cargar el paquete o en ~/.Renviron:

   Sys.setenv(RETICULATE_PYTHON = "/Users/ana/miniforge3/envs/monitoreo/bin/python")
   # Windows: "C:/Users/ana/miniforge3/envs/monitoreo/python.exe"

   Si ya usas rgee con su propio Python, puedes instalar en él el paquete de Python con mo_install().

1. Conectar

library(monitoreo)
mo_init()                        # sin la variable EE_PROJECT: la cuenta compartida de MonitorEO
# mo_init("mi-proyecto-ee")      # o tu cuenta y tu proyecto (o mo_init() con EE_PROJECT); mo_init(""): la compartida
mo_eecu()                       # tus EECU de este mes; con tu proyecto, mo_eecu("Contributor") si es su nivel

Con la cuenta compartida son tus research EECUs (25 EECU-h al mes, solo para investigación): si te queda menos del 10 %, mo_init() y las funciones que calculan dan un warning() (una vez al mes) que pide pasar a tu cuenta, y cuando se acaban la cuenta compartida no calcula hasta el mes siguiente (Conceptos).

Al cargarse, el paquete muestra la descripción y los agradecimientos de MonitorEO; mo_about("es") los da en español (por defecto, en inglés).

2. Variables y regiones

mo_variables(sensor = "Modis", text = "NDVI")[, c("ref", "variable", "sensor", "collection", "resolution_m")]

  ref variable sensor        collection resolution_m
1  25     NDVI  Modis MODIS/061/MOD13Q1          250
2  26     NDVI  Modis MODIS/061/MYD13Q1          250

mo_protected_areas("Doñana")[, c("SITE_ID", "NAME", "DESIG", "IUCN_CAT")]

    SITE_ID                 NAME                                          DESIG IUCN_CAT
1 555588792               Doñana Zona de Importancia Comunitaria ZIC (ZEPA/ZEC)       IV
2 555595662 Doñana Norte y Oeste Zona de Importancia Comunitaria ZIC (ZEPA/ZEC)       IV
3    389012               Doñana                                Parque Nacional       II
4    349268               Doñana                                 Parque Natural        V
...

Las regiones se construyen con funciones:

sn <- mo_wdpa(389011)                          # una o varias áreas protegidas: mo_wdpa(c(389011, 389012))
mo_ecoregion(793)                              # ecorregión por ID, o mo_ecoregion(lon = -3.3, lat = 37.1)
mo_bbox(-3.6, 36.9, -2.6, 37.3)                # oeste, sur, este, norte
mo_point(-3.42, 37.13, buffer = 1500)          # círculo de 1,5 km
zonas <- mo_sf(sf::st_read("zonas.gpkg"), "nombre")   # cualquier objeto sf de polígonos; una serie por polígono

3. Serie temporal

ts <- mo_time_series(sn, "NDVI", "Sentinel-2", "2022-01-01", "2023-12-31", step = "monthly")
head(ts, 3)

        date n_images n_pixels        region     value
1 2022-01-01       30   346520 Sierra Nevada 0.3536522
2 2022-02-01       29   298880 Sierra Nevada 0.4035685
3 2022-03-01       31   224912 Sierra Nevada 0.2237470

Es un data.frame con fechas Date y los atributos title, unit y scale (el píxel usado, aquí 50 m).

mo_plot(ts)                                    # ggplot2: se puede seguir añadiendo capas
write.csv(ts, "ndvi.csv", row.names = FALSE)   # o mo_time_series(..., csv = "ndvi.csv")

Más ejemplos:

mo_time_series(sn, "LST", "Landsat", "2015-01-01", "2023-12-31", step = "annual", metric = "max")
mo_time_series(zonas, "NDVI", "Sentinel-2", "2023-01-01", "2023-12-31", season = c("06-01", "09-30"))
mo_time_series(sn, "NDSI", "Modis", "2022-10-01", "2023-06-30", step = "16-day", quality = list(snow = FALSE))

4. Mapas

Teselas de Earth Engine en leaflet (sobre una ortoimagen):

capa <- mo_map_layer(sn, "NDVI", "Sentinel-2", "2023-06-01", "2023-08-31")
mo_leaflet(capa)

Con rgee, la imagen de Earth Engine y su rampa de colores:

x <- mo_image(sn, "LST", "Modis", "2023-06-01", "2023-08-31", metric = "mean")
rgee::Map$centerObject(x$image, 10)
rgee::Map$addLayer(x$image, x$vis, x$name)

Estadísticas de un mapa en cada polígono:

mo_region_stats(sn, "LST", "Modis", "2023-06-01", "2023-08-31", metric = "mean")
#          region     mean      std      min       p5      p50      p95      max n_pixels
# 1 Sierra Nevada 32.20741 1.976974 27.26778 29.40603 32.08767 34.90077 40.85588     1051

En GeoTIFF, para terra (regiones pequeñas o medianas; las grandes, con mo_export_image() a Google Drive):

f <- mo_download_image(mo_point(-3.42, 37.13, 1500), "NDVI", "Sentinel-2", "2023-01-01", "2023-12-31",
                       path = "ndvi_p95.tif", metric = "p95")
terra::plot(terra::rast(f$file))

5. Tendencias

La tendencia del NDVI de cada mes de 2005 a 2024 (MODIS):

tr <- mo_trend_table(sn, "NDVI", "Modis", "2005-01-01", "2024-12-31", step = "monthly")
tr[, c("label", "n", "tau", "p_value", "significant", "slope")]

       label  n         tau     p_value significant         slope
1    January 20 -0.06315789 0.697031003       FALSE -0.0002659812
2   February 20  0.07368421 0.649670504       FALSE  0.0006624354
3      March 20  0.25263158 0.119393190       FALSE  0.0014753512
4      April 20  0.04210526 0.795207612       FALSE  0.0003653394
5        May 20  0.10526316 0.516412268       FALSE  0.0008568967
6       June 20  0.18947368 0.242809091       FALSE  0.0011871008
7       July 20  0.17894737 0.269981282       FALSE  0.0013375445
8     August 20  0.45263158 0.005267412        TRUE  0.0018529123
9  September 20  0.36842105 0.023140931        TRUE  0.0019901468
10   October 20  0.13684211 0.398920137       FALSE  0.0010093510
11  November 20  0.12631579 0.436178362       FALSE  0.0011457025
12  December 20  0.10526316 0.516412268       FALSE  0.0004527001

Solo agosto y septiembre tienen una tendencia significativa (p < 0,05): unas +0,002 unidades de NDVI por año. attr(tr, "yearly") tiene los valores de cada año y mo_trend_table(..., chart = "tendencias.png") dibuja la gráfica de cada mes con su recta de Sen.

En cada píxel, en GeoTIFF con las bandas slope, tau, z, p_value, intercept y n:

f <- mo_download_trend(sn, "NDVI", "Modis", "2005-01-01", "2024-12-31", path = "tendencia.tif")   # step = "annual"
r <- terra::rast(f$file)
r

class       : SpatRaster
size        : 121, 315, 6  (nrow, ncol, nlyr)
resolution  : 250, 250  (x, y)
coord. ref. : WGS 84 / UTM zone 30N (EPSG:32630)
names       :   slope, tau,         z, p_value, intercept,  n

pendiente <- terra::ifel(r$p_value < 0.05, r$slope, NA)      # solo la tendencia significativa
terra::plot(pendiente, col = hcl.colors(11, "RdYlGn"), range = c(-0.005, 0.005))

Para un mes: step = "monthly", part = 8 (agosto). El mapa de teselas: mo_leaflet(mo_trend_map(...)).

6. Comparar dos periodos

cmp <- mo_compare_series(sn, "NDVI", "Modis", "2005-01-01", "2009-12-31",    # referencia
                         "2020-01-01", "2024-12-31")                         # otro periodo
head(cmp[, c("label", "reference", "other", "difference")], 3)

     label reference     other    difference
1  January 0.4954321 0.4952818 -0.0001502843
2 February 0.4680880 0.4829461  0.0148581591
3    March 0.4335899 0.4606720  0.0270820311

El mapa de la diferencia: mo_compare_map() (teselas) o mo_download_comparison() (GeoTIFF con las bandas difference, reference y other).

7. Severidad de un incendio

s <- mo_burn_severity(mo_wdpa(389013), "2023-05-15", "2023-08-14", "2023-09-15", "2023-11-30",
                      path = "severidad.tif")          # Sentinel-2; sensor = "Landsat" para incendios anteriores
s[s$class >= 1, c("label", "ha")]                     # hectáreas de cada clase de severidad
terra::plot(terra::rast(attr(s, "file"))$severity)
r <- mo_burn_severity(mo_wdpa(389013), "2023-05-15", "2023-08-14", "2023-09-15", "2023-11-30",
                      index = "RBR", path = "rbr.tif")  # el RBR: media, mediana, p10, p90, mín., máx. y ha

8. Exportaciones y recetas

tarea <- mo_export_image(sn, "NDVI", "Sentinel-2", "2023-01-01", "2023-12-31", resolution = 10, folder = "MonitorEO")
mo_task_status(tarea$operation)

mo_download_variables("receta.yaml", "zona.gpkg", "variables", dry_run = TRUE)   # receta: ver la página Python
mo_download_variables("receta.yaml", "zona.gpkg", "variables")
terra::rast(list.files("variables", "\\.tif$", full.names = TRUE))

La receta puede ser un YAML o una lista de R (yaml::read_yaml()); su formato está en Python.

Funciones

Función	Qué hace
mo_init(), mo_install(), mo_about()	Conectar, instalar el paquete de Python, descripción y agradecimientos
mo_eecu()	Tus EECU de Earth Engine gastados este mes (research EECUs con la cuenta compartida)
mo_variables(), mo_describe()	Variables del catálogo
mo_protected_areas()	Buscar áreas protegidas de WDPA
mo_wdpa(), mo_ecoregion(), mo_bbox(), mo_point(), mo_geojson(), mo_sf()	Regiones
mo_time_series(), mo_plot()	Serie temporal y su gráfica
mo_map_layer(), mo_leaflet(), mo_image()	Mapas: teselas, leaflet y rgee
mo_region_stats()	Estadísticas de un mapa por polígono
mo_download_image(), mo_download_variables()	GeoTIFF en el ordenador
mo_trend_table(), mo_trend_map(), mo_download_trend()	Tendencias estacionales
mo_compare_series(), mo_compare_map(), mo_download_comparison()	Comparación de periodos
mo_burn_severity()	Severidad de un incendio (dNBR, clases del USGS; o el RBR con index = "RBR")
mo_export_image(), mo_export_time_series(), mo_task_status()	Exportaciones a Google Drive

Los argumentos son los de las funciones de Python del mismo nombre.