---
title: "R (3/3) ABM Result Analysis: Mapping Vaccinated Agents"
output: html_notebook
---

# Data 
## Agents' vaccination status for selected time steps
FINAL = Simulation results 
Weight: a_100_b_113_c_311_d_131
```{r}
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
```

# Analysis
## MAP1: process map, vax rate, Vaccinated ind. point 
### Grid preparation
```{r}
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
```{r}
## 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()


```


## MAP2: Census block group 
```{r}
# 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"
```

