/
09-usboundaries.R
59 lines (54 loc) · 2.28 KB
/
09-usboundaries.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
# Aim: create animation showing shifting US boundaries
# depends on 17 MB USAboundariesData package
# link to script file that shows changing state boundaries
# install.packages("USAboundaries", repos = "https://ropensci.r-universe.dev")
# install.packages("USAboundariesData", repos = "https://ropensci.r-universe.dev", type = "source")
library(USAboundaries)
library(tidyverse)
library(tmap)
library(sf)
dates = paste(historydata::us_state_populations$year, "01", "01", sep = "-")
dates_unique = unique(dates)
# select all dates earlier than 2000 after this error:
# Error in us_states(map_date, resolution, states) :
# map_date <= as.Date("2000-12-31") is not TRUE
dates_unique = dates_unique[dates_unique <= "2000-12-31"]
usb1 = USAboundaries::us_states(map_date = dates_unique[1])
usb1$year = lubridate::year(dates_unique[1])
plot(usb1$geometry)
usbl = map(dates_unique, ~USAboundaries::us_states(map_date = .))
# usb = do.call(rbind, usbl)
statepop = historydata::us_state_populations |>
select(-GISJOIN) |> rename(name = state)
sel = usb1$name %in% statepop$name
summary(sel)
usb1$name[!sel]
usbj = left_join(usb1, statepop)
plot(usbj["population"])
i = 2
dates_unique[dates_unique > "2000-12-31"] = "2000-12-31"
for(i in 2:length(dates_unique)) {
usbi = USAboundaries::us_states(map_date = dates_unique[i])
#print(st_crs(usbi))
usbi$year = lubridate::year(dates_unique[i])
if (dates_unique[i] == "2000-12-31") usbi$year = 2010
# plot(usbi$geometry)
usbji = left_join(usbi, statepop)
# plot(usbji["population"])
usbj = bind_rows(usbj, usbji)
}
summary(usbj)
usa_contig = usbji[!grepl(pattern = "Alaska|Haw", usbji$name), ]
usbj_contig9311 = st_intersection(usbj, usa_contig[0]) |>
st_transform("EPSG:9311") |>
st_collection_extract("POLYGON")
pal = viridis::viridis(n = 7, direction = -1)
pb = c(0, 1, 2, 5, 10, 20, 30, 40) * 1e6
facet_anim = tm_shape(usbj_contig9311) +
tm_polygons(fill = "population",
fill.scale = tm_scale(values = pal, breaks = pb, label.na = FALSE)) +
tm_facets(by = "year", nrow = 1, ncol = 1, free.coords = FALSE) +
tm_shape(usa_union) + tm_borders(lwd = 2) +
tm_layout(legend.position = c("left", "bottom"), legend.frame = FALSE)
tmap_animation(tm = facet_anim, filename = "09-us_pop.gif", width = 900, height = 600)
browseURL("09-us_pop.gif")