R y GEE
El paquete rgee
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.
Podemos imprimir al estilo R tradicional a al estilo Python de GEE:
print("Hola dataintelligence!")## [1] "Hola dataintelligence!"ee$String("Hola dataintelligence!")$getInfo()## [1] "Hola dataintelligence!"Podemos imprimir con pipes al estilo python:
ee$String("Hola dataintelligence!") %>% ee$String$getInfo()## [1] "Hola dataintelligence!"Ahora extraemos parte de la información almacenada en una imagen:
# a_003 <- ee$Image("LANDSAT/LC08/C01/T1/LC08_044034_20140318") %>% ee$Image$getInfo()
# f <- last(a_003)
# g <- head(f,8)
# saveRDS(g,"g.rds")
g <- readRDS("g.rds")
g## $RADIANCE_MULT_BAND_5
## [1] 0.0061709
##
## $RADIANCE_MULT_BAND_6
## [1] 0.0015346
##
## $RADIANCE_MULT_BAND_3
## [1] 0.011958
##
## $RADIANCE_MULT_BAND_4
## [1] 0.010084
##
## $RADIANCE_MULT_BAND_1
## [1] 0.012673
##
## $RADIANCE_MULT_BAND_2
## [1] 0.012977
##
## $K2_CONSTANT_BAND_11
## [1] 1201.144
##
## $K2_CONSTANT_BAND_10
## [1] 1321.079Vamos a cargar y desplegar una imagen Landsat de la ciudad de Bogotá:
image <- ee$Image("LANDSAT/LC08/C01/T1/LC08_008057_20180317")
Map$centerObject(image)
Map$addLayer(image, name = "Landsat 8 original image")Agregamos como capas los parámetros de visualización:
vizParams <- list(
bands = c("B5", "B4", "B3"),
min = 5000, max = 15000, gamma = 1.3
)
Map$addLayer(image, vizParams, "Landsat 8 False color")Haremos lo mismo para una imagen Landsat del estado de California añadiendo información del estado de la FeatureCollection TIGER:
image <- ee$Image("LANDSAT/LC08/C01/T1/LC08_044034_20140318")
Map$centerObject(image)
Map$addLayer(image, name = "Landsat 8 original image")vizParams <- list(
bands = c("B5", "B4", "B3"),
min = 5000, max = 15000, gamma = 1.3
)
Map$addLayer(image, vizParams, "Landsat 8 False color")counties <- ee$FeatureCollection("TIGER/2016/Counties")
Map$addLayer(
eeObject = counties,
visParams = vizParams,
name = "counties"
)Vamos a cargar un ImageCollection LANDSAT/LC08/C01/T1 con ciertos parámetros filtrado con una cobertura de nubes. Filtraremos nuevamente con los datos geográficos de TIGER.
collection <- ee$ImageCollection("LANDSAT/LC08/C01/T1")
point <- ee$Geometry$Point(-122.262, 37.8719)
start <- ee$Date("2014-06-01")
finish <- ee$Date("2014-10-01")
filteredCollection <- ee$ImageCollection("LANDSAT/LC08/C01/T1")$
filterBounds(point)$
filterDate(start, finish)$
sort("CLOUD_COVER", TRUE)
first <- filteredCollection$first()
# Define visualization parameters in an object literal.
vizParams <- list(
bands = c("B5", "B4", "B3"),
min = 5000,
max = 15000,
gamma = 1.3
)
Map$addLayer(first, vizParams, "Landsat 8 image")featureCollection <- ee$FeatureCollection("TIGER/2016/States")
filteredFC <- featureCollection$filter(ee$Filter$eq("NAME", "California"))
Map$addLayer(
eeObject = filteredFC,
visParams = list(palette = "red"),
name = "California"
)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")Podemos recorrer todos los datos contenidos en una ImageCollection sin necesidad de usar un ciclo for:
library(rgee)
library(kableExtra)## Warning: package 'kableExtra' was built under R version 4.0.3ee_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 |
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"
)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.119Aunque 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"
)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")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"
)# 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"
)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.# 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"
)# # 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))# 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."
)# 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"
)# 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"
)# 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"
)# 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"
)# 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"
)# Load a Landsat 5 image and select the bands we want to unmix.
bands <- c("B1", "B2", "B3", "B4", "B5", "B6", "B7")
image <- ee$Image("LANDSAT/LT05/C01/T1/LT05_044034_20080214")$select(bands)
# Define spectral endmembers.
urban <- c(88, 42, 48, 38, 86, 115, 59)
veg <- c(50, 21, 20, 35, 50, 110, 23)
water <- c(51, 20, 14, 9, 7, 116, 4)
# Unmix the image.
fractions <- image$unmix(list(urban, veg, water))
Map$setCenter(lon = -122.1899, lat = 37.5010) # San Francisco Bay
Map$setZoom(zoom = 10)
Map$addLayer(
eeObject = image,
visParams = list(min = 0, max = 128, bands = c("B4", "B3", "B2")),
name = "image"
) +
Map$addLayer(
eeObject = fractions,
visParams = list(min = 0, max = 2),
name = "unmixed"
)# Load a high-resolution NAIP image.
image = ee$Image('USDA/NAIP/DOQQ/m_3712213_sw_10_1_20140613')
# Zoom to San Francisco, display.
Map$setCenter(lon = -122.466123, lat = 37.769833) # San Francisco Bay
Map$setZoom(zoom = 17)
Map$addLayer(eeObject = image,
visParams = list(min=0, max=255),
name = 'image')# Get the NIR band.
nir = image$select('N')
# Define a neighborhood with a kernel.
square = ee$Kernel$square(radius=4)
# Compute entropy and display.
entropy = nir$entropy(square)
Map$addLayer(eeObject = entropy,
visParams = list(min=1, max=5, palette=c('0000CC', 'CC0000')),
name = 'entropy')# Compute the gray-level co-occurrence matrix (GLCM), get contrast.
glcm = nir$glcmTexture(size=4)
contrast = glcm$select('N_contrast')
Map$addLayer(eeObject = contrast,
visParams = list(min=0, max=1500, palette=c('0000CC', 'CC0000')),
name = 'contrast')# Create a list of weights for a 9x9 kernel.
list = list(1, 1, 1, 1, 1, 1, 1, 1, 1)
# The center of the kernel is zero.
centerList = list(1, 1, 1, 1, 0, 1, 1, 1, 1)
# Assemble a list of lists: the 9x9 kernel weights as a 2-D matrix.
lists = list(list, list, list, list, centerList, list, list, list, list)
# Create the kernel from the weights.
# Non-zero weights represent the spatial neighborhood.
kernel = ee$Kernel$fixed(9, 9, lists, -4, -4, FALSE)
# Convert the neighborhood into multiple bands.
neighs = nir$neighborhoodToBands(kernel)
# Compute local Geary's C, a measure of spatial association.
gearys = nir$subtract(neighs)$pow(2)$reduce(ee$Reducer$sum())$divide(9^2)
Map$addLayer(eeObject = gearys,
visParams = list(min=20, max=2500, palette=c('0000CC', 'CC0000')),
name = "Geary's C")# Make an area of interest geometry centered on San Francisco.
point <- ee$Geometry$Point(-122.1899, 37.5010)
aoi <- point$buffer(10000)
# Import a Landsat 8 image, subset the thermal band, and clip to the
# area of interest.
kelvin <- ee$Image("LANDSAT/LC08/C01/T1_TOA/LC08_044034_20140318")$select(list("B10"), list("kelvin"))$clip(aoi)
# Display the thermal band.
Map$setCenter(lon = -122.1899, lat = 37.5010)
Map$setZoom(zoom = 13)
Map$addLayer(
eeObject = kelvin,
visParams = list(min = 288, max = 305),
name = "Kelvin"
)# Threshold the thermal band to set hot pixels as value 1 and not as 0.
hotspots <- kelvin$gt(303)$selfMask()$rename("hotspots")
# Display the thermal hotspots on the Map.
Map$addLayer(
eeObject = hotspots,
visParams = list(min = 288, max = 305, palette = "FF0000"),
name = "Hotspots"
)# Uniquely label the hotspot image objects.
objectId <- hotspots$connectedComponents(
connectedness = ee$Kernel$plus(1),
maxSize = 128
)
# Display the uniquely ID'ed objects to the Map.
Map$addLayer(
eeObject = objectId$randomVisualizer(),
visParams = list(),
name = "Objects"
)# Compute the number of pixels in each object defined by the "labels" band.
objectSize <- objectId$select("labels")$connectedPixelCount(
maxSize = 128, eightConnected = FALSE
)
# Display object pixel count to the Map.
Map$addLayer(
eeObject = objectSize,
visParams = list(),
name = "Object n pixels"
)# Get a pixel area image.
pixelArea <- ee$Image$pixelArea()
# Multiply pixel area by the number of pixels in an object to calculate
# the object area. The result is an image where each pixel
# of an object relates the area of the object in m^2.
objectArea <- objectSize$multiply(pixelArea)
# Display object area to the Map.
Map$addLayer(
eeObject = objectArea,
visParams = list(),
name = "Object area m^2"
)# Threshold the `objectArea` image to define a mask that will mask out
# objects below a given size (1 hectare in this case).
areaMask <- objectArea$gte(10000)
# Update the mask of the `objectId` layer defined previously using the
# minimum area mask just defined.
objectId <- objectId$updateMask(areaMask)
Map$addLayer(
eeObject = objectId,
visParams = list(),
name = "Large hotspots"
)# Make a suitable image for `reduceConnectedComponents()` by adding a label
# band to the `kelvin` temperature image.
kelvin <- kelvin$addBands(objectId$select("labels"))
# Calculate the mean temperature per object defined by the previously added
# "labels" band.
patchTemp <- kelvin$reduceConnectedComponents(
reducer = ee$Reducer$mean(),
labelBand = "labels"
)
# Display object mean temperature to the Map.
Map$addLayer(
eeObject = patchTemp,
visParams = list(min = 303, max = 304, palette = c("yellow", "red")),
name = "Mean temperature"
)# A rectangle representing Bangui, Central African Republic.
geometry <- ee$Geometry$Rectangle(list(18.5229, 4.3491, 18.5833, 4.4066))
# Create a source image where the geometry is 1, everything else is 0.
sources <- ee$Image()$toByte()$paint(geometry, 1)
# Mask the sources image with itself.
sources <- sources$updateMask(sources)
# The cost data is generated from classes in ESA/GLOBCOVER.
cover <- ee$Image("ESA/GLOBCOVER_L4_200901_200912_V2_3")$select(0)
# Classes 60, 80, 110, 140 have cost 1.
# Classes 40, 90, 120, 130, 170 have cost 2.
# Classes 50, 70, 150, 160 have cost 3.
cost <- cover$eq(60)$Or(cover$eq(80))$Or(cover$eq(110))$Or(cover$eq(140))$
multiply(1)$add(
cover$eq(40)$Or(cover$eq(90))$Or(cover$eq(120))$Or(cover$eq(130))$
Or(cover$eq(170))$
multiply(2)$add(
cover$eq(50)$Or(cover$eq(70))$Or(cover$eq(150))$Or(cover$eq(160))$
multiply(3)
)
)
# Compute the cumulative cost to traverse the lAnd cover.
cumulativeCost <- cost$cumulativeCost(
source = sources,
maxDistance = 80 * 1000 # 80 kilometers
)
# Display the results
Map$setCenter(lon = 18.71, lat = 4.2)
Map$setZoom(zoom = 9)
Map$addLayer(
eeObject = cover,
visParams = list(),
name = "Globcover"
) +
Map$addLayer(
eeObject = cumulativeCost,
visParams = list(min = 0, max = 5e4),
name = "accumulated cost"
) +
Map$addLayer(
eeObject = geometry,
visParams = list(color = "FF0000"),
name = "source geometry"
)# Load the two images to be registered.
image1 <- ee$Image("SKYSAT/GEN-A/PUBLIC/ORTHO/MULTISPECTRAL/s01_20150502T082736Z")
image2 <- ee$Image("SKYSAT/GEN-A/PUBLIC/ORTHO/MULTISPECTRAL/s01_20150305T081019Z")
# Use bicubic resampling during registration.
image1Orig <- image1$resample("bicubic")
image2Orig <- image2$resample("bicubic")
# Choose to register using only the 'R' bAnd.
image1RedBAnd <- image1Orig$select("R")
image2RedBAnd <- image2Orig$select("R")
# Determine the displacement by matching only the 'R' bAnds.
displacement <- image2RedBAnd$displacement(
referenceImage = image1RedBAnd,
maxOffset = 50.0,
patchWidth = 100.0
)
# Compute image offset And direction.
offset <- displacement$select("dx")$hypot(displacement$select("dy"))
angle <- displacement$select("dx")$atan2(displacement$select("dy"))
# Display offset distance And angle.
Map$setCenter(lon = 37.44, lat = 0.58)
Map$setZoom(zoom = 15)
Map$addLayer(
eeObject = offset,
visParams = list(min = 0, max = 20),
name = "offset"
) +
Map$addLayer(
eeObject = angle,
visParams = list(min = -pi, max = pi),
name = "angle"
)# Use the computed displacement to register all Original bAnds.
registered <- image2Orig$displace(displacement)
# Show the results of co-registering the images.
visParams <- list(bands = c("R", "G", "B"), max = 4000)
Map$addLayer(
eeObject = image1Orig,
visParams = visParams,
name = "Reference"
) +
Map$addLayer(
eeObject = image2Orig,
visParams = visParams,
name = "BefOre Registration"
) +
Map$addLayer(
eeObject = registered,
visParams = visParams,
name = "After Registration"
)alsoRegistered <- image2Orig$register(
referenceImage = image1Orig,
maxOffset = 50.0,
patchWidth = 100.0
)
Map$addLayer(
eeObject = alsoRegistered,
visParams = visParams,
name = "Also Registered"
)image <- ee$Image("USDA/NAIP/DOQQ/m_3712213_sw_10_1_20140613")
Map$setCenter(-122.466123, 37.769833, 17)
Map$addLayer(image, list(bands = c("N", "R", "G")), "NAIP")geometry <- image$geometry()
# py_help(image$reduceRegions)
means <- image$reduceRegions(
collection = geometry,
reducer = ee$Reducer$mean()$forEachBand(image),
scale = 10
)
print(means$getInfo())## $type
## [1] "FeatureCollection"
##
## $columns
## $columns$B
## [1] "Float<0.0, 255.0>"
##
## $columns$G
## [1] "Float<0.0, 255.0>"
##
## $columns$N
## [1] "Float<0.0, 255.0>"
##
## $columns$R
## [1] "Float<0.0, 255.0>"
##
## $columns$`system:index`
## [1] "String"
##
##
## $features
## $features[[1]]
## $features[[1]]$type
## [1] "Feature"
##
## $features[[1]]$geometry
## $features[[1]]$geometry$type
## [1] "Polygon"
##
## $features[[1]]$geometry$coordinates
## $features[[1]]$geometry$coordinates[[1]]
## $features[[1]]$geometry$coordinates[[1]][[1]]
## [1] -122.50385 37.81519
##
## $features[[1]]$geometry$coordinates[[1]][[2]]
## [1] -122.5038 37.7473
##
## $features[[1]]$geometry$coordinates[[1]][[3]]
## [1] -122.4336 37.7473
##
## $features[[1]]$geometry$coordinates[[1]][[4]]
## [1] -122.43358 37.81519
##
## $features[[1]]$geometry$coordinates[[1]][[5]]
## [1] -122.50385 37.81519
##
##
##
##
## $features[[1]]$id
## [1] "0"
##
## $features[[1]]$properties
## $features[[1]]$properties$B
## [1] 113.9226
##
## $features[[1]]$properties$G
## [1] 125.4139
##
## $features[[1]]$properties$N
## [1] 86.16627
##
## $features[[1]]$properties$R
## [1] 121.3433# 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
)
# Define visualization parameters in an object literal.
vizParams <- list(
bands = c("B5", "B4", "B3"),
min = 0,
max = 1,
gamma = 1.3
)
Map$centerObject(points, 10)## NOTE: Center obtained from the first element.Map$addLayer(image, vizParams, "Image") +
Map$addLayer(points, list(color = "yellow"), "Training points")ee_print(training)## ---------------------------------------------- Earth Engine FeatureCollection --
## FeatureCollection Metadata:
## - Class : ee$FeatureCollection
## - Number of Features : 98
## - Number of Properties : 1
## Feature Metadata:
## - Number of Properties : 9
## Geometry Metadata (f_index = 0):
## - CRS : WGS 84 (EPSG:4326)
## - proj4string : +proj=longlat +datum=WGS84 +no_defs
## - Geotransform : 1 0 0 0 1 0
## - Geodesic : TRUE
## - Geo Type :
## --------------------------------------------------------------------------------continuará…
# Create arbitrary constant images.
constant1 <- ee$Image(1)
constant2 <- ee$Image(2)
# Create a collection by giving a list to the constructor.
collectionFromConstructor <- ee$ImageCollection(list(constant1, constant2))
collectionFromConstructor$getInfo()## $type
## [1] "ImageCollection"
##
## $bands
## list()
##
## $features
## $features[[1]]
## $features[[1]]$type
## [1] "Image"
##
## $features[[1]]$bands
## $features[[1]]$bands[[1]]
## $features[[1]]$bands[[1]]$id
## [1] "constant"
##
## $features[[1]]$bands[[1]]$data_type
## $features[[1]]$bands[[1]]$data_type$type
## [1] "PixelType"
##
## $features[[1]]$bands[[1]]$data_type$precision
## [1] "int"
##
## $features[[1]]$bands[[1]]$data_type$min
## [1] 1
##
## $features[[1]]$bands[[1]]$data_type$max
## [1] 1
##
##
## $features[[1]]$bands[[1]]$crs
## [1] "EPSG:4326"
##
## $features[[1]]$bands[[1]]$crs_transform
## [1] 1 0 0 0 1 0
##
##
##
## $features[[1]]$properties
## $features[[1]]$properties$`system:index`
## [1] "0"
##
##
##
## $features[[2]]
## $features[[2]]$type
## [1] "Image"
##
## $features[[2]]$bands
## $features[[2]]$bands[[1]]
## $features[[2]]$bands[[1]]$id
## [1] "constant"
##
## $features[[2]]$bands[[1]]$data_type
## $features[[2]]$bands[[1]]$data_type$type
## [1] "PixelType"
##
## $features[[2]]$bands[[1]]$data_type$precision
## [1] "int"
##
## $features[[2]]$bands[[1]]$data_type$min
## [1] 2
##
## $features[[2]]$bands[[1]]$data_type$max
## [1] 2
##
##
## $features[[2]]$bands[[1]]$crs
## [1] "EPSG:4326"
##
## $features[[2]]$bands[[1]]$crs_transform
## [1] 1 0 0 0 1 0
##
##
##
## $features[[2]]$properties
## $features[[2]]$properties$`system:index`
## [1] "1"# Create a collection with fromImages().
collectionFromImages <- ee$ImageCollection$fromImages(
list(
ee$Image(3),
ee$Image(4)
)
)
collectionFromImages$getInfo()## $type
## [1] "ImageCollection"
##
## $bands
## list()
##
## $features
## $features[[1]]
## $features[[1]]$type
## [1] "Image"
##
## $features[[1]]$bands
## $features[[1]]$bands[[1]]
## $features[[1]]$bands[[1]]$id
## [1] "constant"
##
## $features[[1]]$bands[[1]]$data_type
## $features[[1]]$bands[[1]]$data_type$type
## [1] "PixelType"
##
## $features[[1]]$bands[[1]]$data_type$precision
## [1] "int"
##
## $features[[1]]$bands[[1]]$data_type$min
## [1] 3
##
## $features[[1]]$bands[[1]]$data_type$max
## [1] 3
##
##
## $features[[1]]$bands[[1]]$crs
## [1] "EPSG:4326"
##
## $features[[1]]$bands[[1]]$crs_transform
## [1] 1 0 0 0 1 0
##
##
##
## $features[[1]]$properties
## $features[[1]]$properties$`system:index`
## [1] "0"
##
##
##
## $features[[2]]
## $features[[2]]$type
## [1] "Image"
##
## $features[[2]]$bands
## $features[[2]]$bands[[1]]
## $features[[2]]$bands[[1]]$id
## [1] "constant"
##
## $features[[2]]$bands[[1]]$data_type
## $features[[2]]$bands[[1]]$data_type$type
## [1] "PixelType"
##
## $features[[2]]$bands[[1]]$data_type$precision
## [1] "int"
##
## $features[[2]]$bands[[1]]$data_type$min
## [1] 4
##
## $features[[2]]$bands[[1]]$data_type$max
## [1] 4
##
##
## $features[[2]]$bands[[1]]$crs
## [1] "EPSG:4326"
##
## $features[[2]]$bands[[1]]$crs_transform
## [1] 1 0 0 0 1 0
##
##
##
## $features[[2]]$properties
## $features[[2]]$properties$`system:index`
## [1] "1"# Merge two collections.
mergedCollection <- collectionFromConstructor$merge(collectionFromImages)
mergedCollection$getInfo()## $type
## [1] "ImageCollection"
##
## $bands
## list()
##
## $features
## $features[[1]]
## $features[[1]]$type
## [1] "Image"
##
## $features[[1]]$bands
## $features[[1]]$bands[[1]]
## $features[[1]]$bands[[1]]$id
## [1] "constant"
##
## $features[[1]]$bands[[1]]$data_type
## $features[[1]]$bands[[1]]$data_type$type
## [1] "PixelType"
##
## $features[[1]]$bands[[1]]$data_type$precision
## [1] "int"
##
## $features[[1]]$bands[[1]]$data_type$min
## [1] 1
##
## $features[[1]]$bands[[1]]$data_type$max
## [1] 1
##
##
## $features[[1]]$bands[[1]]$crs
## [1] "EPSG:4326"
##
## $features[[1]]$bands[[1]]$crs_transform
## [1] 1 0 0 0 1 0
##
##
##
## $features[[1]]$properties
## $features[[1]]$properties$`system:index`
## [1] "1_0"
##
##
##
## $features[[2]]
## $features[[2]]$type
## [1] "Image"
##
## $features[[2]]$bands
## $features[[2]]$bands[[1]]
## $features[[2]]$bands[[1]]$id
## [1] "constant"
##
## $features[[2]]$bands[[1]]$data_type
## $features[[2]]$bands[[1]]$data_type$type
## [1] "PixelType"
##
## $features[[2]]$bands[[1]]$data_type$precision
## [1] "int"
##
## $features[[2]]$bands[[1]]$data_type$min
## [1] 2
##
## $features[[2]]$bands[[1]]$data_type$max
## [1] 2
##
##
## $features[[2]]$bands[[1]]$crs
## [1] "EPSG:4326"
##
## $features[[2]]$bands[[1]]$crs_transform
## [1] 1 0 0 0 1 0
##
##
##
## $features[[2]]$properties
## $features[[2]]$properties$`system:index`
## [1] "1_1"
##
##
##
## $features[[3]]
## $features[[3]]$type
## [1] "Image"
##
## $features[[3]]$bands
## $features[[3]]$bands[[1]]
## $features[[3]]$bands[[1]]$id
## [1] "constant"
##
## $features[[3]]$bands[[1]]$data_type
## $features[[3]]$bands[[1]]$data_type$type
## [1] "PixelType"
##
## $features[[3]]$bands[[1]]$data_type$precision
## [1] "int"
##
## $features[[3]]$bands[[1]]$data_type$min
## [1] 3
##
## $features[[3]]$bands[[1]]$data_type$max
## [1] 3
##
##
## $features[[3]]$bands[[1]]$crs
## [1] "EPSG:4326"
##
## $features[[3]]$bands[[1]]$crs_transform
## [1] 1 0 0 0 1 0
##
##
##
## $features[[3]]$properties
## $features[[3]]$properties$`system:index`
## [1] "2_0"
##
##
##
## $features[[4]]
## $features[[4]]$type
## [1] "Image"
##
## $features[[4]]$bands
## $features[[4]]$bands[[1]]
## $features[[4]]$bands[[1]]$id
## [1] "constant"
##
## $features[[4]]$bands[[1]]$data_type
## $features[[4]]$bands[[1]]$data_type$type
## [1] "PixelType"
##
## $features[[4]]$bands[[1]]$data_type$precision
## [1] "int"
##
## $features[[4]]$bands[[1]]$data_type$min
## [1] 4
##
## $features[[4]]$bands[[1]]$data_type$max
## [1] 4
##
##
## $features[[4]]$bands[[1]]$crs
## [1] "EPSG:4326"
##
## $features[[4]]$bands[[1]]$crs_transform
## [1] 1 0 0 0 1 0
##
##
##
## $features[[4]]$properties
## $features[[4]]$properties$`system:index`
## [1] "2_1"# # 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()# 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"
)# # 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())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"
)# 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á…
# 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 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