library(tidyverse)
library(RColorBrewer)
library(leaflet)
library(leaflegend)
library(readxl)
library(tidyr)
library(knitr)Project 2 Volcano Mapping
Data 110 Spring 2026
Introduction
My dataset is the Global Volcanism Program dataset (2025), Volcanoes of the World, distributed by the Smithsonian. I chose volcanoes because one of my discord servers always comments on some volcanoes, particularly Kilauea, are active. I joined two of the GVP’s datasets: the Holocene volcano list, and the Eruption search, which between the two of them contained variables such as the volcano name, the coordinates, the country of the volcano, the start year, month, and day of an eruption, sometimes the end year, month, and day, the type of volcano and it’s landform, as well as the Volcano Explosivity Index (VEI).
The VEI is is logarithmic scale from 0-8 that describes the size of an eruption based on magnitude and intensity of the ash cloud. This means that not all eruptions will have a VEI because it was an eruption (such as many of those in Hawai’i) that did not have an ash cloud, or because it was not observed. I set the NAs for VEI to zero, as I used VEI as the basis for my mapping radius.
I created a consolidated Start and End Date columns from the separate start year, start month, start day, end year, end month, and end day columns, which also allowed me to create new variable Duration in Days using a built in R function. Again, some eruptions did not have complete start and end dates, and thus had NAs for duration.
I filtered for volcanoes in Indonesia, as a friend recommended looking into them, and I found a page on the GVP website listing volcanoes with populations living nearby, and the 30km radius list was filled Indonesian entries. I also filtered for eruptions in or after 1900, as the turn of a century seemed like a reasonable starting point to get to under 800 observations.
References
Global Volcanism Program, 2025. [Database] Volcanoes of the World (v. 5.3.5; 31 Mar 2026). Distributed by Smithsonian Institution, compiled by Venzke, E. https://doi.org/10.5479/si.GVP.VOTW5-2025.5.3
https://www.nps.gov/subjects/volcanoes/volcanic-explosivity-index.htm
https://volcano.si.edu/faq/Populations_around_Volcanoes.cfm
https://roh.engineering/posts/2021/05/map-symbols-and-size-legends-for-leaflet/
Load Libraries and Data
setwd("~/Schol Stuff/Montgomery College 2025/Data 110 Data Visualization/Volcano Smithsonian")
volcano_list = read_excel("GVP_Volcano_List_Holocene_202603241259.xlsx")
eruption_list = read_excel("GVP_Eruption_Search_Result.xlsx")head(eruption_list)# A tibble: 6 × 24
`Volcano Number` `Volcano Name` `Eruption Number` `Eruption Category`
<dbl> <chr> <dbl> <chr>
1 300130 Karymsky 22609 Confirmed Eruption
2 273010 Bulusan 22608 Confirmed Eruption
3 334050 Northern EPR at 9.8°N 22613 Confirmed Eruption
4 300260 Klyuchevskoy 22610 Confirmed Eruption
5 371020 Reykjanes 22612 Confirmed Eruption
6 252120 Ulawun 22606 Confirmed Eruption
# ℹ 20 more variables: `Area of Activity` <chr>, VEI <chr>,
# `VEI Modifier` <chr>, `Start Year Modifier` <chr>, `Start Year` <dbl>,
# `Start Year Uncertainty` <chr>, `Start Month` <dbl>,
# `Start Day Modifier` <chr>, `Start Day` <dbl>,
# `Start Day Uncertainty` <chr>, `Evidence Method (dating)` <chr>,
# `End Year Modifier` <chr>, `End Year` <dbl>, `End Year Uncertainty` <lgl>,
# `End Month` <dbl>, `End Day Modifier` <chr>, `End Day` <dbl>, …Data cleaning
remove spaces in variable names. Make VEI and Start and End Day Uncertainty Numeric
#removing spaces in variable names
names(eruption_list) <- gsub(" ","",names(eruption_list))
names(volcano_list) <- gsub(" ","",names(volcano_list))
# making these variables numeric. Absolutely need VEI. Uncertainty might not use because it might not be what I thought it was
eruption_list$VEI<- as.numeric(eruption_list$VEI)
eruption_list$StartDayUncertainty <- as.numeric(eruption_list$StartDayUncertainty)
eruption_list$EndDayUncertainty <- as.numeric(eruption_list$EndDayUncertainty) Create Start and End date columns in eruption list
Create a column that shows the duration between the start and end date of the eruption.
eruption1 <- eruption_list |>
mutate(StartDate = make_date(StartYear, StartMonth, StartDay), .after = StartDay) |>
mutate(EndDate = make_date(EndYear, EndMonth, EndDay), .after = EndDay ) |>
mutate(Duration_days = as.numeric(difftime(EndDate, StartDate, units ="days"))) |>
# need as.numeric to have it not include "days" in the columns
select(-Latitude, -Longitude) # get rid of these so I don't have duplicate columns once joined
head(eruption1)# A tibble: 6 × 25
VolcanoNumber VolcanoName EruptionNumber EruptionCategory AreaofActivity VEI
<dbl> <chr> <dbl> <chr> <chr> <dbl>
1 300130 Karymsky 22609 Confirmed Erupt… <NA> NA
2 273010 Bulusan 22608 Confirmed Erupt… <NA> NA
3 334050 Northern E… 22613 Confirmed Erupt… Tica hydrothe… 0
4 300260 Klyuchevsk… 22610 Confirmed Erupt… <NA> NA
5 371020 Reykjanes 22612 Confirmed Erupt… <NA> NA
6 252120 Ulawun 22606 Confirmed Erupt… <NA> NA
# ℹ 19 more variables: VEIModifier <chr>, StartYearModifier <chr>,
# StartYear <dbl>, StartYearUncertainty <chr>, StartMonth <dbl>,
# StartDayModifier <chr>, StartDay <dbl>, StartDate <date>,
# StartDayUncertainty <dbl>, `EvidenceMethod(dating)` <chr>,
# EndYearModifier <chr>, EndYear <dbl>, EndYearUncertainty <lgl>,
# EndMonth <dbl>, EndDayModifier <chr>, EndDay <dbl>, EndDate <date>,
# EndDayUncertainty <dbl>, Duration_days <dbl>Some eruptions only have a Star Year, if they have a full start Date, don’t have an end date– either way don’t have an end date nor duration.
# some NA values for VEi are becaues the eruption was too small, in that case zero in accurate. Others are because it wasn't observed, in that case it's *not*, but I don't have a way ot separating the two and I want vei as my mapping radius
eruption1$VEI[is.na(eruption1$VEI)] <- 0
head(eruption1)# A tibble: 6 × 25
VolcanoNumber VolcanoName EruptionNumber EruptionCategory AreaofActivity VEI
<dbl> <chr> <dbl> <chr> <chr> <dbl>
1 300130 Karymsky 22609 Confirmed Erupt… <NA> 0
2 273010 Bulusan 22608 Confirmed Erupt… <NA> 0
3 334050 Northern E… 22613 Confirmed Erupt… Tica hydrothe… 0
4 300260 Klyuchevsk… 22610 Confirmed Erupt… <NA> 0
5 371020 Reykjanes 22612 Confirmed Erupt… <NA> 0
6 252120 Ulawun 22606 Confirmed Erupt… <NA> 0
# ℹ 19 more variables: VEIModifier <chr>, StartYearModifier <chr>,
# StartYear <dbl>, StartYearUncertainty <chr>, StartMonth <dbl>,
# StartDayModifier <chr>, StartDay <dbl>, StartDate <date>,
# StartDayUncertainty <dbl>, `EvidenceMethod(dating)` <chr>,
# EndYearModifier <chr>, EndYear <dbl>, EndYearUncertainty <lgl>,
# EndMonth <dbl>, EndDayModifier <chr>, EndDay <dbl>, EndDate <date>,
# EndDayUncertainty <dbl>, Duration_days <dbl>head(volcano_list)# A tibble: 6 × 14
VolcanoNumber VolcanoName Country VolcanicRegionGroup VolcanicRegion
<dbl> <chr> <chr> <chr> <chr>
1 210010 West Eifel Volcanic … Germany European Volcanic … Central Europ…
2 210020 Chaine des Puys France European Volcanic … Western Europ…
3 210030 Olot Volcanic Field Spain European Volcanic … Western Europ…
4 210040 Calatrava Volcanic F… Spain European Volcanic … Western Europ…
5 211004 Colli Albani Italy European Volcanic … Italian Penin…
6 211010 Campi Flegrei Italy European Volcanic … Italian Penin…
# ℹ 9 more variables: VolcanoLandform <chr>, PrimaryVolcanoType <chr>,
# ActivityEvidence <chr>, LastKnownEruption <chr>, Latitude <dbl>,
# Longitude <dbl>, `Elevation(m)` <dbl>, TectonicSetting <chr>,
# DominantRockType <chr>no obvious cleaning needed on volcano list
move ahead to joining
joined <- left_join(eruption1, volcano_list, by = c("VolcanoNumber"))
head(joined)# A tibble: 6 × 38
VolcanoNumber VolcanoName.x EruptionNumber EruptionCategory AreaofActivity
<dbl> <chr> <dbl> <chr> <chr>
1 300130 Karymsky 22609 Confirmed Erupt… <NA>
2 273010 Bulusan 22608 Confirmed Erupt… <NA>
3 334050 Northern EPR at … 22613 Confirmed Erupt… Tica hydrothe…
4 300260 Klyuchevskoy 22610 Confirmed Erupt… <NA>
5 371020 Reykjanes 22612 Confirmed Erupt… <NA>
6 252120 Ulawun 22606 Confirmed Erupt… <NA>
# ℹ 33 more variables: VEI <dbl>, VEIModifier <chr>, StartYearModifier <chr>,
# StartYear <dbl>, StartYearUncertainty <chr>, StartMonth <dbl>,
# StartDayModifier <chr>, StartDay <dbl>, StartDate <date>,
# StartDayUncertainty <dbl>, `EvidenceMethod(dating)` <chr>,
# EndYearModifier <chr>, EndYear <dbl>, EndYearUncertainty <lgl>,
# EndMonth <dbl>, EndDayModifier <chr>, EndDay <dbl>, EndDate <date>,
# EndDayUncertainty <dbl>, Duration_days <dbl>, VolcanoName.y <chr>, …Choosing Indonesia as focus:
Consulted this page on the Global Volcanism Program’s site that shows population within certain radii of volcanoes. Indonesia has many volcanoes on the 30km list.
Indonesia <- joined |> # filter just for Eruptions in Indonesia
filter(Country == "Indonesia")
head(Indonesia)# A tibble: 6 × 38
VolcanoNumber VolcanoName.x EruptionNumber EruptionCategory AreaofActivity
<dbl> <chr> <dbl> <chr> <chr>
1 263340 Raung 22607 Confirmed Erupt… <NA>
2 264230 Lewotolok 22601 Confirmed Erupt… <NA>
3 263200 Dieng Volcanic C… 22600 Confirmed Erupt… <NA>
4 263340 Raung 22593 Confirmed Erupt… <NA>
5 261230 Dempo 22592 Confirmed Erupt… <NA>
6 261170 Kerinci 22568 Confirmed Erupt… <NA>
# ℹ 33 more variables: VEI <dbl>, VEIModifier <chr>, StartYearModifier <chr>,
# StartYear <dbl>, StartYearUncertainty <chr>, StartMonth <dbl>,
# StartDayModifier <chr>, StartDay <dbl>, StartDate <date>,
# StartDayUncertainty <dbl>, `EvidenceMethod(dating)` <chr>,
# EndYearModifier <chr>, EndYear <dbl>, EndYearUncertainty <lgl>,
# EndMonth <dbl>, EndDayModifier <chr>, EndDay <dbl>, EndDate <date>,
# EndDayUncertainty <dbl>, Duration_days <dbl>, VolcanoName.y <chr>, …length(Indonesia$VolcanoNumber) # see how many eruptions I have[1] 1347Filter further for eruptions in or after 1900
Indonesia2 <- Indonesia |>
filter(StartYear >= 1900)
length(Indonesia2$VolcanoNumber)[1] 763Yay, I’m under 800.
Arrange by how many eruptions each volcano has had
most_active <- Indonesia2 |>
group_by(VolcanoName.x, VolcanoLandform) |>
count() |>
arrange(desc(n))
most_active# A tibble: 66 × 3
# Groups: VolcanoName.x, VolcanoLandform [66]
VolcanoName.x VolcanoLandform n
<chr> <chr> <int>
1 Raung Composite 51
2 Karangetang Composite 50
3 Marapi Composite 49
4 Krakatau Caldera 46
5 Soputan Composite 35
6 Kerinci Composite 33
7 Slamet Composite 32
8 Tengger Caldera Composite 32
9 Merapi Composite 31
10 Lokon-Empung Composite 26
# ℹ 56 more rowsThere’s 66 volcanoes and I really want to make a heatmap.
Visualization: Heatmap of most active Indonesian volcanoes by duration of eruptions
top_ten <- Indonesia2 |>
filter(VolcanoName.x %in% c('Raung', 'Karangetang', 'Marapi', 'Krakatau', 'Soputan', 'Kerinci', 'Slamet', 'Tengger Caldera', 'Merapi', 'Lokon-Empung'))
head(top_ten)# A tibble: 6 × 38
VolcanoNumber VolcanoName.x EruptionNumber EruptionCategory AreaofActivity
<dbl> <chr> <dbl> <chr> <chr>
1 263340 Raung 22607 Confirmed Eruption <NA>
2 263340 Raung 22593 Confirmed Eruption <NA>
3 261170 Kerinci 22568 Confirmed Eruption <NA>
4 263310 Tengger Caldera 22555 Confirmed Eruption <NA>
5 261140 Marapi 22550 Confirmed Eruption <NA>
6 267020 Karangetang 22577 Confirmed Eruption Kawah Utama, …
# ℹ 33 more variables: VEI <dbl>, VEIModifier <chr>, StartYearModifier <chr>,
# StartYear <dbl>, StartYearUncertainty <chr>, StartMonth <dbl>,
# StartDayModifier <chr>, StartDay <dbl>, StartDate <date>,
# StartDayUncertainty <dbl>, `EvidenceMethod(dating)` <chr>,
# EndYearModifier <chr>, EndYear <dbl>, EndYearUncertainty <lgl>,
# EndMonth <dbl>, EndDayModifier <chr>, EndDay <dbl>, EndDate <date>,
# EndDayUncertainty <dbl>, Duration_days <dbl>, VolcanoName.y <chr>, …plot1 <- top_ten |>
ggplot(aes(StartYear, VolcanoName.x, fill = Duration_days)) +
geom_tile(color = 'gray') +
scale_fill_gradientn(colors = brewer.pal(9, "OrRd" )) +
labs(title = "Top Ten Most Active Indonesian Volcanos by Duration",
subtitle = "1900 - 2025",
x = "Year",
y = "",
fill = "Duration of Eruption (days)",
caption = "Source: Global Volcanism Program, 2025. [Database] Volcanoes of the World (v. 5.3.5; 31 Mar 2026). Distributed by Smithsonian Institution, compiled by Venzke, E.") +
theme_bw()
plot1 
I’m glad I didn’t drop the NAs for Duration because I’m actually incredibly interested in where they are. I expected many of them to be towards the beginning of the 20th century due to less monitoring methods – or colonials ignoring any system the locals already had for documenting eruptions – but the ones nearing and past the 2000s are a surprise.
I did not expect the range of durations to be so large, which makes it hard to differentiate durations under 1000 days.
Vizualization: Number of Eruptions per volcano by VEI:
vei_counts <- Indonesia2 |>
group_by(VolcanoName.x, VEI) |>
count()
vei_counts# A tibble: 153 × 3
# Groups: VolcanoName.x, VEI [153]
VolcanoName.x VEI n
<chr> <dbl> <int>
1 Agung 3 1
2 Agung 5 1
3 Ambang 1 1
4 Arjuno-Welirang 0 1
5 Arjuno-Welirang 2 1
6 Awu 0 2
7 Awu 1 1
8 Awu 2 3
9 Awu 4 1
10 Banda Api 2 1
# ℹ 143 more rowsplot2 <- vei_counts |>
ggplot(aes(x = VolcanoName.x, y = n, fill = factor(VEI))) +
geom_col(position = "dodge") +
labs(title = 'Indonesia Eruptions per Volcano by VEI',
subtitle = 'since 1900',
x = "",
y = 'Number of Eruptions',
fill = 'Volcano Eruption Index',
caption = "Global Volcanism Program, 2025. [Database] Volcanoes of the World (v. 5.3.5; 31 Mar 2026). Distributed by Smithsonian Institution, compiled by Venzke, E.") +
scale_fill_brewer(palette = "Dark2") +
theme_bw() +
coord_flip() +
geom_text(aes(label = n), hjust = -.4, vjust = -.05, size = 3)
plot2
Pretty much the only thing legible from this plot is that VEI 2 is the most common among Indonesian volcanoes since 1900.
A Version You Can Actually Read
It’s clear from the bar graph of all 66 volcanoes that VEI 2 is the most common, but let’s see that for the ten most active, on a plot that’s actually legible.
plot3 <- top_ten |>
group_by(VolcanoName.x, VEI) |>
count() |>
ggplot(aes(x = VolcanoName.x, y = n, fill = factor(VEI))) +
geom_col(position = "dodge") +
labs(title = 'Indonesia Eruptions per Ten Most Active Volcano by VEI',
subtitle = 'since 1900',
x = "",
y = 'Number of Eruptions',
fill = 'Volcano Eruption Index',
caption = "Global Volcanism Program, 2025. [Database] Volcanoes of the World (v. 5.3.5; 31 Mar 2026). Distributed by Smithsonian Institution, compiled by Venzke, E.") +
scale_fill_brewer(palette = "Dark2") +
theme_bw() +
coord_flip() +
geom_text(aes(label = n), hjust = -.4, vjust = -.01, size = 3)
plot3
This shows that VEI 2 is the most common, with VEI 1 being the next most, The ten volcanoes did not all experience any of the other VEI levels.
Mapping!
Going back to the entire Indonesia2
Making popup first so I can test maps
Central point is approx -2.518722 S (Latitude), 118.015568 E (Longitude)
unique(Indonesia2$VEI)[1] 0 2 1 3 4 5#Created adjusted VEI so I'm not multiplying 0
indonesia3 <- Indonesia2 |>
mutate(vei_adj = VEI +1 )
#get rid of the (es) and (s) at the end of some of the Primary Volcano Type groups, so they don't form additonal groups
indonesia3$PrimaryVolcanoType <- str_remove(indonesia3$PrimaryVolcanoType, "\\(.*")
head(indonesia3)# A tibble: 6 × 39
VolcanoNumber VolcanoName.x EruptionNumber EruptionCategory AreaofActivity
<dbl> <chr> <dbl> <chr> <chr>
1 263340 Raung 22607 Confirmed Erupt… <NA>
2 264230 Lewotolok 22601 Confirmed Erupt… <NA>
3 263200 Dieng Volcanic C… 22600 Confirmed Erupt… <NA>
4 263340 Raung 22593 Confirmed Erupt… <NA>
5 261230 Dempo 22592 Confirmed Erupt… <NA>
6 261170 Kerinci 22568 Confirmed Erupt… <NA>
# ℹ 34 more variables: VEI <dbl>, VEIModifier <chr>, StartYearModifier <chr>,
# StartYear <dbl>, StartYearUncertainty <chr>, StartMonth <dbl>,
# StartDayModifier <chr>, StartDay <dbl>, StartDate <date>,
# StartDayUncertainty <dbl>, `EvidenceMethod(dating)` <chr>,
# EndYearModifier <chr>, EndYear <dbl>, EndYearUncertainty <lgl>,
# EndMonth <dbl>, EndDayModifier <chr>, EndDay <dbl>, EndDate <date>,
# EndDayUncertainty <dbl>, Duration_days <dbl>, VolcanoName.y <chr>, …unique(indonesia3$PrimaryVolcanoType)[1] "Stratovolcano" "Complex" "Caldera" "Lava dome" make the popup after I make the dataset I’m using for the maps!
eruptpop <- paste0(
"<b>Name: </b>", indonesia3$VolcanoName.x,",<b>Year: </b>", indonesia3$StartYear, "<br>",
"<b>Date: </b>", indonesia3$StartDate," <b>Duration: </b>", indonesia3$Duration_days, " Days <br>",
"<b>VEI: </b>", indonesia3$VEI,"<br>",
"<b>Type: </b>", indonesia3$VolcanoLandform,": ", indonesia3$PrimaryVolcanoType, "<br>",
"<b>Elevation: </b>", indonesia3$`Elevation(m)`," m <br>",
"<b>Coordinates: </b>", indonesia3$Latitude," S, ", indonesia3$Longitude," E <br>"
)Option 1: Layering each VEI
so smaller radi are on top and still accessabile
#split dataset into separate dfs per VEI so I can map separate layers
vei0 <- indonesia3 |>
filter(VEI == 0) # also includes the former NA values
vei1 <- indonesia3 |>
filter(VEI == 1)
vei2 <- indonesia3 |>
filter(VEI == 2)
vei3 <- indonesia3 |>
filter(VEI == 3)
vei4 <- indonesia3 |>
filter(VEI == 4)
vei5 <- indonesia3 |>
filter(VEI == 5)
head(vei5)# A tibble: 1 × 39
VolcanoNumber VolcanoName.x EruptionNumber EruptionCategory AreaofActivity
<dbl> <chr> <dbl> <chr> <chr>
1 264020 Agung 16210 Confirmed Eruption <NA>
# ℹ 34 more variables: VEI <dbl>, VEIModifier <chr>, StartYearModifier <chr>,
# StartYear <dbl>, StartYearUncertainty <chr>, StartMonth <dbl>,
# StartDayModifier <chr>, StartDay <dbl>, StartDate <date>,
# StartDayUncertainty <dbl>, `EvidenceMethod(dating)` <chr>,
# EndYearModifier <chr>, EndYear <dbl>, EndYearUncertainty <lgl>,
# EndMonth <dbl>, EndDayModifier <chr>, EndDay <dbl>, EndDate <date>,
# EndDayUncertainty <dbl>, Duration_days <dbl>, VolcanoName.y <chr>, …# setting a color per vei for testing purposes
leaflet() |>
setView(lng = 118.015568, lat = -2.518722, zoom = 5) |>
addProviderTiles("OpenTopoMap") |>
addCircles(
data = vei5,
radius = 1000 * 2**vei5$vei_adj,
color = 'red',
fillColor = 'red',
fillOpacity = 0.5,
popup = eruptpop) |>
addCircles(
data = vei4,
radius = 1000 * 2**vei4$vei_adj,
color = 'orange',
fillColor = 'orange',
fillOpacity = 0.5,
popup = eruptpop) |>
addCircles(
data = vei3,
radius = 1000 * 2**vei3$vei_adj,
color = 'purple',
fillColor = 'purple',
fillOpacity = 0.5,
popup = eruptpop) |>
addCircles(
data = vei2,
radius = 1000 * 2**vei2$vei_adj,
color = 'limegreen',
fillColor = 'limegreen',
fillOpacity = 0.8,
popup = eruptpop) |>
addCircles(
data = vei1,
radius = 1000 * 2**vei1$vei_adj,
color = 'lightblue',
fillColor = 'lightblue',
fillOpacity = 0.8,
popup = eruptpop)Assuming "Longitude" and "Latitude" are longitude and latitude, respectively
Assuming "Longitude" and "Latitude" are longitude and latitude, respectively
Assuming "Longitude" and "Latitude" are longitude and latitude, respectively
Assuming "Longitude" and "Latitude" are longitude and latitude, respectively
Assuming "Longitude" and "Latitude" are longitude and latitude, respectivelyeruptions with the same VEI at the same location are hidden behind each other clickable. This might be fixable by jittering the coordinates in each layer, but I’m going to move on to clustering.
each layer would also also needs it’s own pop up
Final Map: Clustering
pal <- colorFactor(palette = c("red", "orange", "green", "purple"), domain = indonesia3$PrimaryVolcanoType)
leaflet() |>
setView(lng = 118.015568, lat = -2.518722, zoom = 5) |>
addProviderTiles("OpenTopoMap") |>
addCircleMarkers(data = indonesia3,
lng = ~Longitude, lat = ~Latitude,
radius = 2**indonesia3$vei_adj,
color = ~pal(PrimaryVolcanoType),
stroke = TRUE, fillOpacity = 0.8,
clusterOptions = markerClusterOptions(),
popup = eruptpop) |>
addLegend("bottomright",
pal = pal,
values = indonesia3$PrimaryVolcanoType,
title = "Volcano Type") |>
addLegendSize(
values = indonesia3$VEI,
title = 'Volcano Explositivity Index',
shape = 'circle',
color = 'gray',
orientation = 'horizontal',
position = "bottomleft")a map that doesn’t work
Stratovolcano are extremely common, and each volcano stay the same type, so there’s not much variation to the color. Tried to change the color of the circles to gradient based on elevation, but it won’t work
# library(viridis)
# pal <- colorNumeric(palette = "inferno", domain = indonesia3$`Elevation(m)`) |>
#
# leaflet() |>
# setView(lng = 118.015568, lat = -2.518722, zoom = 5) |>
# addProviderTiles("OpenTopoMap") |>
# addCircleMarkers(data = indonesia3,
# lng = ~Longitude, lat = ~Latitude,
# radius = 2**indonesia3$vei_adj,
# color = ~pal(indonesia3$`Elevation(m)`),
# stroke = TRUE, fillOpacity = 0.5,
# clusterOptions = markerClusterOptions(),
# popup = eruptpop)
# I either get “could not find function”pal”” or “simplewarning: some domain values being outside the scale and will be treated as NA” and no map.
Additional Comments
I was not expecting the duration of the eruptions to be so varied, so that was a surprised in the heatmap. My conception of an volcanic eruption was not something that could last 2000 days. A better approach would have been to create an additional categorical variable for bins of durations, and plot based on that. . Unfortunately, I made Duration in the first place to have another numerical variable. I had mostly wanted to make a heatmap because I hadn’t made one yet, and it seemed appropriate with a volcanic dataset. The bar graph by VEI showed just how prevalent VEI 2 is, which while initially surprsiing, makes sense as the scale is logarithmic and it is more massive each increase from each level to the next. Since I had folded the NA values into 0, I had feared that 0 would seem inflated, but it does not look like a problem.
Mapping wise, I tried two approaches. Layering each VEI would have worked, except there are so many of each VEI per volcano that ended in the same space. This probably would have been fixed by jittering the coordinates – an approach that didn’t work on the entire data set as a whole because smaller radii were still blocked under larger ones– but I also would have needed a separate popup per level, and clustering seemed a little easier at that point.
My final map clusters, on the last level, each volcano’s eruptions, which spiral out upon click, to individual circles with radius based on VEI, and colored by the volcano type. I tried adding a fourth variable in the shape of the individual marker, but the documentation for addMarkers and addAwesomeMarkers made it seem like size wasn’t something easily controlled by variable like addCircleMarker radius was. This proved moot, anyway, as the last map shows, the numeric gradient for color of circles didn’t work anyway. (Of course, now that I’m realizing I need to fix my legend, I find coding documentation that shows me what I could have done that might have fixed it.)