Skip to content

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.

Copy link
Owner Author

paleolimbot commented Aug 16, 2018

Output:

pourbaix-anim-1

@Coutelot

This comment has been minimized.

Copy link

Coutelot commented Oct 23, 2019

Hello,
Your code is amazing thank you. It was working fine until a few weeks ago and package update I think, and since then I get this error message, and can't find the answer. Do you have the same problem? Thank you

Error: Can't cast poly$geometry <sfc_MULTIPOLYGON> to poly$geometry <sfc_MULTIPOLYGON>.
Call rlang::last_error() to see a backtrace

rlang::last_error()

message: Can't cast `poly$geometry` to `poly$geometry` . class: `vctrs_error_incompatible_cast` backtrace: 1. base::suppressWarnings(...) 18. vctrs:::vec_cast.default(x = x, to = to, x_arg = x_arg, to_arg = to_arg) 19. vctrs::stop_incompatible_cast(x, to, x_arg = x_arg, to_arg = to_arg) 20. vctrs:::stop_incompatible(...) 21. vctrs:::stop_vctrs(...)
@paleolimbot

This comment has been minimized.

Copy link
Owner Author

paleolimbot commented Oct 23, 2019

The problem is that sfc objects can't be combined because of the way that unnest() works. Off the top of my head, I don't know what the easiest way around this is, except to continue using an old version of tidyr. You could do so using remotes::install_github("cran/tidyr@61a7c3d"). In the next release of the sf package this should be fixed.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.