How to overlay a choropleth map on a base map

R code
1

Describe your issue

  • I would like to “overlay” (not sure if that’s the proper term) my choropleth map (using ggplot2) of attack rates per barangay (village) on a basemap (using OpenStreetMap package). Perhaps adjusting the transparency of the choropleth map to show the underlying basemap.

What steps have you already taken to find an answer?

  • I followed the Epi R Handbook Chapters 28.7 to 28.9 in creating the choropleth map and basemap.
  • I was able to create both maps by following the instructions and modifying it a little bit, however I was unable to “overlay” one map on the other.

Provide an example of your R code

Attached is my code along with a link to the shapefile I used.

# load package

pacman::p_load(rio, here, janitor, tidyverse, reprex, sf, ggrepel, OpenStreetMap)

# create demo data frame

barangay_cases <- data.frame(
  stringsAsFactors = FALSE,
          barangay = c("Buluangan","San Juan",
                       "Barangay VI (Pob.)","Punao","Barangay I (Pob.)",
                       "Barangay IV (Pob.)","Ermita","Barangay V (Pob.)","Palampas",
                       "Rizal","Nataban","Guadalupe","Barangay III (Pob.)",
                       "Prosperidad","Quezon","Barangay II (Pob.)","Codcod",
                       "Bagonbon"),
               pop = c(11168,2714,5465,6197,10821,
                       449,2193,7318,9519,17087,3886,10947,3269,5878,
                       10798,6607,14507,5886),
                 n = c(50,10,20,5,30,5,0,40,10,
                       20,0,0,0,0,0,0,0,0),
       attack_rate = c(4.4,
                       3.6,3.6,0.8,
                       2.7,11.1,0,5.4,
                       1.0,1.1,0,0,0,0,0,0,0,0)
)

# import shapefile
barangay_shp <- read_sf("san_carlos.shp") %>% 
  clean_names() %>% 
  select(-c(region, reg_psgc9, reg_psgc10, province, prov_psgc9, prov_psgc1,
            municity, muni_psgc9, muni_psgc1, brg_psgc9, brg_psgc10))

# join barangay_cases with tmap
barangay_map <- barangay_cases %>% 
  left_join(barangay_shp, by="barangay") %>%
  mutate(barangay = replace(barangay, barangay == "Barangay I (Pob.)", "I"),
         barangay = replace(barangay, barangay == "Barangay II (Pob.)", "II"),
         barangay = replace(barangay, barangay == "Barangay III (Pob.)", "III"),
         barangay = replace(barangay, barangay == "Barangay IV (Pob.)", "IV"),
         barangay = replace(barangay, barangay == "Barangay V (Pob.)", "V"),
         barangay = replace(barangay, barangay == "Barangay VI (Pob.)", "VI")) %>% 
  st_as_sf()

### using ggplot --------------------------------------------------

ggplot(data = barangay_map) +                           
  geom_sf(aes(fill = attack_rate)) +                        
  scale_fill_continuous(high="#54278f", low="#f2f0f7") +    # change color gradient
  theme_void() +
  labs(fill = "Attack Rate \nper 1,000 pop.") +
  # use geom_label_repel() from {ggrepel} to prevent overlap in barangay
  # labels instead of geom_sf_label(aes(label = barangay))
  geom_label_repel(aes(label = barangay, geometry = geometry), 
                   stat = "sf_coordinates")

### add basemap -------------------------------------------------------------

# Fit basemap by range of lat/long coordinates.
map <- OpenStreetMap::openmap(
  upperLeft = c(10.580948, 123.133346),   # limits of basemap tile
  lowerRight = c(10.383456, 123.461219),
  zoom = NULL,
  type = "apple-iphoto")

# Projection WGS84
map_latlon <- openproj(map, projection = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")

# Plot map. Must use "autoplot" in order to work with ggplot
autoplot.OpenStreetMap(map_latlon)

Created on 2023-04-17 with reprex v2.0.2

Session info 
sessionInfo()
#> R version 4.2.2 (2022-10-31 ucrt)
#> Platform: x86_64-w64-mingw32/x64 (64-bit)
#> Running under: Windows 10 x64 (build 22621)
#> 
#> Matrix products: default
#> 
#> locale:
#> [1] LC_COLLATE=English_Philippines.utf8  LC_CTYPE=English_Philippines.utf8   
#> [3] LC_MONETARY=English_Philippines.utf8 LC_NUMERIC=C                        
#> [5] LC_TIME=English_Philippines.utf8    
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#>  [1] OpenStreetMap_0.3.4 ggrepel_0.9.3       sf_1.0-9           
#>  [4] reprex_2.0.2        lubridate_1.9.2     forcats_1.0.0      
#>  [7] stringr_1.5.0       dplyr_1.1.0         purrr_1.0.1        
#> [10] readr_2.1.4         tidyr_1.3.0         tibble_3.2.0       
#> [13] ggplot2_3.4.1       tidyverse_2.0.0     janitor_2.2.0      
#> [16] here_1.0.1          rio_0.5.29         
#> 
#> loaded via a namespace (and not attached):
#>  [1] httr_1.4.5         sp_1.6-0           highr_0.10         cellranger_1.1.0  
#>  [5] yaml_2.3.7         pillar_1.8.1       lattice_0.20-45    glue_1.6.2        
#>  [9] digest_0.6.31      snakecase_0.11.0   colorspace_2.1-0   htmltools_0.5.4   
#> [13] pkgconfig_2.0.3    raster_3.6-20      haven_2.5.2        scales_1.2.1      
#> [17] terra_1.7-18       openxlsx_4.2.5.2   tzdb_0.3.0         timechange_0.2.0  
#> [21] proxy_0.4-27       generics_0.1.3     ellipsis_0.3.2     pacman_0.5.1      
#> [25] withr_2.5.0        cli_3.6.0          magrittr_2.0.3     readxl_1.4.2      
#> [29] mime_0.12          evaluate_0.20      fs_1.6.1           fansi_1.0.4       
#> [33] xml2_1.3.3         foreign_0.8-83     class_7.3-20       tools_4.2.2       
#> [37] data.table_1.14.8  hms_1.1.2          lifecycle_1.0.3    munsell_0.5.0     
#> [41] zip_2.2.2          compiler_4.2.2     e1071_1.7-13       rlang_1.0.6       
#> [45] classInt_0.4-9     units_0.8-1        grid_4.2.2         rstudioapi_0.14   
#> [49] labeling_0.4.2     rmarkdown_2.20     gtable_0.3.1       codetools_0.2-18  
#> [53] DBI_1.1.3          curl_5.0.0         R6_2.5.1           knitr_1.42        
#> [57] rgdal_1.6-5        fastmap_1.1.1      utf8_1.2.3         rprojroot_2.0.3   
#> [61] KernSmooth_2.23-20 rJava_1.0-6        stringi_1.7.12     Rcpp_1.0.10       
#> [65] vctrs_0.5.2        tidyselect_1.2.0   xfun_0.37

shapefile for reference: Error

2 Likes
2

Hi Ian,
@ChrisJ @andysouth @aspina @avaughan or others may be able to help
Neale

1 Like
3

Hi Ian,

You are 99% of the way there, there is just a small addition to your plotting code needed.

The way to use geom_sf() ontop of the OpenStreetMap object is to specify inherit.aes = FALSE in the function. You also need to explicitly state the data too.


# Plot map. Must use "autoplot" in order to work with ggplot
autoplot.OpenStreetMap(map_latlon) +
  geom_sf(data = barangay_map,                                   #Note that I have specified the data here
          mapping = aes(fill = attack_rate), 
          inherit.aes = FALSE) +                                 #inherit.aes = FALSE              
  scale_fill_continuous(high="#54278f", low="#f2f0f7") +    
  theme_void() +
  labs(fill = "Attack Rate \nper 1,000 pop.") +
  geom_label_repel(data = barangay_map,                          #Note that I have specified the data here
                   mapping = aes(label = barangay, geometry = geometry),   
                   stat = "sf_coordinates",
                   inherit.aes = FALSE)                          #inherit.aes = FALSE

osm_sf_map
osm_sf_map1920×1366 119 KB

3 Likes
4

Thank you, @arranh! That worked perfectly. Is there a way to adjust the transparency of the choropleth so I can still see some of the underlying terrain in the basemap?

Regards,
Ian

1 Like
5

Great!

Yes to control transparency you use the argument alpha = , where alpha is a transparency scale from completely transparent (0) to completely opaque (1). This works with pretty much every ggplot() command

# Plot map. Must use "autoplot" in order to work with ggplot
autoplot.OpenStreetMap(map_latlon) +
  geom_sf(data = barangay_map,                #Note that I have specified the data here
          mapping = aes(fill = attack_rate), 
          inherit.aes = FALSE,                           #inherit.aes = FALSE
          alpha = 0.5) +                                      #transparency on a scale from completely transparent (0) to completely opaque (1)        
  scale_fill_continuous(high="#54278f", low="#f2f0f7") +    
  theme_void() +
  labs(fill = "Attack Rate \nper 1,000 pop.") +
  geom_label_repel(data = barangay_map,                          #Note that I have specified the data here
                   mapping = aes(label = barangay, geometry = geometry),   
                   stat = "sf_coordinates",
                   inherit.aes = FALSE)                          #inherit.aes = FALSE
3 Likes
6

I totally forgot about the alpha = argument. Thank you so much for the quick response!

3 Likes