INTRODUCCIÓN

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)

ANÁLISIS DE LA DISTRIBUCIÓN DE LOS GLACIARES SOBRE EL PLANETA

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)

FACTORES QUE EXPLICAN LA DISTRIBUCIÓN DE GLACIARES SOBRE EL PLANETA

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)

Cálculo de correlaciones entre variables

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/

ANÁLISIS REGIONAL DE LOS GLACIARES

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

Procesamiento inicial del regional

Antes de realizar cualquier tipo de análisis, es necesario realizar una serie de procedimientos para obtener y organizar los datos.

Importación del regional

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"

Seleccionar las cordilleras

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"

Transformación en factores de las variables cualitativas (alfanuméricas).

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)

Representación cartográfica

borrame_alpes <- subset(datos, 
                        montana =="Alpes")

borrame_alpes <- borrame_alpes[ , c(5:7)]

Transformación en un objeto espacial.

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:

  • Las columnas que contienen las coordenadas X (oeste-este) e Y (sur- norte).
  • El Sistema de Proyección de coordenadas, que en este caso estará basado en coordenadas geográficas.

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)

Importación de un modelo digital del terreno.

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

Cartografiar el mapa

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

Valores resumen

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

Diferencias entre cordilleras

Análisis de variables cuantitativas: área, zmin, zmediana, pendiente, lmax.

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)

  • Altitud mínima del frente
boxplot(zmin ~ montana,
        datos,
        main = "Altura mínima del frente",
        xlab = "Cordillera",
        ylab = "m",
        outline = FALSE,  
        notch = TRUE)

  • Altitud mediana
boxplot(zmediana ~ montana,
        datos,
        main = "Altura mediana",
        xlab = "Cordillera",
        ylab = "m",
        outline = FALSE,  
        notch = TRUE)

  • Pendiente
boxplot(pendiente ~ montana,
        datos,
        main = "Pendiente media del glaciar",
        xlab = "Cordillera",
        ylab = "º",
        outline = FALSE,  
        notch = TRUE)

  • Longitud máxima del glaciar
boxplot(lmax ~ montana,
        datos,
        main = "Longitud máxima del glaciar",
        xlab = "Cordillera",
        ylab = "m",
        outline = FALSE,  
        notch = TRUE)

Análisis de variables cualitativas: orientación

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)

Significación estadística de las diferencias entre ambas cordilleras

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
  • Altitud mínima (Zmin)
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
  • Altitud mínima (Zmin)
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
  • Pendiente
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
  • Longitud máxima:
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")

Relación entre los parámetros geográficos (latitud y longitud) y los morfológicos (área, zmin, zmediana…).

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?

  • ¿Cuál es la relación entre esos mismos parámetros?

Papel de la orientación en la altitud de la línea de equilibrio.

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)