Índice

2 Nuevos limites glaciares


library(rgee)
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 
## --------------------------------------------------------------------------------

1 Calculo de Areas

(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")
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)
# obtenenos informacion de las medias de una imagen
# mean <- first$reduceRegion(
#   reducer = ee$Reducer$mean(),
#   geometry = box_2,
#   scale = 30,
#   bestEffort= TRUE
# )
# print(mean$getInfo())

construimos un raster con NDGI mayor a 0

Lo que aca hacemos es construir una imagen raster con una columnallamada nd con dos categorias para cada pixel: 0 si es rojo, 1 si es verde

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

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")

Intersectamos con los límites políticos 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 <- ndgi_01$clip(region_0)
Map$setCenter(lon = -73.079, lat = -42.611, zoom = 7)
Map$addLayer(sale_4, ndgiParams, "NDGI")

Calculating area by classified image in Earth Engine
https://gis.stackexchange.com/questions/313288/calculating-area-by-classified-image-in-earth-engine

Calcular el area de los pixeles rojos del mapa anterior. La imagen posee una columna llamada nd con categorias 0 y 1. 0 son los pixeles rojos.

Para determinar el área de cada clase en una región de interés:

1 definir una región como una geometría: region_0

2 poseer una imagen clasificada (con valores de datos categóricos): ndgi_01

Luego se construye una imagen compuesta donde la primera banda es el área del píxel y la segunda banda tiene la información de la clase.

Luego se aplica un reductor agrupado a la imagen compuesta, para sumar todas las áreas de píxeles dentro de cada clase:

Ejemplo gee:

  roi = ee$Geometry$Point(-113.5911, 37.0855)$buffer(100)
  class_image = ee$Image('USGS/NLCD/NLCD2011')$select('landcover')

  var class_areas = ee.Image.pixelArea().addBands(class_image)
  .reduceRegion({
  reducer: ee.Reducer.sum().group({
  groupField: 1,
  groupName: 'code',
  }),
  geometry: roi,
  scale: 1,  // sample the geometry at 1m intervals
  maxPixels: 1e10
  }).get('groups');
class_areas = sale_4$reduceRegion(
    reducer= ee$Reducer$sum(),
    geometry= region_0,
    scale= 30  ,
    maxPixels= 1e9
  );
class_areas$getInfo()
## $nd
## [1] 123899.7

Encontramos el area de las zonas verdes interiores:{#e0011}

(volver al índice)

verdes_in <- sale_4$select('nd')$eq(1)
# Calcula el área de píxeles en kilómetros cuadrados 
verdes_in_km = verdes_in$multiply(ee$Image$pixelArea())$divide(1000*1000)

verdes_in_km_area = verdes_in_km$reduceRegion(
    reducer= ee$Reducer$sum(),
    geometry= region_0,
    scale= 30  ,
    maxPixels= 1e9
  )

verdes_in_km_area$getInfo()
## $nd
## [1] 81.31329

Encontramos el area de las zonas rojas interiores:{#e0012}

(volver al índice)

rojas_in <- sale_4$select('nd')$eq(0)
# Calcula el área de píxeles en kilómetros cuadrados 
rojas_in_km = rojas_in$multiply(ee$Image$pixelArea())$divide(1000*1000)

rojas_in_km_area = rojas_in_km$reduceRegion(
    reducer= ee$Reducer$sum(),
    geometry= region_0,
    scale= 30  ,
    maxPixels= 1e9
  )

rojas_in_km_area$getInfo()
## $nd
## [1] 1282.863

Obtengamos los mismos resultados pero sin correccion a kilometros:

fff <- sale_4$select('nd')$eq(1)

class_areas = fff$reduceRegion(
    reducer= ee$Reducer$sum(),
    geometry= region_0,
    scale= 30  ,
    maxPixels= 1e9
  )

class_areas$getInfo()
## $nd
## [1] 123899.7
fff <- sale_4$select('nd')$eq(0)

class_areas = fff$reduceRegion(
    reducer= ee$Reducer$sum(),
    geometry= region_0,
    scale= 30  ,
    maxPixels= 1e9
  )

class_areas$getInfo()
## $nd
## [1] 1946999

reduceRegion

Statistics of an Image Region
https://developers.google.com/earth-engine/guides/reducers_reduce_region

Ejemplo gee:

  var meanDictionary = image.reduceRegion({
        reducer: ee.Reducer.mean(),
        geometry: region.geometry(),
        scale: 30,
        maxPixels: 1e9
  });

https://spatialthoughts.com/2020/06/19/calculating-area-gee/