MetricT

Modified US State/County maps with US territories

Sep 26th, 2020
457
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  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.  
  54. state_map <-
  55.   read_sf("../Shapefiles/cb_2018_us_state_500k/cb_2018_us_state_500k.shp") %>%
  56.   st_transform(map_crs) %>%
  57.   as(Class = "Spatial")
  58.  
  59. county_map <-
  60.   read_sf("../Shapefiles/cb_2018_us_county_500k/cb_2018_us_county_500k.shp") %>%
  61.   st_transform(map_crs) %>%
  62.   as(Class = "Spatial")
  63.  
  64. state_map@data$id  <- rownames(state_map@data)
  65. county_map@data$id <- rownames(county_map@data)
  66.  
  67. ### We want to strip out states/territories that aren't in the continental US
  68. ### and move them around to make a more convenient map
  69. non_conus_fips <- c("02", "11", "15", "60", "66", "69", "72", "78")
  70.  
  71. ################################################################################
  72. ### Adjust the state-level map
  73. ################################################################################
  74.  
  75. ### Alaska
  76. alaska_state <- state_map[state_map$STATEFP=="02",]
  77. alaska_state <- elide(alaska_state, rotate=-50)
  78. alaska_state <- elide(alaska_state, scale=max(apply(bbox(alaska_state), 1, diff)) / 1.8)
  79. alaska_state <- elide(alaska_state, shift=c(-2400000, -2800000))
  80. proj4string(alaska_state) <- proj4string(state_map)
  81.  
  82. ### Hawaii
  83. hawaii_state <- state_map[state_map$STATEFP=="15",]
  84. hawaii_state <- elide(hawaii_state, rotate=-35)
  85. hawaii_state <- elide(hawaii_state, shift=c(5800000, -1900000))
  86. proj4string(hawaii_state) <- proj4string(state_map)
  87.  
  88. ### Puerto Rico
  89. puertorico_state <- state_map[state_map$STATEFP=="72",]
  90. puertorico_state <- elide(puertorico_state, rotate=+13)
  91. puertorico_state <- elide(puertorico_state, scale=max(apply(bbox(puertorico_state), 1, diff)) / 0.5)
  92. puertorico_state <- elide(puertorico_state, shift=c(+600000, -2600000))
  93. proj4string(puertorico_state) <- proj4string(state_map)
  94.  
  95. ### US Virgin Islands
  96. usvi_state <- state_map[state_map$STATEFP=="78",]
  97. usvi_state <- elide(usvi_state, rotate=+13)
  98. usvi_state <- elide(usvi_state, scale=max(apply(bbox(usvi_state), 1, diff)) / 0.25)
  99. usvi_state <- elide(usvi_state, shift=c(+1500000, -2600000))
  100. proj4string(usvi_state) <- proj4string(state_map)
  101.  
  102. ### Guam
  103. guam_state <- state_map[state_map$STATEFP=="66",]
  104. guam_state <- elide(guam_state, rotate=-65)
  105. guam_state <- elide(guam_state, scale=max(apply(bbox(guam_state), 1, diff)) / 0.15)
  106. guam_state <- elide(guam_state, shift=c(+1200000, -3200000))
  107. proj4string(guam_state) <- proj4string(state_map)
  108.  
  109. ### Northern Mariana Islands
  110. noma_state <- state_map[state_map$STATEFP=="69",]
  111. noma_state <- elide(noma_state, rotate=-55)
  112. noma_state <- elide(noma_state, scale=max(apply(bbox(noma_state), 1, diff)) / 0.85)
  113. noma_state <- elide(noma_state, shift=c(+300000, -3400000))
  114. proj4string(noma_state) <- proj4string(state_map)
  115.  
  116. ### American Samoa
  117. amsam_state <- state_map[state_map$STATEFP=="60",]
  118. amsam_state <- elide(amsam_state, rotate=-55)
  119. amsam_state <- elide(amsam_state, scale=max(apply(bbox(amsam_state), 1, diff)) / 0.25)
  120. amsam_state <- elide(amsam_state, shift=c(-2300000, -3400000))
  121. proj4string(amsam_state) <- proj4string(state_map)
  122.  
  123. ### Add the moved states/territories back to the CONUS map
  124. state_map <- state_map[!state_map$STATEFP %in% non_conus_fips,]
  125. state_map <- rbind(state_map,
  126.                    alaska_state,
  127.                    hawaii_state,
  128.                    puertorico_state,
  129.                    usvi_state,
  130.                    noma_state,
  131.                    guam_state,
  132.                    amsam_state
  133.                    )
  134.  
  135. ################################################################################
  136. ### Adjust the county-level map
  137. ################################################################################
  138.  
  139. ### Alaska
  140. alaska_county <- county_map[county_map$STATEFP=="02",]
  141. alaska_county <- elide(alaska_county, rotate=-50)
  142. alaska_county <- elide(alaska_county, scale=max(apply(bbox(alaska_county), 1, diff)) / 1.8)
  143. alaska_county <- elide(alaska_county, shift=c(-2400000, -2800000))
  144. proj4string(alaska_county) <- proj4string(county_map)
  145.  
  146. ### Hawaii
  147. hawaii_county <- county_map[county_map$STATEFP=="15",]
  148. hawaii_county <- elide(hawaii_county, rotate=-35)
  149. hawaii_county <- elide(hawaii_county, shift=c(5800000, -1900000))
  150. proj4string(hawaii_county) <- proj4string(county_map)
  151.  
  152. ### Puerto Rico
  153. puertorico_county <- county_map[county_map$STATEFP=="72",]
  154. puertorico_county <- elide(puertorico_county, rotate=+13)
  155. puertorico_county <- elide(puertorico_county, scale=max(apply(bbox(puertorico_county), 1, diff)) / 0.5)
  156. puertorico_county <- elide(puertorico_county, shift=c(+600000, -2600000))
  157. proj4string(puertorico_county) <- proj4string(county_map)
  158.  
  159. ### US Virgin Islands
  160. usvi_county <- county_map[county_map$STATEFP=="78",]
  161. usvi_county <- elide(usvi_county, rotate=+13)
  162. usvi_county <- elide(usvi_county, scale=max(apply(bbox(usvi_county), 1, diff)) / 0.25)
  163. usvi_county <- elide(usvi_county, shift=c(+1500000, -2600000))
  164. proj4string(usvi_county) <- proj4string(county_map)
  165.  
  166. ### Guam
  167. guam_county <- county_map[county_map$STATEFP=="66",]
  168. guam_county <- elide(guam_county, rotate=-65)
  169. guam_county <- elide(guam_county, scale=max(apply(bbox(guam_county), 1, diff)) / 0.15)
  170. guam_county <- elide(guam_county, shift=c(+1200000, -3200000))
  171. proj4string(guam_county) <- proj4string(county_map)
  172.  
  173. ### Northern Mariana Islands
  174. noma_county <- county_map[county_map$STATEFP=="69",]
  175. noma_county <- elide(noma_county, rotate=-55)
  176. noma_county <- elide(noma_county, scale=max(apply(bbox(noma_county), 1, diff)) / 0.85)
  177. noma_county <- elide(noma_county, shift=c(+300000, -3400000))
  178. proj4string(noma_county) <- proj4string(county_map)
  179.  
  180. ### American Samoa
  181. amsam_county <- county_map[county_map$STATEFP=="60",]
  182. amsam_county <- elide(amsam_county, rotate=-55)
  183. amsam_county <- elide(amsam_county, scale=max(apply(bbox(amsam_county), 1, diff)) / 0.25)
  184. amsam_county <- elide(amsam_county, shift=c(-2300000, -3400000))
  185. proj4string(amsam_county) <- proj4string(county_map)
  186.  
  187. ### Add the moved states/territories back to the CONUS map
  188. county_map <- county_map[!county_map$STATEFP %in% non_conus_fips,]
  189. county_map <- rbind(county_map,
  190.                    alaska_county,
  191.                    hawaii_county,
  192.                    puertorico_county,
  193.                    usvi_county,
  194.                    noma_county,
  195.                    guam_county,
  196.                    amsam_county
  197. )
  198.  
  199. ################################################################################
  200. ### Save maps, reload, and graph to make sure it worked ok
  201. ###############################################################################
  202.  
  203. ### Save modified maps to new shapefile
  204. state_map  %>% st_as_sf() %>% write_sf(mod_state_map)
  205. county_map %>% st_as_sf() %>% write_sf(mod_county_map)
  206.  
  207. ### Load our newly created maps
  208. state_map  <- read_sf(mod_state_map)
  209. county_map <- read_sf(mod_county_map)
  210.  
  211. ### Overlay the state/county maps so we can verify that the map looks correct
  212. ggplot() +
  213.   theme_void() +
  214.   geom_sf(data = county_map, size = 0.1) +
  215.   geom_sf(data = state_map,  size = 0.6, fill = NA)
  216.  
  217.  
RAW Paste Data

Adblocker detected! Please consider disabling it...

We've detected AdBlock Plus or some other adblocking software preventing Pastebin.com from fully loading.

We don't have any obnoxious sound, or popup ads, we actively block these annoying types of ads!

Please add Pastebin.com to your ad blocker whitelist or disable your adblocking software.

×