29 Day 27: Resources
29.1 Data Source: USGS Geochemical and Mineralogical Resources
st_read(here::here("data/me-counties.json"), stringsAsFactors = FALSE) %>%
st_set_crs(4326) -> me
## Reading layer `cb_2015_maine_county_20m' from data source `/Users/hrbrmstr/books/30-day-map-challenge/data/me-counties.json' using driver `TopoJSON'
## Simple feature collection with 16 features and 10 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: -71.08434 ymin: 43.05975 xmax: -66.9502 ymax: 47.45684
## epsg (SRID): NA
## proj4string: NA# https://pubs.usgs.gov/sir/2017/5118/sir20175118_geo.php
if (!all(file.exists(c(here::here("data/Zn_tif.zip"),
here::here("data/Fe_tif.zip"),
here::here("data/U_tif.zip"),
here::here("data/Cu_tif.zip"))))) {
download.file(
url = c(
"https://pubs.usgs.gov/sir/2017/5118/elements/Zinc/Zn_tif.zip",
"https://pubs.usgs.gov/sir/2017/5118/elements/Iron/Fe_tif.zip",
"https://pubs.usgs.gov/sir/2017/5118/elements/Uranium/U_tif.zip",
"https://pubs.usgs.gov/sir/2017/5118/elements/Copper/Cu_tif.zip"
),
destfile = c(
here::here("data/Zn_tif.zip"),
here::here("data/Fe_tif.zip"),
here::here("data/U_tif.zip"),
here::here("data/Cu_tif.zip")
),
method = "libcurl"
)
unzip(zipfile = here::here("data/Zn_tif.zip"), exdir = here::here("data/usgs"))
unzip(zipfile = here::here("data/Fe_tif.zip"), exdir = here::here("data/usgs"))
unzip(zipfile = here::here("data/U_tif.zip"), exdir = here::here("data/usgs"))
unzip(zipfile = here::here("data/Cu_tif.zip"), exdir = here::here("data/usgs"))
}Each individual mineral transformation takes a while so we’re going to pre-process it all and save off our hard work in the event we want to use it later or we have session issues:
if (!file.exists(here::here("data/all-four.rds"))) {
copper_r <- read_stars(here::here("data/usgs/A_Cu.tif")) # read raster
st_as_sf(copper_r[1], as_points = FALSE, merge = TRUE) %>% # make polygons
lwgeom::st_make_valid() %>% # ensure polygons are valid
st_transform(st_crs(me)) %>% # our CRS
st_intersection(me) %>% # only maine
rename(concentration = A_Cu.tif) %>% # need a common name for the value
mutate(element = "Copper (Cu)") -> copper_me # name the element
zinc_r <- read_stars(here::here("data/usgs/A_Zn.tif"))
st_as_sf(zinc_r[1], as_points = FALSE, merge = TRUE) %>%
lwgeom::st_make_valid() %>%
st_transform(st_crs(me)) %>%
st_intersection(me) %>%
rename(concentration = A_Zn.tif) %>%
mutate(element = "Zinc (Zn)") -> zinc_me
iron_r <- read_stars(here::here("data/usgs/A_Fe.tif"))
st_as_sf(iron_r[1], as_points = FALSE, merge = TRUE) %>%
lwgeom::st_make_valid() %>%
st_transform(st_crs(me)) %>%
st_intersection(me) %>%
rename(concentration = A_Fe.tif) %>%
mutate(element = "Iron (Fe)") -> iron_me
uranium_r <- read_stars(here::here("data/usgs/A_U.tif"))
st_as_sf(uranium_r[1], as_points = FALSE, merge = TRUE) %>%
st_buffer(0) %>%
st_transform(st_crs(me)) %>%
st_intersection(me) %>%
rename(concentration = A_U.tif) %>%
mutate(element = "Uranium (U)") -> uranium_me
all_four <- rbind(copper_me, zinc_me, iron_me, uranium_me)
all_four <- mutate(all_four, quantile = gtools::quantcut(concentration, 5))
saveRDS(all_four, here::here("data/all-four.rds"))
}
all_four <- readRDS(here::here("data/all-four.rds"))29.2 Drawing the Map
ggplot() +
geom_sf(data = all_four, aes(fill = quantile), size = 0.125, color = "white") +
geom_sf(data = me, fill = NA) +
scale_fill_viridis_d(
name = "Concentration Quantile\n(relative to scale of individual mineral)",
option = "magma", labels = 1:5
) +
facet_wrap(~element, ncol = 4) +
guides(
fill = guide_legend(title.position = "top")
) +
labs(
title = "Geochemistry of Maine",
subtitle = "A-horizon (surface) samples of four selected elements",
caption = "Data source: USGS <pubs.usgs.gov/sir/2017/5118/sir20175118_geo.php> • #30DayMapChallenge"
) +
coord_sf(datum = NA) +
theme_ipsum_es(grid="", strip_text_family = font_es_bold) +
theme(legend.position = "bottom") +
theme(legend.direction = "horizontal")