I Región de Arica y Parinacota.

1 Lo primero a hacer es dibujar el shp de la Región de Arica y parinacota

region_arica <- st_read("Regiones_separadas/Region_15.shp", 
        quiet = TRUE) %>% 
        sf_as_ee()
arica <- region_arica$geometry()
Map$setCenter(lon = -69.7522, lat = -18.65665, zoom = 7)
Map$addLayer(arica)

2 lo segundo es asociar una imagen satelital Sentinel filtrafda por fecha, nubosidad y vinculada a la region de estudio. Y establecemos que el parametro sea la media.

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

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

first <- sentinel1$median()

vizParams <- list(
  bands = c("B8","B5" , "B3"),
  #  bands = c("B2", "B3"),
  min = 100,
  max = 1000,
  gamma = 2
)
Map$setCenter(lon = -69.7522, lat = -18.65665, zoom = 7)
Map$addLayer(first, vizParams, "Landsat 8 image")

3 Lo tercero es construir un raster con NDGI mayor a 0 simpre vinculada a la primera imagen satelital que hemos extraído.

Construimos una imagen raster con una columna llamada nd con dos categorías para cada pixel: 0 si es rojo, 1 si es verde

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

ndgi <- getNDGI(first)

## Warning: Ops.ee.image.Image will be deprecated in rgee v.1.1.0. Please install
## rgeeExtra (https://github.com/r-earthengine/rgeeExtra). Deeply sorry for the
## inconveniences.

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

Map$setCenter(lon = -69.7522, lat = -18.65665, zoom = 7)
Map$addLayer(ndgi, ndgiParams, "NDGI")

4 obtenemos los shp que nos indican la presencia de glaciares en la Region de Arica y Parinacota.

glaciar_arica <- st_read("GlaciaresxRegion/Simple_ARICA_Y_PARINACOTA.shp",
        quiet = TRUE) %>%
        sf_as_ee()
glaciar_arica <- glaciar_arica$geometry()
Map$setCenter(lon = -69.7522, lat = -18.65665, zoom = 7)
Map$addLayer(glaciar_arica, ndgiParams, "NDGI")

5 Intersectamos las regiones glaciares de la region de Arica con nuestro ndvi, y de esta forma obtenemos un universo muestral al cual poder extraer nuestro primer set de muestras (pixeles rojos interiores).

ESTA ES LA IMAGEN SOBRE LA QUE EXTRAEREMOS EL PRIMER SET DE MUESTRAS

glaciares_adm_con_ngvi <- ndgi$clip(glaciar_arica)
Map$setCenter(lon = -69.80584, lat = -17.67837, zoom = 15)
Map$addLayer(glaciares_adm_con_ngvi, ndgiParams, "NDGI")

6 Hay que hacer una diferencia (se le resta al shp de la región de Arica los shps de los glaciares) para intersectarla con nuestra imagen ngdi y obtener el segundo universo muestral.

6.1 Primero obtenemos nuestro shapes sin zonas administrativas glaciares:

arica_sin_glacia <- arica$difference(glaciar_arica, ee$ErrorMargin(1))
Map$setCenter(lon = -69.80584, lat = -17.67837, zoom = 7)
Map$addLayer(arica_sin_glacia, ndgiParams, "NDGI")

6.2 Necesitamos intersectar nuestro ndvi con la zona de Arica y Parinacota que excluye los glaciares.

ndvi_arica_sin_glaciares <- ndgi$clip(arica_sin_glacia)
Map$setCenter(lon = -69.80584, lat = -17.67837, zoom = 7)
Map$addLayer(ndvi_arica_sin_glaciares, ndgiParams, "NDGI")

8 Ahora tenemos dos regiones preparadas sobre las cuales podemos extraer nuestras muestras. A efectos de orden nombrémoslas:

8.1. pixeles rojos interiores

8.2. pixeles verdes exteriores

8.1. pixeles rojos interiores

Necesitamos intersectar las regiones administrativas glaciares con nuestro raster ndvi y de ahi extraer las muestras de color rojo que simbolizaran los glaciares.

limites_glaciares_con_rasters <- ndgi$clip(glaciar_arica)
Map$setCenter(lon = -69.80584, lat = -17.67837, zoom = 15)
Map$addLayer(limites_glaciares_con_rasters , ndgiParams, "NDGI")

Los rojos interiores son glaciares, 0 = rojo y son llamados muestras_in_gee_rojas.

muestras_in <- ndgi$sampleRegions(
    collection = ee$Feature(glaciar_arica),
    scale = 100,
    tileScale = 16,
    geometries = TRUE
  )

muestras_in_gee_rojas = muestras_in$filter(ee$Filter$eq('nd', 0))

Map$setCenter(lon = -69.80584, lat = -17.67837, zoom = 7)
Map$addLayer(
  eeObject =  muestras_in_gee_rojas,
  visParams = {},
  name = "puntos glaciares"
)

8.2. Pixeles verdes exteriores

ndvi_arica_sin_glaciares <- ndgi$clip(arica_sin_glacia)
Map$setCenter(lon = -69.80584, lat = -17.67837, zoom = 7)
Map$addLayer(ndvi_arica_sin_glaciares, ndgiParams, "NDGI")

Los verdes exteriores no son glaciares 1 = verde

muestras_in <- ndgi$sampleRegions(
    collection = ee$Feature(ndvi_arica_sin_glaciares),
    scale = 1000,
    tileScale = 16,
    geometries = TRUE
  )

muestras_out_gee_verdes = muestras_in$filter(ee$Filter$eq('nd', 1))

Map$setCenter(lon = -69.80584, lat = -17.67837, zoom = 7)
Map$addLayer(
  eeObject =  muestras_out_gee_verdes,
  visParams = {},
  name = "puntos no glaciares"
)

III debemos unir muestras y generar el random forest

union_total <- muestras_out_gee_verdes$merge(muestras_in_gee_rojas)

3 Aplicamos Random Forest

(volver al índice)

region_arica <- st_read("Regiones_separadas/Region_15.shp", 
        quiet = TRUE) %>% 
        sf_as_ee()
arica <- region_arica$geometry()

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

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

first <- sentinel1$median()

# Use these bands for prediction.
bands <- c("B2", "B3", "B4", "B8")

# 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 = 10
)

# Train a CART classifier with default parameters.
trained <- ee$Classifier$smileRandomForest(10)$train(training, label, bands)

# Classify the image with the same bands used for training.
classified1 <- 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(classified1, viz_class, name = "classification")

18:29 18:30

Le aplicamos la matriz de confusion:

trainAccuracy = trained$confusionMatrix()
trainAccuracy$getInfo()

## [[1]]
## [1] 23  2
## 
## [[2]]
## [1]     3 17752

18:40 15:43

4 Aplicamos CART

(volver al índice)

region_arica <- st_read("Regiones_separadas/Region_15.shp", 
        quiet = TRUE) %>% 
        sf_as_ee()
arica <- region_arica$geometry()

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

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

first <- sentinel1$median()

# Use these bands for prediction.
bands <- c("B2", "B3", "B4", "B8")

# 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 = 10
)

# Train a CART classifier with default parameters.
trained <- ee$Classifier$smileCart(10)$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")

18:47 18:50

trainAccuracy = trained$confusionMatrix()
trainAccuracy$getInfo()

## [[1]]
## [1] 16  9
## 
## [[2]]
## [1]     2 17753

18:51 18 53

4 Aplicamos SVM

(volver al índice)

region_arica <- st_read("Regiones_separadas/Region_15.shp", 
        quiet = TRUE) %>% 
        sf_as_ee()
arica <- region_arica$geometry()

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

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

first <- sentinel1$median()

# Use these bands for prediction.
bands <- c("B2", "B3", "B4", "B8")

# 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 = 10
)
classifier = ee$Classifier$libsvm(
  kernelType = "RBF",
  gamma = 0.5,
  cost = 10
)
# Train a CART classifier with default parameters.
trained = classifier$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")

No se puede obtener la matriz de confusion para libsvm pues excede limites de memoria

# trainAccuracy = trained$confusionMatrix()
# trainAccuracy$getInfo()

19:12

Hay que automatizar el proceso

se utiliza el reduce

2. Una vez obtenido nuestro raster random forest interceptamos por limites glaciares, Y EN CADA UNO DE ELLOS determinamos la cantidad de pixeles rojos

mask_0 <- st_read("GlaciaresxRegion/Simple_ARICA_Y_PARINACOTA.shp")

## Reading layer `Simple_ARICA_Y_PARINACOTA' from data source `C:\Users\usuario\Desktop\ds\ds_rgee\cart\GlaciaresxRegion\Simple_ARICA_Y_PARINACOTA.shp' using driver `ESRI Shapefile'
## Simple feature collection with 327 features and 55 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -69.81277 ymin: -19.03571 xmax: -68.91922 ymax: -17.64658
## Geodetic CRS:  WGS 84

receptaculo <- data.frame()

for (n in 1:1) {
  
  
  reg_15_2 <-  mask_0 %>% 
filter(COD_GLA == mask_0$COD_GLA[n])
varia <- reg_15_2 %>%
        sf_as_ee()
# convertimos el shp en geometria:
region_0 <- varia$geometry()

sale_4 <- classified1$clip(glaciar_arica)

puntos_r <- sale_4$select('classification')$eq(0)

areas_rojas = puntos_r$reduceRegion(
reducer= ee$Reducer$sum(),
geometry= glaciar_arica,
scale= 30 ,
maxPixels= 1e9
)



puntos_v <- sale_4$select('classification')$eq(1)

areas_verdes = puntos_v$reduceRegion(
reducer= ee$Reducer$sum(),
geometry= glaciar_arica,
scale= 30 ,
maxPixels= 1e9
)

rojas <- areas_rojas$getInfo()
verdes <- areas_verdes$getInfo()

texto <- paste("Código de graciar: ",mask_0$COD_GLA[n], " areas rojas", sep = "")
print(texto)
print(rojas)

texto <- paste("Código de graciar: ",mask_0$COD_GLA[n], " areas verdes", sep = "")
print(texto)
print(verdes)

areas_rv <- cbind(mask_0$COD_GLA[n],rojas,verdes,mask_0$AREA_Km2[n])
cbindg <- cbind(areas_rv[[1]],areas_rv[[2]],areas_rv[[3]],areas_rv[[4]])
cbindg <- as.data.frame(cbindg)
receptaculo <- rbind(receptaculo,cbindg)

}

## [1] "Código de graciar: CL101001010 areas rojas"
## $classification
## [1] 1015.745
## 
## [1] "Código de graciar: CL101001010 areas verdes"
## $classification
## [1] 32202.69

colnames(receptaculo) <- c("codigo_g", "area_rojo", "area_verde", "area_g")


receptaculo

##      codigo_g        area_rojo       area_verde          area_g
## 1 CL101001010 1015.74509803922 32202.6862745098 0.0830087413962

uno

CC

25-05-2021

R

I Región de Arica y Parinacota.

1 Lo primero a hacer es dibujar el shp de la Región de Arica y parinacota

2 lo segundo es asociar una imagen satelital Sentinel filtrafda por fecha, nubosidad y vinculada a la region de estudio. Y establecemos que el parametro sea la media.

3 Lo tercero es construir un raster con NDGI mayor a 0 simpre vinculada a la primera imagen satelital que hemos extraído.

4 obtenemos los shp que nos indican la presencia de glaciares en la Region de Arica y Parinacota.

5 Intersectamos las regiones glaciares de la region de Arica con nuestro ndvi, y de esta forma obtenemos un universo muestral al cual poder extraer nuestro primer set de muestras (pixeles rojos interiores).

6 Hay que hacer una diferencia (se le resta al shp de la región de Arica los shps de los glaciares) para intersectarla con nuestra imagen ngdi y obtener el segundo universo muestral.

6.1 Primero obtenemos nuestro shapes sin zonas administrativas glaciares:

6.2 Necesitamos intersectar nuestro ndvi con la zona de Arica y Parinacota que excluye los glaciares.

8 Ahora tenemos dos regiones preparadas sobre las cuales podemos extraer nuestras muestras. A efectos de orden nombrémoslas:

8.1. pixeles rojos interiores

8.2. Pixeles verdes exteriores

3 Aplicamos Random Forest

4 Aplicamos CART

4 Aplicamos SVM

Hay que automatizar el proceso

2. Una vez obtenido nuestro raster random forest interceptamos por limites glaciares, Y EN CADA UNO DE ELLOS determinamos la cantidad de pixeles rojos