-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathgenerate_maps
141 lines (114 loc) · 4.02 KB
/
generate_maps
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
library(ggplot2)
library(ggmap)
library(ggpubr)
library(maps)
library(mapproj)
library(mapdata)
library(sf)
library(tidyverse)
library(geojsonio)
library(broom)
library(viridis)
setwd("../PHI_scaleit/SCALEIT_2.0")
load('ZCTAroche.RData')
load('ZCTAvitros.RData')
samples_collected<-read_csv("reshaped_redux.csv")
counts<-samples_collected %>% count(zcta)
spdf <- geojson_read("../Bay Area ZIP Codes.geojson", what = "sp")
spdf<-spdf[substr(spdf$zip, 1,3) %in% c("941") , ]
spdf<-spdf[spdf$zip != "94128" , ]
# fortify the data
spdf_fortified <- tidy(spdf, region = "zip")
spdf_fortified$id<-as.numeric(spdf_fortified$id)
# spdf_fortified_totalsamp<-spdf_fortified %>%
# left_join(. , as.data.frame(counts), by=c("id"="zcta"))
#
spdf_fortified_r = spdf_fortified %>%
left_join(. , as.data.frame(zcta_Roche[[4]]), by=c("id"="strat_cat"))
spdf_fortified_v = spdf_fortified %>%
left_join(. , as.data.frame(zcta_vitros[[4]]), by=c("id"="strat_cat"))
#spdf_fortified$mean[ is.na(spdf_fortified$mean )] = 0.0
# now we plot
ggplot() +
geom_polygon(data = spdf_fortified_totalsamp, aes(fill=n, x = long, y = lat, group = group)) +
theme_void() +
coord_map()+
scale_colour_brewer(palette = "RdBu")
spdf_fortified_v %>% filter(id == "94134")
seroprev_r<-
ggplot() +
geom_polygon(data = spdf_fortified_r, aes(fill = mean, x = long, y = lat, group = group)) +
theme_void() +
coord_map()+
#theme(legend.position = "none")+
labs_pubr()+
rremove("axis.text")+
rremove("xlab")+
rremove("ylab")+
scale_fill_distiller(palette ="RdBu")+
#scale_fill_viridis(limits = range(spdf_fortified_r$mean,spdf_fortified_v$mean, na.rm = TRUE))+
# scale_fill_viridis()+
ggtitle("Seroprevalence (Roche)")
seroprev_v<-
ggplot() +
geom_polygon(data = spdf_fortified_v, aes(fill = mean, x = long, y = lat, group = group)) +
theme_void() +
coord_map()+
#theme(legend.position = "none")+
theme(legend.title = element_blank())+
labs_pubr()+
rremove("axis.text")+
rremove("xlab")+
rremove("ylab")+
scale_fill_distiller(palette ="PuOr")+
#scale_fill_viridis(limits = range(spdf_fortified_r$mean,spdf_fortified_v$mean, na.rm = TRUE))+
#scale_fill_viridis()+
ggtitle("Seroprevalence (Vitros)")
ci_r<-
ggplot() +
geom_polygon(data = spdf_fortified_r, aes(fill = (hici-lowci), x = long, y = lat, group = group)) +
theme_void() +
coord_map()+
# theme(legend.position = "none")+
scale_fill_viridis(limits = range(spdf_fortified_r$mean,spdf_fortified_v$mean, na.rm = TRUE))+
#scale_fill_viridis()+
labs_pubr()+
rremove("axis.text")+
rremove("xlab")+
rremove("ylab")+
ggtitle("Credible interval range (Roche)")
ci_v<-ggplot() +
geom_polygon(data = spdf_fortified_v, aes(fill = (hici-lowci), x = long, y = lat, group = group)) +
theme_void() +
#theme(legend.position = "none")+
coord_map()+
scale_fill_viridis(limits = range(spdf_fortified_r$mean,spdf_fortified_v$mean, na.rm = TRUE))+
#scale_fill_viridis()+
labs_pubr()+
rremove("axis.text")+
rremove("xlab")+
rremove("ylab")+
ggtitle("Credible interval range (Vitros)")
spdf_fortified_v$diff<-spdf_fortified_v$mean - spdf_fortified_r$mean
spdf_fortified_v$diff
r_minus_vitros<-ggplot() +
geom_polygon(data = spdf_fortified_v, aes(fill = diff, x = long, y = lat, group = group)) +
theme_void() +
#theme(legend.position = "none")+
coord_map()+
labs_pubr()+
rremove("axis.text")+
rremove("xlab")+
rremove("ylab")+
scale_fill_distiller(palette = "BrBG")+
#scale_fill_viridis(limits = range(spdf_fortified_r$mean,spdf_fortified_v$mean, na.rm = TRUE))+
#scale_fill_viridis()+
ggtitle("Vitros - Roche")
ggarrange(seroprev_r, ci_r, seroprev_v,ci_v,
labels = c("a", "b", "c","d"),
ncol = 2, nrow = 2, align="hv")
ggsave("map-notsamescale-withci.tiff", units="in", width=8, height=5, dpi=300, compression = 'lzw')
ggarrange(seroprev_r, seroprev_v, r_minus_vitros,
labels = c("a", "b", "c"),
ncol = 3, nrow = 1, align="hv")
ggsave("map-notsamescale.tiff", units="in", width=9, height=4, dpi=300, compression = 'lzw')