Mapping_230926
Heeyoung
2023-09-27
setwd("C:/Users/User/OneDrive/R FILE/Spatial Statistics/mapping with r")
library(tigris)## To enable caching of data, set `options(tigris_use_cache = TRUE)`
## in your R script or .Rprofile.library(tidycensus)
library(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.3 ✔ readr 2.1.4
## ✔ forcats 1.0.0 ✔ stringr 1.5.0
## ✔ ggplot2 3.4.3 ✔ tibble 3.2.1
## ✔ lubridate 1.9.2 ✔ tidyr 1.3.0
## ✔ purrr 1.0.2## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorslibrary(tmap)## The legacy packages maptools, rgdal, and rgeos, underpinning the sp package,
## which was just loaded, will retire in October 2023.
## Please refer to R-spatial evolution reports for details, especially
## https://r-spatial.org/r/2023/05/15/evolution4.html.
## It may be desirable to make the sf package available;
## package maintainers should consider adding sf to Suggests:.
## The sp package is now running under evolution status 2
## (status 2 uses the sf package in place of rgdal)
## Breaking News: tmap 3.x is retiring. Please test v4, e.g. with
## remotes::install_github('r-tmap/tmap')Set the dataset
First, let’s suppose that we have county level dataset
data <- readRDS("./sampledata.rds")
glimpse(data)## Rows: 3,220
## Columns: 64
## $ fips <chr> "01001", "01003", "01005", "01007", "01009", …
## $ county <chr> "Autauga County, Alabama", "Baldwin County, A…
## $ year <chr> "2019", "2019", "2019", "2019", "2019", "2019…
## $ m_pop <dbl> 26934, 103496, 13421, 12150, 28495, 5571, 908…
## $ f_pop <dbl> 28446, 109334, 11940, 10343, 29186, 4677, 107…
## $ m_b_pop <dbl> 4989, 9315, 6409, 3247, 535, 3969, 3987, 1094…
## $ f_b_pop <dbl> 5549, 10403, 5657, 1767, 393, 3697, 5014, 128…
## $ asian_pop <dbl> 573, 1969, 134, 27, 212, 50, 69, 1014, 376, 6…
## $ m_w_pop <dbl> 20220, 85703, 6100, 8536, 24576, 1409, 4872, …
## $ f_w_pop <dbl> 21095, 91100, 5518, 8227, 25545, 782, 5360, 4…
## $ m_h_pop <dbl> 856, 5115, 610, 270, 2882, 115, 103, 2317, 41…
## $ f_h_pop <dbl> 709, 4596, 495, 309, 2460, 147, 174, 2035, 38…
## $ med_age_blk <dbl> 35.4, 33.4, 35.5, 36.6, 35.6, 37.6, 35.8, 33.…
## $ med_age_asi <dbl> 34.6, 39.3, 40.4, NA, 27.9, 6.1, 27.8, 41.4, …
## $ med_age_whi <dbl> 40.5, 45.6, 48.5, 42.1, 43.4, 53.6, 48.0, 43.…
## $ med_age_his <dbl> 23.0, 28.5, 25.6, 28.2, 24.8, 40.1, 32.6, 23.…
## $ total_pop <dbl> 55380, 212830, 25361, 22493, 57681, 10248, 19…
## $ us_born <dbl> 54081, 204906, 24672, 22153, 55062, 10176, 19…
## $ fore_born <dbl> 1299, 7924, 689, 340, 2619, 72, 153, 3051, 61…
## $ fore_born_ci <dbl> 750, 3650, 335, 140, 966, 72, 78, 1523, 135, …
## $ fore_born_nonci <dbl> 118, 980, 25, 19, 82, 14, 58, 344, 31, 0, 46,…
## $ fore_born_origin_asia <dbl> 496, 1395, 114, 41, 153, 9, 63, 799, 305, 64,…
## $ fore_born_origin_e_asia <dbl> 187, 725, 28, 0, 72, 9, 30, 283, 213, 29, 9, …
## $ fore_born_origin_c_asia <dbl> 0, 395, 74, 41, 38, 0, 20, 166, 56, 0, 0, 5, …
## $ fore_born_origin_se_asia <dbl> 309, 254, 12, 0, 41, 0, 13, 350, 36, 35, 127,…
## $ fore_born_origin_w_asia <dbl> 0, 21, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ pop_over_m_25 <dbl> 17680, 72164, 9562, 8683, 19273, 3843, 6087, …
## $ pop_educ_m_ba <dbl> 3097, 15507, 717, 474, 1370, 266, 475, 3681, …
## $ pop_educ_m_ma <dbl> 1495, 4883, 177, 173, 544, 55, 165, 1509, 172…
## $ pop_educ_m_pr <dbl> 153, 1690, 89, 51, 119, 43, 86, 589, 121, 235…
## $ pop_educ_m_phd <dbl> 289, 821, 32, 47, 91, 42, 11, 444, 9, 62, 25,…
## $ pop_over_f_25 <dbl> 19687, 78948, 8402, 7485, 20518, 3258, 7732, …
## $ pop_educ_f_ba <dbl> 2922, 16294, 650, 569, 2062, 261, 868, 4630, …
## $ pop_educ_f_ma <dbl> 1380, 6929, 318, 315, 849, 157, 561, 2975, 46…
## $ pop_educ_f_pr <dbl> 346, 1222, 39, 29, 147, 26, 35, 398, 88, 134,…
## $ pop_educ_f_phd <dbl> 247, 802, 58, 20, 28, 6, 29, 394, 51, 45, 65,…
## $ pov_pop <dbl> 8340, 21704, 6875, 3740, 7739, 2825, 4397, 19…
## $ pov_b_pop <dbl> 3324, 4684, 4579, 1127, 72, 2711, 2931, 7171,…
## $ pov_a_pop <dbl> 100, 155, 0, 0, 15, 9, 20, 174, 51, 25, 0, 0,…
## $ pov_w_pop <dbl> 4558, 14680, 1551, 2613, 6126, 69, 1192, 1066…
## $ pov_h_pop <dbl> 68, 1756, 607, 0, 1359, 35, 191, 1116, 441, 5…
## $ hhincome <dbl> 58731, 58320, 32525, 47542, 49358, 37785, 406…
## $ hhincome_b <dbl> 28808, 36616, 22357, 21964, 77518, 28889, 295…
## $ hhincome_a <dbl> 2499, 47269, 52841, NA, 101071, NA, 92188, 37…
## $ hhincome_w <dbl> 65992, 61872, 47175, 52543, 49529, 65795, 510…
## $ hhincome_h <dbl> 86220, 41851, 30563, 46103, 47529, NA, NA, 36…
## $ pov_pop_rate <dbl> 0.1505959, 0.1019781, 0.2710855, 0.1662740, 0…
## $ pov_w_propor <dbl> 0.54652278, 0.67637302, 0.22560000, 0.6986631…
## $ pov_b_propor <dbl> 0.398561151, 0.215812753, 0.666036364, 0.3013…
## $ pov_h_propor <dbl> 0.008153477, 0.080906745, 0.088290909, 0.0000…
## $ pov_a_propor <dbl> 0.011990408, 0.007141541, 0.000000000, 0.0000…
## $ pov_b_to_w <dbl> 0.72926722, 0.31907357, 2.95228885, 0.4313050…
## $ pov_h_to_w <dbl> 0.01491882, 0.11961853, 0.39136041, 0.0000000…
## $ hhincome_b_to_w <dbl> 0.4365378, 0.5918024, 0.4739163, 0.4180195, 1…
## $ hhincome_h_to_w <dbl> 1.3065220, 0.6764126, 0.6478643, 0.8774337, 0…
## $ educ_over_ba_m <dbl> 5034, 22901, 1015, 745, 2124, 406, 737, 6223,…
## $ educ_over_ba_f <dbl> 4895, 25247, 1065, 933, 3086, 450, 1493, 8397…
## $ educ_m_to_f <dbl> 0.9723878, 1.1024409, 1.0492611, 1.2523490, 1…
## $ asia_to_total <dbl> 0.0103466956, 0.0092515153, 0.0052837033, 0.0…
## $ fore_asia_to_fore_born <dbl> 0.38183218, 0.17604745, 0.16545718, 0.1205882…
## $ fore_e_asia_to_fore_born <dbl> 0.143956890, 0.091494195, 0.040638607, 0.0000…
## $ fore_born_nonci_to_ci <dbl> 0.090839107, 0.123674912, 0.036284470, 0.0558…
## $ asian_ratio <dbl> 0.0103466956, 0.0092515153, 0.0052837033, 0.0…
## $ fore_born_nonci_to_total <dbl> 0.0021307331, 0.0046046140, 0.0009857655, 0.0…Â
Next, let’s make a spatial data using tigris()
# state
us_states <- tigris::states(cb = TRUE, resolution = "20m") %>%
shift_geometry()## Retrieving data for the year 2021##
|
| | 0%
|
|====== | 8%
|
|============ | 17%
|
|============= | 19%
|
|======================== | 34%
|
|============================== | 43%
|
|======================================= | 56%
|
|================================================ | 69%
|
|====================================================== | 77%
|
|=============================================================== | 90%
|
|======================================================================| 100%# county
us_counties <- tigris::counties(cb = TRUE, resolution = "20m") %>%
shift_geometry()## Retrieving data for the year 2021##
|
| | 0%
|
|= | 2%
|
|=== | 4%
|
|===== | 7%
|
|======== | 11%
|
|========== | 14%
|
|============= | 19%
|
|=============== | 21%
|
|================== | 26%
|
|=================== | 27%
|
|==================== | 28%
|
|===================== | 30%
|
|======================= | 33%
|
|========================= | 36%
|
|============================ | 40%
|
|============================== | 43%
|
|=============================== | 44%
|
|================================= | 47%
|
|=================================== | 50%
|
|====================================== | 54%
|
|======================================== | 57%
|
|=========================================== | 61%
|
|============================================= | 64%
|
|============================================== | 66%
|
|================================================ | 69%
|
|================================================== | 71%
|
|===================================================== | 76%
|
|======================================================= | 78%
|
|========================================================== | 83%
|
|============================================================ | 85%
|
|=============================================================== | 90%
|
|================================================================= | 93%
|
|==================================================================== | 97%
|
|======================================================================| 100%Â
Finally, simply use left_join
# ORDER: spatial object must be on the left-hand side
combined <- us_counties %>% left_join(data,
by = c("GEOID"="fips"))Now, I’m ready!
Â
Â
ggplot2
combined %>%
ggplot(aes(fill = asia_to_total)) +
geom_sf() +
scale_fill_distiller(palette = "GrandBudapest",
direction = 1) +
labs(title = "Total Population, 2019",
caption = "Data source: 2019 5-year ACS, US Census Bureau",
fill = "ACS estimate") +
theme_void()## Warning in pal_name(palette, type): Unknown palette GrandBudapest
library(wesanderson)## Warning: package 'wesanderson' was built under R version 4.3.1combined <- combined %>%
mutate(asia_to_total2 = asia_to_total*100)
summary(combined$asia_to_total2)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 0.0000 0.2708 0.6132 1.3546 1.2896 42.6568 2combined <- mutate(combined, asia_to_total3 = cut(asia_to_total2,
breaks = c(-Inf, 5, 10, 20, 30, Inf),
labels = c("0-5", "5-10", "10-20", "20-30", "30+")))
combined %>%
filter(!is.na(asia_to_total3)) %>%
ggplot(aes(fill = asia_to_total3)) +
geom_sf() +
scale_fill_brewer(palette = "PuRd", direction = 1) +
labs(title = "Total Population, 2019",
caption = "Data source: 2019 5-year ACS, US Census Bureau",
fill = "% of Asian Population") +
theme_void() +
theme(plot.title = element_text(face = "bold", size = 20, hjust = 0.5),
plot.subtitle = element_text(size = 14, hjust = 0.5),
plot.caption = element_text(size = 10, hjust = 1),
legend.title = element_text(face = "bold", size = 12),
legend.text = element_text(size = 10))
tmap
tm_shape(combined) +
tm_polygons()
tm_shape(combined) +
tm_polygons(col = "asia_to_total2")
tm_shape(combined) +
tm_polygons(col = "asia_to_total2",
style = "quantile",
n = 4,
palette = "Purples",
title = "Asian Population")
# style = "jenks"
# palette = "magma"
tm_shape(combined) +
tm_polygons(col = "asia_to_total2",
style = "jenks",
n = 4,
palette = "Reds",
title = "Asian Population")
combined %>%
filter(!is.na(asia_to_total2)) %>%
tm_shape() +
tm_polygons(col = "asia_to_total2",
style = "jenks",
n = 4,
palette = "Reds",
title = "Percentage (%)",
legend.hist = F,
border.col = "grey",
border.alpha = 0.7) +
tm_layout(title = "",
frame = F,
legend.outside = F,
legend.outside.position = "left",
legend.position = c("left", "top"))
texas <- combined %>% filter(STATE_NAME=="Texas")
tm_shape(texas) +
tm_polygons() +
tm_bubbles(size = "asia_to_total2",
alpha = 0.5,
col = "navy",
title.size = "% of Asian") +
tm_layout(title = "Asian Population \n by County",
frame = T,
legend.outside = TRUE,
legend.outside.position = "bottom",
bg.color = "white") 
texas <- combined %>% filter(STATE_NAME=="Texas")
california <- combined %>% filter(STATE_NAME=="California")
newyork <- combined %>% filter(STATE_NAME=="New York")
washington <- combined %>% filter(STATE_NAME=="Washington")
t <- tm_shape(texas) +
tm_polygons() +
tm_bubbles(size = "asia_to_total2",
alpha = 0.5,
col = "darkorange2",
title.size = "Percentage (%)") +
tm_layout(title = "Texas",
frame = F,
legend.outside = TRUE,
legend.outside.position = "left",
legend.width = 1.5,
bg.color = "white")
c <- tm_shape(california) +
tm_polygons() +
tm_bubbles(size = "asia_to_total2",
alpha = 0.5,
col = "darkorange2",
title.size = "Percentage (%)") +
tm_layout(title = "California",
frame = F,
legend.outside = TRUE,
legend.outside.position = "left",
legend.width = 1.5,
bg.color = "white")
n <- tm_shape(newyork) +
tm_polygons() +
tm_bubbles(size = "asia_to_total2",
alpha = 0.5,
col = "darkorange2",
title.size = "Percentage (%)") +
tm_layout(title = "New York",
frame = F,
legend.outside = TRUE,
legend.outside.position = "left",
legend.width = 1.5,
bg.color = "white")
w <- tm_shape(washington) +
tm_polygons() +
tm_bubbles(size = "asia_to_total2",
alpha = 0.5,
col = "darkorange2",
title.size = "Percentage (%)") +
tm_layout(title = "Washington",
frame = F,
legend.outside = TRUE,
legend.outside.position = "left",
legend.width = 1.5,
bg.color = "white")
tmap_arrange(w, n, c, t,
widths = c(0.5,0.5))## Warning in pre_process_gt(x, interactive = interactive, orig_crs =
## gm$shape.orig_crs): legend.width controls the width of the legend within a map.
## Please use legend.outside.size to control the width of the outside legend
## Warning in pre_process_gt(x, interactive = interactive, orig_crs =
## gm$shape.orig_crs): legend.width controls the width of the legend within a map.
## Please use legend.outside.size to control the width of the outside legend
## Warning in pre_process_gt(x, interactive = interactive, orig_crs =
## gm$shape.orig_crs): legend.width controls the width of the legend within a map.
## Please use legend.outside.size to control the width of the outside legend
## Warning in pre_process_gt(x, interactive = interactive, orig_crs =
## gm$shape.orig_crs): legend.width controls the width of the legend within a map.
## Please use legend.outside.size to control the width of the outside legend
## Warning in pre_process_gt(x, interactive = interactive, orig_crs =
## gm$shape.orig_crs): legend.width controls the width of the legend within a map.
## Please use legend.outside.size to control the width of the outside legend
## Warning in pre_process_gt(x, interactive = interactive, orig_crs =
## gm$shape.orig_crs): legend.width controls the width of the legend within a map.
## Please use legend.outside.size to control the width of the outside legend
## Warning in pre_process_gt(x, interactive = interactive, orig_crs =
## gm$shape.orig_crs): legend.width controls the width of the legend within a map.
## Please use legend.outside.size to control the width of the outside legend
## Warning in pre_process_gt(x, interactive = interactive, orig_crs =
## gm$shape.orig_crs): legend.width controls the width of the legend within a map.
## Please use legend.outside.size to control the width of the outside legend## The legend is too narrow to place all symbol sizes.
combined <- combined %>% mutate(fore_born_nonci_to_total2 = fore_born_nonci_to_total*100)
combined %>%
filter(!is.na(fore_born_nonci_to_total2)) %>%
tm_shape() +
tm_polygons(col = "fore_born_nonci_to_total2",
style = "jenks",
n = 4,
palette = "Greens",
title = "Percentage (%)",
legend.hist = F,
border.col = "grey",
border.alpha = 0.7) +
tm_layout(title = "",
frame = F,
legend.outside = F,
legend.outside.position = "left",
legend.position = c("left", "top"))
texas <- combined %>% filter(STATE_NAME=="Texas")
california <- combined %>% filter(STATE_NAME=="California")
newyork <- combined %>% filter(STATE_NAME=="New York")
washington <- combined %>% filter(STATE_NAME=="Washington")
t <- tm_shape(texas) +
tm_polygons() +
tm_bubbles(size = "fore_born_nonci_to_total2",
alpha = 0.6,
col = "darkgreen",
title.size = "Percentage (%)") +
tm_layout(title = "Texas",
frame = F,
legend.outside = TRUE,
legend.outside.position = "left",
legend.width = 1.5,
bg.color = "white")
c <- tm_shape(california) +
tm_polygons() +
tm_bubbles(size = "fore_born_nonci_to_total2",
alpha = 0.6,
col = "darkgreen",
title.size = "Percentage (%)") +
tm_layout(title = "California",
frame = F,
legend.outside = TRUE,
legend.outside.position = "left",
legend.width = 1.5,
bg.color = "white")
n <- tm_shape(newyork) +
tm_polygons() +
tm_bubbles(size = "fore_born_nonci_to_total2",
alpha = 0.6,
col = "darkgreen",
title.size = "Percentage (%)") +
tm_layout(title = "New York",
frame = F,
legend.outside = TRUE,
legend.outside.position = "left",
legend.width = 1.5,
bg.color = "white")
w <- tm_shape(washington) +
tm_polygons() +
tm_bubbles(size = "fore_born_nonci_to_total2",
alpha = 0.6,
col = "darkgreen",
title.size = "Percentage (%)") +
tm_layout(title = "Washington",
frame = F,
legend.outside = TRUE,
legend.outside.position = "left",
legend.width = 1.5,
bg.color = "white")
tmap_arrange(w, n, c, t,
widths = c(0.5,0.5))## Warning in pre_process_gt(x, interactive = interactive, orig_crs =
## gm$shape.orig_crs): legend.width controls the width of the legend within a map.
## Please use legend.outside.size to control the width of the outside legend
## Warning in pre_process_gt(x, interactive = interactive, orig_crs =
## gm$shape.orig_crs): legend.width controls the width of the legend within a map.
## Please use legend.outside.size to control the width of the outside legend
## Warning in pre_process_gt(x, interactive = interactive, orig_crs =
## gm$shape.orig_crs): legend.width controls the width of the legend within a map.
## Please use legend.outside.size to control the width of the outside legend
## Warning in pre_process_gt(x, interactive = interactive, orig_crs =
## gm$shape.orig_crs): legend.width controls the width of the legend within a map.
## Please use legend.outside.size to control the width of the outside legend
## Warning in pre_process_gt(x, interactive = interactive, orig_crs =
## gm$shape.orig_crs): legend.width controls the width of the legend within a map.
## Please use legend.outside.size to control the width of the outside legend## Legend labels were too wide. Therefore, legend.text.size has been set to 0.62. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.## The legend is too narrow to place all symbol sizes.## Warning in pre_process_gt(x, interactive = interactive, orig_crs =
## gm$shape.orig_crs): legend.width controls the width of the legend within a map.
## Please use legend.outside.size to control the width of the outside legend
## Warning in pre_process_gt(x, interactive = interactive, orig_crs =
## gm$shape.orig_crs): legend.width controls the width of the legend within a map.
## Please use legend.outside.size to control the width of the outside legend
## Warning in pre_process_gt(x, interactive = interactive, orig_crs =
## gm$shape.orig_crs): legend.width controls the width of the legend within a map.
## Please use legend.outside.size to control the width of the outside legend
reference https://walker-data.com/census-r/mapping-census-data-with-r.html https://bookdown.org/nicohahn/making_maps_with_r5/docs/tmap.html https://thinking-spatial.org/courses/angewandte_geodatenverarbeitung/kurs04/ https://r-tmap.github.io/tmap-book/layout.html https://workshop.mhermans.net/thematic-maps-r/04_plot.html