################################################################################

### This script takes state/county shapefiles from the Census Bureau, and moves

### Alaska, Hawaii, Puerto Rico, US Virgin Islands, Guam, Samoa, and the Mariana

### Islands to be just under the continental US.   This format is very useful for

### mapping data.

###

### I originally used the "state_laea" and "county_laea" maps from the 

### "tidycensus" package, but kept hitting various issues:

###

### 1) The Tidycensus maps are low-resolution, so counties don't look good if 

###    you zoom in on a particular state

###

### 2) The "Tidycensus" maps are old and require various FIPS mods every time I 

###    use them, such as changing fips 46113 -> 46102 to handle Oglala Co SD, 

###    changing fips 02270 -> 02158 to handle Kusilvak Census Area, AK

###

### 3) The Tidycensus map lacked the smaller US territories.

###

### I took inspiration for this script from the URL below:

###

### https://rud.is/b/2014/11/16/moving-the-earth-well-alaska-hawaii-with-r/

###

### If you find any bugs or issues with the script or with the maps it generates

### please let me know - /u/MetricT

###

### KNOWN BUGS:

###

### * Sourcing the script will cause errors at these lines near the bottom:

###

###    > state_map  %>% st_as_sf() %>% write_sf(mod_state_map)

###    There were 50 or more warnings (use warnings() to see the first 50)

###    > county_map %>% st_as_sf() %>% write_sf(mod_county_map)

###    There were 50 or more warnings (use warnings() to see the first 50)

###

###   Just run them manually at that point and it should work fine.

################################################################################

library(tidyverse)

library(maptools)

library(mapproj)

library(rgeos)

library(rgdal)

library(RColorBrewer)

library(ggplot2)

library(sf)

### Path/filename to save modified shapefiles at when we're done

mod_state_map  <- "../Shapefiles/us_state/us_state.shp"

mod_county_map <- "../Shapefiles/us_county/us_county.shp"

### We want our map to use the Albers projection

map_crs <- "+proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs"

#map_crs <- "epsg:5070"

state_map <- 

  read_sf("../Shapefiles/cb_2018_us_state_500k/cb_2018_us_state_500k.shp") %>%

  st_transform(map_crs) %>%

  as(Class = "Spatial")

county_map <- 

  read_sf("../Shapefiles/cb_2018_us_county_500k/cb_2018_us_county_500k.shp") %>%

  st_transform(map_crs) %>%

  as(Class = "Spatial")

state_map@data$id  <- rownames(state_map@data)

county_map@data$id <- rownames(county_map@data)

### We want to strip out states/territories that aren't in the continental US

### and move them around to make a more convenient map

non_conus_fips <- c("02", "11", "15", "60", "66", "69", "72", "78")

################################################################################

### Adjust the state-level map

################################################################################

### Alaska

alaska_state <- state_map[state_map$STATEFP=="02",]

alaska_state <- elide(alaska_state, rotate=-50)

alaska_state <- elide(alaska_state, scale=max(apply(bbox(alaska_state), 1, diff)) / 1.8)

alaska_state <- elide(alaska_state, shift=c(-2400000, -2800000))

proj4string(alaska_state) <- proj4string(state_map)

### Hawaii

hawaii_state <- state_map[state_map$STATEFP=="15",]

hawaii_state <- elide(hawaii_state, rotate=-35)

hawaii_state <- elide(hawaii_state, shift=c(5800000, -1900000))

proj4string(hawaii_state) <- proj4string(state_map)

### Puerto Rico

puertorico_state <- state_map[state_map$STATEFP=="72",]

puertorico_state <- elide(puertorico_state, rotate=+13)

puertorico_state <- elide(puertorico_state, scale=max(apply(bbox(puertorico_state), 1, diff)) / 0.5)

puertorico_state <- elide(puertorico_state, shift=c(+600000, -2600000))

proj4string(puertorico_state) <- proj4string(state_map)

### US Virgin Islands

usvi_state <- state_map[state_map$STATEFP=="78",]

usvi_state <- elide(usvi_state, rotate=+13)

usvi_state <- elide(usvi_state, scale=max(apply(bbox(usvi_state), 1, diff)) / 0.25)

usvi_state <- elide(usvi_state, shift=c(+1500000, -2600000))

proj4string(usvi_state) <- proj4string(state_map)

### Guam

guam_state <- state_map[state_map$STATEFP=="66",]

guam_state <- elide(guam_state, rotate=-65)

guam_state <- elide(guam_state, scale=max(apply(bbox(guam_state), 1, diff)) / 0.15)

guam_state <- elide(guam_state, shift=c(+1200000, -3200000))

proj4string(guam_state) <- proj4string(state_map)

### Northern Mariana Islands

noma_state <- state_map[state_map$STATEFP=="69",]

noma_state <- elide(noma_state, rotate=-55)

noma_state <- elide(noma_state, scale=max(apply(bbox(noma_state), 1, diff)) / 0.85)

noma_state <- elide(noma_state, shift=c(+300000, -3400000))

proj4string(noma_state) <- proj4string(state_map)

### American Samoa

amsam_state <- state_map[state_map$STATEFP=="60",]

amsam_state <- elide(amsam_state, rotate=-55)

amsam_state <- elide(amsam_state, scale=max(apply(bbox(amsam_state), 1, diff)) / 0.25)

amsam_state <- elide(amsam_state, shift=c(-2300000, -3400000))

proj4string(amsam_state) <- proj4string(state_map)

### Add the moved states/territories back to the CONUS map

state_map <- state_map[!state_map$STATEFP %in% non_conus_fips,]

state_map <- rbind(state_map, 

                   alaska_state, 

                   hawaii_state, 

                   puertorico_state,

                   usvi_state, 

                   noma_state, 

                   guam_state, 

                   amsam_state

                   )

################################################################################

### Adjust the county-level map

################################################################################

### Alaska

alaska_county <- county_map[county_map$STATEFP=="02",]

alaska_county <- elide(alaska_county, rotate=-50)

alaska_county <- elide(alaska_county, scale=max(apply(bbox(alaska_county), 1, diff)) / 1.8)

alaska_county <- elide(alaska_county, shift=c(-2400000, -2800000))

proj4string(alaska_county) <- proj4string(county_map)

### Hawaii

hawaii_county <- county_map[county_map$STATEFP=="15",]

hawaii_county <- elide(hawaii_county, rotate=-35)

hawaii_county <- elide(hawaii_county, shift=c(5800000, -1900000))

proj4string(hawaii_county) <- proj4string(county_map)

### Puerto Rico

puertorico_county <- county_map[county_map$STATEFP=="72",]

puertorico_county <- elide(puertorico_county, rotate=+13)

puertorico_county <- elide(puertorico_county, scale=max(apply(bbox(puertorico_county), 1, diff)) / 0.5)

puertorico_county <- elide(puertorico_county, shift=c(+600000, -2600000))

proj4string(puertorico_county) <- proj4string(county_map)

### US Virgin Islands

usvi_county <- county_map[county_map$STATEFP=="78",]

usvi_county <- elide(usvi_county, rotate=+13)

usvi_county <- elide(usvi_county, scale=max(apply(bbox(usvi_county), 1, diff)) / 0.25)

usvi_county <- elide(usvi_county, shift=c(+1500000, -2600000))

proj4string(usvi_county) <- proj4string(county_map)

### Guam

guam_county <- county_map[county_map$STATEFP=="66",]

guam_county <- elide(guam_county, rotate=-65)

guam_county <- elide(guam_county, scale=max(apply(bbox(guam_county), 1, diff)) / 0.15)

guam_county <- elide(guam_county, shift=c(+1200000, -3200000))

proj4string(guam_county) <- proj4string(county_map)

### Northern Mariana Islands

noma_county <- county_map[county_map$STATEFP=="69",]

noma_county <- elide(noma_county, rotate=-55)

noma_county <- elide(noma_county, scale=max(apply(bbox(noma_county), 1, diff)) / 0.85)

noma_county <- elide(noma_county, shift=c(+300000, -3400000))

proj4string(noma_county) <- proj4string(county_map)

### American Samoa

amsam_county <- county_map[county_map$STATEFP=="60",]

amsam_county <- elide(amsam_county, rotate=-55)

amsam_county <- elide(amsam_county, scale=max(apply(bbox(amsam_county), 1, diff)) / 0.25)

amsam_county <- elide(amsam_county, shift=c(-2300000, -3400000))

proj4string(amsam_county) <- proj4string(county_map)

### Add the moved states/territories back to the CONUS map

county_map <- county_map[!county_map$STATEFP %in% non_conus_fips,]

county_map <- rbind(county_map, 

                   alaska_county, 

                   hawaii_county, 

                   puertorico_county,

                   usvi_county, 

                   noma_county, 

                   guam_county, 

                   amsam_county

)

################################################################################

### Save maps, reload, and graph to make sure it worked ok

###############################################################################

### Save modified maps to new shapefile

state_map  %>% st_as_sf() %>% write_sf(mod_state_map)

county_map %>% st_as_sf() %>% write_sf(mod_county_map)

### Load our newly created maps

state_map  <- read_sf(mod_state_map)

county_map <- read_sf(mod_county_map)

### Overlay the state/county maps so we can verify that the map looks correct

ggplot() +

  theme_void() +

  geom_sf(data = county_map, size = 0.1) +

  geom_sf(data = state_map,  size = 0.6, fill = NA)