-
Notifications
You must be signed in to change notification settings - Fork 2
/
make_map.R
152 lines (118 loc) · 5.06 KB
/
make_map.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
library(tidyverse)
library(toxEval)
library(maptools)
library(maps)
library(nhdplusTools)
library(sf)
path_to_data <- Sys.getenv("PASSIVE_PATH")
options(dplyr.summarise.inform = FALSE)
source(file = "read_chemicalSummary.R")
source(file = "R/report/combo_plot2.R")
source(file = "R/report/stack_plots.R")
site_info <- tox_list$chem_site
site_info <- site_info %>%
bind_rows(filter(site_info, SiteID == "04231600"))
site_info$SiteID[duplicated(site_info$SiteID)] <- "04232007"
lake <- st_read("Great_Lakes-shp/Great_Lakes.shp")
glri_states <- map_data("state") %>%
filter(region %in% c("wisconsin", "illinois", "michigan",
"ohio","indiana", "minnesota",
"pennsylvania", "new york"))
site_info <- site_info %>%
filter(!is.na(dec_lat),
!is.na(dec_lon))
p <- sf::st_as_sf(data.frame(lon = site_info$dec_lon,
lat = site_info$dec_lat),
coords = c("lon", "lat"), crs = 4326)
# bbox <- sf::st_bbox(c(xmin = min(site_info$dec_lon, na.rm = TRUE)-3, ymin = min(site_info$dec_lat, na.rm = TRUE)-3,
# xmax = max(site_info$dec_lon, na.rm = TRUE)+3, ymax = max(site_info$dec_lat, na.rm = TRUE)+3),
# crs = "+proj=longlat +datum=WGS84 +no_defs")
bbox <- sf::st_bbox(c(xmin = min(glri_states$long, na.rm = TRUE), ymin = min(glri_states$lat, na.rm = TRUE),
xmax = max(glri_states$long, na.rm = TRUE), ymax = max(glri_states$lat, na.rm = TRUE)),
crs = "+proj=longlat +datum=WGS84 +no_defs")
data <- plot_nhdplus(bbox = bbox, streamorder = 4, actually_plot = FALSE)
# saveRDS(data, file = "plot_nhd_out_3.rds")
data <- readRDS("plot_nhd_out_2.rds")
flowlines <- data$flowline
indexes <- get_flowline_index(flowlines,
p,
search_radius = 0.1,
max_matches = 1)
indexes <- left_join(sf::st_sf(id = c(1:nrow(p)),
geom = sf::st_geometry(p)),
indexes, by = "id")
flowlines_small <- flowlines %>%
filter(REACHCODE %in% indexes$REACHCODE)
flowlines_med <- flowlines %>%
filter(TerminalPa %in% unique(flowlines_small$TerminalPa))
site_info <- prep_site_list(site_info)
graphData <- graph_chem_data(chemical_summary = chemicalSummary)
gd <- graphData %>%
group_by(site) %>%
summarize(maxmax = sum(meanEAR)) %>%
left_join(site_info, by = c("site"="SiteID")) %>%
arrange(`Short Name`)
gdplot <- ggplot() +
geom_polygon(aes(x = long, y = lat, group = group),
data = glri_states, fill = "grey90",
color = "white") +
geom_sf(data = lake, fill = "lightblue") +
geom_point(data = gd,
aes(x = dec_lon, y = dec_lat,
fill = maxmax),
size = 3, pch = 21, color = "black") +
scale_fill_gradient(low = "grey", high = "red") +
# scale_fill_gradient(low = "white", high = "red") +
geom_sf(data = flowlines_med , color = "blue") +
theme_minimal() +
theme(panel.grid = element_blank(),
axis.text = element_blank(),
axis.title = element_blank(),
panel.border = element_blank(),
legend.position="none",
plot.background = element_rect(fill = "transparent",colour = NA)
)
gdplot
ggsave(gdplot, file = "map3.png", bg = "transparent")
library(grid)
color_map <- class_colors(tox_list)
gd_stack <- graphData %>%
left_join(site_info, by = c("site"="SiteID"))
top4 <- gd_stack %>%
group_by(Class) %>%
summarise(sumEARs = sum(meanEAR)) %>%
arrange(desc(sumEARs))
gd_stack$Class <- as.character(gd_stack$Class)
gd_stack$Class[!gd_stack$Class %in% as.character(top4$Class[1:5])] <- "Other (14)"
gd_stack$Class <- factor(gd_stack$Class,
levels = c(as.character(top4$Class[1:5]), "Other (14)"))
names(color_map)[names(color_map) == "Other"] <- "Other (14)"
stack1 <- ggplot() +
geom_col(data = gd_stack,
aes(x=site, y=meanEAR, fill = Class),
position = position_stack()) +
geom_hline(yintercept = 0) +
theme_minimal() +
ylab("Exposure-Activity Ratio (EAR)") +
facet_grid(. ~ site_grouping , switch = "x",
scales="free", space="free") +
scale_y_continuous(breaks = scales::pretty_breaks(n = 5)) +
scale_fill_manual(values = color_map, drop=TRUE) +
guides(fill = guide_legend(nrow = 1)) +
theme(strip.background = element_blank(),
strip.text.x = element_text(size = 9),
panel.grid.minor = element_blank(),
panel.grid.major.x = element_blank(),
panel.grid.major.y = element_blank(),
legend.position = "bottom",
legend.title = element_blank(),
legend.text = element_text(size = 8),
# legend.position = c(.95, .95),
# legend.justification = c("right", "top"),
# legend.box.just = "right",
# legend.margin = margin(6, 6, 6, 6),
axis.text = element_blank(),
axis.title.x = element_blank(),
plot.background = element_rect(fill = "transparent",colour = NA))
stack1
ggsave(stack1, file = "stack1.png", bg = "transparent")