What if it was your city?

hiroshima.R

# load libraries 
library(tidyverse)
library(osmdata)
library(sf)
library(ggtext)
library(exactextractr)
library(fasterize)
library(ggtext)
library(tidygeocoder)
library(glue)
library(purrr)

epsg <- 32632
country <- "germany"

aoi_citys <- c("Freiburg im Breisgau, Germany",
               "München, Germany",
               "Berlin, Germany",
               "Erfurt, Germany",
               "Dresden, Germany",
               "Göttingen, Germany",
               "Kiel, Germany",
               "Hannover, Germany",
               "Essen, Germany",
               "Koblenz, Germany",
               "Greifswald, Germany",
               "Trier, Germany",
               "Kaiserslautern, Germany",
               "Nürnberg, Germany", 
               "Köln, Germany",
               "Regensburg, Germany",
               "Heidelberg, Germany",
               "Rothenburg ob der Tauber, Germany",
               "Bamberg, Germany",
               "Bremen, Germany",
               "Frankfurt, Germany",
               "Schwerin, Germany",
               "Potsdam, Germany",
               "Bonn, Germany",
               "Bacharach, Germany",
               "Tübingen, Germany",
               "Marburg, Germany",
               "Baden-Baden, Germany",
               "Stuttgart, Germany",
               "Sankt Goar, Germany",
               "Fritzlar, Germany",
               "Mittenwald, Germany",
               "Goslar, Germany",
               "Dortmund, Germany",
               "Düsseldorf, Germany",
               "Leipzig, Germany",
               "Duisburg, Germany",
               "Düsseldorf, Germany",
               "Dortmund, Germany",
               "Kassel, Germany",
               "Erfenbach, Germany",
               "Hamburg, Germany")

walk(aoi_citys, function(aoi_city){
  city <- tidygeocoder::geo(aoi_city, method = "google")  %>% 
    st_as_sf(coords = c("long", "lat"), crs = 4326) %>% 
    st_transform(32632)
  
  # streets data
  big_streets <-  getbb(aoi_city)%>%
    opq()%>%
    add_osm_feature(key = "highway", 
                    value = c("motorway", "primary", "motorway_link", "primary_link")) %>%
    osmdata_sf()
  big_streets <- big_streets$osm_lines %>% 
    st_transform(32632)
  
  med_streets <-  getbb(aoi_city)%>%
    opq()%>%
    add_osm_feature(key = "highway", 
                    value = c("secondary", "tertiary", "secondary_link", "tertiary_link")) %>%
    osmdata_sf()
  med_streets <- med_streets$osm_lines %>% 
    st_transform(32632)
  
  small_streets <- getbb(aoi_city)%>%
    opq()%>%
    add_osm_feature(key = "highway", 
                    value = c("residential", "living_street",
                              "unclassified",
                              "service", "footway"
                    )) %>%
    osmdata_sf()
  small_streets <- small_streets$osm_lines %>% 
    st_transform(32632)
  
  river <- getbb(aoi_city)%>%
    opq()%>%
    add_osm_feature(key = "waterway", value = "river") %>%
    osmdata_sf()
  river <- river$osm_lines %>% 
    st_transform(32632)
  
  radius_thermal_radiation <- st_buffer(city, 2100)
  radius_air_blast <- st_buffer(city, 1720)
  radius_radiation <- st_buffer(city, 1310)
  radius_air_blast_inner <- st_buffer(city, 760)
  radius_fireball <- st_buffer(city, 200)
  crop_buffer <- st_buffer(city, 2100)
  
  river_cropped <- st_crop(river, radius_thermal_radiation)
  river_cropped <- st_intersection(river_cropped, crop_buffer)
  
  small_streets <- st_crop(small_streets, radius_thermal_radiation)
  small_streets_2 <- st_intersection(small_streets, radius_thermal_radiation)  
  med_streets <- st_crop(med_streets, radius_thermal_radiation)
  med_streets <- st_intersection(med_streets, radius_thermal_radiation)  
  big_streets <- st_crop(big_streets, radius_thermal_radiation)
  big_streets <- st_intersection(big_streets, radius_thermal_radiation)
  
  border <- radius_thermal_radiation %>% st_cast("LINESTRING")
  
  legend_fire <- "**FIREBALL RADIUS** <br><br><br> 0.2km (0.13km<sup>2</sup>)"
  legend_airi <- "**AIR BLAST RADIUS** <br>(inner)<br><br> 0.76km (1.82km<sup>2</sup>)"
  legend_radi <- "**RADIATION RADIUS** <br><br><br> 1.31km (5.42km<sup>2</sup>)"
  legend_airb <- "**AIR BLAST RADIUS** <br>(outer)<br><br> 1.72km (9.26km<sup>2</sup>)"
  legend_ther <- "**THERMAL** <br> **RADIATION RADIUS** <br><br> 2.21km (15.4km<sup>2</sup>)"
  
  ggplot(radius_thermal_radiation) +
    geom_sf(fill = "#F1EFED") +
    geom_sf(data = small_streets, color = "#171717") +
     geom_sf(data = med_streets, color = "#171717", size = 1.4) +
    # geom_sf(data = big_streets, color = "#171717", size = 1.9) +
    geom_sf(data = river_cropped, inherit.aes = FALSE, color = "steelblue", size = 1.5) +
    geom_sf(data = radius_air_blast, aes(fill = "4"), alpha = 0.4, show.legend = "polygon") +
    geom_sf(data = radius_radiation, aes(fill = "3"), alpha = 0.4, show.legend = "polygon") +
    geom_sf(data = radius_air_blast_inner, aes(fill = "2"), alpha = 0.4, show.legend = "polygon") +
    geom_sf(data = radius_fireball, aes(fill = "1"), alpha = 0.4, show.legend = "polygon") +
    geom_sf(data = border, color = "#F1EFED", linetype = "dotted", show.legend = F) +
    labs(title = "What if it was your City?",
         subtitle = glue("A hypothetical atomic bomb hit on <span style = 'margin-bottom: 500px; font-size:22pt'>**_{aoi_city}_**</span>. <br>
                          The simulated explosion is based on the bomb Fat Man hitting Nagasaki in 1945."),
         caption = "Data: OSM | Visualisation: Marco Sciaini <br><br>
                  Original Visualisation by Alberto Lucas López for the South China Morning Post") +
    scale_fill_manual("", 
                      guide = "legend",
                      values=c("1"="#ff3333", "2"="#ff794d","3"="#c6c6c6","4"="#ffffb1", "5" = "#F1EFED"),
                      labels = c(legend_fire, legend_airi, legend_radi, legend_airb, legend_ther)) +
    theme_void() +
    guides(fill = guide_legend(keyheight = unit(1, units = "mm"),  
                               keywidth = unit(5, units = "mm"),
                               direction = "horizontal",
                               nrow = 1,
                               ticks.colour = "white",
                               label.position = "bottom",
                               title.position = "bottom",
                               title.hjust = 0.5)) +
    theme(legend.direction="horizontal",
          legend.position="bottom", 
          legend.box = "vertical",
          text=element_text(family="West"),
          panel.spacing.y = unit(1, "lines"),
          plot.title.position = "plot",
          plot.title = element_markdown(size = 32, 
                                        face = "bold", 
                                        colour = '#F1EFED',
                                        hjust = 0,
                                        margin = margin(.3,0,.5,-2, "cm")),
          plot.subtitle = element_markdown(size = 20, 
                                           face = "bold", 
                                           colour = '#F1EFED',
                                           hjust = 0,
                                           margin = margin(.3,0,.5,-2, "cm"),
                                           lineheight = 1.3,),
          legend.text=element_markdown(color='#F1EFED',size=14),
          legend.key = element_rect(color = NULL, linetype = NULL),
          plot.caption = element_markdown(size = 12,
                                          colour = '#F1EFED',
                                          margin = margin(.3,-2,.3,0, "cm")),
          plot.margin = margin(0, 150, 0, 150),
          plot.background = element_rect(fill = "#171717"))
  
  ggsave(glue("plots/germany/{aoi_city}.png"), width = 12, height = 12, dpi = 600)        
})