FINAL = Simulation results Weight: a_100_b_113_c_311_d_131
data <- syn_ind_urb_empID_sch_smedia_teen_coord
lt_step <- c(12,15,20,49)
lt_run <- 1:5
str_weight <- "a_100_b_113_c_311_d_131"
abm_df_vax <- list()
for (i in 1:length(lt_step)) {
df <- data
for (j in 1:length(lt_run)) {
path_name = sprintf(
"plot/11_model_0409_verification_validation/finl_%s/agent_step%s_weight_%s_7_days (%s).csv",
str_weight, lt_step[i], str_weight, lt_run[j])
print(path_name)
tt <- read.csv(path_name) %>% select(ind_id, vax_status) %>%
setnames(old="vax_status", new = sprintf("step%s_r%s",lt_step[i],lt_run[j]))
df <- left_join(df, tt, by=c("ind_id"="ind_id"))
}
abm_df_vax[[i]] <- df
}
names(abm_df_vax) <- sprintf("step_%01d",lt_step)
i=5
lt_run <- 1:5
df <- data
for (j in lt_run) {
path_name = sprintf(
"plot/11_model_0409_verification_validation/finl_%s/agent_simplified_weight_%s_7_days (1%s).csv",
str_weight, str_weight, lt_run[j])
print(path_name)
tt <- read.csv(path_name) %>% select(ind_id, vax_status) %>%
setnames(old="vax_status", new = sprintf("step%s_r%s",73,lt_run[j]))
df <- left_join(df, tt, by=c("ind_id"="ind_id"))
}
abm_df_vax[["step_73"]] <- df
cell_size = 2500
# create grid covers studyarea
c_grid <- c_county %>% st_make_grid(cellsize = cell_size) %>% st_as_sf() ## unit = m; 52*41=2132 CELLS
c_grid <- c_grid[lengths(st_intersects(c_county, c_grid))>0, ]
c_grid["idx"] <- sprintf("idx_%04d", 1:nrow(c_grid))
# individual agent location and count ind. in each grid
data <- syn_ind_urb_empID_sch_smedia_teen_coord[c("long","lat")] %>%
st_as_sf(., coords = c("long","lat"), crs=4326) %>%
st_transform(., crs = st_crs(c_grid))
c_grid["n_pop"] <- lengths(st_intersects(c_grid, data))
# count vaccined ind in each step
lt_step <- c(12,15,20,49,73)
lt_label_1 <- sprintf("avg_vax_s%s", lt_step)
lt_label_2 <- sprintf("p_avg_vax_s%s", lt_step)
n_run = 5
for (i in 1:length(lt_step)) {
df <- c_grid[c("idx","n_pop")]
x <- abm_df_vax[[i]]
for (j in 1:n_run) {
temp <- sprintf("vax_s%s_r%s", lt_step[i], j)
tt <- x[c(18,19,19+j)] %>% setnames(., new=c("long","lat","vaxed")) %>% filter(vaxed == 1) %>%
st_as_sf(., coords = c("long","lat"), crs=4326) %>%
st_transform(., crs = st_crs(c_grid))
df[temp] <- lengths(st_intersects(df, tt))
}
st_geometry(df) <- NULL
df[lt_label_1[i]] <- round(apply(df[,3:7], 1, mean))
df[lt_label_2[i]] <- round(100 * df[,8] / df["n_pop"],2)
df[is.na(df)] <- 0
c_grid <- left_join(c_grid, df[,c("idx",lt_label_1[i], lt_label_2[i])], by=c("idx"="idx"))
}
# save result
c_grid_pvax <- c_grid
## PLOT ## Vaccination Rate
n=40
p_alpha = 0.6
p_palette = "GnBu"
brks = seq(0, 100, length = n+1)
pdf(sprintf("plot/05_grid_vaxrate_process_%s.pdf", cell_size), width = 10, height = 5)
par(mfrow=c(1,4), oma=c(5,4,0,0)+0.1, mar=c(0,0,1,1)+0.1, bty="n")
i=1
p1 <- tm_shape(c_grid_pvax) +
tm_polygons(col = as.character(lt_label_2[[i]]), border.alpha = 0,
palette = p_palette, alpha = p_alpha, breaks = brks) +
tm_shape(c_county) + tm_polygons(alpha = 0) +
# tm_shape(c_urban) + tm_polygons(alpha = 0) +
tm_layout(main.title =lt_label_2[[i]], legend.show = F, frame = F)
i=2
p2 <- tm_shape(c_grid_pvax) +
tm_polygons(col = as.character(lt_label_2[[i]]), border.alpha = 0,
palette = p_palette, alpha = p_alpha, breaks = brks) +
tm_shape(c_county) + tm_polygons(alpha = 0) +
# tm_shape(c_urban) + tm_polygons(alpha = 0) +
tm_layout(main.title =lt_label_2[[i]], legend.show = F, frame = F)
i=3
p3 <- tm_shape(c_grid_pvax) +
tm_polygons(col = as.character(lt_label_2[[i]]), border.alpha = 0,
palette = p_palette, alpha = p_alpha, breaks = brks) +
tm_shape(c_county) + tm_polygons(alpha = 0) +
# tm_shape(c_urban) + tm_polygons(alpha = 0) +
tm_layout(main.title =lt_label_2[[i]], legend.show = F, frame = F)
i=4
p4 <- tm_shape(c_grid_pvax) +
tm_polygons(col = as.character(lt_label_2[[i]]), border.alpha = 0,
palette = p_palette, alpha = p_alpha, breaks = brks) +
tm_shape(c_county) + tm_polygons(alpha = 0) +
# tm_shape(c_urban) + tm_polygons(alpha = 0) +
tm_layout(main.title =lt_label_2[[i]], legend.show = F, frame = F)
i=5 ## legend
p5 <- tm_shape(c_grid_pvax) +
tm_polygons(col = as.character(lt_label_2[[i]]), border.alpha = 0,
palette = p_palette, alpha = p_alpha, breaks = brks) +
tm_shape(c_county) + tm_polygons(alpha = 0) +
tm_layout(main.title =lt_label_2[[i]], legend.show = T, frame = F, legend.outside = T)
tmap_arrange(p1, p2, p3, p4, p5, widths = c(0.5,0.5), nrow = 1)
dev.off()
## PLOT ## n of Vaccinated Individuals
n=50
p_alpha = 0.6
cols = mako(n=(n),alpha =1, begin = 1, end = 0)
brks = seq(0, 7000, length = n+1)
pdf(sprintf("plot/06_grid_vaxabs_process_%s.pdf", cell_size), width = 10, height = 5)
i=1
data <- c_grid_pvax[lt_label_1[i]] %>% filter(get(lt_label_1[[i]])>0)
p1 <- tm_shape(c_county) + tm_polygons(alpha = 0) +
tm_shape(data) +
tm_polygons(col = lt_label_1[i], breaks = brks, palette = cols, border.alpha = 0) +
tm_shape(c_county) + tm_polygons(alpha = 0) +
tm_shape(c_urban) + tm_polygons(alpha = 0) +
tm_layout(main.title =lt_label_1[[i]], legend.show = F, frame = F, title.size = 0.5)
i=2
data <- c_grid_pvax[lt_label_1[i]] %>% filter(get(lt_label_1[[i]])>0)
p2 <- tm_shape(c_county) + tm_polygons(alpha = 0) +
tm_shape(data) +
tm_polygons(col = lt_label_1[i], breaks = brks, palette = cols, border.alpha = 0) +
tm_shape(c_county) + tm_polygons(alpha = 0) +
tm_shape(c_urban) + tm_polygons(alpha = 0) +
tm_layout(main.title =lt_label_1[[i]], legend.show = F, frame = F, title.size = 0.5)
i=3
data <- c_grid_pvax[lt_label_1[i]] %>% filter(get(lt_label_1[[i]])>0)
p3 <- tm_shape(c_county) + tm_polygons(alpha = 0) +
tm_shape(data) +
tm_polygons(col = lt_label_1[i], breaks = brks, palette = cols, border.alpha = 0) +
tm_shape(c_county) + tm_polygons(alpha = 0) +
tm_shape(c_urban) + tm_polygons(alpha = 0) +
tm_layout(main.title =lt_label_1[[i]], legend.show = F, frame = F, title.size = 0.5)
i=4
data <- c_grid_pvax[lt_label_1[i]] %>% filter(get(lt_label_1[[i]])>0)
p4 <- tm_shape(c_county) + tm_polygons(alpha = 0) +
tm_shape(data) +
tm_polygons(col = lt_label_1[i], breaks = brks, palette = cols, border.alpha = 0) +
tm_shape(c_county) + tm_polygons(alpha = 0) +
tm_shape(c_urban) + tm_polygons(alpha = 0) +
tm_layout(main.title =lt_label_1[[i]], legend.show = F, frame = F, title.size = 0.5)
i=5 ## legend
data <- c_grid_pvax[lt_label_1[i]] %>% filter(get(lt_label_1[[i]])>0)
p5 <- tm_shape(c_county) + tm_polygons(alpha = 0) +
tm_shape(data) +
tm_polygons(col = lt_label_1[i], breaks = brks, palette = cols, border.alpha = 0) +
tm_shape(c_county) + tm_polygons(alpha = 0) +
tm_shape(c_urban) + tm_polygons(alpha = 0) +
tm_layout(main.title =lt_label_1[[i]], legend.show = T, frame = F, title.size = 0.5, legend.outside = T)
tmap_arrange(p1, p2, p3, p4, p5, widths = c(0.5,0.5), nrow = 1)
dev.off()
# prepare averaging data
data <- abm_df_vax[["step_73"]] %>%
st_as_sf(., coords = c("long","lat"), crs=4326) %>%
st_transform(., st_crs(c_county))
df <- c_bgroup
# all individual/agent (n=127584)
df["n_pop"] <- lengths(st_intersects(df, data))
lt_label <- sprintf("step73_r%01d",1:5)
lt_label_1 <- sprintf("step73_nvax_r%01d",1:5)
for (i in 1:length(lt_label)) {
tt <- data %>% filter(get(lt_label[i])==1)
df[as.character(lt_label_1[i])] <- lengths(st_intersects(df, tt))
}
ttt <- df[c("GEOID",lt_label_1)]
st_geometry(ttt) <- NULL
ttt["step73_nvax_mean"] <- round(apply(ttt[lt_label_1], 1, mean))
df <- left_join(df, ttt[c("GEOID","step73_nvax_mean")], by=c("GEOID"="GEOID"))
df["step73_pvax_mean"] <- round(100 * df$step73_nvax_mean / df$n_pop, 2)
df <- na.omit(df)
pdf("plot/07_cbgrou_pred_vaxrate_a_100_b_113_c_311_d_131.pdf")
tm_shape(c_county) + tm_polygons(col=NA, alpha = 0) +
tm_shape(df) + tm_polygons(col = "step73_pvax_mean", palette = "GnBu", alpha = 0.7,
border.alpha = 0, lwd = 0.1) +
tm_shape(c_urban) + tm_polygons(alpha = 0, border.col = "black", lwd = 0.1) +
tm_shape(c_county) + tm_polygons(alpha = 0, border.col = "black", lwd = 0.5) +
tm_layout(legend.show = T, frame = FALSE, legend.outside = T) +
tm_compass(position = c("left", "top"), north = 0) + tm_scale_bar(position = c("left", "top"), width = 0.15)
dev.off()
## test
data <- df[df$step73_pvax_mean>69, ]
tm_shape(data) + tm_polygons(col = "step73_pvax_mean", palette = "GnBu", alpha = 0.7,
border.alpha = 0.1, lwd = 0.1)
#"GnBu"