Extracción de datos de GEE 2

1. Cargamos y analizamos el shapefile de glaciares

Leemos el shp:

mask_1 <- st_read("Lim_glaciares.shp")

## Reading layer `Lim_glaciares' from data source `C:\Users\usuario\Desktop\ds\ds_rgee\puntos\Lim_glaciares.shp' using driver `ESRI Shapefile'
## Simple feature collection with 2085 features and 1 field
## Geometry type: POLYGON
## Dimension:     XYZ
## Bounding box:  xmin: -72.95502 ymin: -44.14786 xmax: -71.1558 ymax: -41.39342
## z_range:       zmin: 0 zmax: 0
## Geodetic CRS:  WGS 84

Vemos lo que hay adentro y cuantas filas tiene:

mask_1

## Simple feature collection with 2085 features and 1 field
## Geometry type: POLYGON
## Dimension:     XYZ
## Bounding box:  xmin: -72.95502 ymin: -44.14786 xmax: -71.1558 ymax: -41.39342
## z_range:       zmin: 0 zmax: 0
## Geodetic CRS:  WGS 84
## First 10 features:
##    ID                       geometry
## 1   1 POLYGON Z ((-72.70323 -44.1...
## 2   2 POLYGON Z ((-72.70229 -44.1...
## 3   3 POLYGON Z ((-72.69812 -44.1...
## 4   4 POLYGON Z ((-71.84016 -44.1...
## 5   5 POLYGON Z ((-72.69504 -44.1...
## 6   6 POLYGON Z ((-71.84184 -44.1...
## 7   7 POLYGON Z ((-72.70229 -44.1...
## 8   8 POLYGON Z ((-72.07368 -44.1...
## 9   9 POLYGON Z ((-72.05738 -44.1...
## 10 10 POLYGON Z ((-71.73046 -44.1...

nrow(mask_1)

## [1] 2085

Extraemos los primeros 4 poligonos

mask_100 <- mask_1[c(1:4),]
nrow(mask_100)

## [1] 4

mask_100

## Simple feature collection with 4 features and 1 field
## Geometry type: POLYGON
## Dimension:     XYZ
## Bounding box:  xmin: -72.7038 ymin: -44.14786 xmax: -71.83893 ymax: -44.14276
## z_range:       zmin: 0 zmax: 0
## Geodetic CRS:  WGS 84
##   ID                       geometry
## 1  1 POLYGON Z ((-72.70323 -44.1...
## 2  2 POLYGON Z ((-72.70229 -44.1...
## 3  3 POLYGON Z ((-72.69812 -44.1...
## 4  4 POLYGON Z ((-71.84016 -44.1...

Extraemos el 4 del grupo de 4:

mmm es el cuarto poligono

mmm <- mask_100$geometry[4]
mmm

## Geometry set for 1 feature 
## Geometry type: POLYGON
## Dimension:     XYZ
## Bounding box:  xmin: -71.84109 ymin: -44.14395 xmax: -71.83893 ymax: -44.14276
## z_range:       zmin: 0 zmax: 0
## Geodetic CRS:  WGS 84

## POLYGON Z ((-71.84016 -44.14287 0, -71.83998 -4...

Graficamos: mmm es el cuarto poligono

polig <- ggplot() + geom_sf(data = mmm, size = 0.1, color = "red", fill = "cyan1") + ggtitle("poligono de glaciar") + coord_sf()
polig

2. Intersectemos:

Guardemos una roi de respaldo:

# habitats = ee$FeatureCollection(
#     list(
#       ee$Feature(ee$Geometry$Polygon(
#      list(
#           c(-72.7032312220809, -44.14785800381241),
#           c(-72.7037943383522, -44.1478591101313),
#           c(-72.703796, -44.147822),
#           c(-72.70323299999992, -44.14781),
#           c(-72.7032312220809, -44.1478580038124)
#         )
#       ))
#     )
#   )
# habitats$geometry

El objeto k contendra las geometrias del poligono mmm

mmm_l <- as.list.data.frame(mmm)
lista <- mmm_l[[1]]
lista2 <- as.list.data.frame(lista)
k <- as.list(as.data.frame(lista2))
k

## $X1
##  [1] -71.84016 -71.83998 -71.83991 -71.83984 -71.83981 -71.83974 -71.83967
##  [8] -71.83958 -71.83959 -71.83946 -71.83915 -71.83907 -71.83901 -71.83899
## [15] -71.83905 -71.83901 -71.83893 -71.83893 -71.83894 -71.83905 -71.83924
## [22] -71.83938 -71.83954 -71.83984 -71.84012 -71.84033 -71.84045 -71.84073
## [29] -71.84103 -71.84109 -71.84101 -71.84094 -71.84053 -71.84025 -71.84019
## [36] -71.84016
## 
## $X2
##  [1] -44.14287 -44.14291 -44.14286 -44.14282 -44.14281 -44.14283 -44.14286
##  [8] -44.14293 -44.14306 -44.14320 -44.14320 -44.14327 -44.14339 -44.14341
## [15] -44.14345 -44.14357 -44.14362 -44.14367 -44.14377 -44.14376 -44.14390
## [22] -44.14391 -44.14392 -44.14395 -44.14394 -44.14391 -44.14378 -44.14366
## [29] -44.14349 -44.14330 -44.14312 -44.14285 -44.14279 -44.14276 -44.14280
## [36] -44.14287
## 
## $X3
##  [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

Construimos el objeto habitats con los 8 primeros puntos del poligono

habitats = ee$FeatureCollection(
    list(
      ee$Feature(ee$Geometry$Polygon(
     list(
          c(k$X1[1], k$X2[1]),
          c(k$X1[2], k$X2[2]),
          c(k$X1[3], k$X2[3]),
                    c(k$X1[4], k$X2[4]),
                    c(k$X1[5], k$X2[5]),
                    c(k$X1[6], k$X2[6]),
                    c(k$X1[7], k$X2[7]),         
          c(k$X1[8], k$X2[8])
        )
      ))
    )
  )
habitats$geometry

## <bound method ApiFunction.importApi.<locals>.MakeBoundFunction.<locals>.<lambda> of <ee.featurecollection.FeatureCollection>>

Creamos una coleccion de imagenes satelitales S2

sentinel1 = ee$ImageCollection('COPERNICUS/S2')$filterDate("2019-03-01", "2020-03-01")
sentinel1_mean = sentinel1$median()
  collection <- ee$ImageCollection$fromImages(list(sentinel1_mean))$toBands()
collection

## EarthEngine Object: Image

Lo intersectamos

  mySamples <- collection$sampleRegions(
  collection = ee$Feature(habitats$first()),
    scale= 10,
    geometries=TRUE
  )

y para el conjunto de 8 puntos obtenemos 4 filas de datos de bandas

  mySamples_sf_batch <-  ee_as_sf(mySamples) 
  mySamples_sf_batch <- as.data.frame((mySamples_sf_batch))
  mySamples_sf_batch

##    X0_B1 X0_B10 X0_B11 X0_B12  X0_B2    X0_B3   X0_B4 X0_B5  X0_B6  X0_B7
## 1 6418.0  126.0   2133   1730 5989.0 5513.500 5833.50  6131 6197.5 6100.5
## 2 6418.0  126.0   2133   1730 6109.0 5633.000 5982.75  6131 6197.5 6100.5
## 3 6407.5  103.5   2092   1628 6157.5 5789.667 6079.00  6237 6334.0 6422.0
## 4 6407.5  103.5   2092   1628 6297.5 5833.250 6293.00  6237 6334.0 6422.0
##    X0_B8  X0_B8A X0_B9 X0_QA10 X0_QA20 X0_QA60                   geometry
## 1 5825.0 5966.00  3133       0       0    1024 POINT (-71.84002 -44.1429)
## 2 5813.5 5966.00  3133       0       0    1024 POINT (-71.83985 -44.1429)
## 3 5853.5 6187.75  3044       0       0    1024 POINT (-71.83976 -44.1429)
## 4 6013.5 6187.75  3044       0       0    1024 POINT (-71.83967 -44.1429)

  mySamples_sf_batch$geometry

## Geometry set for 4 features 
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -71.84002 ymin: -44.1429 xmax: -71.83967 ymax: -44.1429
## Geodetic CRS:  WGS 84

## POINT (-71.84002 -44.1429)

## POINT (-71.83985 -44.1429)

## POINT (-71.83976 -44.1429)

## POINT (-71.83967 -44.1429)

lo anterior son valores de un conjunto de registros de variables para el poligono:

polig

3. Calculo del NDGI

calculo_NDGI <-  ( mySamples_sf_batch$X0_B3 - mySamples_sf_batch$X0_B4) / ( mySamples_sf_batch$X0_B3 + mySamples_sf_batch$X0_B4)

calculo_NDGI

## [1] -0.02820129 -0.03010998 -0.02437791 -0.03791362

Intento de automatizacion de la lista:

# h<-""
# secuencilla <- seq(1,36)
# f <- paste("c(",k$X1[secuencilla] , ",", k$X2[secuencilla],")")
# for (i in 1 : 36) {
# h <- paste(f[i],",",h)
# }
# h

# habitats = ee$FeatureCollection(
#     list(
#       ee$Feature(ee$Geometry$Polygon(
#      list(
#           h
#         )
#       ))
#     )
#   )
# 
# 
# habitats$geometry

##################### ##################### #####################

4 Lo que sigue

Generar un indice a partir de 0 < (B3-B4)/ B3+B4 < 0.5

El criterio de clasificacion sería:

El índice a partir de (B3-B4)/ B3+B4 > 0.5 NDGI y dentro de límites del shp de glaciares deben ser los criterios para establecer los márgenes para extraer los pixeles de muestra a entrenar

No glaciar es todo lo que este fuera y < -0.5

y sacar una muestra de 0.01% de pixeles o 5000 registros que correspondan a glaciar? esta muestra es con la que vamos a entrenar el clasificador.

CART: /ramdom forest aca entran las bandas de B1 a B8 como input

4.1 Que extraer?

Necesitamos seleccionar (datos) raster con columnas informativas de las imágenes satelitales para cruzar con otra data.

¿A quién le podría ser de interés esta información?

a los productores que requieren ciertos parámetros para cosechar en forma óptima sus productos, por ejemplo:

1 un índice de interés sería la precipitación para quienes cultivan arroz.

Variables que incidan en la compra y venta de un bien inmueble, por ejemplo:

2 sería de interés medir la temperatura en las áreas urbanas para agregarlo como parámetro en la decisión de comprar una casa.

Hacer esto con datos climáticos y datos categóricos:

1 Índices de vegetación: Data Agro
2 Propiedades del suelo: Omar.
2 Temperaturas, precipitación y evapotranspiración

4.2 Cómo procesar la información extraída?

Extraer estadísticas asociadas a un punto: min max media sd.

Extraer promedios de valores incluídos dentro de un campo (las mismas métricas dentro de un polígono)