Erupciones Volcanicas en Indonesia

Indonesia, es un archipiélago en Oceanía con más de 13.000 islas, está situado en el “Anillo de Fuego del Pacífico”, una serie de fallas sísmicas alrededor del océano Pacífico. Hay 127 volcanes activos repartidos por todo el país, desde el norte de Sumatra en el oeste hasta la isla de Banda en el este. Con una población de más de 270 millones de personas, este evento natural es de alto riesgo para las personas que viven en el país.

Este informe hablará sobre los datos de las erupciones mortales desde 1300 hasta la fecha, para visualizarlos y entender mejor este fenómeno.

Primero, cargaremos las librerias necesarias:

library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.4.4     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(dplyr)
library(lubridate)
library(janitor)

## 
## Attaching package: 'janitor'
## 
## The following objects are masked from 'package:stats':
## 
##     chisq.test, fisher.test

library(skimr)
library(cowplot)

## 
## Attaching package: 'cowplot'
## 
## The following object is masked from 'package:lubridate':
## 
##     stamp

options(warn=-1, message=FALSE)

Ahora, crearemos una función para la función para ver la gráfica:

fig <- function(width, heigth){
     options(repr.plot.width = width, repr.plot.height = heigth)
}

##Cargamos los datos. Ahora, cargamos los datos necesarios para el analisis:

volcano_events <- read_csv("volcano-events-1300-2021.csv", show_col_type=FALSE)

Colocamos nombres a las columnas

v_events <- volcano_events %>% clean_names() 
# glimpse(v_events)
# skim_without_charts(v_events)

Un dato faltante importante es el VEI; según Newhall y Self, 1982, a las erupciones históricas que fueron definitivamente explosivas, pero que no contienen otra información descriptiva, se les asigna un VEI predeterminado de 2.

v_events <- v_events %>% mutate(vei = if_else(is.na(vei), 2, vei),
                               total_deaths = if_else(is.na(total_deaths), 0, total_deaths))
# glimpse(v_events)

Observación de los Datos

Ahora, miremos los datos de diferentes ángulos ¿Qué nos revelará cada ángulo?

Localización de las Montañas

Solamente en la isla Java, hay 114 erupciones, seguida por Nusatenggara con 24 islas. Más de 150 personas viven en la isla Java..

fig(12, 7)
n_eruptbyisland <- v_events %>% 
  group_by(island) %>% 
  summarize(cnt_mount=n())

ggplot(data=n_eruptbyisland, aes(x=reorder(island, -cnt_mount), y=cnt_mount, fill=island)) +
  geom_bar(stat="identity", width=0.9, color="white") +
  geom_text(aes(label=cnt_mount), position=position_stack(vjust=0.5), size=5) +
  labs(title="Número de Erupciones por isla") +
  xlab("Islas") +
  ylab("Número de Montañas") +
  theme(legend.position = "none", text = element_text(size=18))

Top 10 erupciones de montañas

muestra el top diez de erupciones.

fig(12, 7)
n_eruptbymountain <- v_events %>% 
  group_by(name) %>% 
  summarise(cnt_bymount=n()) %>% 
  arrange(desc(cnt_bymount))

top10_eruptbymountain <- n_eruptbymountain %>% top_n(10) %>% arrange(desc(cnt_bymount))

## Selecting by cnt_bymount

ggplot(data=top10_eruptbymountain, aes(x=reorder(name, -cnt_bymount), y=cnt_bymount, fill=name)) +
  geom_bar(stat="identity", width=0.9, color="white") +
  scale_x_discrete(guide = guide_axis(angle=45)) +
  geom_text(aes(label=cnt_bymount), position=position_stack(vjust=0.5), size=5) +
  labs(title="Número de Erupciones por Montaña") +
  xlab("Nombre de la Montaña") +
  ylab("Erupciones") +
  theme(legend.position = "none", text = element_text(size=18))

Indice de Explosividad volcánica (VEI)

Un esquema de clasificación ampliamente utilizado para describir el tamaño de las erupciones explosivas. Se basa principalmente en la masa o volumen de la erupción de un depósito Newhall y Self, 1982. A las erupciones históricas que fueron definitivamente explosivas, pero que no contienen otra información descriptiva, se les asigna un VEI predeterminado de 2. La puntuación del VEI es de 0 a 8:

VEI	Descripción General	Altura de la columna de nubes (km)	Volumen (m3)	Descripción Cualitativa	Clasificación	Con que frecuenia	Ejemplo
0	no explosivo	<0.1	1x104	Gentil	Hawaiian	Diario	Kilauea
1	Pequeño	0.1-1	1x106	Efusivo	Haw/Strombolian	Diario	Stromboli
2	Moderado	1-5	1x107	Explosivo	Strom/Vulcanian	Semanal	Galeras, 1992
3	Moderado-Grande	3-15	1x108	Explosivo	Vulcanian	Anual	Ruiz, 1985
4	Grande	10-25	1x109	Explosivo	Vulc/Plinian	Decenio	Galunggung, 1982
5	Muy Grande	>25	1x1010	Cataclismico	Plinian	Siglo	St. Helens, 1981
6	>25 km	1x1011	Paroxistico	Plin/Ultra-Plinian	Siglo	Krakatau, 1883
7	>25 km	1x1012	Colosal	Ultra-Plinian	1000’s of years	Tambora, 1815
8	>25 km	>1x1012	Colosal	Ultra-Plinian	10,000’s of years	Yellowstone, 2 Ma

vei_eruptions <- v_events %>% group_by(vei) %>% summarise(cnt_vei=n())
ggplot(data=vei_eruptions, aes(x=reorder(vei, -cnt_vei), y=cnt_vei, fill=vei)) +
  geom_bar(stat="identity", width=0.8, color="white", position = position_dodge(width = 0.1)) +
  geom_text(aes(label=cnt_vei), position=position_stack(vjust=0.5), size=5) +
  labs(title="Puntaje VEI") +
  xlab("Puntaje VEI") +
  ylab("Erupciones") +
  theme(legend.position = "none", text = element_text(size=18))

## Erupciones por Tipo de Montaña Muestra el número de erupciones por tipo de montaña.

n_eruptbytype <- v_events %>% 
  group_by(type) %>% 
  summarise(cnt_bytype=n())
ggplot(data=n_eruptbytype, aes(x=reorder(type, -cnt_bytype), y=cnt_bytype, fill=type)) +
      geom_bar(stat="identity", width=0.8, color="white", position = position_dodge(width = 0.1)) +
      geom_text(aes(label=cnt_bytype), position=position_stack(vjust=0.5), size=5) +
      scale_x_discrete(guide = guide_axis(angle=45)) +
      labs(title="Tipo de Montañas") +
      xlab("Tipo de Montañas") +
      ylab("Número de Montañas") +
      theme(legend.position = "none", text = element_text(size=18))

Por Elevación de Montaña

fig(12,7)
e_mounts <- v_events %>% select(name, elevation_m) %>% distinct()
ggplot(data=e_mounts, aes(x=reorder(name, -elevation_m), y=elevation_m, fill=name)) +
  geom_bar(stat="identity", width=0.8, color="white", position = position_dodge(width = 0.1)) +
  scale_x_discrete(guide = guide_axis(angle=45)) +
  # geom_text(aes(label=elevation_m), position=position_stack(vjust=0.5), size=5) +
  labs(title="Mountain Elevation") +
  xlab("Nombre de Montañas") +
  ylab("Elevación") +
  theme(legend.position = "none", text = element_text(size=15))

Erupciones Mortales

Aunque las erupciones extremas con gran explosión y que causen muertes son raras, es necesario observar los eventos.

n_deathspereruptions <- v_events %>% select(name, erupt_date, total_deaths) %>% arrange(desc(total_deaths)) %>% top_n(10)

## Selecting by total_deaths

ggplot(data=n_deathspereruptions, aes(x=reorder(name, -total_deaths), y=total_deaths, fill=name)) +
  geom_bar(stat="identity", width=0.9, color="white") +
  scale_x_discrete(guide = guide_axis(angle=45)) +
  geom_text(aes(label=total_deaths), position=position_stack(vjust=0.5), size=5) +
  labs(title="Total de Muertes por Montaña") +
  xlab("Nombre de Montañas") +
  ylab("Muertes") +
  theme(legend.position = "none", text = element_text(size=15))

La erupción más mortífera fue la del Monte Tambora el 8 de octubre de 1822, donde mató a 60.000 personas, seguida por el Monte Krakatau el 27 de agosto de 1883, y luego Kelut entró en erupción dos veces en el Año Nuevo de 1586 y Mei, el 19 de agosto de 1919, mató a 10.000 y 5.000 respectivamente.

Ahora echemos un vistazo a las erupciones mortales que ocurrieron en los últimos 30 años.

Analisis de la Distribución Normal

El gráfico de histograma muestra la visualización de erupciones basadas en erupciones anuales.

n_mounts <- v_events %>% select(year) %>% filter(year>1900) 

ggplot(n_mounts, aes(year)) +
  geom_histogram(aes(x=year, y=..density..), bins = 50, fill="#d3d3d3", color="black") +
  stat_function(fun=dnorm, args = list(mean=mean(n_mounts$year), sd=sd(n_mounts$year)), color="red") +
  geom_density(color="blue") +
  ggtitle("Histograma de Erupciones") +
  xlab("Año") +
  ylab("Densidad") +
  theme_bw()

Gráfico de distribución normal basado en la elevación de las montañas.

e_mounts <- v_events %>% select(elevation_m)
ggplot(e_mounts, aes(elevation_m)) +
  geom_histogram(aes(x=elevation_m, y=..density..), bins = 50, fill="#d3d3d3", color="black") +
  stat_function(fun=dnorm, args = list(mean=mean(e_mounts$elevation_m), sd=sd(e_mounts$elevation_m)), color="red") +
  geom_density(color="blue") +
  ggtitle("Histograma de Erupciones") +
  xlab("Elevación (m)") +
  ylab("Densidad") +
  theme_bw()

Este gráfico muestra montañas con una elevación entre 2000 y 3000 (m) que son los volcanes más activos.

##Analisis QQ Un gráfico QQ visualiza datos basados en las cantidades de la variable proporcionada frente a los cuantiles que existirían si los datos estuvieran distribuidos normalmente. Los datos que siguen la distribución normal deben estar en línea con una pendiente establecida. Incluir stat_qq() genera un gráfico QQ. El siguiente ejemplo genera un gráfico QQ de la variable erupciones (n_erupts).

y     <- quantile(n_mounts$year, c(0.25, 0.75)) # Encuentra el 1er y 3er cuartil
x     <- qnorm( c(0.25, 0.75))           #Encuentra los valores normales correspondientes en el eje x
slope <- diff(y) / diff(x)               # Computa la pendiente
int   <- y[1] - slope * x[1]             # Comput la intercepción

ggplot(n_mounts, aes(sample = year)) +
  stat_qq() + 
  geom_abline(intercept = int, slope = slope, color = "blue", linewidth = 2, alpha = 0.5 )

y     <- quantile(e_mounts$elevation_m, c(0.25, 0.75)) # Encuentra el 1er y 3er cuartil
x     <- qnorm( c(0.25, 0.75))           # Encuentra los valores normales correspondientes del en el eje x
slope <- diff(y) / diff(x)               # Computa la pendiente
int   <- y[1] - slope * x[1]             # Computa la intercepción

ggplot(e_mounts, aes(sample = elevation_m)) +
  stat_qq() + 
  geom_abline(intercept = int, slope = slope, color = "blue", size = 2, alpha = 0.5 )

Plot de Caja

Ahora veamos en el plot de caja por tipo de montaña vs VEI.

p <- ggplot(v_events, aes(island, y=elevation_m))
p + geom_boxplot(fill = "white", colour = "#3366FF")

##Correlación de Elevación Este gráfico de dispersión indica la tendencia de las erupciones. El eje X es el año, los círculos son el registro de erupciones, el tamaño indica el número de heridos y el tamaño del círculo es el número de muertos. Aunque todavía es discutible el argumento de cómo se registran y documentan los datos, especialmente para registros muy antiguos de hace más de 100 años. El gráfico muestra que las erupciones son cada vez mayores en número e intensidad.

# Cargamos librerias para que le diagrama se vea mejor
library(colorspace) # para darken()
library(ggrepel)    # para geom_text_repel()

fig(12, 7)
island_cols <- c("#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#999999", "#889EAF", "#D4B499")
e_events <- v_events 
ggplot(data=e_events, aes(x=year, y=elevation_m)) +
  geom_point(aes(color = island, fill = island, size = total_deaths), alpha = 0.5, shape = 21) +
  geom_smooth(
        aes(color = "y ~ log(x)", fill = "y ~ log(x)"),
        method = 'loess', formula = y~log(x), se = FALSE, fullrange = TRUE) +
  scale_size_continuous(range = c(3, 10)) +
  scale_color_manual(values = darken(island_cols, 0.5)) +
  scale_fill_manual(values = island_cols) +
  labs(title="Total de Erupciones en los Ultimos 30 años")

Eruptions After 1970, zoom the data to show only from 1970 records. Mountains become more active, since the year 1800, however the impact is getting lower, both in the number of deaths as well as injuries.

fig(18, 12)
e_events_1970 <- v_events %>% filter(year>1970 & elevation_m>0)
# e_events_top10 <- e_events_1970 %>% filter(vei>2)
ggplot(data=e_events_1970, aes(x=year, y=elevation_m)) +
  geom_point(aes(color = island, fill = island, size = vei), alpha = 0.5, shape = 21) +
  geom_smooth(aes(color = "y ~ log(x)", fill = "y ~ log(x)", alpha=.1), size=3,
    method = 'loess', formula = y~log(x), se = FALSE, fullrange = TRUE) +
  geom_text_repel(
    aes(label = name),
    color = "black",
    size = 11/.pt, # tamaño de fuente 9pt
    point.padding = 0.1, 
    box.padding = .6,
    min.segment.length = 0,
    max.overlaps = 1000,
    seed = 7654
  ) +
  scale_size_continuous(range = c(3, 15)) +
  scale_color_manual(values = darken(island_cols, 0.5)) +
  scale_fill_manual(values = island_cols) +
  scale_x_continuous(
    name = "Eventos de Erupciones en los años 1970-2021",
  ) +
  scale_y_continuous(
    name = "Elevación (m)",
    limits = c(10, 4100),
    breaks = c(1, 900, 1800, 2700, 3600, 3900),
    expand = c(0, 0)
  ) +
  labs(title="Total de Erupciones en los Ultimos 30 años") +
  theme_minimal_hgrid(12, rel_small = 1) +
  theme(
    legend.position = "bottom",
    legend.justification = "right",
    legend.text = element_text(size = 9),
    legend.box.spacing = unit(0, "pt")
  )

fig(12, 7)
island_cols <- c("#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#999999", "#889EAF", "#D4B499")
ve_data <- v_events %>% filter(elevation_m>0) %>% select(island, type, name, vei, elevation_m, total_deaths)
ggplot(ve_data, aes(vei, elevation_m)) +
  geom_smooth(aes(color = "y ~ log(x)", fill = "y ~ log(x)"),
        method = 'lm', formula = y~log(x), se = FALSE, fullrange = TRUE) +
  geom_point(aes(color = island, fill = island, size = total_deaths), alpha = 0.5, shape = 21) +
  labs(title="Elevación (m) vs VEI") +
  xlab("VEI") +
  ylab("Elevación (m)") +
  scale_color_manual(values = darken(island_cols, 0.3)) +
  scale_fill_manual(values = island_cols) +
  scale_size_continuous(range = c(3, 21)) +
  guides(size=guide_legend("Total de Muertes")) +
  theme_minimal_hgrid(12, rel_small = 1) + 
  theme(
    legend.position = "right",
    legend.justification = "right",
    legend.text = element_text(size = 9),
    legend.box.spacing = unit(3, "pt")
  )

fig(12, 7)
island_cols <- c("#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#999999", "#889EAF", "#D4B499")
ve_data <- v_events %>% filter(elevation_m>0 & total_deaths<10000) %>% select(island, type, name, vei, elevation_m, total_deaths)
ggplot(ve_data, aes(total_deaths, elevation_m)) +
  geom_smooth(aes(color = "y ~ log(x)", fill = "y ~ log(x)"),
        method = 'lm', formula = y~log(x), se = FALSE, fullrange = TRUE) +
  geom_point(aes(color = island, fill = island, size = total_deaths), alpha = 0.5, shape = 21) +
  labs(title="Elevación (m) vs Total de Muertes") +
  xlab("Total de Muertes") +
  ylab("Elevación (m)") +
  scale_color_manual(values = darken(island_cols, 0.3)) +
  scale_fill_manual(values = island_cols) +
  scale_size_continuous(range = c(3, 21)) +
  guides(size=guide_legend("Total de Muertes")) +
  guides(colour=guide_legend("Isla")) +
  theme_minimal_hgrid(12, rel_small = 1)

# instale el paquete ("sf")
# Instale el paquete ("rworldmap")
# instale el paquete ("geos")
library(geos)
library(sf)        # para manipulacion de objectos con caracteristicas simples

## Linking to GEOS 3.11.2, GDAL 3.8.2, PROJ 9.3.1; sf_use_s2() is TRUE

library(rworldmap) # para getMap()

## Loading required package: sp

## ### Welcome to rworldmap ###

## For a short introduction type :   vignette('rworldmap')

Mapa de Volcanes

La visualización de datos en un mapa es muy importante para mostrar dónde están los objetos que se comentan y al mismo tiempo visualizar más datos que se ajusten al mapa.

require("maps")

## Loading required package: maps

## 
## Attaching package: 'maps'

## The following object is masked from 'package:purrr':
## 
##     map

fig(20, 9)
island_cols <- c("#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#999999", "#889EAF", "#D4B499")
type_cols <- c("#6F69AC", "#95DAC1", "#FFEBA1", "#FD6F96")
mydata <- e_events %>% 
    dplyr::group_by(island, type, name, longitude, latitude) %>%
    dplyr::summarise(sum_erpt = n(), m_vei = mean(vei), m_elev = max(elevation_m), m_deaths = max(total_deaths), .groups="drop") %>%
    select(island, type, name, m_vei, m_elev, longitude, latitude, m_deaths)
global <- map_data("world") # Datos de Latitud y Longitud del mundo 
ggplot() + 
    geom_polygon(data = global, aes(x=long, y = lat, group = group), 
        fill = "gray85", color = "gray80") + 
    coord_fixed(1.3) +
    xlim(94,142) + ylim(-11,7.5) +
    geom_point(data = mydata, aes(x = longitude, y = latitude, color = type, fill = type, size = m_deaths), 
        alpha = 0.8, show.legend = F, shape=24) +
    geom_text_repel(data = mydata,
        aes(x = longitude, y = latitude, label = name),
        color = "grey30",
        size = 11/.pt, # Tamaño de fuente 9 pt
        point.padding = 0.1, 
        box.padding = .6,
        min.segment.length = 0,
        max.overlaps = 1000,
        seed = 7654
      ) +
    ggtitle ("Erupciones de Volcanes") + 
    scale_color_manual(values = darken(type_cols, 0.7)) +
    scale_fill_manual(values = type_cols) +
    scale_size_continuous(range = c(3, 9)) +
    guides(size=guide_legend("Total de Muertes")) +
    guides(color=guide_legend("Tipo")) +
    theme_minimal_hgrid(12, rel_small = 1) + 
    theme(
      legend.position = "bottom",
      legend.justification = "right",
      legend.text = element_text(size = 9),
      legend.box.spacing = unit(3, "pt"))

Erupciones Volcanicas en Indonesia

Santiago Tinoco Román

2024-02-14

Erupciones Volcanicas en Indonesia

Observación de los Datos

Localización de las Montañas

Top 10 erupciones de montañas

Indice de Explosividad volcánica (VEI)

Por Elevación de Montaña

Erupciones Mortales

Analisis de la Distribución Normal

Plot de Caja

Mapa de Volcanes

Conclusión