Instantly share code, notes, and snippets.

Embed
What would you like to do?
Manganeese 3D pourbaix diagram using PHREEQC and tidyphreeqc
library(tidyverse)
library(gganimate)
library(tidyphreeqc)
phr_use_db_minteq()
theme_set(theme_bw())
library(ggspatial)
result <- phr_run(
phr_solution_list(
pH = seq(3, 13, 0.1),
pe = seq(-10, 20, 0.1),
Mn = 10^2, # mmol / kg water
temp = seq(0, 70, 1)
),
phr_selected_output(
activities = c("MnOH+", "Mn+2", "Mn+3", "MnO4-2", "MnO4-"),
saturation_indices = c(
"Birnessite", "Bixbyite", "Hausmannite", "Manganite",
"Nsutite", "Pyrocroite", "Pyrolusite"
),
totals = "Mn",
temp = TRUE
)
)
result_df <- result %>%
as_tibble() %>%
select(pH, pe, temp = `temp(C)`, Mn = starts_with("Mn"), starts_with("la_"), starts_with("si_")) %>%
# these values aren't quite equal each time
mutate(Mn = 10^round(log10(Mn)))
result_activity <- result_df %>%
select(pH, pe, temp, Mn, starts_with("la_")) %>%
gather(species, log_activity, starts_with("la_")) %>%
mutate(species = str_remove(species, "^la_")) %>%
group_by(pH, pe, temp, Mn) %>%
summarise(species = species[which.max(log_activity)], log_activity = max(log_activity)) %>%
ungroup()
result_si <- result_df %>%
select(pH, pe, temp, Mn, starts_with("si_")) %>%
gather(species, saturation_index, starts_with("si_")) %>%
mutate(species = str_remove(species, "^si_")) %>%
group_by(pH, pe, temp, Mn) %>%
summarise(species = species[which.max(saturation_index)], saturation_index = max(saturation_index)) %>%
ungroup()
activity_polygons <- suppressWarnings(
result_activity %>%
mutate(species = factor(species), species_int = as.integer(species)) %>%
group_by(Mn, temp) %>%
nest() %>%
mutate(
rast = map(data, ~raster::rasterFromXYZ(tibble(.$pH, .$pe, .$species_int))),
poly = map2(
rast, data,
~sf::st_as_sf(raster::rasterToPolygons(.x, dissolve = TRUE)) %>%
mutate(species = levels(.y$species)[..species_int])
)
) %>%
unnest(poly) %>%
sf::st_as_sf() %>%
sf::st_simplify(preserveTopology = TRUE)
)
si_polygons <- suppressWarnings(
result_si %>%
mutate(species = factor(species), species_int = as.integer(species)) %>%
group_by(Mn, temp) %>%
nest() %>%
mutate(
rast = map(data, ~raster::rasterFromXYZ(tibble(.$pH, .$pe, .$species_int))),
poly = map2(
rast, data,
~sf::st_as_sf(raster::rasterToPolygons(.x, dissolve = TRUE)) %>%
mutate(species = levels(.y$species)[..species_int])
)
) %>%
unnest(poly) %>%
sf::st_as_sf() %>%
sf::st_simplify(preserveTopology = TRUE)
)
activity_labels <- activity_polygons %>%
sf::st_cast("MULTIPOLYGON") %>%
sf::st_cast("POLYGON", warn = FALSE) %>%
mutate(area = sf::st_area(.)) %>%
filter(area > 0.2) %>%
sf::st_centroid() %>%
df_spatial()
si_labels <- si_polygons %>%
sf::st_cast("MULTIPOLYGON") %>%
sf::st_cast("POLYGON", warn = FALSE) %>%
mutate(area = sf::st_area(.)) %>%
filter(area > 0.2) %>%
sf::st_centroid() %>%
df_spatial()
anim <- ggplot(NULL, aes(x, y)) +
geom_polygon(
aes(group = piece_id, fill = species),
data = si_polygons %>% df_spatial(),
col = NA
) +
geom_polygon(
aes(group = piece_id),
data = activity_polygons %>% df_spatial(),
col = "black", fill = NA
) +
# geom_label(
# aes(label = species),
# data = si_labels,
# col = "grey50", alpha = NA
# ) +
geom_label(
aes(label = species),
data = activity_labels,
col = "black", alpha = NA
) +
scale_fill_discrete(
breaks = c("Birnessite", "Bixbyite", "Hausmannite", "Manganite", "Nsutite", "Pyrocroite", "Pyrolusite"),
labels = c(
"Birnessite (MnO2)", "Bixbyite (Mn2O3)", "Hausmannite (Mn3O4)",
"Manganite (MnOOH)", "Nsutite (MnO2)", "Pyrocroite (Mn(OH)2)",
"Pyrolusite (MnO2)"
)
) +
labs(
title = "Dominant species of Mn at 0.1 mol/kg water, {round(frame_time)}°C",
caption = "As modelled by PHREEQC/minteq database using tidyphreeqc",
fill = "Most saturated species",
alpha = "Log Activity",
x = "pH", y = "pe"
) +
transition_time(temp)
animate(anim, nframes = 300, fps = 10, width = 600, height = 600)
anim_save("mn_animation.gif")
@paleolimbot

This comment has been minimized.

Owner

paleolimbot commented Aug 16, 2018

Output:

pourbaix-anim-1

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment