API tools faq
paste
MetricT

Modified US State/County maps with US territories

MetricT
Sep 26th, 2020 (edited)
988
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
R 9.04 KB | None | 0 0
raw download clone embed print report
  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)
Add Comment
Please, Sign In to add comment
Public Pastes
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!