Glaciares en Chile

III Parte: La muestra y el mecanismo de Random Forest

DataIntelligence

Abstract

Podemos extraer muestras, como FeatureCollections de rasters, filtrarlas y agruparlas. Podemos llevar la escala del proceso de muestreo con sampleRegions a un mínimo de 400. Sabemos que para Sentinel 100 no está mal.

Se nos pide poder hacerlo a escala 10 y extraer muestras de ésta selección.

Índice

1 Primer conjunto de muestras

2 Segundo conjunto de muestras

3 Random Forest

Primer conjunto de muestras

(volver al índice)

Extraemos la imagen satelital Sentinel filtrada por la region de los lagos

mask <- st_read("region_los_lagos.shp", 
        quiet = TRUE) %>% 
        sf_as_ee()

# convertimos el shp en geometria:
region <- mask$geometry()

start <- ee$Date("2019-03-01")
finish <- ee$Date("2020-03-01")
cc <- 20

sentinel1 = ee$ImageCollection('COPERNICUS/S2')$filterDate(start, finish)$filterBounds(region)$
  filter(ee$Filter$lt("CLOUDY_PIXEL_PERCENTAGE", cc))

first <- sentinel1$median()

# definimos los parametros de visualizacion:
vizParams <- list(
  bands = c("B8","B5" , "B3"),
  #  bands = c("B2", "B3"),
  min = 100,
  max = 1000,
  gamma = 2
)
Map$setCenter(lon = -73.079, lat = -42.611, zoom = 6)
Map$addLayer(first, vizParams, "Landsat 8 image")

construimos un raster con NDGI mayor a 0

getNDGI <- function(image) 
{
  #image$normalizedDifference(c("B3", "B4")) 
   image$normalizedDifference(c("B3", "B4")) > 0
}

ndgi_01 <- getNDGI(first)

ndgiParams <- list(palette = c(
  "#d73027", "#f46d43", "#fdae61",
  "#fee08b", "#d9ef8b", "#a6d96a",
  "#66bd63", "#1a9850"
))

Map$setCenter(lon = -73.079, lat = -42.611, zoom = 6)
Map$addLayer(ndgi_01, ndgiParams, "NDGI")

construimos un cuadrado con el que intersectaremos nuestro raster

box_1 <- ee$Geometry$Rectangle(coords = c(-74.83774,-44.06706,-71.58093,-40.23878),
                             proj = "EPSG:4326",
                             geodesic = FALSE)
Map$addLayer(box_1)

intersectamos

sale_2 <- ndgi_01$clip(box_1)
Map$addLayer(sale_2, ndgiParams, "NDGI")

intersectamos con los limites politicos de los glaciares

mask_0 <- st_read("Lim_glaciares_v2_Simplify.shp",
        quiet = TRUE) %>%
        sf_as_ee()

# convertimos el shp en geometria:
region_0 <- mask_0$geometry()

sale_4 <- sale_2$clip(region_0)
Map$setCenter(lon = -73.079, lat = -42.611, zoom = 7)
Map$addLayer(sale_4, ndgiParams, "NDGI")

Extraemos muestras a una escala de 150 mts

 mySamples <- sale_4$sampleRegions(
    collection = ee$Feature(region_0),
    scale = 400,
    tileScale = 16,
    geometries = TRUE
  )

filtramos para seleccionar solo los pixeles en rojo:

largeSheds_1 = mySamples$filter(ee$Filter$eq('nd', 0))

#ee_help(ee$Filter)

2 Segundo conjunto de muestras

(volver al índice)

creamos un rectangulo al que le vamos a extraer los limites politicos glaciares

box_2 <- ee$Geometry$Rectangle(coords = c(-73,-43,-72.2,-42),
                             proj = "EPSG:4326",
                             geodesic = FALSE)
Map$setCenter(lon = -73.079, lat = -42.611, zoom = 7)
Map$addLayer(box_2)

diff1 <- box_2$difference(region_0, ee$ErrorMargin(1))
Map$addLayer(diff1, list(color = "FFFF00"), "diff1")

intersectaremos nuestro raster ndgi

sale_3 <- sale_2$clip(diff1)
Map$setCenter(lon = -73.079, lat = -42.611, zoom = 7)
Map$addLayer(sale_3, ndgiParams, "NDGI")

extraemos el segundo set de muestras

 mySamples_2 <- sale_3$sampleRegions(
    collection = ee$Feature(diff1),
    scale= 400,
    tileScale=16,
    geometries = TRUE
  )

Filtramos para obtener solo los pixeles verdes:

largeSheds_2 = mySamples_2$filter(ee$Filter$eq('nd', 1))

Unimos ambos sets:

union_total <- largeSheds_1$merge(largeSheds_2)

3 Aplicamos Random Forest

(volver al índice)

mask <- st_read("region_los_lagos.shp", 
        quiet = TRUE) %>% 
        sf_as_ee()

# convertimos el shp en geometria:
region <- mask$geometry()

start <- ee$Date("2019-03-01")
finish <- ee$Date("2020-03-01")
cc <- 20

sentinel1 = ee$ImageCollection('COPERNICUS/S2')$filterDate(start, finish)$filterBounds(region)$
  filter(ee$Filter$lt("CLOUDY_PIXEL_PERCENTAGE", cc))

first <- sentinel1$median()

# Use these bands for prediction.
bands <- c("B2", "B3", "B4", "B5", "B6", "B7", "B10", "B11")

# This property of the table stores the land cover labels.
label <- "nd"

# Overlay the points on the imagery to get training.
training <- first$select(bands)$sampleRegions(
  collection = union_total,
  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 <- first$select(bands)$classify(trained)

# Viz parameters.
# viz_img <- list(bands = c("B4", "B3", "B2"), max = 0.4)
vizParams <- list(
  bands = c("B8","B5" , "B3"),
  #  bands = c("B2", "B3"),
  min = 100,
  max = 1000,
  gamma = 2
)

viz_class <- list(palette = c("red", "green", "blue"), min = 0, max = 2)

Map$addLayer(first, vizParams, name = "image") +
Map$addLayer(classified, viz_class, name = "classification")