-
-
Save paleolimbot/e854c098c6d95330c361014276b5a8db to your computer and use it in GitHub Desktop.
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") |
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.
I want to make Pourbaix diagrams that reflect a database I created. Is there a way to load my own database into tidyphreeqc?
Yes! You can do tidyphreeqc::phr_use_db_llnl(readLines("path_to_file.dat"))
Whoops...just tidyphreeqc::phr_use_db(readLines("path_to_file.dat"))
.
The 3D pourbaix diagram is beautiful and elegant but is a little complex. Will you consider creating a simple 2-D pourbaix diagram (for example, in 25 Celsius degree only), which is more often used. I guess it need convert regions into polygons as well, right?
Please I would like to inquire about how to calculate saturation index (SI) with tidyphreeqc package
I have a set of 50 wells analyzes
What does data entry look like?
How to do it?
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
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(...)poly$geometry
<sfc_MULTIPOLYGON> topoly$geometry
<sfc_MULTIPOLYGON>.Call
rlang::last_error()
to see a backtrace