View difference between Paste ID: <a href="/mi1krugp">mi1krugp</a> and <a href="/P7CrXugY">P7CrXugY</a>

################################################################################
1		################################################################################
2		### This script takes state/county shapefiles from the Census Bureau, and moves
3		### Alaska, Hawaii, Puerto Rico, US Virgin Islands, Guam, Samoa, and the Mariana
4		### Islands to be just under the continental US. This format is very useful for
5		### mapping data.
6		###
7		### I originally used the "state_laea" and "county_laea" maps from the
8		### "tidycensus" package, but kept hitting various issues:
9		###
10		### 1) The Tidycensus maps are low-resolution, so counties don't look good if
11		### you zoom in on a particular state
12		###
13		### 2) The "Tidycensus" maps are old and require various FIPS mods every time I
14		### use them, such as changing fips 46113 -> 46102 to handle Oglala Co SD,
15		### changing fips 02270 -> 02158 to handle Kusilvak Census Area, AK
16		###
17		### 3) The Tidycensus map lacked the smaller US territories.
18		###
19		### I took inspiration for this script from the URL below:
20		###
21		### https://rud.is/b/2014/11/16/moving-the-earth-well-alaska-hawaii-with-r/
22		###
23		### If you find any bugs or issues with the script or with the maps it generates
24		### please let me know - /u/MetricT
25		###
26		### KNOWN BUGS:
27		###
28		### * Sourcing the script will cause errors at these lines near the bottom:
29		###
30		### > state_map %>% st_as_sf() %>% write_sf(mod_state_map)
31		### There were 50 or more warnings (use warnings() to see the first 50)
32		### > county_map %>% st_as_sf() %>% write_sf(mod_county_map)
33		### There were 50 or more warnings (use warnings() to see the first 50)
34		###
35		### Just run them manually at that point and it should work fine.
36		################################################################################
37
38		library(tidyverse)
39		library(maptools)
40		library(mapproj)
41		library(rgeos)
42		library(rgdal)
43		library(RColorBrewer)
44		library(ggplot2)
45		library(sf)
46
47		### Path/filename to save modified shapefiles at when we're done
48		mod_state_map <- "../Shapefiles/us_state/us_state.shp"
49		mod_county_map <- "../Shapefiles/us_county/us_county.shp"
50
51		### We want our map to use the Albers projection
52		map_crs <- "+proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs"
53		#map_crs <- "epsg:5070"
54
55		state_map <-
56		read_sf("../Shapefiles/cb_2018_us_state_500k/cb_2018_us_state_500k.shp") %>%
57		st_transform(map_crs) %>%
58		as(Class = "Spatial")
59
60		county_map <-
61		read_sf("../Shapefiles/cb_2018_us_county_500k/cb_2018_us_county_500k.shp") %>%
62		st_transform(map_crs) %>%
63		as(Class = "Spatial")
64
65		state_map@data$id <- rownames(state_map@data)
66		county_map@data$id <- rownames(county_map@data)
67
68		### We want to strip out states/territories that aren't in the continental US
69		### and move them around to make a more convenient map
70		non_conus_fips <- c("02", "11", "15", "60", "66", "69", "72", "78")
71
72		################################################################################
73		### Adjust the state-level map
74		################################################################################
75
76		### Alaska
77		alaska_state <- state_map[state_map$STATEFP=="02",]
78		alaska_state <- elide(alaska_state, rotate=-50)
79		alaska_state <- elide(alaska_state, scale=max(apply(bbox(alaska_state), 1, diff)) / 1.8)
80		alaska_state <- elide(alaska_state, shift=c(-2400000, -2800000))
81		proj4string(alaska_state) <- proj4string(state_map)
82
83		### Hawaii
84		hawaii_state <- state_map[state_map$STATEFP=="15",]
85		hawaii_state <- elide(hawaii_state, rotate=-35)
86		hawaii_state <- elide(hawaii_state, shift=c(5800000, -1900000))
87		proj4string(hawaii_state) <- proj4string(state_map)
88
89		### Puerto Rico
90		puertorico_state <- state_map[state_map$STATEFP=="72",]
91		puertorico_state <- elide(puertorico_state, rotate=+13)
92		puertorico_state <- elide(puertorico_state, scale=max(apply(bbox(puertorico_state), 1, diff)) / 0.5)
93		puertorico_state <- elide(puertorico_state, shift=c(+600000, -2600000))
94		proj4string(puertorico_state) <- proj4string(state_map)
95
96		### US Virgin Islands
97		usvi_state <- state_map[state_map$STATEFP=="78",]
98		usvi_state <- elide(usvi_state, rotate=+13)
99		usvi_state <- elide(usvi_state, scale=max(apply(bbox(usvi_state), 1, diff)) / 0.25)
100		usvi_state <- elide(usvi_state, shift=c(+1500000, -2600000))
101		proj4string(usvi_state) <- proj4string(state_map)
102
103		### Guam
104		guam_state <- state_map[state_map$STATEFP=="66",]
105		guam_state <- elide(guam_state, rotate=-65)
106		guam_state <- elide(guam_state, scale=max(apply(bbox(guam_state), 1, diff)) / 0.15)
107		guam_state <- elide(guam_state, shift=c(+1200000, -3200000))
108		proj4string(guam_state) <- proj4string(state_map)
109
110		### Northern Mariana Islands
111		noma_state <- state_map[state_map$STATEFP=="69",]
112		noma_state <- elide(noma_state, rotate=-55)
113		noma_state <- elide(noma_state, scale=max(apply(bbox(noma_state), 1, diff)) / 0.85)
114		noma_state <- elide(noma_state, shift=c(+300000, -3400000))
115		proj4string(noma_state) <- proj4string(state_map)
116
117		### American Samoa
118		amsam_state <- state_map[state_map$STATEFP=="60",]
119		amsam_state <- elide(amsam_state, rotate=-55)
120		amsam_state <- elide(amsam_state, scale=max(apply(bbox(amsam_state), 1, diff)) / 0.25)
121		amsam_state <- elide(amsam_state, shift=c(-2300000, -3400000))
122		proj4string(amsam_state) <- proj4string(state_map)
123
124		### Add the moved states/territories back to the CONUS map
125		state_map <- state_map[!state_map$STATEFP %in% non_conus_fips,]
126		state_map <- rbind(state_map,
127		alaska_state,
128		hawaii_state,
129		puertorico_state,
130		usvi_state,
131		noma_state,
132		guam_state,
133		amsam_state
134		)
135
136		################################################################################
137		### Adjust the county-level map
138		################################################################################
139
140		### Alaska
141		alaska_county <- county_map[county_map$STATEFP=="02",]
142		alaska_county <- elide(alaska_county, rotate=-50)
143		alaska_county <- elide(alaska_county, scale=max(apply(bbox(alaska_county), 1, diff)) / 1.8)
144		alaska_county <- elide(alaska_county, shift=c(-2400000, -2800000))
145		proj4string(alaska_county) <- proj4string(county_map)
146
147		### Hawaii
148		hawaii_county <- county_map[county_map$STATEFP=="15",]
149		hawaii_county <- elide(hawaii_county, rotate=-35)
150		hawaii_county <- elide(hawaii_county, shift=c(5800000, -1900000))
151		proj4string(hawaii_county) <- proj4string(county_map)
152
153		### Puerto Rico
154		puertorico_county <- county_map[county_map$STATEFP=="72",]
155		puertorico_county <- elide(puertorico_county, rotate=+13)
156		puertorico_county <- elide(puertorico_county, scale=max(apply(bbox(puertorico_county), 1, diff)) / 0.5)
157		puertorico_county <- elide(puertorico_county, shift=c(+600000, -2600000))
158		proj4string(puertorico_county) <- proj4string(county_map)
159
160		### US Virgin Islands
161		usvi_county <- county_map[county_map$STATEFP=="78",]
162		usvi_county <- elide(usvi_county, rotate=+13)
163		usvi_county <- elide(usvi_county, scale=max(apply(bbox(usvi_county), 1, diff)) / 0.25)
164		usvi_county <- elide(usvi_county, shift=c(+1500000, -2600000))
165		proj4string(usvi_county) <- proj4string(county_map)
166
167		### Guam
168		guam_county <- county_map[county_map$STATEFP=="66",]
169		guam_county <- elide(guam_county, rotate=-65)
170		guam_county <- elide(guam_county, scale=max(apply(bbox(guam_county), 1, diff)) / 0.15)
171		guam_county <- elide(guam_county, shift=c(+1200000, -3200000))
172		proj4string(guam_county) <- proj4string(county_map)
173
174		### Northern Mariana Islands
175		noma_county <- county_map[county_map$STATEFP=="69",]
176		noma_county <- elide(noma_county, rotate=-55)
177		noma_county <- elide(noma_county, scale=max(apply(bbox(noma_county), 1, diff)) / 0.85)
178		noma_county <- elide(noma_county, shift=c(+300000, -3400000))
179		proj4string(noma_county) <- proj4string(county_map)
180
181		### American Samoa
182		amsam_county <- county_map[county_map$STATEFP=="60",]
183		amsam_county <- elide(amsam_county, rotate=-55)
184		amsam_county <- elide(amsam_county, scale=max(apply(bbox(amsam_county), 1, diff)) / 0.25)
185		amsam_county <- elide(amsam_county, shift=c(-2300000, -3400000))
186		proj4string(amsam_county) <- proj4string(county_map)
187
188		### Add the moved states/territories back to the CONUS map
189		county_map <- county_map[!county_map$STATEFP %in% non_conus_fips,]
190		county_map <- rbind(county_map,
191		alaska_county,
192		hawaii_county,
193		puertorico_county,
194		usvi_county,
195		noma_county,
196		guam_county,
197		amsam_county
198		)
199
200		################################################################################
201		### Save maps, reload, and graph to make sure it worked ok
202		###############################################################################
203
204		### Save modified maps to new shapefile
205		state_map %>% st_as_sf() %>% write_sf(mod_state_map)
206		county_map %>% st_as_sf() %>% write_sf(mod_county_map)
207
208		### Load our newly created maps
209		state_map <- read_sf(mod_state_map)
210		county_map <- read_sf(mod_county_map)
211
212		### Overlay the state/county maps so we can verify that the map looks correct
213		ggplot() +
214		theme_void() +
215		geom_sf(data = county_map, size = 0.1) +
216		geom_sf(data = state_map, size = 0.6, fill = NA)