Glaciares en Chile
Areas glaciares en el tiempo y su medida.
Abstract
Se calcula la superficie glaciar para los primeros cuatrimestres del 2018 y 2019, obteniéndose un valor de 493.7959 el 2018 y de 384.6978 el 2019 en la región de los Lagos de Chile. Se consigue una metodología para capturar glaciar a glaciar, y se ejemplifica para Arica y Parinacota. Se obtiene un método para sumar pixeles rojos y verdes de los polígonos glaciares clasificados con RF. Resta automatizar los ciclos for.
del 1 al 566 no tocar!!!
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
## --------------------------------------------------------------------------------
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-01-01")
finish <- ee$Date("2019-04-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")
Creamos un cuadrado arbitrario
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 columna llamada 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
}
ndgi_01 <- 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 = -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")
Los rojos interiores son glaciares 0 = rojo
muestras_in <- sale_4$sampleRegions(
collection = ee$Feature(region_0),
scale = 1500,
tileScale = 16,
geometries = TRUE
)
muestras_in_gee_rojas = muestras_in$filter(ee$Filter$eq('nd', 0))
#muestras_in_gee_verdes = muestras_in$filter(ee$Filter$eq('nd', 1))
#Obtencion de puntos glaciares:
Map$setCenter(-72.7, -43.4, 6)
Map$addLayer(
eeObject = muestras_in_gee_rojas,
visParams = {},
name = "puntos glaciares"
)
Obtencion de puntos no 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")
sale_3 <- ndgi_01$clip(diff1)
Map$setCenter(lon = -73.079, lat = -42.611, zoom = 7)
Map$addLayer(sale_3, ndgiParams, "NDGI")
muestras_out <- sale_3$sampleRegions(
collection = ee$Feature(diff1),
scale = 1500,
tileScale = 16,
geometries = TRUE
)
muestras_in_gee_verdes = muestras_in$filter(ee$Filter$eq('nd', 1))
#muestras_in_gee_rojas = muestras_in$filter(ee$Filter$eq('nd', 0))
#Obtencion de puntos glaciares:
Map$setCenter(-72.7, -43.4, 6)
Map$addLayer(
eeObject = muestras_in_gee_verdes,
visParams = {},
name = "puntos glaciares"
)
union_total <- muestras_in_gee_verdes$merge(muestras_in_gee_rojas)
mask <- st_read("region_los_lagos.shp",
quiet = TRUE) %>%
sf_as_ee()
# convertimos el shp en geometria:
region <- mask$geometry()
start <- ee$Date("2019-01-01")
finish <- ee$Date("2019-04-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", "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.
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")
verdes_in <- classified$select('classification')$eq(0)
# 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= box_2,
scale= 100 ,
maxPixels= 1e9
)
verdes_in_km_area$getInfo()
## $classification
## [1] 346.3606
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("2018-01-01")
finish <- ee$Date("2018-04-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
}
ndgi_01 <- 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 = -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")
Los pixeles rojos son glaciares
muestras_in <- sale_4$sampleRegions(
collection = ee$Feature(region_0),
scale = 1500,
tileScale = 16,
geometries = TRUE
)
muestras_in_gee_rojas = muestras_in$filter(ee$Filter$eq('nd', 0))
#muestras_in_gee_verdes = muestras_in$filter(ee$Filter$eq('nd', 1))
#Obtencion de puntos glaciares:
Map$setCenter(-72.7, -43.4, 6)
Map$addLayer(
eeObject = muestras_in_gee_rojas,
visParams = {},
name = "puntos glaciares"
)
Obtencion de puntos no 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")
sale_3 <- ndgi_01$clip(diff1)
Map$setCenter(lon = -73.079, lat = -42.611, zoom = 7)
Map$addLayer(sale_3, ndgiParams, "NDGI")
muestras_out <- sale_3$sampleRegions(
collection = ee$Feature(diff1),
scale = 1500,
tileScale = 16,
geometries = TRUE
)
muestras_in_gee_verdes = muestras_in$filter(ee$Filter$eq('nd', 1))
#muestras_in_gee_rojas = muestras_in$filter(ee$Filter$eq('nd', 0))
#Obtencion de puntos glaciares:
Map$setCenter(-72.7, -43.4, 6)
Map$addLayer(
eeObject = muestras_in_gee_rojas,
visParams = {},
name = "puntos glaciares"
)
union_total <- muestras_in_gee_verdes$merge(muestras_in_gee_rojas)
mask <- st_read("region_los_lagos.shp",
quiet = TRUE) %>%
sf_as_ee()
# convertimos el shp en geometria:
region <- mask$geometry()
start <- ee$Date("2018-01-01")
finish <- ee$Date("2018-04-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", "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.
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")
verdes_in <- classified$select('classification')$eq(0)
# 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= box_2,
scale= 100 ,
maxPixels= 1e9
)
verdes_in_km_area$getInfo()
## $classification
## [1] 368.7491
Ahora lo que necesitamos es la obtencion clasificada de glaciar por glaciar:
Lo primero que podriamos hacer es obtener muestras dentro de los limites glciares que sean rojos en la region de arica y tarapaca
mask <- st_read("regiones_separadas/Region_15.shp",
quiet = TRUE) %>%
sf_as_ee()
# convertimos el shp en geometria:
region <- mask$geometry()
Map$setCenter(lon = -69.0055, lat = -18.5774, zoom = 7)
Map$addLayer(region)
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 = -69.0055, lat = -18.5774, zoom = 7)
Map$addLayer(first,vizParams, "Landsat 8 image")
getNDGI <- function(image)
{
#image$normalizedDifference(c("B3", "B4"))
image$normalizedDifference(c("B3", "B4")) < 0
}
ndgi_01 <- 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.0055, lat = -18.5774, zoom = 7)
Map$addLayer(ndgi_01, ndgiParams, "NDGI")
Intersectamos con los límites políticos de los glaciares
glaciar_01 <- st_read("GlaciaresxRegion/Simple_ARICA_Y_PARINACOTA.shp",
quiet = TRUE) %>%
sf_as_ee()
# convertimos el shp en geometria:
region_0 <- glaciar_01$geometry()
sale_4 <- ndgi_01$clip(region_0)
Map$setCenter(lon = -69.0055, lat = -18.5774, zoom = 7)
Map$addLayer(sale_4, ndgiParams, "NDGI")
Los rojos interiores son glaciares 0 = rojo
muestras_in <- sale_4$sampleRegions(
collection = ee$Feature(region_0),
scale = 100,
tileScale = 16,
geometries = TRUE
)
###Los pixels rojos son los 0
muestras_in_gee_rojas = muestras_in$filter(ee$Filter$eq('nd', 0))
#Obtencion de puntos glaciares:
Map$setCenter(lon = -69.0055, lat = -18.5774, zoom = 7)
Map$addLayer(
eeObject = muestras_in_gee_rojas,
visParams = {},
name = "puntos glaciares"
)
mask <- st_read("regiones_separadas/Region_15.shp",
quiet = TRUE) %>%
sf_as_ee()
# convertimos el shp en geometria:
region <- mask$geometry()
Map$setCenter(lon = -69.0055, lat = -18.5774, zoom = 7)
Map$addLayer(region)
diff1 <- region$difference(region_0, ee$ErrorMargin(1))
Map$addLayer(diff1, list(color = "FFFF00"), "diff1")
sale_3 <- ndgi_01$clip(diff1)
Map$setCenter(lon = -69.0055, lat = -18.5774, zoom = 7)
Map$addLayer(sale_3, ndgiParams, "NDGI")
muestras_in <- sale_3$sampleRegions(
collection = ee$Feature(diff1),
scale = 1000,
tileScale = 16,
geometries = TRUE
)
muestras_in_gee_verdes = muestras_in$filter(ee$Filter$eq('nd', 1))
#muestras_in_gee_verdes = muestras_in$filter(ee$Filter$eq('nd', 1))
#Obtencion de puntos glaciares:
Map$setCenter(lon = -69.0055, lat = -18.5774, zoom = 7)
Map$addLayer(eeObject = muestras_in_gee_verdes,visParams = {}, name = "puntos glaciares")
union_total <- muestras_in_gee_verdes$merge(muestras_in_gee_rojas)
mask <- st_read("regiones_separadas/Region_15.shp",
quiet = TRUE) %>%
sf_as_ee()
# convertimos el shp en geometria:
region <- mask$geometry()
start <- ee$Date("2018-01-01")
finish <- ee$Date("2018-04-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", "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.
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")
Este paso en muy, muy importante, y es el area poligonal oficial donde nosotros estudiamos la variacion glaciar. Ahora debemos cortar por los limites graciares y obtener los shps. No nos interesan las zonas no glaciares!!!
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
reg_15_2 <- mask_0 %>%
filter(COD_GLA =="CL101010001")
varia <- reg_15_2 %>%
sf_as_ee()
# convertimos el shp en geometria:
region_0 <- varia$geometry()
sale_4 <- classified$clip(region_0)
Map$setCenter(lon = -69.0055, lat = -18.5774, zoom = 7)
Map$addLayer(sale_4, viz_class, "NDGI")
fff <- sale_4$select('classification')$eq(0)
class_areas = fff$reduceRegion(
reducer= ee$Reducer$sum(),
geometry= region_0,
scale= 30 ,
maxPixels= 1e9
)
class_areas$getInfo()
## $classification
## [1] 282.3686
Anexo.
glaciar_001 <- read_sf("GlaciaresxRegion/Simple_ARICA_Y_PARINACOTA.shp")
glaciar_001
## 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
## Warning: `...` is not empty.
##
## We detected these problematic arguments:
## * `needs_dots`
##
## These dots only exist to allow future extensions and should be empty.
## Did you misspecify an argument?
## # A tibble: 327 x 56
## Id COD_GLA NOMBRE CLASIFICA REGION COMUNA DATUM HUSO ESTE NORTE
## <int> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <int> <int>
## 1 11 CL1010~ S/N GLACIAR ~ ARICA~ PUTRE WGS ~ 19SUR 490679 8.01e6
## 2 1 CL1010~ VOLCA~ GLACIAR ~ ARICA~ PUTRE WGS ~ 19SUR 484875 7.99e6
## 3 8 CL1010~ S/N GLACIAR ~ ARICA~ PUTRE WGS ~ 19SUR 494376 7.97e6
## 4 11 CL1010~ S/N GLACIAR ~ ARICA~ PUTRE WGS ~ 19SUR 494388 7.97e6
## 5 14 CL1010~ S/N GLACIAR ~ ARICA~ PUTRE WGS ~ 19SUR 494057 7.97e6
## 6 18 CL1010~ S/N GLACIARE~ ARICA~ PUTRE WGS ~ 19SUR 494908 7.97e6
## 7 2 CL1010~ S/N GLACIARE~ ARICA~ PUTRE WGS ~ 19SUR 493272 7.98e6
## 8 37 CL1010~ VOLCA~ GLACIAR ~ ARICA~ PUTRE WGS ~ 19SUR 484474 7.99e6
## 9 36 CL1010~ VOLCA~ GLACIAR ~ ARICA~ PUTRE WGS ~ 19SUR 484792 7.99e6
## 10 12 CL1010~ S/N GLACIARE~ ARICA~ PUTRE WGS ~ 19SUR 490628 8.01e6
## # ... with 317 more rows, and 46 more variables: FUENTE_DIG <chr>,
## # FUENTE_FEC <chr>, INVENT_FEC <chr>, COD_BNA <chr>, AREA_Km2 <dbl>,
## # COD_CUEN <chr>, NOMB_CUEN <chr>, COD_REGION <chr>, COD_SCUEN <chr>,
## # COD_SSCUEN <chr>, N_CUENDRE <chr>, PRECM <chr>, LARGO_PROM <dbl>,
## # LMAXTOTAL <dbl>, L_MAX_ACUM <dbl>, L_MAX_EXP <dbl>, LMAXABLAC <dbl>,
## # ANCHO_PROM <dbl>, ESP_MED <dbl>, AREA_EXP <dbl>, AREA_CUB <dbl>,
## # AREA_ABLAC <dbl>, AREA_ACUM <dbl>, VOL_M3 <dbl>, HMAX <dbl>, HMEDIA <dbl>,
## # HMINTOTAL <dbl>, HMEDABLAC <dbl>, HMINEXP <dbl>, LATITUD <dbl>,
## # LONGITUD <dbl>, ORIENACUM <chr>, ORIENABLAC <chr>, ORIENTA <chr>,
## # CLAS_WGI <chr>, CLAS_2_CUB <chr>, ZONA_GLACI <chr>, HMEDIANA <dbl>,
## # PENDIENTE <dbl>, ERROR_KM2 <dbl>, ERROR_PORC <dbl>, InPoly_FID <dbl>,
## # SimPgnFlag <int>, MaxSimpTol <dbl>, MinSimpTol <dbl>,
## # geometry <MULTIPOLYGON [°]>