Glaciares en Chile
Glaciares
VE-CC
DataIntelligence23-04-2021
Abstract
Extraemos una roi, que será la región de Los Lagos en forma de un shp, la cual convertimos en un raster con data filtrada del Sentinel. Construímos imágenes con el despliegue de ciertos valores NDGI y procedemos a intersectarlos con otro shape que delimitará los márgenes de los glaciares. Continuamos con el proceso de Random Forest tomando muestras de pixeles, asignándoles las etiquetas correspondientes.
1. Gráfica del shape de la región de los lagos
Tenemos una función que grafica inmediatamente la forma de un shape:
# leemos el shp:
mask <- st_read("region_los_lagos.shp",
quiet = TRUE) %>%
sf_as_ee()# convertimos el shp en geometria:
region <- mask$geometry()Map$setCenter(lon = -73.079, lat = -42.611, zoom = 7)
Map$addLayer(region)Cruzamos con data Sentinel filtrada:
start <- ee$Date("2014-06-01")
finish <- ee$Date("2014-10-01")
cc <- 20
sentinel1 = ee$ImageCollection('COPERNICUS/S2')$filterDate("2019-03-01", "2020-03-01")$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")2. Matemáticas de mapas
Cálculo del Normalized Difference Glacier Index (NDGI):
\[ NDGI = {B_3-B_4 \over B_3+B_4} \]
normalizedDifference calcula la diferencia normalizada entre dos bandas. Es:
\[{B1 - B2 } \over { B1 + B2}\]
Los valores de entrada negativos se fuerzan a 0 para que el resultado se limite al rango (-1, 1).
3. Obtengamos los píxeles del NDGI para tres rangos:
El Índice Glaciar Diferencial Normalizado (NDGI) se utiliza para ayudar a detectar y monitorear glaciares utilizando las bandas espectrales verde y roja. Esta ecuación se utiliza comúnmente en la detección de glaciares y en aplicaciones de monitoreo de glaciares.
- NDGI
getNDGI <- function(image)
{
#image$normalizedDifference(c("B3", "B4"))
image$normalizedDifference(c("B3", "B4"))
}
ndgi <- 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, ndgiParams, "NDGI")- NDGI < 0.4
getNDGI <- function(image)
{
#image$normalizedDifference(c("B3", "B4"))
image$normalizedDifference(c("B3", "B4")) < 0.4
}
ndgi_04 <- 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_04, ndgiParams, "NDGI")- NDGI < 0.1
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")4. Cortamos con geometrías geográficas de 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()
Map$setCenter(lon = -73.079, lat = -42.611, zoom = 7)
Map$addLayer(region_0)# mask_0 <- st_read("cl_lagos_geo/cl_lagos_geoPolygon.shp",
# quiet = TRUE) %>%
# sf_as_ee()
#
# # convertimos el shp en geometria:
# region_0 <- mask_0$geometry()
#
# Map$setCenter(lon = -73.079, lat = -42.611, zoom = 7)
# Map$addLayer(region_0)Tomemos el raster generado con NDGI < 0.1: ngdi_01 en una caja coords = c(-73,-43,-72,-42.5)
box <- ee$Geometry$Rectangle(coords = c(-73,-43,-72,-42.5),
## WGS 84
proj = "EPSG:4326",
geodesic = FALSE)sale <- ndgi_01$clip(box)Map$setCenter(lon = -73.079, lat = -42.611, zoom = 7)
Map$addLayer(sale, ndgiParams, "NDGI")el raster sale contiene el corte de ndgi_01 con la cajita.
tomemos una cajita mas amplia:
box_1 <- ee$Geometry$Rectangle(coords = c(-74.83774,-44.06706,-71.58093,-40.23878),
proj = "EPSG:4326",
geodesic = FALSE)sale_2 <- ndgi_01$clip(box_1)Map$setCenter(lon = -73.079, lat = -42.611, zoom = 7)
Map$addLayer(sale_2, ndgiParams, "NDGI")y ahora intersectamos con el shp de los glaciares de Chile:
sale_3 <- sale_2$clip(region_0)
Map$setCenter(lon = -73.079, lat = -42.611, zoom = 7)
Map$addLayer(sale_3, ndgiParams, "NDGI")Creemos haber cumplido con los pasos 1-6:
- Cargar el shapefile de límite de la zona de estudio https://1drv.ms/u/s!AipQ_fMuljs8jP1boyYhdGH-QUbZQw?e=bgq2lT
- Cargar el shapefile de glaciares https://1drv.ms/u/s!AipQ_fMuljs8jP1c0--utfgc4tCvyQ?e=3zzImb
- Acceder a la colección de imágenes de Sentinel 2 ee.ImageCollection(“COPERNICUS/S2”)
https://developers.google.com/earth-engine/datasets/catalog/COPERNICUS_S2
- Filtrar la colección (colección filtrada):
Filtrar para la zona de estudio: .filterBounds(Limite)
Filtrar por nivel de nubosidad: .filterMetadata(‘CLOUD_COVERAGE_ASSESSMENT’, ‘less_than’, 20)
Filtrar por fecha: ee.Filter.date(‘2019-01-01’, ‘2019-12-31’)
Obtener una sola imagen a partir de la mediana de los píxeles seleccionados cortada con la zona de estudio. colección filtrara. median().clip(Limite)
Calcular el Normalized Difference Glacier Index (NDGI):
NDWI = (B3 – B4) / (B3 + B4)
vamos aqui
5 Seleccionar una muestra de entrenamiento
- NDGI > 0.5 que coinciden con el shapefile de límite de glaciares
“pixeles de referencia presencia”, asignándoles el valor de 1.
- NDGI < -0.5 que NO coinciden con el shapefile de límite de glaciares,
“pixeles de referencia de ausencia”, asignándoles el valor de 0.
6. Obtener una muestra del 0.01% de los píxeles de referencia de presencia y también de los de ausencia, esta será la “muestra de entrenamiento”.
7. Aplicar CART clasificación utilizando la muestra de entrenamiento y las bandas de Sentinel de la 1 a la 8.
CART
Ejercicio ejemplo:
# 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")
points## EarthEngine Object: FeatureCollection# 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")Ejemplo 1: grafica de una caja
box <- ee$Geometry$Rectangle(coords = c(-73,-43,-72.462,-42.811),
## WGS 84
proj = "EPSG:4326",
geodesic = FALSE)
Map$setCenter(lon = -73.079, lat = -42.611, zoom = 7)
Map$addLayer(box)Ejemplo 2: estadisticas de un shp:
vias<-readOGR(".","region_los_lagos")## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\usuario\Desktop\ds\ejercicio", layer: "region_los_lagos"
## with 1 features
## It has 17 fields
## Integer64 fields read as strings: FID_1 TOTAL_VIVI PARTICULAR COLECTIVAS TOTAL_PERS HOMBRES MUJERESsummary(vias)## Object of class SpatialPolygonsDataFrame
## Coordinates:
## min max
## x -74.83774 -71.58093
## y -44.06706 -40.23878
## Is projected: FALSE
## proj4string : [+proj=longlat +datum=WGS84 +no_defs]
## Data attributes:
## FID_1 REGION NOM_REGION TOTAL_VIVI
## Length:1 Length:1 Length:1 Length:1
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
##
##
##
## PARTICULAR COLECTIVAS TOTAL_PERS HOMBRES
## Length:1 Length:1 Length:1 Length:1
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
##
##
##
## MUJERES DENSIDAD INDICE_MAS INDICE_DEP
## Length:1 Min. :17.11 Min. :97.64 Min. :47.03
## Class :character 1st Qu.:17.11 1st Qu.:97.64 1st Qu.:47.03
## Mode :character Median :17.11 Median :97.64 Median :47.03
## Mean :17.11 Mean :97.64 Mean :47.03
## 3rd Qu.:17.11 3rd Qu.:97.64 3rd Qu.:47.03
## Max. :17.11 Max. :97.64 Max. :47.03
## IND_DEP_JU IND_DEP_VE SHAPE_Leng Shape_Le_1
## Min. :30.55 Min. :16.48 Min. :7760164 Min. :60.65
## 1st Qu.:30.55 1st Qu.:16.48 1st Qu.:7760164 1st Qu.:60.65
## Median :30.55 Median :16.48 Median :7760164 Median :60.65
## Mean :30.55 Mean :16.48 Mean :7760164 Mean :60.65
## 3rd Qu.:30.55 3rd Qu.:16.48 3rd Qu.:7760164 3rd Qu.:60.65
## Max. :30.55 Max. :16.48 Max. :7760164 Max. :60.65
## Shape_Area
## Min. :5.259
## 1st Qu.:5.259
## Median :5.259
## Mean :5.259
## 3rd Qu.:5.259
## Max. :5.2591. Cruce entre el shp de Los Lagos y la mediana de una colección de imágenes filtrada Sentinel.
Tenemos problema para asir el shape de glaciares
Extraigamos su informacion
mask_1 <- st_read("Lim_glaciares.shp") ## Reading layer `Lim_glaciares' from data source `C:\Users\usuario\Desktop\ds\ejercicio\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# st_geometry_type(mask_1)st_crs(mask_1)## Coordinate Reference System:
## User input: WGS 84
## wkt:
## GEOGCRS["WGS 84",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["latitude",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["longitude",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4326]]st_bbox(mask_1)## xmin ymin xmax ymax
## -72.95502 -44.14786 -71.15580 -41.39342mask_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...summary(mask_1)## ID geometry
## Min. : 1 POLYGON Z :2085
## 1st Qu.: 522 epsg:4326 : 0
## Median :1043 +proj=long...: 0
## Mean :1043
## 3rd Qu.:1564
## Max. :2085r <- mask_1$geometry
r[4]## 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...mmm <- mask_1$geometry[6]
polig <- ggplot() + geom_sf(data = mmm, size = 0.1, color = "red", fill = "cyan1") + ggtitle("poligono de glaciar") + coord_sf()
polig
intentamos graficar otros shapes:
No podemos graficar geometria con los limites de glaciares:
mask_1 <- st_read("Lim_glaciares.shp") ## Reading layer `Lim_glaciares' from data source `C:\Users\usuario\Desktop\ds\ejercicio\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# mask_1summary(mask_1)## ID geometry
## Min. : 1 POLYGON Z :2085
## 1st Qu.: 522 epsg:4326 : 0
## Median :1043 +proj=long...: 0
## Mean :1043
## 3rd Qu.:1564
## Max. :2085mmm <- mask_1$geometry[6]# mmm$geometry()polig <- ggplot() + geom_sf(data = mmm, size = 0.1, color = "red", fill = "cyan1") + ggtitle("poligono de glaciar") + coord_sf()
polig
hay que llevar mask_1 a geometria gee
# ndgi# box <- ee$Geometry$Rectangle(coords = c(-73,-43,-72,-42.5),
box <- ee$Geometry$Rectangle(coords = c(-74.83774,-44.06706,-71.58093,-40.23878),
## WGS 84
proj = "EPSG:4326",
geodesic = FALSE)intentar funciones
ee\(Image\)clip
https://developers.google.com/earth-engine/apidocs/ee-image-clip
Recorta una imagen a una geometría o entidad.
Las bandas de salida corresponden exactamente a las bandas de entrada, excepto que los datos no cubiertos por la geometría están enmascarados. La imagen de salida conserva los metadatos de la imagen de entrada.
Utilice clipToCollection para recortar una imagen a FeatureCollection.
Devuelve la imagen recortada.
To clip the data to a region mask_1
# mmm <- mask_1$geometry[6]
# mmm# img_02 <- ee_as_raster(
# image = ndgi,
# region = mmm
# )Referencias:
https://drive.google.com/drive/folders/1yDB9oZBS6ZSZ-U1Rf7WsNlLdaHdHl-3z
https://rstudio-pubs-static.s3.amazonaws.com/643255_63554e7e1f9f466dad4f9e62ff977a88.html
https://drive.google.com/drive/folders/1yDB9oZBS6ZSZ-U1Rf7WsNlLdaHdHl-3z?usp=sharing
https://rstudio-pubs-static.s3.amazonaws.com/639598_f8b124e23b7949a49250693dc3c5a6a7.html
https://github.com/csaybar/rgee/blob/examples//GetStarted/03_finding_images.R
https://rpubs.com/daniballari/raster
https://rpubs.com/ohfrancom/618462
https://www.ide.cl/index.php/limites-y-fronteras/item/1528-division-politica-administrativa-2020
Catastro de Lagos https://www.ide.cl/index.php/aguas-continentales/item/1508-catastro-de-lagos
Lagos de Chile
http://datos.cedeus.cl/layers/geonode:cl_lagos_geo
# glaciar_ndgi <- ndgi$clip(mask_1)
# # glaciar_ndgi
# spplot(glaciar_ndgi)# Map$setCenter(lon = -73.079, lat = -42.611, zoom = 6)
# Map$addLayer(glaciar_ndgi , ndgiParams, "NDGI")mmm <- mask_1$geometry[6]
# mmmpolig <- ggplot() + geom_sf(data = mmm, size = 0.1, color = "red", fill = "cyan1") + ggtitle("poligono de glaciar") + coord_sf()
polig
# ndvi1