Análisis de Erupciones Volcanicas en Indonesia

Indonesia se ubica entre Oceania y el Sudeste Asiatico. Se compone de 17000 islas y tiene 1,904,569 kilometros cuadrados y 278 millones de habitantes. Su geografía se compone principalmente de vocanes formados por la subducción entre la placa Eurosiatica e Indoastrualiana. Los 127 volcanes forman parte del cinturon de Fuego del Pacifico.

Este informe pretende mostrar datos recolectados desde 1300 hasta el año 2021 con el proposito de explorar y comprender mejor este fenomeno natural y así poder tener una mejor gestión de riesgo en el futuro.

Primero se debe crear una funcion para generar las graficas

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

Ahora se deben agregar las librerias que se necesitan

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 se debe cargar los datos que se deben cargar para el informe

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

Se deben nombrar las columnas de la tabla que se hará

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

Un dato faltante importante es el VEI; según Newhal 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)

Se visuliza número de erupciones por islas

Solo en Java hay 114 de erupciones, más de 150 personas viven en esta isla.

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("Isla") +
  ylab("Número de montañas") +
  theme(legend.position = "none", text = element_text(size=18))

Se visualiza el Top 10 de las erupciones por montaña

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="Number of Eruption by Mountain") +
  xlab("Mountains Name") +
  ylab("Eruptions") +
  theme(legend.position = "none", text = element_text(size=18))

Clasificacíon de explosividad de los volcanes (VEI)

Newhall y Self (1982) proponen que un sistema de clasificacíon adecuado se basa en la masa o el volumen del deposito. Como se explicó antes las explosiones que fueron explosivas pero no se conoce mas información tiene una calificación de 2. La escala va de 0 a 8

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

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="Erupciones por tipo de montaña") +
      xlab("Tipo de montaña") +
      ylab("Número montañas") +
      theme(legend.position = "none", text = element_text(size=18))

Erupciones por elevación de las montañas

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="Elevación de las montañas") +
  xlab("Nombre de la montaña") +
  ylab("Elevación") +
  theme(legend.position = "none", text = element_text(size=15))

Erupciones mortales

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 muertes por montaña") +
  xlab("Nombre de las montañas") +
  ylab("Muertes") +
  theme(legend.position = "none", text = element_text(size=15))

La erupcíon del volcán Tambora sucedió el 10 de abril de 1815, este alcanzó una magnitud 7 en el indice de explosividad volcánica (VEI) y provocó la muerte de 60.000 personas y desencadenó efectos importantes en el clima de Europa

Histograma de las erupciones por año con sus respectivas densidades

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

Tabla de distribución con información de elevación de 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()

Se denota la relación que hay entre la altura de la montaña y su actividad volcanica.

Análisis 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. Se genera un gráfico QQ de la variable “Erupciones”, con respecto al año y a la elvación

y     <- quantile(n_mounts$year, c(0.25, 0.75)) # Encontrar el primer y tercer cuartil
x     <- qnorm( c(0.25, 0.75))           # Encuentre los valores normales coincidentes en el eje x
slope <- diff(y) / diff(x)               # Calcular la pendiente de la línea
int   <- y[1] - slope * x[1]             # Clacular el intercepto de la linea

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

y     <- quantile(e_mounts$elevation_m, c(0.25, 0.75)) 
x     <- qnorm( c(0.25, 0.75))           
slope <- diff(y) / diff(x)               
int   <- y[1] - slope * x[1]            

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

Visualización de un diagrama de caja (Tipo vs VEI)

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

Grafico de dispersión que indica la tendencia de las erupciones

Donde el eje x= año, círculos= registro de erupciones, el tamaño indica el numero de heridos y el tamaño de los circulos indica las muertes.

library(colorspace) 
library(ggrepel)    
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 erupciones en los ultimos 30 años")

Como se muestra en el gráfico, la erupciones han aumentado en número.

Erupciones despues de 1970

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, # font size 9 pt
    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 = "Eruption Events Year 1970-2021",
  ) +
  scale_y_continuous(
    name = "Elevation (m)",
    limits = c(10, 4100),
    breaks = c(1, 900, 1800, 2700, 3600, 3900),
    expand = c(0, 0)
  ) +
  labs(title="Total 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")
  )

Al ampliar el gráfico se evidencia que aunque la actividad volcánica haya aumentado, el impacto en la sociedad ha disminuido pues hay un menor número de muertos y heridos.

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ónn (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 Deaths")) +
  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 Deaths") +
  ylab("Elevation (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) 

library(sf)        
## Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1; sf_use_s2() is TRUE
library(rworldmap)
## Loading required package: sp
## ### Welcome to rworldmap ###
## For a short introduction type :   vignette('rworldmap')

Visualización en mapa de las erupciones

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") # World longitude and latitude data 
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, # font size 9 pt
        point.padding = 0.1, 
        box.padding = .6,
        min.segment.length = 0,
        max.overlaps = 1000,
        seed = 7654
      ) +
    ggtitle ("Erupciones volcanicas") + 
    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 Deaths")) +
    guides(color=guide_legend("Type")) +
    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"))

Conclusiones

Apartir de la visualización y analisis de las gráficas se pueden concluir multiples cosas:

  • Las erupciones no estan uniformemente distribuidas en las islas.

  • Java es la isla más volcanica de la región y teniendo en cuenta su alta densidad poblacional es vital tener una buena gestión y mitigación del riesgo.

  • La escala VEI permite comprender que la mayoría de las erupciones tienen una clasificación baja, lo que sugiere una menor explosividad, pocos heridos y/o muertes

  • Como se menciono en la respectiva gráfica, a mayor altura de la montaña (principlamente entre 2000 y 3000 m), los volcanes tienen mayor actividad