R y GEE
El paquete rgee
VE-MS-CC
DataIntelligence07-04-2021
Abstract
Éste informe, aún en construcción, se basa en la página https://csaybar.github.io/rgee-examples/ que contiene un poco más de 250 ejemplos y revisa las diferentes funcionalidades que ofrece el paquete rgee para desarrollar análisis sobre las imágenes satelitales de Google Earth Engine. La intención es explorar un poco el código R para identificar las similitudes con el JavaScript nativo de GEE y así adaptarlo a nuestros propios requerimientos. El trabajo se divide en 5 partes. La primera entrega una introducción a las características básicas, la segunda, ejercicios sobre algunos modelos de aprendizaje automático (CART, SVM y KMeans) para pasar luego a revisar el tema de imágenes. La cuarta parte aborda el tema de las colecciones (ImageCollection) y una quinta, Geometry, Feature y FeatureCollection.
Rgee es un híbrido entre Python y R. Su ventaja es el dar la posibilidad de utilizar los estilos python o R a conveniencia.
Catalogo global de imagenes landsat:
https://developers.google.com/earth-engine/datasets/catalog
Índice
1. La información contenida dentro de las imágenes landsat
1.1 Estructura de la colección Landsat:
El la encuesta geológica de los EE.UU (USGS) produce datos en 3 niveles (categorías) para cada satélite:
- Nivel 1 (T1): datos que cumplen con los requisitos de calidad geométricos y radiométricos
- Nivel 2 (T2): datos que no cumplen con los requisitos del Nivel 1
- Tiempo real (RT): datos que aún no se han evaluado (se necesitan hasta un mes).
LANDSAT/LC08/C01/T1_RT Landsat 8, Colección 1, Nivel 1 + Tiempo real
LANDSAT/LC08/C01/T1 Landsat 8, Colección 1, Nivel 1 únicamente
LANDSAT/LC08/C01/T2 Landsat 8, Colección 1, Nivel 2 únicamente.
Vamos a analizar la información extraíble de una imagen LANDSAT/LC08/C01/T1.
1.1.1 Descripción general de los datos de una imagen.
Despleguemos las primeras 30 columnas de la base de datos:
a_001 <- ee$Image("LANDSAT/LC08/C01/T1/LC08_044034_20140403") %>% ee$Image$getInfo()
write.csv2(a_001, "a_001.csv")
# saveRDS(g,"g.rds")
a_001 <- as.data.frame(a_001)
head(colnames(a_001),30)## [1] "type" "bands.id"
## [3] "bands.data_type.type" "bands.data_type.precision"
## [5] "bands.data_type.min" "bands.data_type.max"
## [7] "bands.dimensions" "bands.crs"
## [9] "bands.crs_transform" "bands.id.1"
## [11] "bands.data_type.type.1" "bands.data_type.precision.1"
## [13] "bands.data_type.min.1" "bands.data_type.max.1"
## [15] "bands.dimensions.1" "bands.crs.1"
## [17] "bands.crs_transform.1" "bands.id.2"
## [19] "bands.data_type.type.2" "bands.data_type.precision.2"
## [21] "bands.data_type.min.2" "bands.data_type.max.2"
## [23] "bands.dimensions.2" "bands.crs.2"
## [25] "bands.crs_transform.2" "bands.id.3"
## [27] "bands.data_type.type.3" "bands.data_type.precision.3"
## [29] "bands.data_type.min.3" "bands.data_type.max.3"Que se corresponden con las propiedades de las imagenes contenidas en: USGS Landsat 8 Collection 1 Tier 1 Raw Scenes https://developers.google.com/earth-engine/datasets/catalog/LANDSAT_LC08_C01_T1#image-properties
1.1.2 Los datos asociados a una colección de imágenes filtrada: información general
Necesitamos el conjunto de escenas landsat que cubran Chile
El Path y el Row de Landsat. Son parámetros numéricos que permiten identificar una imagen satelital Landsat de forma análoga a los valores de longitud y latitud. Landsat mapea de manera constante las mismas superficies a través de una malla de cuadrículas. Cada cuadrícula presenta de manera constante el mismo Path y Row, por lo que, buscando por estos parámetros, se consiguen todas las imágenes históricas de Landsat para el mismo lugar.
Ahora, cuales son el Path y el Row para Santiago? Encontramos un convertidor automatico en linea: https://landsat.usgs.gov/landsat_acq#convertPathRow
# Realiza una llamada a alguna de las colecciones de Landsat
object <- ee$ImageCollection("LANDSAT/LC08/C01/T1")$
filter(ee$Filter()$eq("WRS_PATH", 233))$
filter(ee$Filter()$eq("WRS_ROW", 83))$
filterMetadata ('CLOUD_COVER', 'Less_Than', 20)$
filterDate("2020-03-01", "2020-04-01")$aside(ee_print)## ------------------------------------------------ Earth Engine ImageCollection --
## ImageCollection Metadata:
## - Class : ee$ImageCollection
## - Number of Images : 2
## - Number of Properties : 38
## - Number of Pixels* : 1415544744
## - Approximate size* : 28.86 GB
## Image Metadata (img_index = 0):
## - ID : LANDSAT/LC08/C01/T1/LC08_233083_20200314
## - system:time_start : 2020-03-14 14:33:32
## - system:time_end : 2020-03-14 14:33:32
## - Number of Bands : 12
## - Bands names : B1 B2 B3 B4 B5 B6 B7 B8 B9 B10 B11 BQA
## - Number of Properties : 119
## - Number of Pixels* : 707772372
## - Approximate size* : 14.43 GB
## Band Metadata (img_band = 'B1'):
## - EPSG (SRID) : WGS 84 / UTM zone 19N (EPSG:32619)
## - proj4string : +proj=utm +zone=19 +datum=WGS84 +units=m +no_defs
## - Geotransform : 30 0 227985 0 -30 -3556485
## - Nominal scale (meters) : 30
## - Dimensions : 7643 7717
## - Number of Pixels : 58981031
## - Data type : INT
## - Approximate size : 1.20 GB
## --------------------------------------------------------------------------------
## NOTE: (*) Properties calculated considering a constant geotransform and data type.# 1.1.2 Los datos asociados a una colección de imágenes filtrada: información de las bandas.
La tabla rescatada de una imagen posee una serie de 11 conjuntos de datos de 7 columnas cada una que contiene información relativa a las bandas. Son 11 bandas.
library(kableExtra)## Warning: package 'kableExtra' was built under R version 4.0.3addNDVI <- function(image) {
return(image$addBands(image$normalizedDifference(c("B5", "B4"))))
}
ndviCollection <- object$map(addNDVI)
first <- ndviCollection$first()
a_001 <- first$getInfo()
a_001 <- as.data.frame(a_001)
kbl(a_001[,c(2:9)]) %>%
kable_styling(bootstrap_options = c("striped", "hover")) %>%
kable_paper() %>%
scroll_box(width = "100%", height = "300px")| bands.id | bands.data_type.type | bands.data_type.precision | bands.data_type.min | bands.data_type.max | bands.dimensions | bands.crs | bands.crs_transform |
|---|---|---|---|---|---|---|---|
| B1 | PixelType | int | 0 | 65535 | 7671 | EPSG:32619 | 30 |
| B1 | PixelType | int | 0 | 65535 | 7731 | EPSG:32619 | 0 |
| B1 | PixelType | int | 0 | 65535 | 7671 | EPSG:32619 | 227985 |
| B1 | PixelType | int | 0 | 65535 | 7731 | EPSG:32619 | 0 |
| B1 | PixelType | int | 0 | 65535 | 7671 | EPSG:32619 | -30 |
| B1 | PixelType | int | 0 | 65535 | 7731 | EPSG:32619 | -3556485 |
kbl(a_001[,c(10:17)])%>%
kable_styling(bootstrap_options = c("striped", "hover")) %>%
kable_paper() %>%
scroll_box(width = "100%", height = "300px")| bands.id.1 | bands.data_type.type.1 | bands.data_type.precision.1 | bands.data_type.min.1 | bands.data_type.max.1 | bands.dimensions.1 | bands.crs.1 | bands.crs_transform.1 |
|---|---|---|---|---|---|---|---|
| B2 | PixelType | int | 0 | 65535 | 7671 | EPSG:32619 | 30 |
| B2 | PixelType | int | 0 | 65535 | 7731 | EPSG:32619 | 0 |
| B2 | PixelType | int | 0 | 65535 | 7671 | EPSG:32619 | 227985 |
| B2 | PixelType | int | 0 | 65535 | 7731 | EPSG:32619 | 0 |
| B2 | PixelType | int | 0 | 65535 | 7671 | EPSG:32619 | -30 |
| B2 | PixelType | int | 0 | 65535 | 7731 | EPSG:32619 | -3556485 |
kbl(a_001[,c(18:25)])%>%
kable_styling(bootstrap_options = c("striped", "hover")) %>%
kable_paper() %>%
scroll_box(width = "100%", height = "300px")| bands.id.2 | bands.data_type.type.2 | bands.data_type.precision.2 | bands.data_type.min.2 | bands.data_type.max.2 | bands.dimensions.2 | bands.crs.2 | bands.crs_transform.2 |
|---|---|---|---|---|---|---|---|
| B3 | PixelType | int | 0 | 65535 | 7671 | EPSG:32619 | 30 |
| B3 | PixelType | int | 0 | 65535 | 7731 | EPSG:32619 | 0 |
| B3 | PixelType | int | 0 | 65535 | 7671 | EPSG:32619 | 227985 |
| B3 | PixelType | int | 0 | 65535 | 7731 | EPSG:32619 | 0 |
| B3 | PixelType | int | 0 | 65535 | 7671 | EPSG:32619 | -30 |
| B3 | PixelType | int | 0 | 65535 | 7731 | EPSG:32619 | -3556485 |
kbl(a_001[,c(26:33)])%>%
kable_styling(bootstrap_options = c("striped", "hover")) %>%
kable_paper() %>%
scroll_box(width = "100%", height = "300px")| bands.id.3 | bands.data_type.type.3 | bands.data_type.precision.3 | bands.data_type.min.3 | bands.data_type.max.3 | bands.dimensions.3 | bands.crs.3 | bands.crs_transform.3 |
|---|---|---|---|---|---|---|---|
| B4 | PixelType | int | 0 | 65535 | 7671 | EPSG:32619 | 30 |
| B4 | PixelType | int | 0 | 65535 | 7731 | EPSG:32619 | 0 |
| B4 | PixelType | int | 0 | 65535 | 7671 | EPSG:32619 | 227985 |
| B4 | PixelType | int | 0 | 65535 | 7731 | EPSG:32619 | 0 |
| B4 | PixelType | int | 0 | 65535 | 7671 | EPSG:32619 | -30 |
| B4 | PixelType | int | 0 | 65535 | 7731 | EPSG:32619 | -3556485 |
kbl(a_001[,c(34:41)])%>%
kable_styling(bootstrap_options = c("striped", "hover")) %>%
kable_paper() %>%
scroll_box(width = "100%", height = "300px")| bands.id.4 | bands.data_type.type.4 | bands.data_type.precision.4 | bands.data_type.min.4 | bands.data_type.max.4 | bands.dimensions.4 | bands.crs.4 | bands.crs_transform.4 |
|---|---|---|---|---|---|---|---|
| B5 | PixelType | int | 0 | 65535 | 7671 | EPSG:32619 | 30 |
| B5 | PixelType | int | 0 | 65535 | 7731 | EPSG:32619 | 0 |
| B5 | PixelType | int | 0 | 65535 | 7671 | EPSG:32619 | 227985 |
| B5 | PixelType | int | 0 | 65535 | 7731 | EPSG:32619 | 0 |
| B5 | PixelType | int | 0 | 65535 | 7671 | EPSG:32619 | -30 |
| B5 | PixelType | int | 0 | 65535 | 7731 | EPSG:32619 | -3556485 |
kbl(a_001[,c(42:49)])%>%
kable_styling(bootstrap_options = c("striped", "hover")) %>%
kable_paper() %>%
scroll_box(width = "100%", height = "300px")| bands.id.5 | bands.data_type.type.5 | bands.data_type.precision.5 | bands.data_type.min.5 | bands.data_type.max.5 | bands.dimensions.5 | bands.crs.5 | bands.crs_transform.5 |
|---|---|---|---|---|---|---|---|
| B6 | PixelType | int | 0 | 65535 | 7671 | EPSG:32619 | 30 |
| B6 | PixelType | int | 0 | 65535 | 7731 | EPSG:32619 | 0 |
| B6 | PixelType | int | 0 | 65535 | 7671 | EPSG:32619 | 227985 |
| B6 | PixelType | int | 0 | 65535 | 7731 | EPSG:32619 | 0 |
| B6 | PixelType | int | 0 | 65535 | 7671 | EPSG:32619 | -30 |
| B6 | PixelType | int | 0 | 65535 | 7731 | EPSG:32619 | -3556485 |
kbl(a_001[,c(50:57)])%>%
kable_styling(bootstrap_options = c("striped", "hover")) %>%
kable_paper() %>%
scroll_box(width = "100%", height = "300px")| bands.id.6 | bands.data_type.type.6 | bands.data_type.precision.6 | bands.data_type.min.6 | bands.data_type.max.6 | bands.dimensions.6 | bands.crs.6 | bands.crs_transform.6 |
|---|---|---|---|---|---|---|---|
| B7 | PixelType | int | 0 | 65535 | 7671 | EPSG:32619 | 30 |
| B7 | PixelType | int | 0 | 65535 | 7731 | EPSG:32619 | 0 |
| B7 | PixelType | int | 0 | 65535 | 7671 | EPSG:32619 | 227985 |
| B7 | PixelType | int | 0 | 65535 | 7731 | EPSG:32619 | 0 |
| B7 | PixelType | int | 0 | 65535 | 7671 | EPSG:32619 | -30 |
| B7 | PixelType | int | 0 | 65535 | 7731 | EPSG:32619 | -3556485 |
kbl(a_001[,c(58:70)]) %>%
kable_styling(bootstrap_options = c("striped", "hover")) %>%
kable_paper() %>%
scroll_box(width = "100%", height = "300px")| bands.id.7 | bands.data_type.type.7 | bands.data_type.precision.7 | bands.data_type.min.7 | bands.data_type.max.7 | bands.dimensions.7 | bands.crs.7 | bands.crs_transform.15L | bands.crs_transform.0L | bands.crs_transform.227992.5 | bands.crs_transform.0L.1 | bands.crs_transform..15L | bands.crs_transform..3556492.5 |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| B8 | PixelType | int | 0 | 65535 | 15341 | EPSG:32619 | 15 | 0 | 227992.5 | 0 | -15 | -3556493 |
| B8 | PixelType | int | 0 | 65535 | 15461 | EPSG:32619 | 15 | 0 | 227992.5 | 0 | -15 | -3556493 |
| B8 | PixelType | int | 0 | 65535 | 15341 | EPSG:32619 | 15 | 0 | 227992.5 | 0 | -15 | -3556493 |
| B8 | PixelType | int | 0 | 65535 | 15461 | EPSG:32619 | 15 | 0 | 227992.5 | 0 | -15 | -3556493 |
| B8 | PixelType | int | 0 | 65535 | 15341 | EPSG:32619 | 15 | 0 | 227992.5 | 0 | -15 | -3556493 |
| B8 | PixelType | int | 0 | 65535 | 15461 | EPSG:32619 | 15 | 0 | 227992.5 | 0 | -15 | -3556493 |
kbl(a_001[,c(71:78)]) %>%
kable_styling(bootstrap_options = c("striped", "hover")) %>%
kable_paper() %>%
scroll_box(width = "100%", height = "300px")| bands.id.8 | bands.data_type.type.8 | bands.data_type.precision.8 | bands.data_type.min.8 | bands.data_type.max.8 | bands.dimensions.8 | bands.crs.8 | bands.crs_transform.7 |
|---|---|---|---|---|---|---|---|
| B9 | PixelType | int | 0 | 65535 | 7671 | EPSG:32619 | 30 |
| B9 | PixelType | int | 0 | 65535 | 7731 | EPSG:32619 | 0 |
| B9 | PixelType | int | 0 | 65535 | 7671 | EPSG:32619 | 227985 |
| B9 | PixelType | int | 0 | 65535 | 7731 | EPSG:32619 | 0 |
| B9 | PixelType | int | 0 | 65535 | 7671 | EPSG:32619 | -30 |
| B9 | PixelType | int | 0 | 65535 | 7731 | EPSG:32619 | -3556485 |
kbl(a_001[,c(79:86)]) %>%
kable_styling(bootstrap_options = c("striped", "hover")) %>%
kable_paper() %>%
scroll_box(width = "100%", height = "300px")| bands.id.9 | bands.data_type.type.9 | bands.data_type.precision.9 | bands.data_type.min.9 | bands.data_type.max.9 | bands.dimensions.9 | bands.crs.9 | bands.crs_transform.8 |
|---|---|---|---|---|---|---|---|
| B10 | PixelType | int | 0 | 65535 | 7671 | EPSG:32619 | 30 |
| B10 | PixelType | int | 0 | 65535 | 7731 | EPSG:32619 | 0 |
| B10 | PixelType | int | 0 | 65535 | 7671 | EPSG:32619 | 227985 |
| B10 | PixelType | int | 0 | 65535 | 7731 | EPSG:32619 | 0 |
| B10 | PixelType | int | 0 | 65535 | 7671 | EPSG:32619 | -30 |
| B10 | PixelType | int | 0 | 65535 | 7731 | EPSG:32619 | -3556485 |
kbl(a_001[,c(87:110)]) %>%
kable_styling(bootstrap_options = c("striped", "hover")) %>%
kable_paper() %>%
scroll_box(width = "100%", height = "300px")| bands.id.10 | bands.data_type.type.10 | bands.data_type.precision.10 | bands.data_type.min.10 | bands.data_type.max.10 | bands.dimensions.10 | bands.crs.10 | bands.crs_transform.9 | bands.id.11 | bands.data_type.type.11 | bands.data_type.precision.11 | bands.data_type.min.11 | bands.data_type.max.11 | bands.dimensions.11 | bands.crs.11 | bands.crs_transform.10 | bands.id.12 | bands.data_type.type.12 | bands.data_type.precision.12 | bands.data_type.min.12 | bands.data_type.max.12 | bands.dimensions.12 | bands.crs.12 | bands.crs_transform.11 |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| B11 | PixelType | int | 0 | 65535 | 7671 | EPSG:32619 | 30 | BQA | PixelType | int | 0 | 65535 | 7671 | EPSG:32619 | 30 | nd | PixelType | float | -1 | 1 | 7671 | EPSG:32619 | 30 |
| B11 | PixelType | int | 0 | 65535 | 7731 | EPSG:32619 | 0 | BQA | PixelType | int | 0 | 65535 | 7731 | EPSG:32619 | 0 | nd | PixelType | float | -1 | 1 | 7731 | EPSG:32619 | 0 |
| B11 | PixelType | int | 0 | 65535 | 7671 | EPSG:32619 | 227985 | BQA | PixelType | int | 0 | 65535 | 7671 | EPSG:32619 | 227985 | nd | PixelType | float | -1 | 1 | 7671 | EPSG:32619 | 227985 |
| B11 | PixelType | int | 0 | 65535 | 7731 | EPSG:32619 | 0 | BQA | PixelType | int | 0 | 65535 | 7731 | EPSG:32619 | 0 | nd | PixelType | float | -1 | 1 | 7731 | EPSG:32619 | 0 |
| B11 | PixelType | int | 0 | 65535 | 7671 | EPSG:32619 | -30 | BQA | PixelType | int | 0 | 65535 | 7671 | EPSG:32619 | -30 | nd | PixelType | float | -1 | 1 | 7671 | EPSG:32619 | -30 |
| B11 | PixelType | int | 0 | 65535 | 7731 | EPSG:32619 | -3556485 | BQA | PixelType | int | 0 | 65535 | 7731 | EPSG:32619 | -3556485 | nd | PixelType | float | -1 | 1 | 7731 | EPSG:32619 | -3556485 |
La columna 111 del dataset nos da la identificación de la scene:
kbl(a_001[,c(111)]) %>%
kable_styling(bootstrap_options = c("striped", "hover")) %>%
kable_paper() %>%
scroll_box(width = "100%", height = "300px")| x |
|---|
| LANDSAT/LC08/C01/T1/LC08_233083_20200314 |
| LANDSAT/LC08/C01/T1/LC08_233083_20200314 |
| LANDSAT/LC08/C01/T1/LC08_233083_20200314 |
| LANDSAT/LC08/C01/T1/LC08_233083_20200314 |
| LANDSAT/LC08/C01/T1/LC08_233083_20200314 |
| LANDSAT/LC08/C01/T1/LC08_233083_20200314 |
object <- ee$ImageCollection("LANDSAT/LC08/C01/T1")$
filter(ee$Filter()$eq("WRS_PATH", 233))$
filter(ee$Filter()$eq("WRS_ROW", 83))$
filterMetadata ('CLOUD_COVER', 'Less_Than', 20)$
filterDate("2020-03-01", "2020-04-01")
# Representa la imagen a traves de una composicion RGB
first <- object$first()
vizParams <- list(
bands = c("B5", "B4", "B3"),
min = 5000,
max = 15000,
gamma = 1.3
)
Map$setCenter(lon = -70.64724, lat = -33.47269, zoom = 8)
Map$addLayer(first, vizParams, "Landsat 8")library(rgee)
#ee_install()
#ee_check()
ee_Initialize()## -- rgee 1.0.9 --------------------------------------- earthengine-api 0.1.259 --
## v email: not_defined
## v Initializing Google Earth Engine:
v Initializing Google Earth Engine: DONE!
##
v Earth Engine user: users/tarredwall
## --------------------------------------------------------------------------------l8 <- ee$ImageCollection("LANDSAT/LC08/C01/T1")
image_2015 <- ee$Algorithms$Landsat$simpleComposite(
collection = l8$filterDate("2015-01-01", "2015-12-31"),
asFloat = TRUE
)
bands <- c( "B5", "B6")
getNDBI <- function(image) {
image$normalizedDifference(c("B5", "B4"))
}
ndbi1 <- getNDBI(image_2015)
ndbiParams <- list(palette = c(
"#000000", "#ff0000", "#ffffff","#0000ff"
))
landMask <- ee$Image("CGIAR/SRTM90_V4")$mask()
maskedDifference <- ndbi1$updateMask(landMask)
Map$setCenter(lon = -70.64724, lat = -33.47269, zoom = 10)
Map$addLayer(maskedDifference, ndbiParams, "NDBI")Descarga de precipitación desde GEE https://rpubs.com/Bonilla-Escoto/720385
library(rgee)
library(sf)## Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1require(tidyverse)## Loading required package: tidyverse## -- Attaching packages --------------------------------------- tidyverse 1.3.0 --## v ggplot2 3.3.2 v purrr 0.3.4
## v tibble 3.0.1 v dplyr 1.0.2
## v tidyr 1.1.0 v stringr 1.4.0
## v readr 1.3.1 v forcats 0.5.0## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x readr::col_factor() masks scales::col_factor()
## x purrr::discard() masks scales::discard()
## x tidyr::extract() masks raster::extract()
## x dplyr::filter() masks stats::filter()
## x dplyr::group_rows() masks kableExtra::group_rows()
## x dplyr::lag() masks stats::lag()
## x dplyr::select() masks raster::select()punto <- st_point(c(617175,1428076), dim="XY")
puntoUTM <- st_sfc(punto,crs=32616)
puntoGeo <- st_transform(puntoUTM,4326)
punto_ee <- sf_as_ee(puntoGeo)
gee1 <- ee$ImageCollection("ECMWF/ERA5/MONTHLY")$select("total_precipitation")$
filterDate("2005-01-01", "2020-01-01")
datos_gee1 <- ee_extract(x=gee1,y=punto_ee,scale = 250,sf=T)
st_geometry(datos_gee1) <- NULL
datos_t1 <- as.data.frame(t(datos_gee1[,-1]))
head(datos_t1)## 1
## X200502_total_precipitation 0.0029145
## X200503_total_precipitation 0.0158447
## X200504_total_precipitation 0.0500956
## X200505_total_precipitation 0.1628986
## X200506_total_precipitation 0.2862051
## X200507_total_precipitation 0.1587650sal <- rownames_to_column(datos_t1, var="fecha_apilada")
colnames(sal) <- c("Fecha_apilada","Valores")
Fecha <- seq.Date(from = as.Date("01-01-2005", format="%d-%m-%Y", sep="-"),
to=as.Date("01-12-2019", format="%d-%m-%Y",sep="-"),
by="month")
#sal$Fecha <- Fecha
sal$Valores_mm <- sal$Valores*1000
mensual <- as.data.frame(matrix(sal$Valores_mm, ncol = 12, byrow = T))## Warning in matrix(sal$Valores_mm, ncol = 12, byrow = T): la longitud de los
## datos [179] no es un submúltiplo o múltiplo del número de filas [15] en la
## matrizcolnames(mensual) <- c("Ene","Feb","Mar","Abr","May","Jun","Jul","Ago","Sep","Oct","Nov","Dic")
mensual$Anio <- seq.int(from=2005, to=2019,by=1)
head(mensual)## Ene Feb Mar Abr May Jun Jul Ago
## 1 2.9145 15.8447 50.0956 162.8986 286.2051 158.7650 156.5607 130.7581
## 2 14.5964 7.5814 17.5079 92.0442 180.7454 95.7995 77.7299 99.9864
## 3 5.1926 17.1376 55.7835 193.5994 126.5605 97.5545 176.3386 203.9109
## 4 13.0933 9.3654 16.8935 255.3760 156.1875 162.6278 152.4938 219.4597
## 5 9.9886 7.9723 12.8181 257.8217 223.0067 73.8023 82.5352 109.9282
## 6 24.8376 17.5173 131.0137 431.3926 170.1693 179.8437 423.4943 269.9235
## Sep Oct Nov Dic Anio
## 1 195.7654 101.2689 23.9303 22.9142 2005
## 2 155.8129 43.0326 46.5991 16.9111 2006
## 3 321.9736 49.9877 18.0156 18.9402 2007
## 4 287.7165 33.1901 15.6414 10.3437 2008
## 5 108.7499 42.6289 16.7708 9.9737 2009
## 6 117.7062 61.2735 5.1110 14.4748 2010p<-ggplot(data=mensual, aes(x=Anio, y=Dic)) +
geom_bar(stat="identity", fill="steelblue")+
geom_text(aes(label=Dic), vjust=1.6, color="black", size=3.5)+
theme_minimal()
p
Construimos un histograma del Sentinel
Quiero muestrear píxeles aleatoriamente dentro de una imagen que tiene tres bandas, en un conjunto de regiones diferentes. Con esos datos muestreados, quiero crear un gráfico que contenga histogramas para cada región de muestra, para una banda determinada. Entonces tendré tres parcelas (una para cada banda).
En mi caso, una región representa un tipo específico de hábitat dentro de la imagen, por lo que tengo un polígono para cada tipo de hábitat.
library(ggplot2)
###############################################################
# Create a geometry for each habitat type that I want to sample
###############################################################
unimproved_grazing = ee$Geometry$Polygon(
list(
c(-1.345864517292552, 59.97338078318311),
c(-1.345864517292552, 59.97237144577356),
c(-1.3438904114575911, 59.97237144577356),
c(-1.3438904114575911, 59.97338078318311)
)
)
improved_cut = ee$Geometry$Polygon(
list(
c(-1.329604032461742, 59.92455155105425),
c(-1.329604032461742, 59.923922477208166),
c(-1.3280268935609851, 59.923922477208166),
c(-1.3280268935609851, 59.92455155105425)
)
)
heath = ee$Geometry$Polygon(
list(
c(-1.351189448199126, 59.94125887709115),
c(-1.351189448199126, 59.940151826756974),
c(-1.3490865963314502, 59.940151826756974),
c(-1.3490865963314502, 59.94125887709115)
)
)
arable = ee$Geometry$Polygon(
list(
c(-1.3263107736431534, 59.969782903189916),
c(-1.3263107736431534, 59.96911710113049),
c(-1.325023313316005, 59.96911710113049),
c(-1.325023313316005, 59.969782903189916)
)
)
##################################################
# Create feature collection of all habitat regions
##################################################
habitats = ee$FeatureCollection(
list(
ee$Feature(unimproved_grazing, list(name = "unimproved")),
ee$Feature(improved_cut, list(name = "improved")),
ee$Feature(heath, list(name = "heath")),
ee$Feature(arable, list(name = "arable"))
)
)
sentinel1 = ee$ImageCollection('COPERNICUS/S1_GRD')$
# Filter date range - one year of images
filterDate("2019-03-01", "2020-03-01")
# Filter by metadata properties
vvvh = sentinel1$
filter(ee$Filter$listContains('transmitterReceiverPolarisation', 'VV'))$
filter(ee$Filter$listContains('transmitterReceiverPolarisation', 'VH'))$
filter(ee$Filter$eq('instrumentMode', 'IW'))$
filter(ee$Filter$eq('resolution_meters', 10))
# Filter to get images from different look angles.
vhAscending = vvvh$filter(ee$Filter$eq('orbitProperties_pass', 'ASCENDING'))
vhDescending = vvvh$filter(ee$Filter$eq('orbitProperties_pass', 'DESCENDING'))
# Create a composite from means at different polarizations and look angles.
# Calculate mean of VH ascending polarisation
VH_Ascending_mean <- vhAscending$select("VH")$mean()
# Calculate mean of VH descending polarisation
VH_Descending_mean <- vhDescending$select("VH")$mean()
# Merge VV and VH and calculate mean
VV_Ascending_Descending_mean <- vhAscending$select("VV") %>%
ee$ImageCollection$merge(vhDescending$select("VV")) %>%
ee$ImageCollection$mean()
# Create single image containing bands of interest
collection <- ee$ImageCollection$fromImages(list(
VH_Ascending_mean,
VV_Ascending_Descending_mean,
VH_Descending_mean
))$toBands()
######################################
# Sample pixels in each habitat region
######################################
mySamples = collection$sampleRegions(
collection = habitats,
scale= 10
)
# ee_help(collection$sampleRegions)
for (index in 0:3) {
habitat_f <- habitats %>% ee_get(index)
mySamples <- collection$sampleRegions(
collection = ee$Feature(habitat_f$first()),
scale= 10,
geometries=TRUE
)
cat(sprintf("Processing geom: %02d \n", index + 1))
mySamples_sf <- ee_as_sf(mySamples)
if (index == 0) {
mySamples_sf_batch <- mySamples_sf
} else {
mySamples_sf_batch <- rbind(mySamples_sf_batch, mySamples_sf)
}
}## Processing geom: 01
## Processing geom: 02
## Processing geom: 03
## Processing geom: 04ggplot(mySamples_sf_batch, aes(X1_VV)) +
geom_histogram(color="darkblue", fill="lightblue") +
facet_grid(~name)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
1.4 La matemática de las bandas
Haremos una sustracción entre los Índices de Vegetación de Diferencia Normalizada (NDVI) de dos imágenes.
getNDVI <- function(image) {
return(image$normalizedDifference(c("B4", "B3")))
}
image1 <- ee$Image("LANDSAT/LT05/C01/T1_TOA/LT05_044034_19900604")
image2 <- ee$Image("LANDSAT/LT05/C01/T1_TOA/LT05_044034_20100611")
ndvi1 <- getNDVI(image1)
ndvi2 <- getNDVI(image2)
ndviDifference <- ndvi2$subtract(ndvi1)
ndviParams <- list(palette = c(
"#d73027", "#f46d43", "#fdae61", "#fee08b",
"#d9ef8b", "#a6d96a", "#66bd63", "#1a9850"
))
ndwiParams <- list(min = -0.5, max = 0.5, palette = c("FF0000", "FFFFFF", "0000FF"))
Map$centerObject(ndvi1)
Map$addLayer(ndvi1, ndviParams, "NDVI 1") +
Map$addLayer(ndvi2, ndviParams, "NDVI 2") +
Map$addLayer(ndviDifference, ndwiParams, "NDVI difference")1.5 La función de mapeo
Podemos recorrer todos los datos contenidos en una ImageCollection sin necesidad de usar un ciclo for:
library(rgee)
library(kableExtra)
ee_Initialize()## -- rgee 1.0.9 --------------------------------------- earthengine-api 0.1.259 --
## v email: not_defined
## v Initializing Google Earth Engine:
v Initializing Google Earth Engine: DONE!
##
v Earth Engine user: users/tarredwall
## --------------------------------------------------------------------------------addNDVI <- function(image) {
return(image$addBands(image$normalizedDifference(c("B5", "B4"))))
}
collection <- ee$ImageCollection("LANDSAT/LC08/C01/T1")$
filterBounds(ee$Geometry$Point(-122.262, 37.8719))$
filterDate("2014-06-01", "2014-10-01")
ndviCollection <- collection$map(addNDVI)
first <- ndviCollection$first()
a_004 <- first$getInfo()
f_004 <- tail(a_004)
g_004 <- head(f_004,2)
write.csv2(g_004,"g_004.csv")
s <- read.csv2("g_004.csv")
s[,c(1:20)]%>%
kbl(caption = "Datos contenidos en LANDSAT/LC08/C01/T1") %>%
kable_classic(full_width = F, html_font = "Cambria")%>%
scroll_box(width = "100%", height = "200px")| X | type | bands.id | bands.data_type.type | bands.data_type.precision | bands.data_type.min | bands.data_type.max | bands.dimensions | bands.crs | bands.crs_transform | bands.id.1 | bands.data_type.type.1 | bands.data_type.precision.1 | bands.data_type.min.1 | bands.data_type.max.1 | bands.dimensions.1 | bands.crs.1 | bands.crs_transform.1 | bands.id.2 | bands.data_type.type.2 |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | Image | B1 | PixelType | int | 0 | 65535 | 7671 | EPSG:32610 | 30 | B2 | PixelType | int | 0 | 65535 | 7671 | EPSG:32610 | 30 | B3 | PixelType |
| 2 | Image | B1 | PixelType | int | 0 | 65535 | 7801 | EPSG:32610 | 0 | B2 | PixelType | int | 0 | 65535 | 7801 | EPSG:32610 | 0 | B3 | PixelType |
| 3 | Image | B1 | PixelType | int | 0 | 65535 | 7671 | EPSG:32610 | 463785 | B2 | PixelType | int | 0 | 65535 | 7671 | EPSG:32610 | 463785 | B3 | PixelType |
| 4 | Image | B1 | PixelType | int | 0 | 65535 | 7801 | EPSG:32610 | 0 | B2 | PixelType | int | 0 | 65535 | 7801 | EPSG:32610 | 0 | B3 | PixelType |
| 5 | Image | B1 | PixelType | int | 0 | 65535 | 7671 | EPSG:32610 | -30 | B2 | PixelType | int | 0 | 65535 | 7671 | EPSG:32610 | -30 | B3 | PixelType |
| 6 | Image | B1 | PixelType | int | 0 | 65535 | 7801 | EPSG:32610 | 4264515 | B2 | PixelType | int | 0 | 65535 | 7801 | EPSG:32610 | 4264515 | B3 | PixelType |
1.6 La función Reduce.
Podemos calcular la mediana de cada banda de las última 5 escenas de nubes:
collection <- ee$ImageCollection("LANDSAT/LC08/C01/T1")$
filterBounds(ee$Geometry$Point(-122.262, 37.8719))$
filterDate("2014-01-01", "2014-12-31")$
sort("CLOUD_COVER")
median <- collection$limit(5)$reduce(ee$Reducer$median())
vizParams <- list(
bands = c("B5_median", "B4_median", "B3_median"),
min = 5000, max = 15000, gamma = 1.3
)
Map$addLayer(
eeObject = median,
visParams = vizParams,
name = "Median image"
)1.7 Estadísticas de las imágenes
Cargamos y desplegamos una imagen de la parte superior de la atmósfera calibrada (TOA) del Landsat 8:
image <- ee$Image("LANDSAT/LC08/C01/T1_TOA/LC08_044034_20140318")
Map$addLayer(
eeObject = image,
visParams = list(bands = c("B4", "B3", "B2"), max = 30000),
name = "Landsat 8"
)Creamos un rectángulo arbitrario como una región y lo desplegamos:
region <- ee$Geometry$Rectangle(-122.2806, 37.1209, -122.0554, 37.2413)
Map$addLayer(
eeObject = region,
name = "Region"
)Obtenemos un diccionario de las medias en la región:
mean <- image$reduceRegion(
reducer = ee$Reducer$mean(),
geometry = region,
scale = 30
)
print(mean$getInfo())## $B1
## [1] 0.1015666
##
## $B10
## [1] 288.3895
##
## $B11
## [1] 287.4494
##
## $B2
## [1] 0.07771911
##
## $B3
## [1] 0.0553003
##
## $B4
## [1] 0.03670828
##
## $B5
## [1] 0.2213975
##
## $B6
## [1] 0.07888247
##
## $B7
## [1] 0.03508945
##
## $B8
## [1] 0.04710009
##
## $B9
## [1] 0.00156097
##
## $BQA
## [1] 2720.1191.8 Enmascaramiento
Aunque la composición con la mediana es una mejora, es posible enmascarar partes de la imagen. Enmascarar píxeles de una imagen hace que esos píxeles sean transparentes y los excluye del análisis. Cada píxel de cada banda de una imagen tiene una máscara. Aquellos con un valor de máscara de 0 o menos serán transparentes. Aquellos con una máscara de cualquier valor por encima de 0 serán renderizados.
En el siguiente código creamos una función que obtiene una diferencia normalizada entre dos bandas, cargamos dos imágenes del Landsat 5, aplicamos la función, calculamos su diferencia, cargamos la máscara como archivos de demostración de videojuego (DEM) desde la Misión de topografía de radar Shuttle (SRTM), actualizamos la máscara de diferencia NDVI con la landMask y visualizamos el resultado enmascarado.
getNDVI <- function(image) {
return(image$normalizedDifference(c("B4", "B3")))
}
image1 <- ee$Image("LANDSAT/LT05/C01/T1_TOA/LT05_044034_19900604")
image2 <- ee$Image("LANDSAT/LT05/C01/T1_TOA/LT05_044034_20100611")
ndvi1 <- getNDVI(image1)
ndvi2 <- getNDVI(image2)
ndviDifference <- ndvi2$subtract(ndvi1)
landMask <- ee$Image("CGIAR/SRTM90_V4")$mask()
maskedDifference <- ndviDifference$updateMask(landMask)
vizParams <- list(
min = -0.5,
max = 0.5,
palette = c("FF0000", "FFFFFF", "0000FF")
)
Map$addLayer(
eeObject = maskedDifference,
visParams = vizParams,
name = "NDVI difference"
)2. Aprendizaje automático
2.1 Clasificación supervisada
2.1.1 Árboles de clasificación y regresión (CART)
Tuvimos que hacer un cambio en el código contenido en nuestra página madre. cart da el siguiente error:
cart: Error in py_call_impl(callable, dots\(args, dots\)keywords) : EEException: Classifier.cart: This classifier has been removed. For more information see: http://goo.gle/deprecated-classifiers.
Debió ser reemplazada por: smileCart
Referencia: https://developers.google.com/earth-engine/transitions#new-style-classifiers
Las imágenes de entrada son un compuesto Landsat 8 sin nubes, cuyas bandas utilizamos para la predicción. Cargamos los puntos de entrenamiento. La propiedad numérica ‘clase’ almacena etiquetas conocidas. Ésta propiedad de la tabla almacena las etiquetas de cobertura terrestre. Superponemos los puntos en las imágenes para entrenar. Entrenamos un clasificador CART con parámetros predeterminados. Clasificamos la imagen con las mismas bandas utilizadas para el entrenamiento. Establecemos los parámetros de visualización y visualizamos las entradas y los resultados.
# Input imagery is a cloud-free Landsat 8 composite.
l8 <- ee$ImageCollection("LANDSAT/LC08/C01/T1")
image <- ee$Algorithms$Landsat$simpleComposite(
collection = l8$filterDate("2018-01-01", "2018-12-31"),
asFloat = TRUE
)
# Use these bands for prediction.
bands <- c("B2", "B3", "B4", "B5", "B6", "B7", "B10", "B11")
# Load training points. The numeric property 'class' stores known labels.
points <- ee$FeatureCollection("GOOGLE/EE/DEMOS/demo_landcover_labels")
# This property of the table stores the land cover labels.
label <- "landcover"
# Overlay the points on the imagery to get training.
training <- image$select(bands)$sampleRegions(
collection = points,
properties = list(label),
scale = 30
)
# Train a CART classifier with default parameters.
trained <- ee$Classifier$smileCart()$train(training, label, bands)
# Classify the image with the same bands used for training.
classified <- image$select(bands)$classify(trained)
# Viz parameters.
viz_img <- list(bands = c("B4", "B3", "B2"), max = 0.4)
viz_class <- list(palette = c("red", "green", "blue"), min = 0, max = 2)
# Display the inputs and the results.
Map$centerObject(points)## NOTE: Center obtained from the first element.Map$addLayer(image, viz_img, name = "image") +
Map$addLayer(classified, viz_class, name = "classification")2.1.2 Máquinas de soporte vectorial (SVM)
svm: Error in py_call_impl(callable, dots\(args, dots\)keywords) : EEException: Classifier.svm: This classifier has been replaced. For more information see: http://goo.gle/deprecated-classifiers.
tuvo que ser reemplazado por libsvm
Referencia: https://developers.google.com/earth-engine/guides/classification
# Input imagery is a cloud-free Landsat 8 composite.
l8 = ee$ImageCollection("LANDSAT/LC08/C01/T1")
image = ee$Algorithms$Landsat$simpleComposite(
collection = l8$filterDate("2018-01-01", "2018-12-31"),
asFloat = TRUE
)
# Use these bands for prediction.
bands = c("B2", "B3", "B4", "B5", "B6", "B7", "B10", "B11")
# Manually created polygons.
forest1 = ee$Geometry$Rectangle(-63.0187, -9.3958, -62.9793, -9.3443)
forest2 = ee$Geometry$Rectangle(-62.8145, -9.206, -62.7688, -9.1735)
nonForest1 = ee$Geometry$Rectangle(-62.8161, -9.5001, -62.7921, -9.4486)
nonForest2 = ee$Geometry$Rectangle(-62.6788, -9.044, -62.6459, -8.9986)
# Make a FeatureCollection from the hand-made geometries.
polygons = ee$FeatureCollection(c(
ee$Feature(nonForest1, list(class = 0)),
ee$Feature(nonForest2, list(class = 0)),
ee$Feature(forest1, list(class = 1)),
ee$Feature(forest2, list(class = 1))
))
# Get the values for all pixels in each polygon in the training.
training = image$sampleRegions(
# Get the sample from the polygons FeatureCollection.
collection = polygons,
# Keep this list of properties from the polygons.
properties = list("class"),
# Set the scale to get Landsat pixels in the polygons.
scale = 30
)
# Create an SVM classifier with custom parameters.
classifier = ee$Classifier$libsvm(
kernelType = "RBF",
gamma = 0.5,
cost = 10
)
# Train the classifier.
trained = classifier$train(training, "class", bands)
# Classify the image.
classified = image$classify(trained)
# Display the classification result and the input image.
geoviz_image = list(bands = c("B4", "B3", "B2"), max = 0.5, gamma = 2)
geoviz_class = list(min = 0, max = 1, palette = c("red", "green"))
Map$setCenter(-62.836, -9.2399, 9)
Map$addLayer(
eeObject = image,
visParams = geoviz_image,
name = "image"
) +
Map$addLayer(
eeObject = classified,
visParams = geoviz_class,
name = "deforestation"
) +
Map$addLayer(
eeObject = polygons,
name = "training polygons"
)2.2 Clasificación no supervisada
2.2.1 Clustering por Kmeans
# Load a pre-computed Landsat composite for input.
input <- ee$Image("LANDSAT/LE7_TOA_1YEAR/2001")
# Define a region in which to generate a sample of the input.
region <- ee$Geometry$Rectangle(29.7, 30, 32.5, 31.7)
# Display the sample region.
Map$setCenter(31.5, 31.0, 7)
Map$addLayer(
eeObject = ee$Image()$paint(region, 0, 2),
name = "region"
)# Make the training dataset.
training <- input$sample(
region = region,
scale = 30,
numPixels = 5000
)
# Instantiate the clusterer and train it.
clusterer <- ee$Clusterer$wekaKMeans(15)$train(training)
# Cluster the input using the trained clusterer.
result <- input$cluster(clusterer)
# Display the clusters with random colors.
Map$centerObject(region)
Map$addLayer(
eeObject = result$randomVisualizer(),
name = "clusters"
)3. Imágenes
3.1 Introducción
Crea una imagen constante con un valor de píxel de 1 Crea una imagen constante con un valor de píxel de 1 Concatenar imágenes Cambie la opción de impresión por: “simplemente”, “json”, “ee_print” Cambiar el nombre de las imágenes usando ee $ Image $ select
# Create a constant Image with a pixel value of 1
image1 <- ee$Image(1)
print(image1, type = "json")## ee.Image({
## "functionInvocationValue": {
## "functionName": "Image.constant",
## "arguments": {
## "value": {
## "constantValue": 1.0
## }
## }
## }
## })print(image1, type = "simply")## EarthEngine Object: Imageprint(image1, type = "ee_print")## ---------------------------------------------------------- Earth Engine Image --
## Image Metadata:
## - Class : ee$Image
## - ID : no_id
## - Number of Bands : 1
## - Bands names : constant
## - Number of Properties : 0
## - Number of Pixels* : 64800
## - Approximate size* : 101.25 KB
## Band Metadata (img_band = constant):
## - EPSG (SRID) : WGS 84 (EPSG:4326)
## - proj4string : +proj=longlat +datum=WGS84 +no_defs
## - Geotransform : 1 0 0 0 1 0
## - Nominal scale (meters) : 111319.5
## - Dimensions : 360 180
## - Number of Pixels : 64800
## - Data type : INT
## - Approximate size : 101.25 KB
## --------------------------------------------------------------------------------
## NOTE: (*) Properties calculated considering a constant geotransform and data type.# Create a constant Image with a pixel value of 1
image2 <- ee$Image(2)
# Concatenate Images
image3 <- ee$Image$cat(c(image1, image2))
ee_print(image3, clean = TRUE)## ---------------------------------------------------------- Earth Engine Image --
## Image Metadata:
## - Class : ee$Image
## - ID : no_id
## - Number of Bands : 2
## - Bands names : constant constant_1
## - Number of Properties : 0
## - Number of Pixels* : 129600
## - Approximate size* : 202.50 KB
## Band Metadata (img_band = constant):
## - EPSG (SRID) : WGS 84 (EPSG:4326)
## - proj4string : +proj=longlat +datum=WGS84 +no_defs
## - Geotransform : 1 0 0 0 1 0
## - Nominal scale (meters) : 111319.5
## - Dimensions : 360 180
## - Number of Pixels : 64800
## - Data type : INT
## - Approximate size : 101.25 KB
## --------------------------------------------------------------------------------
## NOTE: (*) Properties calculated considering a constant geotransform and data type.# Change print option by: "simply","json","ee_print"
options(rgee.print.option = "simply")
multiband <- ee$Image(c(1, 2, 3))
print(multiband)## EarthEngine Object: Image# Rename images using ee$Image$select
renamed <- multiband$select(
opt_selectors = c("constant", "constant_1", "constant_2"),
opt_names = c("band1", "band2", "band3")
)
ee_print(renamed)## ---------------------------------------------------------- Earth Engine Image --
## Image Metadata:
## - Class : ee$Image
## - ID : no_id
## - Number of Bands : 3
## - Bands names : band1 band2 band3
## - Number of Properties : 0
## - Number of Pixels* : 194400
## - Approximate size* : 303.75 KB
## Band Metadata (img_band = band1):
## - EPSG (SRID) : WGS 84 (EPSG:4326)
## - proj4string : +proj=longlat +datum=WGS84 +no_defs
## - Geotransform : 1 0 0 0 1 0
## - Nominal scale (meters) : 111319.5
## - Dimensions : 360 180
## - Number of Pixels : 64800
## - Data type : INT
## - Approximate size : 101.25 KB
## --------------------------------------------------------------------------------
## NOTE: (*) Properties calculated considering a constant geotransform and data type.3.2 Visualización de imágenes
# Simple RGB Vizualization -Landsat
image <- ee$Image("LANDSAT/LC8_L1T_TOA/LC80440342014077LGN00")
vizParams <- list(
min = 0,
max = 0.5,
bands = c("B5", "B4", "B3"),
gamma = c(0.95, 1.1, 1)
)
Map$setCenter(-122.1899, 37.5010, 5)
Map$addLayer(image, vizParams)# Simple Single-Band Vizualization
ndwi <- image$normalizedDifference(c("B3", "B5"))
ndwiViz <- list(
min = 0.5,
max = 1,
palette = c("00FFFF", "0000FF")
)
## Masking
ndwiMasked <- ndwi$updateMask(ndwi$gte(0.4))
## Produce an RGB layers.
imageRGB <- image$visualize(
list(
bands = c("B5", "B4", "B3"),
max = 0.5
)
)
ndwiRGB <- ndwiMasked$visualize(
list(
min = 0.5,
max = 1,
palette = c("00FFFF", "0000FF")
)
)
# Mosaic the visualization layers and display( or export).
roi <- ee$Geometry$Point(c(-122.4481, 37.7599))$buffer(20000)
Map$centerObject(image$clip(roi))
Map$addLayer(
eeObject = image$clip(roi),
visParams = vizParams,
name = "Landsat 8"
) +
Map$addLayer(
eeObject = ndwiMasked$clip(roi),
visParams = ndwiViz,
name = "NDWI"
)3.3 Información sobre las imágenes y su metadata
# # Load an Landsat8 Image - San Francisco, California.
# image <- ee$Image("LANDSAT/LC8_L1T/LC80440342014077LGN00")
#
# # Get Band names of an ee$Image
# bandNames <- image$bandNames()
# #ee_help(bandNames)
# cat("Band names: ", paste(bandNames$getInfo(), collapse = " "))
#
# # Get Band names of an ee$Image
# b1proj <- image$select("B1")$projection()$getInfo()
# cat("Band 1 projection: ")
# cat("type: ", b1proj$type)
# cat("crs: ", b1proj$crs)
# cat("geotransform: ", paste0(b1proj$transform, " "))
#
# b1scale <- image$select("B1")$projection()$nominalScale()
# cat("Band 1 scale: ", b1scale$getInfo())
#
# b8scale <- image$select("B8")$projection()$nominalScale()
# cat("Band 8 scale: ", b8scale$getInfo())
#
# properties <- image$propertyNames()$getInfo()
# cat("Metadata properties: \n-", paste0(properties, collapse = "\n- "))
#
# cloudiness <- image$get("CLOUD_COVER")$getInfo()
# cat("CLOUD_COVER: ", cloudiness)
#
# iso_date <- eedate_to_rdate(image$get("system:time_start"))
# iso_timestamp <- eedate_to_rdate(
# ee_date = image$get("system:time_start"),
# js = TRUE
# )
# cat("ISO Date: ", as.character(iso_date))
# cat("Timestamp : ", format(iso_timestamp, scientific = FALSE))3.4 Operaciones matemáticas
# Load two 5-year Landsat 7 composites.
landsat1999 <- ee$Image("LANDSAT/LE7_TOA_5YEAR/1999_2003")
landsat2008 <- ee$Image("LANDSAT/LE7_TOA_5YEAR/2008_2012")
# Compute NDVI the hard way.
ndvi1999 <- landsat1999$select("B4")$subtract(landsat1999$select("B3"))$divide(landsat1999$select("B4")$add(landsat1999$select("B3")))
# Compute NDVI the easy way.
ndvi2008 <- landsat2008$normalizedDifference(c("B4", "B3"))
# Compute the multi-band difference image.
diff <- landsat2008$subtract(landsat1999)
# Compute the squared difference in each band.
squaredDifference <- diff$pow(2)
Map$setZoom(zoom = 3)
Map$addLayer(
eeObject = diff,
visParams = list(bands = c("B4", "B3", "B2"), min = -32, max = 32),
name = "difference"
) +
Map$addLayer(
eeObject = squaredDifference,
visParams = list(bands = c("B4", "B3", "B2"), max = 1000),
name = "squared diff."
)3.5 Relational, conditional and Boolean operations
# Load a Landsat 8 image.
image <- ee$Image("LANDSAT/LC08/C01/T1_TOA/LC08_044034_20140318")
# Create NDVI and NDWI spectral indices.
ndvi <- image$normalizedDifference(c("B5", "B4"))
ndwi <- image$normalizedDifference(c("B3", "B5"))
# Create a binary layer using logical operations.
bare <- ndvi$lt(0.2)$And(ndwi$lt(0))
# Mask and display the binary layer.
Map$setCenter(lon = -122.3578, lat = 37.7726)
Map$setZoom(zoom = 10)
Map$addLayer(
eeObject = bare$updateMask(bare),
visParams = list(min = 0, max = 20),
name = "bare"
)3.6 Convoluciones
# Load and display an image.
image <- ee$Image("LANDSAT/LC08/C01/T1_TOA/LC08_044034_20140318")
# Define a boxcar or low-pass kernel.
# boxcar = ee$Kernel$square(list(
# radius = 7, units = 'pixels', 'normalize' = T
# ))
boxcar <- ee$Kernel$square(7, "pixels", T)
# Smooth the image by convolving with the boxcar kernel.
smooth <- image$convolve(boxcar)
# Define a Laplacian, or edge-detection kernel.
laplacian <- ee$Kernel$laplacian8(1, F)
# Apply the edge-detection kernel.
edgy <- image$convolve(laplacian)
Map$setCenter(lon = -121.9785, lat = 37.8694)
Map$setZoom(zoom = 11)
Map$addLayer(
eeObject = image,
visParams = list(bands = c("B5", "B4", "B3"), max = 0.5),
name = "input image"
) +
Map$addLayer(
eeObject = smooth,
visParams = list(bands = c("B5", "B4", "B3"), max = 0.5),
name = "smoothed"
) +
Map$addLayer(
eeObject = edgy,
visParams = list(bands = c("B5", "B4", "B3"), max = 0.5),
name = "edges"
)3.7 Operaciones morfológicas
# Load a Landsat 8 image, select the NIR band, threshold, display.
image <- ee$Image("LANDSAT/LC08/C01/T1_TOA/LC08_044034_20140318")$select(4)$gt(0.2)
# Define a kernel.
kernel <- ee$Kernel$circle(radius = 1)
# Perform an erosion followed by a dilation, display.
opened <- image$focal_min(kernel = kernel, iterations = 2)$focal_max(kernel = kernel, iterations = 2)
Map$setCenter(lon = -122.1899, lat = 37.5010)
Map$setZoom(zoom = 13)
Map$addLayer(
eeObject = image,
visParams = list(min = 0, max = 1),
name = "NIR threshold"
) +
Map$addLayer(
eeObject = opened,
visParams = list(min = 0, max = 1),
name = "opened"
)3.8 Gradientes
# Load a Landsat 8 image and select the panchromatic band.
image <- ee$Image("LANDSAT/LC08/C01/T1/LC08_044034_20140318")$select("B8")
# Compute the image gradient in the X and Y directions.
xyGrad <- image$gradient()
# Compute the magnitude of the gradient.
gradient <- xyGrad$select("x")$pow(2)$add(xyGrad$select("y")$pow(2))$sqrt()
# Compute the direction of the gradient.
direction <- xyGrad$select("y")$atan2(xyGrad$select("x"))
# Display the results.
Map$setCenter(lon = -122.054, lat = 37.7295)
Map$setZoom(zoom = 10)
Map$addLayer(
eeObject = direction,
visParams = list(min = -2, max = 2),
name = "direction"
) +
Map$addLayer(
eeObject = gradient,
visParams = list(min = -7, max = 7),
name = "opened"
)3.9 Detección de edges
# Load a Landsat 8 image, select the panchromatic band$
image <- ee$Image("LANDSAT/LC08/C01/T1/LC08_044034_20140318")$select("B8")
# Perform Canny edge detection and display the result$
canny <- ee$Algorithms$CannyEdgeDetector(image, threshold = 10, sigma = 1)
# Perform Hough transform of the Canny result and display$
hough <- ee$Algorithms$HoughTransform(canny, 256, 600, 100)
# Load a Landsat 8 image, select the panchromatic band$
image <- ee$Image("LANDSAT/LC08/C01/T1/LC08_044034_20140318")$select("B8")
# Define a "fat" Gaussian kernel$
fat <- ee$Kernel$gaussian(
radius = 3,
sigma = 3,
units = "pixels",
normalize = T,
magnitude = -1
)
# Define a "skinny" Gaussian kernel.
skinny <- ee$Kernel$gaussian(
radius = 3,
sigma = 1,
units = "pixels",
normalize = T
)
# Compute a difference-of-Gaussians (DOG) kernel.
dog <- fat$add(skinny)
# Compute the zero crossings of the second derivative, display.
zeroXings <- image$convolve(dog)$zeroCrossing()
Map$setCenter(lon = -122.054, lat = 37.7295)
Map$setZoom(zoom = 10)
Map$addLayer(
eeObject = canny,
visParams = list(min = 0, max = 1),
name = "canny"
) +
Map$addLayer(
eeObject = hough,
visParams = list(min = 0, max = 1),
name = "hough"
) +
Map$addLayer(
eeObject = image,
visParams = list(min = 0, max = 12000),
name = "L_B8"
) +
Map$addLayer(
eeObject = zeroXings$updateMask(zeroXings),
visParams = list(palette = "FF0000"),
name = "zero crossings"
)continuará…
4. ImageCollection
4.2 ImageCollection Information and Metadata
# # Load a Landsat 8 ImageCollection for a single path-row.
#
# collection <- ee$ImageCollection("LANDSAT/LC08/C01/T1_TOA")$
# filter(ee$Filter$eq("WRS_PATH", 44))$
# filter(ee$Filter$eq("WRS_ROW", 34))$
# filterDate("2014-03-01", "2014-09-01")
# ee_print(collection)
#
# # Get the number of images.
#
# count <- collection$size()
# cat("Count: ", count$getInfo())
#
# # Get the date range of images in the collection.
#
# range <- collection$reduceColumns(
# ee$Reducer$minMax(),
# list("system:time_start")
# )
#
# col_min <- eedate_to_rdate(range$get("min"))
# col_max <- eedate_to_rdate(range$get("max"))
# cat("Date range: ", as.character(col_min), as.character(col_max))
#
# # Get statistics for a property of the images in the collection.
#
# sunStats <- collection$aggregate_stats("SUN_ELEVATION")
# cat("Sun elevation statistics: ")
# sunStats$getInfo()
#
# # Sort by a cloud cover property, get the least cloudy image.
#
# image <- ee$Image(collection$sort("CLOUD_COVER")$first())
# cat("Least cloudy image: ")
# image$getInfo()
#
# # Limit the collection to the 10 most recent images.
#
# recent <- collection$sort("system:time_start", FALSE)$limit(10)
# cat("Recent images: ")
# recent$getInfo()4.3 Filtering an ImageCollection
{r} # # Load Landsat 5 data, filter by date and bounds. # collection <- ee$ImageCollection("LANDSAT/LT05/C01/T2")$ # filterDate("1987-01-01", "1990-05-01")$ # filterBounds(ee$Geometry$Point(25.8544, -18.08874)) # # # Also filter the collection by the IMAGE_QUALITY property. # filtered <- collection$filterMetadata("IMAGE_QUALITY", "equals", 9) # # # Create two composites to check the effect of filtering by IMAGE_QUALITY. # badComposite <- ee$Algorithms$Landsat$simpleComposite(collection, 75, 3) # goodComposite <- ee$Algorithms$Landsat$simpleComposite(filtered, 75, 3) # # # Display the composites. # Map$setCenter(lon = 25.8544, lat = -18.08874) # Map$setZoom(zoom = 13) # # Map$addLayer( # eeObject = badComposite, # visParams = list(bands = c("B3", "B2", "B1"), gain = 3.5), # name = "bad composite" # ) + # Map$addLayer( # eeObject = goodComposite, # visParams = list(bands = c("B3", "B2", "B1"), gain = 3.5), # name = "good composite" # ) #
4.4 Mapping over an ImageCollection
# # This function adds a band representing the image timestamp.
# addTime <- function(image) {
# return(image$addBands(image$metadata("system:time_start")))
# }
#
# conditional <- function(image) {
# return(ee$Algorithms$If(
# ee$Number(image$get("SUN_ELEVATION"))$gt(40),
# image,
# ee$Image(0)
# ))
# }
#
# # Load a Landsat 8 collection for a single path-row.
# collection <- ee$ImageCollection("LANDSAT/LC08/C01/T1_TOA")$
# filter(ee$Filter$eq("WRS_PATH", 44))$
# filter(ee$Filter$eq("WRS_ROW", 34))
#
# # Map the function over the collection and display the result.
# print(collection$map(addTime)$getInfo())
#
#
# # Load a Landsat 8 collection for a single path-row.
# collection <- ee$ImageCollection("LANDSAT/LC8_L1T_TOA")$
# filter(ee$Filter$eq("WRS_PATH", 44))$
# filter(ee$Filter$eq("WRS_ROW", 34))
#
# # This function uses a conditional statement to return the image if
# # the solar elevation > 40 degrees. Otherwise it returns a zero image.
# # conditional = function(image) {
# # return ee.Algorithms.If(ee.Number(image.get('SUN_ELEVATION')).gt(40),
# # image,
# # ee.Image(0))
# # }
#
# # Map the function over the collection, convert to a List and print the result.
# print(collection$map(conditional)$getInfo())4.5 Reducing an ImageCollection
```{r x_002}
addTime <- function(image) {
image\(addBands(image\)metadata(“system:time_start”)\(divide(1000 * 60 * 60 * 24 * 365)) # } # # # Load a Landsat 8 collection for a single path-row. # collection <- ee\)ImageCollection(“LANDSAT/LC08/C01/T1_TOA”)$
filter(ee\(Filter\)eq(“WRS_PATH”, 44))$
filter(ee\(Filter\)eq(“WRS_ROW”, 34))$
filterDate(“2014-01-01”, “2015-01-01”)
# Compute a median image and display.
median <- collection\(median() # Map\)setCenter(lon = -122.3578, lat = 37.7726)
Map\(setZoom(zoom = 12) # # Map\)addLayer(
eeObject = median,
visParams = list(bands = c(“B4”, “B3”, “B2”), max = 0.3),
name = “median”
)
# Reduce the collection with a median reducer.
median <- collection\(reduce(ee\)Reducer\(median()) # # # Display the median image. # Map\)addLayer(
eeObject = median,
visParams = list(bands = c(“B4_median”, “B3_median”, “B2_median”), max = 0.3),
name = “also median”
)
# This function adds a band representing the image timestamp.
addTime = function(image) {
return image.addBands(image.metadata(‘system:time_start’)
# Convert milliseconds from epoch to years to aid in
# interpretation of the following trend calculation.
.divide(1000 * 60 * 60 * 24 * 365))
}
Load a MODIS collection, filter to several years of 16 day mosaics,
and map the time band function over it.
collection <- ee\(ImageCollection("MODIS/006/MYD13A1")\) filterDate(“2004-01-01”, “2010-10-31”)$ map(addTime)
Select the bands to model with the independent variable first.
trend <- collection\(select(list("system:time_start", "EVI"))\) reduce(ee\(Reducer\)linearFit())
Display the trend with increasing slopes in green, decreasing in red.
Map\(setCenter(lon = -96.943, lat = 39.436) Map\)setZoom(zoom = 5)
Map$addLayer( eeObject = trend, visParams = list(bands = c(“scale”, “scale”, “offset”), min = 0, max = c(-100, 100, 10000)), name = “EVI trend” )
## 4.6 Compositing and Mosaicking
```r
# Load three NAIP quarter quads in the same location, different times.
naip2004_2012 <- ee$ImageCollection("USDA/NAIP/DOQQ")$
filterBounds(ee$Geometry$Point(-71.08841, 42.39823))$
filterDate("2004-07-01", "2012-12-31")$
select(c("R", "G", "B"))
# Temporally composite the images with a maximum value function.
composite <- naip2004_2012$max()
Map$setCenter(lon = -71.12532, lat = 42.3712)
Map$setZoom(zoom = 12)
Map$addLayer(
eeObject = composite,
visParams = list(),
name = "max value composite"
)# Load four 2012 NAIP quarter quads, different locations.
naip2012 <- ee$ImageCollection("USDA/NAIP/DOQQ")$
filterBounds(ee$Geometry$Rectangle(-71.17965, 42.35125, -71.08824, 42.40584))$
filterDate("2012-01-01", "2012-12-31")
# Spatially mosaic the images in the collection and display.
mosaic <- naip2012$mosaic()
Map$addLayer(
eeObject = mosaic,
visParams = list(),
name = "spatial mosaic"
)# Load a NAIP quarter quad, display.
naip <- ee$Image("USDA/NAIP/DOQQ/m_4207148_nw_19_1_20120710")
Map$setCenter(lon = -71.0915, lat = 42.3443)
Map$setZoom(zoom = 14)
Map$addLayer(
eeObject = naip,
visParams = list(),
name = "NAIP DOQQ"
)# Create the NDVI and NDWI spectral indices.
ndvi <- naip$normalizedDifference(c("N", "R"))
ndwi <- naip$normalizedDifference(c("G", "N"))
# Create some binary images from thresholds on the indices.
# This threshold is designed to detect bare land.
bare1 <- ndvi$lt(0.2)$And(ndwi$lt(0.3))
# This detects bare land with lower sensitivity. It also detects shadows.
bare2 <- ndvi$lt(0.2)$And(ndwi$lt(0.8))
# Define visualization parameters for the spectral indices.
ndviViz <- list(min = -1, max = 1, palette = c("FF0000", "00FF00"))
ndwiViz <- list(min = 0.5, max = 1, palette = c("00FFFF", "0000FF"))
# Mask and mosaic visualization images. The last layer is on top.
mosaic <- ee$ImageCollection(list(
# NDWI > 0.5 is water. Visualize it with a blue palette.
ndwi$updateMask(ndwi$gte(0.5))$visualize(ndwiViz),
# NDVI > 0.2 is vegetation. Visualize it with a green palette.
ndvi$updateMask(ndvi$gte(0.2))$visualize(ndviViz),
# Visualize bare areas with shadow (bare2 but not bare1) as gray.
bare2$updateMask(bare2$And(bare1$Not()))$visualize(list(palette = c("AAAAAA"))),
# Visualize the other bare areas as white.
bare1$updateMask(bare1)$visualize(list(palette = c("FFFFFF")))
))$mosaic()
Map$addLayer(
eeObject = ndwi,
visParams = list(),
name = "Visualization mosaic"
)Continuará…
5. Geometry, Feature y FeatureCollection
5.2 Geodesic vs. Planar Geometries
# Create a geodesic polygon.
polygon <- ee$Geometry$Polygon(
list(
c(-5, 40), c(65, 40), c(65, 60), c(-5, 60), c(-5, 60)
)
)
# Create a planar polygon.
planarPolygon <- ee$Geometry(polygon, {}, FALSE)
polygon <- ee$FeatureCollection(polygon)
planarPolygon <- ee$FeatureCollection(planarPolygon)
# Display the polygons by adding them to the map.
Map$centerObject(polygon, zoom = 2)## NOTE: Center obtained from the first element.Map$addLayer(polygon, list(color = "FF0000"), "geodesic polygon") +
Map$addLayer(planarPolygon, list(color = "000000"), "planar polygon")5.3 Geometry Visualization and Information
# Create a geodesic polygon.
polygon <- ee$Geometry$Polygon(
list(
c(-5, 40), c(65, 40), c(65, 60), c(-5, 60), c(-5, 60)
)
)
# Create a planar polygon.
planarPolygon <- ee$Geometry(polygon, {}, FALSE)
polygon <- ee$FeatureCollection(polygon)
planarPolygon <- ee$FeatureCollection(planarPolygon)
# Display the polygons by adding them to the map.
Map$centerObject(polygon, zoom = 2)## NOTE: Center obtained from the first element.Map$addLayer(polygon, list(color = "FF0000"), "geodesic polygon") +
Map$addLayer(planarPolygon, list(color = "000000"), "planar polygon")# Create two circular geometries.
poly1 <- ee$Geometry$Point(c(-50, 30))$buffer(1e6)
poly2 <- ee$Geometry$Point(c(-40, 30))$buffer(1e6)
# Display polygon 1 in red and polygon 2 in blue.
Map$setCenter(-45, 30)
Map$addLayer(poly1, list(color = "FF0000"), "poly1") +
Map$addLayer(poly2, list(color = "0000FF"), "poly2")# Compute the intersection, display it in blue.
intersection <- poly1$intersection(poly2, ee$ErrorMargin(1))
Map$addLayer(intersection, list(color = "00FF00"), "intersection")# Compute the union, display it in magenta.
union <- poly1$union(poly2, ee$ErrorMargin(1))
Map$addLayer(union, list(color = "FF00FF"), "union")# Compute the difference, display in yellow.
diff1 <- poly1$difference(poly2, ee$ErrorMargin(1))
Map$addLayer(diff1, list(color = "FFFF00"), "diff1")# Compute symmetric difference, display in black.
symDiff <- poly1$symmetricDifference(poly2, ee$ErrorMargin(1))
Map$addLayer(symDiff, list(color = "000000"), "symmetric difference")# Create a geodesic polygon.
polygon <- ee$Geometry$Polygon(
list(
c(-5, 40), c(65, 40), c(65, 60), c(-5, 60), c(-5, 60)
)
)
# Create a planar polygon.
planarPolygon <- ee$Geometry(polygon, {}, FALSE)
polygon <- ee$FeatureCollection(polygon)
planarPolygon <- ee$FeatureCollection(planarPolygon)
# Display the polygons by adding them to the map.
Map$centerObject(polygon, zoom = 2)## NOTE: Center obtained from the first element.Map$addLayer(polygon, list(color = "FF0000"), "geodesic polygon") +
Map$addLayer(planarPolygon, list(color = "000000"), "planar polygon")library(rgee)
# ee_reattach() # reattach ee as a reserved word
ee_Initialize()## -- rgee 1.0.9 --------------------------------------- earthengine-api 0.1.259 --
## v email: not_defined
## v Initializing Google Earth Engine:
v Initializing Google Earth Engine: DONE!
##
v Earth Engine user: users/tarredwall
## --------------------------------------------------------------------------------# This function gets NDVI from Landsat 8 imagery.
addNDVI <- function(image) {
return(image$addBands(image$normalizedDifference(c("B5", "B4"))))
}
# This function masks cloudy pixels.
cloudMask <- function(image) {
clouds <- ee$Algorithms$Landsat$simpleCloudScore(image)$select("cloud")
return(image$updateMask(clouds$lt(10)))
}
# Load a Landsat collection, map the NDVI and cloud masking functions over it.
collection <- ee$ImageCollection("LANDSAT/LC08/C01/T1_TOA")$
filterBounds(ee$Geometry$Point(c(-122.262, 37.8719)))$
filterDate("2014-03-01", "2014-05-31")$
map(addNDVI)$
map(cloudMask)
# Reduce the collection to the mean of each pixel and display.
meanImage <- collection$reduce(ee$Reducer$mean())
vizParams <- list(
bands = c("B5_mean", "B4_mean", "B3_mean"),
min = 0,
max = 0.5
)
Map$addLayer(
eeObject = meanImage,
visParams = vizParams,
name = "mean"
)# Load a region in which to compute the mean and display it.
counties <- ee$FeatureCollection("TIGER/2016/Counties")
santaClara <- ee$Feature(counties$filter(
ee$Filter$eq("NAME", "Santa Clara")
)$first())
Map$addLayer(
eeObject = santaClara,
visParams = list(palette = "yellow"),
name = "Santa Clara"
)# Get the mean of NDVI in the region.
mean <- meanImage$select("nd_mean")$reduceRegion(
reducer = ee$Reducer$mean(),
geometry = santaClara$geometry(),
scale = 30
)
# Print mean NDVI for the region.
cat("Santa Clara spring mean NDVI:", mean$get("nd_mean")$getInfo())## Santa Clara spring mean NDVI: 0.4650797