El hielo glaciar cubre actualmente casi 16 millones de \(km^2\) de la superficie de la Tierra, la mayor parte del cual en la Antártida (aproximadamente 13,5 millones de \(km^2\)) y Groenlandia (1,74 millones de \(km^2\)). El 3% restante, unos 500.000 \(km^2\), aparece como casquetes polares, campos de hielo y glaciares de montaña distribuidos por todos los continentes. A escala global, los factores que controlan la aparición de glaciares son la latitud, la altitud y la proximidad a las fuentes de humedad. A continuación, se realizará una actividad guiada para que los alumnos comprendan el papel de cada uno de estos factores.
Para ello, se importará el regional analisis_global.csv
,
que incluye un resumen, por regiones, de la base de datos Randolph Glacier Inventory – A Dataset
of Global Glacier Outlines: Version 6.0. El regional puede
descargarse del siguiente link
global <- read.csv2("https://personales.unican.es/rasillad/docencia/g171/Glaciar/analisis_global.csv", na.strings = -9999)
Su estructura interna es la siguiente:
str(global)
## 'data.frame': 87 obs. of 13 variables:
## $ region : chr "Alaska" "Alaska" "Alaska" "Alaska" ...
## $ montana : chr "Alaska Pena (Aleutians)" "Alaska Ra (Wrangell/Kilbuck)" "N Alaska" "N Coast Ranges" ...
## $ lon : num -158 -150 -147 -133 -141 ...
## $ lat : num 56.7 61.8 68.7 57.8 60.6 ...
## $ area_total : num 1912 16278 346 22963 33174 ...
## $ area_promedio : num 2.193 2.816 0.561 2.176 6.594 ...
## $ zmin : num 761 1503 1717 1297 1532 ...
## $ zmax : num 1345 2016 2129 1769 2024 ...
## $ zmediana : num 1013 1748 1905 1546 1776 ...
## $ lmax : num 1875 1651 1016 1404 1796 ...
## $ zona : int 2 2 2 2 2 2 5 5 5 5 ...
## $ continentalidad_eurasia: int 0 0 0 0 0 0 0 0 0 0 ...
## $ latitudinalidad_america: int 1 1 1 1 1 1 1 0 0 1 ...
o region: nombre de la región.
o montaña: nombre de la cordillera.
o lon: centroide de la longitud de las cordilleras (grados centesimales).
o lat: centroide de la latitud del las cordilleras (grados centesimales).
o area_total: área total de los glaciares (en \(km^2\))
o area_promedio: superficie media, producto de dividir el área total entre el número de glaciares (en \(km^2\))
o zmin: promedio de la altitud mínima de los glaciares de la cordillera (m).
o zmax: altitud máxima de la cordillera (m)
o zmed: altitud mediana de cada glaciar (m). Este parámetro permite el cálculo de la altitud de la línea de equilibrio (ELA), que marca el área o zona de un glaciar donde la acumulación está en equilibrio con la ablación. Cuando el balance de masa neto anual es negativo, la ELA asciende, y cuando el balance de masa neto anual es positivo, la ELA desciende de alitud. Las fluctuaciones en el ELA constituyen un buen indicador de la respuesta de los glaciares a los cambios climáticos (precipitación en la temporada de acumulación, temperatura en la temporada de ablación y direcciones predominantes del viento con nieve). Hay que tener en cuenta que estudios realizados previamente señalan que este método funciona bien para glaciares pequeños con distribuciones uniformes de área/altitud, pero tiende a sobreestimar la ELA, ya que no incluye las variaciones en la morfología del valle, que afectan al cálculo del área y de la elevación de un glaciar.
o zona: zona latitudinal a la que pertenece cada región (m)
o continentalidad_eurasia: número para analizar el papel de la continentalidad.
o latitudinalidad_america: número para analizar el papel de la latitud.
En caso de querer conocer las etiquetas correspondientes a alguna
variable importada como character
se puede aplicar función
unique()
. En el caso de la variable region
unique(global$region)
## [1] "Alaska" "Antarctic and Subantarctica"
## [3] "Arctic Canada (North)" "Arctic Canada (South)"
## [5] "Caucasus and Middle East" "Central Asia"
## [7] "Central Europe" "Greenland"
## [9] "Iceland" "Low Latitudes"
## [11] "New Zealand" "North Asia"
## [13] "Russian Arctic" "Scandinavia"
## [15] "South Asia (East)" "South Asia (West)"
## [17] "Southern Andes" "Svalbard and Jan Mayen"
## [19] "Western Canada and USA"
En el caso de la variable montana
unique(global$montana)
## [1] "Alaska Pena (Aleutians)" "Alaska Ra (Wrangell/Kilbuck)"
## [3] "N Alaska" "N Coast Ranges"
## [5] "St Elias Mtns" "W Chugach Mtns (Talkeetna)"
## [7] "Alexander Island 7H2" "Amery Ice Shelf 7B"
## [9] "Balleny Islands" "Bellingshausen Sea 7H1"
## [11] "E Queen Maud Land 7A" "Marie Byrd Land 7F"
## [13] "NE Antarctic Pena 7I2" "Pine Island Bay 7G"
## [15] "Ross Ice Shelf 7E" "SE Antarctic Pena 7I3"
## [17] "South Shetlands and South Orkneys" "Subantarctic (Atlantic)"
## [19] "Subantarctic (Indian)" "Subantarctic (Pacific)"
## [21] "Victoria Land 7D" "W Antarctic Pena 7I1"
## [23] "W Queen Maud Land 7K" "Wilkes Land 7C"
## [25] "Axel Heiberg and Meighen Is" "Devon Island"
## [27] "Melville Island" "N Ellesmere Island"
## [29] "S Ellesmere Island (NW Devon)" "SC Ellesmere Island"
## [31] "Bylot Island" "Cumberland Sound"
## [33] "EC Baffin Island" "Frobisher Bay"
## [35] "Labrador" "N Baffin Island"
## [37] "NE Baffin Island" "SE Baffin Island"
## [39] "W Baffin Island" "Greater Caucasus"
## [41] "Middle East" "E Kun Lun (Altyn Tagh)"
## [43] "E Tien Shan (Dzhungaria)" "Hissar Alay"
## [45] "Inner Tibet" "Pamir (Safed Khirs/W Tarim)"
## [47] "Qilian Shan" "S and E Tibet"
## [49] "W Kun Lun" "W Tien Shan"
## [51] "Alps" "Southern and Eastern Europe"
## [53] "Greenland periphery" "Iceland"
## [55] "E Africa" "Equatorial Andes"
## [57] "Mexico" "New Guinea"
## [59] "Peru-Bolivian Andes" "New Zealand"
## [61] "Altay and Sayan" "Central Siberia"
## [63] "Cherskiy and Suntar Khayata Ranges" "E Chukotka"
## [65] "NE Russia" "Ural Mountains"
## [67] "Franz Josef Land" "Novaya Zemlya"
## [69] "Severnaya Zemlya" "N Scandinavia"
## [71] "SE Scandinavia" "SW Scandinavia"
## [73] "C Himalaya" "E Himalaya"
## [75] "Hengduan Shan" "Hindu Kush"
## [77] "Karakoram" "W Himalaya"
## [79] "C Andes" "Patagonia"
## [81] "Jan Mayen" "Svalbard"
## [83] "Cascade Ra and Sa Nevada" "Mackenzie and Selwyn Mtns"
## [85] "N Rocky Mtns" "S Coast Ranges"
## [87] "S Rocky Mtns"
Como es habitual, algunas variables deben transformarse en factores.
global$zona <- factor(global$zona,
levels = c(1,2,3,4,5),
labels = c("Ártico", "HN templado", "Tropical", "HS templado", "Antártica"))
global$continentalidad_eurasia <- as.factor(global$continentalidad_eurasia)
global$latitudinalidad_america <- as.factor(global$latitudinalidad_america)
En este apartado se analizará la distribución espacial de glaciares en todo el planeta. En primer lugar, se representarán sobre un mapa de todo el planeta, utilizando las funciones del paquete ggplot.
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.3.2
library(RColorBrewer)
A continuación se crea un mapa de las regiones glaciares en función de su tipología.
ggplot() + borders("world", fill = "white") + theme_bw() +
geom_point(data = global,
aes(x = lon, y = lat, color = zona, size = area_total)) +
scale_color_brewer(palette="Dark2")
En el mapa puede observarse claramente que los glaciares se concentran en las latitudes altas o en regiones montañosas de todo el mundo, preferentemente en el hemisferio norte, como los Alpes y el Himalaya, así como alrededor del Océano Ártico. En el hemisferio sur se aprecian sobre los Andes, Nueva Zelanda e islas subantárticas como Georgia del Sur, así como en algunos márgenes de la Antártida. El número de glaciares en latitudes intertropicales es escaso (Andes, Kilimanjaro, Nueva Guinea)
La distribución de los glaciares refleja la interacción entre la precipitación, la temperatura y la topografía. Estos factores varían sistemáticamente en todo el mundo, especialmente con la latitud, la altitud y la distancia de una fuente de humedad. A continuación se evalúa la importancia de estos factores. En primer lugar, se creará una nueva variable con el valor absoluto de la latitud, puesto que esta variable adopta valores positivos (Hemisferio Norte) y negativos (Hemisferio Sur)
global$lat_abs <- abs(global$lat)
Para obtener una matriz de correlaciones sólo entre las variables cuantitativas, se crea un nuevo objeto que sólo incluye éstas últimas:
regional_corr <- global[ , c(3,14,5:10)]
str(regional_corr)
## 'data.frame': 87 obs. of 8 variables:
## $ lon : num -158 -150 -147 -133 -141 ...
## $ lat_abs : num 56.7 61.8 68.7 57.8 60.6 ...
## $ area_total : num 1912 16278 346 22963 33174 ...
## $ area_promedio: num 2.193 2.816 0.561 2.176 6.594 ...
## $ zmin : num 761 1503 1717 1297 1532 ...
## $ zmax : num 1345 2016 2129 1769 2024 ...
## $ zmediana : num 1013 1748 1905 1546 1776 ...
## $ lmax : num 1875 1651 1016 1404 1796 ...
A continuación, se cargará el paquete Hmisc
.
# install.packages("Hmisc")
library("Hmisc")
## Warning: package 'Hmisc' was built under R version 4.3.3
##
## Attaching package: 'Hmisc'
## The following objects are masked from 'package:base':
##
## format.pval, units
Cálculo del coeficiente de correlación de Pearson. Devuelve los coeficientes de correlación y p-value de los pares de columnas.
matriz.corr <- rcorr(as.matrix(regional_corr))
matriz.corr
## lon lat_abs area_total area_promedio zmin zmax zmediana lmax
## lon 1.00 -0.28 -0.16 -0.11 0.35 0.33 0.37 -0.09
## lat_abs -0.28 1.00 0.24 0.27 -0.89 -0.88 -0.89 0.40
## area_total -0.16 0.24 1.00 0.17 -0.16 -0.13 -0.16 0.24
## area_promedio -0.11 0.27 0.17 1.00 -0.30 -0.30 -0.28 0.88
## zmin 0.35 -0.89 -0.16 -0.30 1.00 1.00 0.99 -0.39
## zmax 0.33 -0.88 -0.13 -0.30 1.00 1.00 0.99 -0.37
## zmediana 0.37 -0.89 -0.16 -0.28 0.99 0.99 1.00 -0.38
## lmax -0.09 0.40 0.24 0.88 -0.39 -0.37 -0.38 1.00
##
## n= 87
##
##
## P
## lon lat_abs area_total area_promedio zmin zmax zmediana
## lon 0.0092 0.1421 0.3074 0.0009 0.0016 0.0004
## lat_abs 0.0092 0.0270 0.0106 0.0000 0.0000 0.0000
## area_total 0.1421 0.0270 0.1178 0.1285 0.2213 0.1423
## area_promedio 0.3074 0.0106 0.1178 0.0047 0.0053 0.0079
## zmin 0.0009 0.0000 0.1285 0.0047 0.0000 0.0000
## zmax 0.0016 0.0000 0.2213 0.0053 0.0000 0.0000
## zmediana 0.0004 0.0000 0.1423 0.0079 0.0000 0.0000
## lmax 0.4302 0.0001 0.0227 0.0000 0.0002 0.0004 0.0003
## lmax
## lon 0.4302
## lat_abs 0.0001
## area_total 0.0227
## area_promedio 0.0000
## zmin 0.0002
## zmax 0.0004
## zmediana 0.0003
## lmax
ACTIVIDAD: ¿Cuáles son las variables mejor relacionadas con las dimensiones y altitud de los glaciares?
A continuación se elabora un gráfico con las correlaciones entre las citadas variables, pero incluyendo la localización latitudinal como variable adicional.
grupos <- global[, 11] # Variable convertida en factor y denominada grupos
pairs(regional_corr,
labels = colnames(regional_corr), # Nombres de las variables
pch = 21, # Símbolo pch
bg = rainbow(5)[grupos], # Color de fondo del símbolo (pch 21 a 25)
col = rainbow(5)[grupos], # Color de borde del símbolo
main = "Glaciares", # Título del gráfico
row1attop = TRUE, # Si FALSE, cambia la dirección de la diagonal
gap = 1, # Distancia entre subplots
cex.labels = NULL, # Tamaño del texto de la diagonal
font.labels = 1)
Un gráfico similar se obtiene con el paquete
PerformanceAnalytics
# install.packages("PerformanceAnalytics")
library(PerformanceAnalytics)
## Loading required package: xts
## Warning: package 'xts' was built under R version 4.3.2
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
##
## Attaching package: 'PerformanceAnalytics'
## The following object is masked from 'package:graphics':
##
## legend
chart.Correlation(regional_corr, # ATENCION: dataframe original, no la matriz de corrrelación
method="spearman",
histogram=TRUE,
pch=16)
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in cor.test.default(as.numeric(x), as.numeric(y), method = method):
## Cannot compute exact p-value with ties
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in cor.test.default(as.numeric(x), as.numeric(y), method = method):
## Cannot compute exact p-value with ties
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in cor.test.default(as.numeric(x), as.numeric(y), method = method):
## Cannot compute exact p-value with ties
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in cor.test.default(as.numeric(x), as.numeric(y), method = method):
## Cannot compute exact p-value with ties
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in cor.test.default(as.numeric(x), as.numeric(y), method = method):
## Cannot compute exact p-value with ties
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in cor.test.default(as.numeric(x), as.numeric(y), method = method):
## Cannot compute exact p-value with ties
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in cor.test.default(as.numeric(x), as.numeric(y), method = method):
## Cannot compute exact p-value with ties
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
Dado que la pareja de variables mejor correlacionada es la latitud (absoluta) y la altura mínima que alcanza el glaciar, se analizará con más detalle su relación.
plot(global$lat_abs ~ global$zmin, col = global$zona) # dependiente ~ independiente
legend(x = "bottomleft",
legend = c("Ártico", "HN templado", "Trópicos", "HS templado", "Antártida"),
title = "Zona",
fill = c("black", "red", "green", "blue", "lightblue"))
A igualdad de otros factores, el tamaño de los glaciares aumenta hacia los polos porque el bajo ángulo solar significa que hay menos energía disponible para fundir el hielo. Al mismo tiempo, en una misma latitud, las montañas más elevadas tendrán una mayor superficie glaciada porque la temperatura disminuye con la altitud.
La interacción de factores latitudinales y altitudinales crea un amplio patrón global de glaciación en el que los glaciares solo existen en regiones ecuatoriales a grandes altitudes y ocupan progresivamente altitudes cada vez más bajas hacia los polos. Por ejemplo, los glaciares solo existen por encima de 4500-5000 m en las montañas tropicales, pero alcanzan el nivel del mar en latitudes altas como el Alto Ártico canadiense y Groenlandia.
papel_latitud <- subset(global,
latitudinalidad_america == 1)
papel_latitud <- papel_latitud[order(papel_latitud$lat), ]
dotchart(papel_latitud$zmin,
labels = papel_latitud$montana, pch = 21, bg = "green", pt.cex = 1.0)
Este patrón emerge claramente en la figura anterior, que muestra la
variable zmin
en un transecto norte-sur a través de las
cordilleras de América. Esa variable se eleva desde ambos polos hacia el
ecuador, alcanzando un máximo en las tierras altas subtropicales
mexicanas y los Andes colombianos y ecuatorianos, aunque la relación
lineal se complica por las variaciones en la precipitación, por lo que,
en algunos casos, los glaciares alcanzan altitudes más bajas en áreas
con altas precipitaciones. Este es el caso de la vertiente sur
fuertemente glaciar de la cordillera de Alaska y los Campos de Hielo
Patagónicos, donde las nevadas pueden superar los 7 m por año. Del mismo
modo, un estudio más detallado mostraría que las zmin sobre
relativamente más bajas en Colombia y Ecuador frente a Perú y Bolivia,
al estar afectadas aquéllas por la Zona de Convergencia
Intertropical.
A pesar de que una zona puede ser suficientemente fría para la supervivencia de un glaciar, también puede estar a una distancia considerable de la fuente de humedad, por lo que las cantidades de precipitación recibidas serán muy bajas. Este es el caso del interior de los continentes, que a pesar de experimentar temperaturas extremadamente bajas durante el invierno (e inviernos prolongados), las precipitaciones recibidas pueden ser insuficientes para mantener glaciares. Este efecto se observa en el interior de Eurasia.
papel_continentalidad <- subset(global,
continentalidad_eurasia == 1)
papel_continentalidad <- papel_continentalidad[order(papel_continentalidad$lon), ]
dotchart(papel_continentalidad$zmin,
labels = papel_continentalidad$montana, pch = 21, bg = "green", pt.cex = 1.0)
En este caso, se observa un progresivo incremento de la altura mínima de los glaciares a medida que nos internamos en las latitudes medias de Eurasia. Esto supone una menor precipitación, ya que ésta procede fundamentalmente de los vientos del W, que cargan humedad en el océano Atlántico, y en el caso de las cordilleras más orientales, del monzón. No obstante, existen algunas anomalías que merece la pena analizar.
A partir de la base de datos anterior, selecciona la
region
“Central Asia” y analiza la relación entre las
diferentes variables que caracterizan cada región. Aplica las
consideraciones anteriores acerca del papel de la altitud y del clima
usando como material de apoyo los climodiagramas que podrás utilizar en
la siguiente página web http://www.klimadiagramme.de/
Además de la latitud, la aparición de glaciares y sus dimensiones está condicionada por la altitud de las cordilleras. El Himalaya, los Alpes europeos, las Rocosas y muchas otras cordilleras poseen suficiente altitud como para albergar un considerable número de glaciares.
A escala regional, tabmién juega un papel importante en la ubicación concreta de los aparatos glaciares, especialmente en áreas marginales, ya que en esas zonas la ELA se encuentra muy cerca de la línea de cubres.
A escalas locales, otras variables, como la orientación, la pendiente o la exposición a los vientos húmedos pueden modificar la ELA, debido a la recepción diferencial de radiación solar y precipitación. Normalmente, las ELAs más bajas en muchas cordilleras del hemisferio norte se encuentran en cuencas orientadas al noreste, debido a que reciben neos radiación solar (sobre todo en verano) y están a sotavento de los vientos predominantes del suroeste, que barren la nieve hacia esa orientación. Por otro lado, las laderas con mayor pendiente no suelen facilitar la acumulación de nieve, ya que las avalanchas trasladan esa nieve a altitudes más bajas.
Para analizar estos aspectos, se importará el regional
analisis_regional.csv
, que es una versión simplificada, de
la base de datos Randolph Glacier
Inventory – A Dataset of Global Glacier Outlines: Version 6.0. El
regional puede descargarse del siguiente link
Antes de realizar cualquier tipo de análisis, es necesario realizar una serie de procedimientos para obtener y organizar los datos.
regional <- read.csv2("https://personales.unican.es/rasillad/docencia/g171/Glaciar/analisis_regional.csv", na.strings = -9999)
Su estructura interna es la siguiente:
str(regional)
## 'data.frame': 87250 obs. of 13 variables:
## $ region : chr "Western Canada and USA" "Western Canada and USA" "Western Canada and USA" "Western Canada and USA" ...
## $ cdg_region : int 2 2 2 2 2 2 2 2 2 2 ...
## $ montana : chr "S Coast Ranges" "S Coast Ranges" "S Coast Ranges" "S Coast Ranges" ...
## $ cdg_montana: int 2 2 2 2 2 2 2 2 2 2 ...
## $ lon : num -128 -129 -128 -129 -128 ...
## $ lat : num 56.1 56 56 56 56 ...
## $ area : num 1.138 0.1 0.249 0.412 1.055 ...
## $ zmin : int 1396 1780 1689 1630 1521 1500 1561 1597 1702 1495 ...
## $ zmax : int 1978 1883 1833 1898 1988 1905 1980 1934 1917 1649 ...
## $ zmediana : int 1875 1825 1775 1775 1925 1725 1825 1825 1775 1625 ...
## $ pendiente : num 23.6 18.6 15.2 22.6 24.5 18.5 20.2 20.3 26.6 26.8 ...
## $ aspecto : int 26 97 6 353 357 7 56 9 108 3 ...
## $ lmax : int 689 313 511 588 1035 1294 850 878 372 362 ...
o reg: nombre de la región.
o cdg_reg: código de la región.
o montana: nombre de la cordillera.
o cdg_montana: código de la cordillera.
o lon: longitud del glaciar (grados centesimales).
o lat: latitud del glaciar (grados centesimales).
o area: superficie del glaciar (en \(km^2\))
o zmin: altitud mínima del glaciar (m).
o zmax: altitud máxima del glaciar (m)
o zmediana: altitud mediana del glaciar (m)
o pendiente: pendiente del glaciar (grados centesimales)
o aspecto: orientación del glaciar (grados)
o lmax: longitud del glaciar (m)
Comprobación de las etiquetas de las variables
names(regional)
## [1] "region" "cdg_region" "montana" "cdg_montana" "lon"
## [6] "lat" "area" "zmin" "zmax" "zmediana"
## [11] "pendiente" "aspecto" "lmax"
unique(regional$montana)
## [1] "S Coast Ranges" "N Rocky Mtns"
## [3] "Cascade Ra and Sa Nevada" "S Rocky Mtns"
## [5] "N Scandinavia" "SW Scandinavia"
## [7] "Altay and Sayan" "Alps"
## [9] "Greater Caucasus" "Pamir (Safed Khirs/W Tarim)"
## [11] "W Tien Shan" "E Tien Shan (Dzhungaria)"
## [13] "W Kun Lun" "E Kun Lun (Altyn Tagh)"
## [15] "Qilian Shan" "Inner Tibet"
## [17] "S and E Tibet" "E Himalaya"
## [19] "Peru-Bolivian Andes"
datos <- subset(regional,
montana =="Alps" | montana == "Greater Caucasus")
Cambiamos el nombre del inglés al español
datos$montana[datos$montana == "Alps"] <- "Alpes"
datos$montana[datos$montana == "Greater Caucasus"] <- "Cáucaso"
datos$orientacion[datos$aspecto <= 45] <- "N"
datos$orientacion[datos$aspecto > 45 & datos$aspecto <= 135] <- "E"
datos$orientacion[datos$aspecto > 135 & datos$aspecto <= 225] <- "S"
datos$orientacion[datos$aspecto > 225 & datos$aspecto <= 315] <- "W"
datos$orientacion[datos$aspecto > 315] <- "N"
datos <- subset(datos, select =-aspecto)
datos$orientacion <- ordered(datos$orientacion,
levels = c("N", "E", "S", "W"))
Reclasificación de la variable área
en categorías
datos$area_categoria[datos$area <= 0.03] <- 1
datos$area_categoria[datos$area > 0.03 & datos$area <= 0.10] <- 2
datos$area_categoria[datos$area > 0.10 & datos$area <= 0.25] <- 3
datos$area_categoria[datos$area > 0.25 & datos$area <= 0.79] <- 4
datos$area_categoria[datos$area > 0.79 & datos$area <= 6.06] <- 5
datos$area_categoria[datos$area > 6.06] <- 6
datos$area_categoria <- ordered(datos$area_categoria,
levels = c(1, 2, 3, 4, 5, 6),
labels = c("Diminuto", "Pequeño", "Mediano", "Grande", "Muy_grande", "Enorme"))
Reclasificación de la variable lmax
en categorías
datos$lmax_categoria[datos$lmax <= 157] <- 1
datos$lmax_categoria[datos$lmax > 157 & datos$lmax <= 392] <- 2
datos$lmax_categoria[datos$lmax > 392 & datos$lmax <= 705] <- 3
datos$lmax_categoria[datos$lmax > 705 & datos$lmax <= 1408] <- 4
datos$lmax_categoria[datos$lmax > 1408] <- 5
datos$lmax_categoria <- ordered(datos$lmax_categoria,
levels = c(1, 2, 3, 4, 5),
labels = c("Muy_corto", "Corto", "Mediano", "Largo", "Muy_largo"))
A partir de este momento, trabajaremos únicamente con el dataframe
datos
. La función attach
lo permite.
attach(datos)
borrame_alpes <- subset(datos,
montana =="Alpes")
borrame_alpes <- borrame_alpes[ , c(5:7)]
El primer paso para el estudio de los movimientos de ladera registrados en un país es ubicarlos geográficamente. Para representarlos gráficamente sobre un mapa, debe convertirse el dataframe en un objeto sf. Para ello, es necesario especificar:
Para realizar esa conversión, se usa la función
sf::st_as_sf()
.
library(sf)
## Warning: package 'sf' was built under R version 4.3.2
## Linking to GEOS 3.11.2, GDAL 3.7.2, PROJ 9.3.0; sf_use_s2() is TRUE
alpes_sf <- st_as_sf(borrame_alpes,
coords = c("lon", "lat"),
crs = "+proj=longlat +datum=WGS84")
Conviene que la ventana espacial asociada a este subconjunto de datos sea un poco mayor que el propio ámbito administrativo (recuérdese que el recorte de un objeto ráster recortado por un objeto vectorial deja algunos píxeles fuera).
st_bbox(alpes_sf)
## xmin ymin xmax ymax
## 5.93000 44.11700 14.55300 47.57464
bbox_new <- st_bbox(alpes_sf)
xrange <- bbox_new$xmax - bbox_new$xmin
yrange <- bbox_new$ymax - bbox_new$ymin
bbox_new[1] <- bbox_new[1] - 1
bbox_new[3] <- bbox_new[3] + 1
bbox_new[2] <- bbox_new[2] - 1
bbox_new[4] <- bbox_new[4] + 1
bbox_new <- bbox_new %>%
st_as_sfc() %>%
st_sf
Los siguientes objetos no se van a volver a utilizar, eliminándose
para no acumularse en la ventana del
Global Environment
.
rm(borrame_alpes, xrange, yrange)
Para analizar las condiciones topográficas favorables a la génesis de
los movimientos de ladera, se utilizará la información contenida en un
modelo digital del terreno. Para ello se utilizará el paquete
elevatr
, a través del que se puede
descargar directamente un mdt con diferentes niveles de resolución.
Este paquete puede descargarse directamente desde desde el repositorio CRAN.
library(raster)
## Warning: package 'raster' was built under R version 4.3.2
## Loading required package: sp
## Warning: package 'sp' was built under R version 4.3.2
#install.packages("elevatr")
library(elevatr)
## Warning: package 'elevatr' was built under R version 4.3.3
## elevatr v0.99.0 NOTE: Version 0.99.0 of 'elevatr' uses 'sf' and 'terra'. Use
## of the 'sp', 'raster', and underlying 'rgdal' packages by 'elevatr' is being
## deprecated; however, get_elev_raster continues to return a RasterLayer. This
## will be dropped in future versions, so please plan accordingly.
La función elevatr::get_elev_raster()
descarga un mdt en
formato ráster desde el sitio web “AWS Open Data Terrain Tiles”. Para
hacerlo necesitamos varios argumentos:
locations es el marco espacial del ámbito analizado.
z el nivel de resolución del mdt veáse
clip objeto espacial para recortar
mdt_alpes <- get_elev_raster(locations = bbox_new,
z = 7,
clip = "locations")
## Mosaicing & Projecting
## Clipping DEM to locations
## Note: Elevation units are in meters.
Actualización de los atributos del Raster, especialmente los valores máximo y mínimo, que se almacenan en el objeto ráster creado.
mdt_alpes <- setMinMax(mdt_alpes)
mdt_alpes # Comprobación en los atributos de la inclusión del nuevo parámetro
## class : RasterLayer
## dimensions : 1070, 2082, 2227740 (nrow, ncol, ncell)
## resolution : 0.005102136, 0.005102136 (x, y)
## extent : 4.929887, 15.55253, 43.11627, 48.57555 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : memory
## names : file124ce40797f
## values : -2633, 4622 (min, max)
Si existen valores inferiores a 0 (topografía submarina) deben eliminarse.
mdt_alpes <- calc(mdt_alpes, function(x){x[x < 0] <- NA;return(x)})
Para usar ggplot, es necesario convertir el objeto ráster en un
dataframe, mediante la función as.data.frame()
,
especificando las dimensiones x e y.
mdt_alpes_df <- as.data.frame(mdt_alpes,
xy = TRUE)
colnames(mdt_alpes_df)[3] <- "altitud"
Eliminamos las filas con uno o más valores ausentes (“NA”), usando la
función complete.cases()
.
mdt_alpes_df <- mdt_alpes_df[complete.cases(mdt_alpes_df), ]
mdt_alpes_df <- subset(mdt_alpes_df,
altitud >= 0)
A continuación ya se pueden representar gráficamente un mapa topográfico de fondo.
library(ggplot2)
#install.packages("marmap")
library(marmap)
## Registered S3 methods overwritten by 'adehabitatMA':
## method from
## print.SpatialPixelsDataFrame sp
## print.SpatialPixels sp
##
## Attaching package: 'marmap'
## The following object is masked from 'package:raster':
##
## as.raster
## The following object is masked from 'package:grDevices':
##
## as.raster
ggplot() +
geom_raster(data = mdt_alpes_df, aes(x = x, y = y, fill = altitud), show.legend = FALSE) +
geom_sf(data = alpes_sf, color = "red", fill = NA) +
coord_sf() +
scale_fill_etopo() +
labs(title = "Mapa topográfico de los Alpes", x = "Longitud", y = "Latitud", fill = "Altitud (metros)")
Esta función también permite aplicar una función a una matriz, array o data frame calculando estadísticos, pero no sólo por filas, sino también por columas.
promedios <- aggregate(cbind(zmin, zmax, zmediana, pendiente, lmax) ~ montana,
FUN = mean,
data = datos,
na.rm=TRUE)
datos["numero"] <- 1
promedios2 <- aggregate(cbind(area, numero) ~ montana,
FUN = sum,
data = datos,
na.rm=TRUE)
promedios <- cbind(promedios, promedios2)
View(promedios)
promedios <- promedios[ , c(1:6, 8:9)] # Se elimina la columan -variable- duplicada.
Por último, se grabará como datos excel
library(xlsx)
write.xlsx(promedios, # Data frame a ser exportado
"D:/G171_Procesos_Geomorfologicos_2022/5_Procesos_Glaciares/Cambio_global_2_Database_glaciares/resultados.xlsx", # Ruta completa
sheetName = "promedios", # Nombre de la hoja de Excel
col.names = TRUE, # Incluir los nombres de las columnas (TRUE) o no (FALSE)
row.names = TRUE, # Incluir los nombres de las filas (TRUE) o no (FALSE)
append = TRUE, # En este caso, sí agregar al archivo existente "resultados.xlsx"
showNA = TRUE) # Si TRUE, los NA serán celdas vacías
Para ello, podemos utilizar un gráfico de caja y bigotes incorporando una variable cuantitativa (nombre de la cordillera). A continuación se representa la superficie de los glaciares
boxplot(area ~ montana,
datos,
ylab = "km^2",
outline = FALSE, # sin outliers
notch = TRUE)
boxplot(zmin ~ montana,
datos,
main = "Altura mínima del frente",
xlab = "Cordillera",
ylab = "m",
outline = FALSE,
notch = TRUE)
boxplot(zmediana ~ montana,
datos,
main = "Altura mediana",
xlab = "Cordillera",
ylab = "m",
outline = FALSE,
notch = TRUE)
boxplot(pendiente ~ montana,
datos,
main = "Pendiente media del glaciar",
xlab = "Cordillera",
ylab = "º",
outline = FALSE,
notch = TRUE)
boxplot(lmax ~ montana,
datos,
main = "Longitud máxima del glaciar",
xlab = "Cordillera",
ylab = "m",
outline = FALSE,
notch = TRUE)
Con los datos anteriores se puede construir una tabla de contingencia, utilizando la función table(), ya vista en el capítulo anterior. En este ejemplo, se va a analizar la posible relación entre las variables garage y número de baños La tabla así configurada es la siguiente.
tabla_orientacion_ni <- table(datos$montana, datos$orientacion)
tabla_orientacion_ni
##
## N E S W
## Alpes 1533 980 605 774
## Cáucaso 582 440 253 362
La tabla de frecuencias relativas se calcula anidando la función prop.table() a la tabla calculada previamente, ¡no al vector original!. Su estructura es análoga a la tabla anterior. Las frecuencias relativas según filas (según cordillera) son:
tabla_orientacion_fi <- prop.table(tabla_orientacion_ni, margin= 1) * 100
tabla_orientacion_fi
##
## N E S W
## Alpes 39.38849 25.17986 15.54471 19.88695
## Cáucaso 35.55284 26.87844 15.45510 22.11362
Si convertimos esas tablas en vectores, podemos unirlos en una única tabla con la función cbind(), que puede visualizarse con View().
tabla_orientacion <- rbind(tabla_orientacion_ni, tabla_orientacion_fi)
View(tabla_orientacion)
Por último, se grabará en formato excel:
write.xlsx(tabla_orientacion, # Data frame a ser exportado
"D:/G171_Procesos_Geomorfologicos_2022/5_Procesos_Glaciares/Cambio_global_2_Database_glaciares/resultados.xlsx", # Ruta completa
sheetName = "tabla_orientacion", # Nombre de la hoja de Excel
col.names = TRUE, # Incluir los nombres de las columnas (TRUE) o no (FALSE)
row.names = TRUE, # Incluir los nombres de las filas (TRUE) o no (FALSE)
append = TRUE, # En este caso, sí agregar al archivo existente "resultados.xlsx"
showNA = TRUE) # Si TRUE, los NA serán celdas vacías
Para facilitar el trabajo, se va a representar la orientación con un gráfico de barras adyacentes
barplot(tabla_orientacion_fi,
main = "Frecuencia de glaciares según su orientación",
ylim = c(0,40),
xlab = "Orientación",
ylab = "Frecuencia (%)",
axes = TRUE,
beside = TRUE,
legend.text = rownames(tabla_orientacion_fi)) # Leyenda)
A continuación, la superficie de los glaciares según la categoría
tabla_area_categoria_ni <- table(datos$montana, datos$area_categoria)
tabla_area_categoria_ni
##
## Diminuto Pequeño Mediano Grande Muy_grande Enorme
## Alpes 1119 1033 650 606 423 61
## Cáucaso 53 288 535 425 298 38
tabla_area_categoria_fi <- prop.table(tabla_area_categoria_ni, margin= 1) * 100
tabla_area_categoria_fi
##
## Diminuto Pequeño Mediano Grande Muy_grande Enorme
## Alpes 28.751285 26.541624 16.700925 15.570401 10.868448 1.567318
## Cáucaso 3.237630 17.593158 32.681735 25.962126 18.204032 2.321319
tabla_area_categoria <- rbind(tabla_area_categoria_ni, tabla_area_categoria_fi)
View(tabla_area_categoria)
barplot(tabla_area_categoria_fi,
main = "Frecuencia de glaciares según su superficie",
ylim = c(0,40),
xlab = "Superficie (ha)",
ylab = "Frecuencia (%)",
axes = TRUE,
beside = TRUE,
legend.text = rownames(tabla_area_categoria_fi)) # Leyenda)
A continuación, la superficie de los glaciares según la categoría
tabla_lmax_categoria_ni <- table(datos$montana, datos$lmax_categoria)
tabla_lmax_categoria_ni
##
## Muy_corto Corto Mediano Largo Muy_largo
## Alpes 569 1403 758 674 488
## Cáucaso 24 302 511 437 363
tabla_lmax_categoria_fi <- prop.table(tabla_lmax_categoria_ni, margin= 1) * 100
tabla_lmax_categoria_fi
##
## Muy_corto Corto Mediano Largo Muy_largo
## Alpes 14.619733 36.048304 19.475848 17.317575 12.538541
## Cáucaso 1.466097 18.448381 31.215638 26.695174 22.174710
tabla_lmax_categoria <- rbind(tabla_lmax_categoria_ni, tabla_lmax_categoria_fi)
View(tabla_lmax_categoria)
barplot(tabla_lmax_categoria_fi,
main = "Frecuencia de glaciares según su longitud",
ylim = c(0,40),
xlab = "Longitud (m)",
ylab = "Frecuencia (%)",
axes = TRUE,
beside = TRUE,
legend.text = rownames(tabla_lmax_categoria_fi)) # Leyenda)
sig_est_area <- t.test(area ~ montana,
data = datos,
var.equal = TRUE)
sig_est_area
##
## Two Sample t-test
##
## data: area by montana
## t = -3.557, df = 5527, p-value = 0.0003782
## alternative hypothesis: true difference in means between group Alpes and group Cáucaso is not equal to 0
## 95 percent confidence interval:
## -0.3572045 -0.1033671
## sample estimates:
## mean in group Alpes mean in group Cáucaso
## 0.5367325 0.7670183
sig_est_zmin <- t.test(zmin ~ montana,
data = datos,
var.equal = TRUE)
sig_est_zmin
##
## Two Sample t-test
##
## data: zmin by montana
## t = -28.254, df = 5527, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group Alpes and group Cáucaso is not equal to 0
## 95 percent confidence interval:
## -311.3723 -270.9674
## sample estimates:
## mean in group Alpes mean in group Cáucaso
## 2798.545 3089.715
sig_est_zmed <- t.test(zmediana ~ montana,
data = datos,
var.equal = TRUE)
sig_est_zmed
##
## Two Sample t-test
##
## data: zmediana by montana
## t = 4.3489, df = 5527, p-value = 1.393e-05
## alternative hypothesis: true difference in means between group Alpes and group Cáucaso is not equal to 0
## 95 percent confidence interval:
## 42.48848 112.23436
## sample estimates:
## mean in group Alpes mean in group Cáucaso
## 2959.173 2881.811
sig_est_pendiente <- t.test(pendiente ~ montana,
data = datos,
var.equal = TRUE)
sig_est_pendiente
##
## Two Sample t-test
##
## data: pendiente by montana
## t = -5.4578, df = 5527, p-value = 5.031e-08
## alternative hypothesis: true difference in means between group Alpes and group Cáucaso is not equal to 0
## 95 percent confidence interval:
## -1.8189039 -0.8575463
## sample estimates:
## mean in group Alpes mean in group Cáucaso
## 27.72775 29.06597
sig_est_lmax <- t.test(lmax ~ montana,
data = datos,
var.equal = TRUE)
sig_est_lmax
##
## Two Sample t-test
##
## data: lmax by montana
## t = -10.66, df = 5527, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group Alpes and group Cáucaso is not equal to 0
## 95 percent confidence interval:
## -462.7354 -318.9779
## sample estimates:
## mean in group Alpes mean in group Cáucaso
## 758.9234 1149.7801
sig_est <- as.data.frame(cbind(sig_est_area$p.value, sig_est_zmin$p.value, sig_est_zmed$p.value,
sig_est_pendiente$p.value, sig_est_lmax$p.value))
names(sig_est)[1:5] <- c("Área", "Zmin", "Zmed","Pendiente", "Lmax")
La función `require() sirve para indicar si un paquete está cargado o no en memoria.
require("Hmisc")
Como en el caso anterior, extraemos las variables cuatitativas en un
dataframe para proceder al cálculo del coeficiente de correlación de
Pearson. La función rcorr()
del paquete Hmisc
devuelve esos coeficientes más el p-value.
Primero se realizará el cálculo para los Alpes.
datos_corr_alpes <- subset(datos, montana =="Alpes")
datos_corr_alpes <- datos_corr_alpes[ , c(5: 12)]
datos_corr_alpes_matriz_corr <- rcorr(as.matrix(datos_corr_alpes))
datos_corr_alpes_matriz_corr
## lon lat area zmin zmax zmediana pendiente lmax
## lon 1.00 0.79 -0.03 -0.14 -0.25 -0.21 -0.21 -0.06
## lat 0.79 1.00 0.00 -0.16 -0.24 -0.20 -0.21 -0.03
## area -0.03 0.00 1.00 -0.29 0.31 0.07 -0.18 0.86
## zmin -0.14 -0.16 -0.29 1.00 0.50 0.80 0.19 -0.40
## zmax -0.25 -0.24 0.31 0.50 1.00 0.84 0.10 0.45
## zmediana -0.21 -0.20 0.07 0.80 0.84 1.00 0.14 0.08
## pendiente -0.21 -0.21 -0.18 0.19 0.10 0.14 1.00 -0.28
## lmax -0.06 -0.03 0.86 -0.40 0.45 0.08 -0.28 1.00
##
## n= 3892
##
##
## P
## lon lat area zmin zmax zmediana pendiente lmax
## lon 0.0000 0.0319 0.0000 0.0000 0.0000 0.0000 0.0003
## lat 0.0000 0.7792 0.0000 0.0000 0.0000 0.0000 0.1158
## area 0.0319 0.7792 0.0000 0.0000 0.0000 0.0000 0.0000
## zmin 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
## zmax 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
## zmediana 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
## pendiente 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
## lmax 0.0003 0.1158 0.0000 0.0000 0.0000 0.0000 0.0000
A continuación se realiza el mismo procedimiento para el Cáucaso.
datos_corr_caucaso <- subset(datos, montana =="Cáucaso")
datos_corr_caucaso <- datos_corr_caucaso[ , c(5: 12)]
datos_corr_caucaso_matriz_corr <- rcorr(as.matrix(datos_corr_caucaso))
datos_corr_caucaso_matriz_corr
## lon lat area zmin zmax zmediana pendiente lmax
## lon 1.00 -0.95 -0.02 0.44 0.34 -0.42 0.07 0.01
## lat -0.95 1.00 0.04 -0.32 -0.21 0.45 -0.05 0.02
## area -0.02 0.04 1.00 -0.29 0.38 0.08 -0.25 0.91
## zmin 0.44 -0.32 -0.29 1.00 0.60 0.30 0.16 -0.36
## zmax 0.34 -0.21 0.38 0.60 1.00 0.41 0.06 0.43
## zmediana -0.42 0.45 0.08 0.30 0.41 1.00 -0.07 0.05
## pendiente 0.07 -0.05 -0.25 0.16 0.06 -0.07 1.00 -0.31
## lmax 0.01 0.02 0.91 -0.36 0.43 0.05 -0.31 1.00
##
## n= 1637
##
##
## P
## lon lat area zmin zmax zmediana pendiente lmax
## lon 0.0000 0.3542 0.0000 0.0000 0.0000 0.0058 0.8275
## lat 0.0000 0.0825 0.0000 0.0000 0.0000 0.0438 0.3927
## area 0.3542 0.0825 0.0000 0.0000 0.0012 0.0000 0.0000
## zmin 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
## zmax 0.0000 0.0000 0.0000 0.0000 0.0000 0.0089 0.0000
## zmediana 0.0000 0.0000 0.0012 0.0000 0.0000 0.0055 0.0368
## pendiente 0.0058 0.0438 0.0000 0.0000 0.0089 0.0055 0.0000
## lmax 0.8275 0.3927 0.0000 0.0000 0.0000 0.0368 0.0000
ACTIVIDAD: - ¿Cuál es el papel de las variables geográficas en las características y dimensiones de los glaciares de ambas montañas?
Para determinar qué papel tiene la orientación en la altitud de la
línea de equilibrio regional, se seleccionarán los glaciares más
pequeños (categorías Diminuto y Pequeño) y se cuantificará su valor
mediante el análisis de la variable zmediana
.
ela <- subset(datos, area_categoria == "Diminuto" | area_categoria == "Pequeño")
ela_promedio <- aggregate(zmediana ~ montana + orientacion,
FUN = mean,
data = ela,
na.rm=TRUE)
ela_promedio
## montana orientacion zmediana
## 1 Alpes N 2845.524
## 2 Cáucaso N 3370.082
## 3 Alpes E 2929.107
## 4 Cáucaso E 3376.807
## 5 Alpes S 3084.214
## 6 Cáucaso S 3463.136
## 7 Alpes W 2991.893
## 8 Cáucaso W 3432.792
Con etos datos se puede realizar también un diagrama de caja y bigotes conjunto.
boxplot(zmediana ~ orientacion * montana,
data = ela,
notch = TRUE,
main = "Altura de la ELA según la orientación de los glaciares",
xlab = "",
ylab = "Altura [m]",
xaxt = "n",
col = c("Blue", "Green", "Yellow", "Red"))
axis(side = 1, at = c(2.5, 6.5), labels = c("Alpes", "Cáucaso"))
legend("bottomright", title = "Orientacion",
legend = c("N", "E", "S", "W"),
fill = c("Blue", "Green", "Yellow", "Red"),
cex = 0.75) # Cambiar el tamaño)