Merge pull request #5 from mansueto-institute/updates-prepub4

updates to repo

Author: nmarchio

complexity-analysis.R

@@ -11,9 +11,7 @@ library(readxl)
 library(readr)
 library(writexl)
 library(units)
-library(ggsn)
 library(arrow)
-#library(geoarrow)
 library(sfarrow)
 library(osmdata)
 library(viridis)
@@ -23,6 +21,9 @@ library(DescTools)
 options(scipen = 9999)
 options(lifecycle_verbosity = "warning")
 
+library(ggsn)
+# devtools::install_github('mansueto-institute/ggsn') # fork of 'oswaldosantos/ggsn'
+
 # Load aggregation function
 setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
 getwd()
@@ -134,17 +135,6 @@ write_xlsx(list('dictionary_labels' = file_labels, 'dictionary_levels' = file_gr
 
 # Aggregation arguments ---------------------------------------------------
 
-log_10 <- function(x) {
-  x = replace_na(na_if(na_if(na_if(log10(x), NaN), Inf), -Inf), 0)
-  x = ifelse(x < 1, 1, x)
-  return(x)
-}
-
-share <- function(x) {
-  x = replace_na(na_if(na_if(na_if(x / sum(x), NaN), Inf), -Inf), 0)
-  return(x)
-}
-
 sum_cols = c("k_complexity_weighted_landscan_un", "k_complexity_weighted_worldpop_un", 
              "block_count", "block_area_m2", "block_hectares", "block_area_km2", "block_perimeter_meters", "building_area_m2", 
              "building_count", "parcel_count", "landscan_population", "landscan_population_un", "worldpop_population", "worldpop_population_un", 
@@ -1045,12 +1035,14 @@ if (!file.exists(paste0("data/africa_geodata.parquet"))) {
   print('africa_geodata.parquet already downloaded.')
 }
 
-area_data <- arrow::open_dataset(paste0('data/africa_geodata.parquet')) %>% 
+#area_data <- arrow::open_dataset(paste0('data/africa_geodata.parquet')) %>% 
+area_data <- arrow::open_dataset('data/lusaka_blocks.parquet') %>%
   filter(urban_id %in% c('ghsl_3798','periurban_925')) %>% 
   read_sf_dataset() %>%
   st_set_crs(4326) %>%
   st_make_valid()
 
+
 # Zoom map ----------------------------------------------------------------
 
 sites <- data.frame(place=c('nga','sle','zmb'),
@@ -1109,8 +1101,8 @@ area_data_framed <- area_data %>%
     theme_void() + inset_element(zoom_inset, left = -.5, bottom = .12, right = 1, top = .55, align_to = 'full')
   )
 
-df_combined_prep %>% group_by(country_code, country_name) %>%
-  tally() %>% print(n = 100)
+#df_combined_prep %>% group_by(country_code, country_name) %>%
+#  tally() %>% print(n = 100)
 
 # buildings <- arrow::open_dataset('buildings_polygons_ZMB.parquet')) %>%
 #   filter(gadm_code %in% unique(area_data_framed$gadm_code))%>% 
@@ -1162,7 +1154,7 @@ colorhexes <- colorRampPalette(c('#93328E','#CF3F80','#F7626B','#FF925A','#FFC55
             fill = 'white', alpha = 0, linewidth = 1.7) +
     geom_sf(data = sites %>% filter(place == place_var), color = '#434343', 
             fill = '#434343', alpha = 0, linewidth = .5) +
-    scale_fill_manual(expand = c(0,0), values = c(grey2, colorhexes), name = 'Block complexity') + 
+    scale_fill_manual(expand = c(0,0), values = c(grey2, colorhexes), name = 'Block\ncomplexity') + 
     theme_void() + theme(text = element_text(color = "#333333"),
                          legend.position = 'none',
                          legend.spacing.x = unit(1, 'pt'),
@@ -1176,10 +1168,10 @@ colorhexes <- colorRampPalette(c('#93328E','#CF3F80','#F7626B','#FF925A','#FFC55
                          legend.title = element_blank(),
                          plot.caption = element_text(size = 11, hjust = .5, vjust = 5), #margin=margin(0,0,0,0)),
                          axis.text = element_blank()) +
-    ggsn::scalebar(y.min = st_bbox(area_data_framed)$ymin - (height_decdegs*.03), 
-                   x.min = st_bbox(area_data_framed)$xmin, 
-                   y.max = st_bbox(area_data_framed)$ymax, 
-                   x.max = st_bbox(area_data_framed)$xmax, 
+    ggsn::scalebar(y.min = st_bbox(area_data_framed)$ymin - (height_decdegs*.03),
+                   x.min = st_bbox(area_data_framed)$xmin,
+                   y.max = st_bbox(area_data_framed)$ymax,
+                   x.max = st_bbox(area_data_framed)$xmax,
                    location = 'bottomleft',
                    height = .01, box.fill = c('#333333','#ffffff'),
                    border.size = .4, st.color = '#333333', st.size = 2.5, box.color = '#333333',
@@ -1193,7 +1185,7 @@ colorhexes <- colorRampPalette(c('#93328E','#CF3F80','#F7626B','#FF925A','#FFC55
     geom_sf(data = streets, color = 'white', fill = 'white', linewidth = 1) +
     geom_sf(data = sites %>% filter(place == place_var), color = 'white', fill = 'white', alpha = 0, linewidth = 1.5) +
     scale_fill_manual(expand = c(0,0), breaks =  c('1', '2', '3', '4', '5', '6', '7', '8', '9', '10+', 'Off\nnetwork'), 
-                      values = c(grey2, colorhexes), name = 'Block complexity') + 
+                      values = c(grey2, colorhexes), name = 'Block\ncomplexity') + 
     theme_void() +
     theme(legend.position = 'none'))
 
@@ -1216,7 +1208,7 @@ colorhexes <- colorRampPalette(c('#93328E','#CF3F80','#F7626B','#FF925A','#FFC55
     geom_sf(data = layers %>% filter(block_property == 'on-network-streets',
                                      block_id %in% c("ZMB.5.6_2_9371")), 
             fill = 'yellow', color = 'white', alpha = .8, linewidth = .7, lineend = 'round') + # linetype = "dashed", 
-    scale_fill_manual(name = 'Block complexity', values = c(grey2, colorhexes), breaks = c('1', '2', '3', '4', '5', '6', '7', '8', '9', '10+', 'Off\nnetwork')) +
+    scale_fill_manual(name = 'Block\ncomplexity', values = c(grey2, colorhexes), breaks = c('1', '2', '3', '4', '5', '6', '7', '8', '9', '10+', 'Off\nnetwork')) +
     labs(subtitle = '') +
     theme_void() + theme(legend.title = element_text(size = 14, face = 'bold'),
                          legend.text = element_text(size = 14),
@@ -1237,7 +1229,7 @@ colorhexes <- colorRampPalette(c('#93328E','#CF3F80','#F7626B','#FF925A','#FFC55
     geom_sf(data = layers %>% filter(block_property == 'on-network-streets',
                                      block_id %in% c("ZMB.5.6_2_9339", "ZMB.5.6_2_9334", "ZMB.5.6_2_9335")), 
             fill = 'yellow', color = 'white', alpha = .8, linewidth = .7, lineend = 'round') + # linetype = "dashed", 
-    scale_fill_manual(name = 'Block complexity', values = c(grey2, colorhexes), breaks = c('1', '2', '3', '4', '5', '6', '7', '8', '9', '10+', 'Off\nnetwork')) +
+    scale_fill_manual(name = 'Block\ncomplexity', values = c(grey2, colorhexes), breaks = c('1', '2', '3', '4', '5', '6', '7', '8', '9', '10+', 'Off\nnetwork')) +
     labs(subtitle = '') +
     theme_void() + theme(legend.position = 'none',
                          plot.subtitle = element_text(size = 10),
@@ -1246,11 +1238,16 @@ colorhexes <- colorRampPalette(c('#93328E','#CF3F80','#F7626B','#FF925A','#FFC55
 design <- "
   AAAABBBBBB
   AAAABBBBBB
-  CCCCDDDDEE
-  CCCCDDDDEE
+  AAAABBBBBB
+  AAAABBBBBB
+  AAAABBBBBB
+  CCCDDDEEFF
+  CCCDDDEEFF
+  CCCDDDEEFF
+  CCCDDDEEFF
 "
 
-(layers_viz <- africa_zoom +  plot_k_discrete + map_block + l1 + l2 + 
+(layers_viz <- africa_zoom +  plot_k_discrete + l1 + map_block + l2 + guide_area() +
     #plot_layout(design = design, guides = "collect", tag_level = 'new') + plot_annotation( tag_levels = c('A')) &
     plot_layout(design = design, guides = "collect") + plot_annotation(tag_levels = list(c("A", "B", "C", "D", "E", "F"))) & 
   theme(plot.tag = element_text(size = 12)))
@@ -1393,7 +1390,7 @@ colorhexes <- colorRampPalette(c('#93328E','#CF3F80','#F7626B','#FF925A','#FFC55
 
 (plot_k_discrete <- ggplot() +
     geom_sf(data = area_data, aes(fill = as.factor(k_labels)), color = '#ffffff', linewidth = .0075) +   
-    scale_fill_manual(expand = c(0,0), values = c(grey2, colorhexes), name = 'Block complexity') + 
+    scale_fill_manual(expand = c(0,0), values = c(grey2, colorhexes), name = 'Block\ncomplexity') + 
     labs(caption = paste0('Population-weighted average block complexity: ', area_data %>% st_drop_geometry() %>% 
                              summarise(wm_var = weighted.mean(as.integer(k_complexity), landscan_population_un)) %>% pull() %>% round(.,2))) +
     #guides(color = guide_legend(nrow = 1, label.position = "bottom", keywidth = 2, keyheight = 1), 

dhs-analysis.R

@@ -6,6 +6,7 @@ library(scales)
 library(sf)
 library(rmapshaper)
 library(patchwork)
+library(ggrepel)
 library(arrow)
 library(sfarrow)
 library(readxl)
@@ -39,7 +40,7 @@ load_saved = TRUE
 
 # RData file
 if (load_saved == TRUE) {
-  load("data/dhs_data.RData")
+  load(paste0(wd_path,"/data/dhs_data.RData"))
   # Requires authorization from DHS: https://uchicago.box.com/s/anw4fxc376dgtgol9tt87tuozipnk0kj
 
 } else {
@@ -55,10 +56,10 @@ if (load_saved == TRUE) {
 # Load staging files
 # Box.com link with access to staging files:
 # https://uchicago.box.com/s/anw4fxc376dgtgol9tt87tuozipnk0kj
-# subnat_indicators <- read_csv(paste0(wd_path,'/dhs_download.csv'))
-# subnat_all <- st_read(paste0(wd_path,'/dhs_geographies.geojson'))
-# subnat_all_wide <- st_read(paste0(wd_path,'/dhs_all_wide.geojson'))
-# subnat_to_blocks <- read_parquet(paste0(wd_path,'/blocks_to_dhs.parquet'))
+# subnat_indicators <- read_csv(paste0(wd_path,'/data/dhs_download.csv'))
+# subnat_all <- st_read(paste0(wd_path,'/data/dhs_geographies.geojson'))
+# subnat_all_wide <- st_read(paste0(wd_path,'/data/dhs_all_wide.geojson'))
+subnat_to_blocks <- read_parquet(paste0(wd_path,'/data/blocks_to_dhs.parquet'))
 
 # CHANGE TO DHS USER LOGIN CREDENTIAL
 # Follow instructions here: # https://docs.ropensci.org/rdhs/articles/introduction.html
@@ -258,7 +259,7 @@ if (load_saved == FALSE) {
       #   st_make_valid() %>% mutate(is_valid = st_is_valid(geometry)) %>%
       #   filter(is_valid == TRUE) %>% select(-any_of('is_valid'))
 
-      iso_blocks <- arrow::open_dataset(paste0(wd_input,'/africa_geodata.parquet')) %>% 
+      iso_blocks <- arrow::open_dataset(paste0('data/africa_geodata.parquet')) %>% 
         select(block_id, gadm_code, country_code, geometry) %>%
         filter(country_code %in% c(country_iso)) %>% 
         read_sf_dataset() %>%
@@ -412,7 +413,7 @@ if (load_saved == FALSE) {
   un_slums_k <- un_slums %>%
     left_join(., urban_k, by = c('country_code'='country_code'))
 
-  city_k <- read_parquet(paste0(wd_input,'/africa_data.parquet')) %>%
+  city_k <- read_parquet(paste0('data/africa_data.parquet')) %>%
     mutate(population_k_1 = case_when(k_labels == '1' ~ landscan_population_un, TRUE ~ as.numeric(0)),
            population_k_2 = case_when(k_labels == '2' ~ landscan_population_un, TRUE ~ as.numeric(0)),
            population_k_3 = case_when(k_labels == '3' ~ landscan_population_un, TRUE ~ as.numeric(0)),
@@ -1477,7 +1478,7 @@ k_country <- read_parquet(paste0('data/africa_data.parquet')) %>%
   ungroup() %>%
   mutate(k_complexity_wt = k_complexity_wt/landscan_population_un)
 
-street_validation <- read_csv('data/streets_external_validation.csv') %>%
+street_validation <- read_csv('data/streets_validation.csv') %>%
   mutate(osm_cia_ratio = osm_total_streets_km / cia_roadways_km,
          osm_ecopia_ratio = osm_total_streets_km / ecopia_ml_roads_km,
          osm_irf_ratio = osm_total_streets_km / irf_all_roads_km,
@@ -1764,7 +1765,7 @@ k_thresh_corr <- k_levels_vs_pc1 %>%
     geom_line(group = 1) + 
     scale_y_continuous(expand = c(0,0), limits = c(-.75,.7)) + 
     scale_size_continuous( labels = label_percent(), name = 'Share of population') +
-    labs(subtitle = 'Block complexity correlation with PC1 begins to slightly weaken\nwhen k = 15+ and for off-network blocks across DHS regions',
+    labs(#subtitle = 'Block complexity correlation with PC1 begins to slightly weaken\nwhen k = 15+ and for off-network blocks across DHS regions',
          y = 'Spearman correlation with PC1 at each block complexity level', x = 'Population share at each block complexity level correlated with PC1') + 
     theme_classic() +
     theme(plot.subtitle = element_text(hjust = .5),
@@ -1847,16 +1848,17 @@ k_15_ratios <- street_validation %>%
     #geom_text(x = 4, y = log10(15), label = 'OSM reports more vehicular roadway than CIA') +
     #scale_x_continuous(labels = scales::percent, breaks = c(0, .05, .1, .15, .2, .25, .3)) +
     scale_y_continuous(breaks = c( log10(.5), log10(1), log10(10), log10(100)), labels = c( .5, 1, 10, 100)) +
-    labs(subtitle = 'South Sudan and Eritrea have relatively high block complexity\nand report less roadway in OSM than the CIA World Factbook',
+    labs(#subtitle = 'South Sudan and Eritrea have relatively high block complexity\nand report less roadway in OSM than the CIA World Factbook',
          #x = "Ratio of population share in k = 15+ / off-network blocks to non-urban population share",
          x = "Average block complexity (weighted to population)",
          y = "OSM completeness\n(ratio of OSM vehicular street distance to CIA roadway distance)") +
     theme_classic() +
     theme(legend.position = 'none', 
           plot.subtitle = element_text(hjust = .5)))
 
+line_correl_k_thresh  + scatter_15plus_outliers
 (k_thresh_charts <- line_correl_k_thresh  + scatter_15plus_outliers & 
-  plot_layout(ncol = 2) &
+  #plot_layout(ncol = 2) &
   plot_annotation(tag_levels = list(c("A", "B"))) & 
   theme(plot.tag = element_text(size = 13)))