Map projections for the USA

projections.R

# Makes sure the required packs are loaed, if not, intalles them.
pkgs <-c ('sp','ggplot2','rgdal','broom','maptools','tigris', 'blscrapeR')
for(p in pkgs) if(p %in% rownames(installed.packages()) == FALSE) {install.packages(p, repos='http://cran.us.r-project.org')}
for(p in pkgs) suppressPackageStartupMessages(library(p, quietly=TRUE, character.only=TRUE))

# Before we get into the projections, let's download some unemployment data to map.
county_dat <- get_bls_county()

# We also need a shapefile with the appropriate FIPS codes.
state <- counties(cb = TRUE, year = 2015)

# Finally, remove anything that's not the continental US>
# Remove Alaska (02), Hawaii (15), Puerto Rico (72), Guam (66), Virgin Islands (78), 
#American Samoa (60) Mariana Islands (69)
# Micronesia (64), Marshall Islands (68), Palau (70), Minor Islands (74)
state <- state[!state$STATEFP %in% c("02", "15", "72", "66", "78", "60", "69",
                                     "64", "68", "70", "74"),]
# Make sure other outling islands are removed.
state <- state[!state$STATEFP %in% c("81", "84", "86", "87", "89", "71", "76",
                                     "95", "79"),]

# Default projection of the map can be found by running the following on the SpatialPolyGonsDataFrame.
raster::crs(county)

########################################################################
 
 # Chose ONE of the following map projections. 
 # Running them all in succession will only give you the last projection.
  
########################################################################

# Mercator projection.
us.map <- spTransform(state, CRS("+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0  
                                 +x_0=0.0 +y_0=0 +k=1.0 +units=m +no_defs"))

# Lambert Azimuthal Equal Area
us.map  <- spTransform(state, CRS("+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 
                          +ellps=GRS80 +units=m +no_defs")) 

# Lambert Conformal Conic
us.map <- spTransform(state, CRS("+proj=lcc +lat_1=33 +lat_2=45 +lat_0=39 +lon_0=-96 
                                 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs")) 

# Albers Equal Area Conic
us.map <- spTransform(state, CRS("+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=37.5 +lon_0=-96 
                          +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs")) 

# Equidistant Conic
us.map <- spTransform(state, CRS("+proj=eqdc +lat_0=39 +lon_0=-96 +lat_1=33 +lat_2=45 +x_0=0 
                                 +y_0=0 +datum=NAD83 +units=m +no_defs")) 

# Please only select one of the us.map methodes above and execute the following.
# Render spatial object to data frame.
county_map <- broom::tidy(us.map, region="GEOID")


# "Roll your own" ggplot.
ggplot() +
    geom_map(data=county_map, map=county_map,
             aes(x=long, y=lat, map_id=id, group=group),
             fill="#ffffff", color="#0e0e0e", size=0.15) +
    geom_map(data=county_dat, map=county_map, aes_string(map_id="fips", fill="unemployed_rate"),
             color="#0e0e0e", size=0.15) +
    scale_fill_gradientn(colors = c("green", "red")) +
    coord_equal() +
    labs(title="Albers Equal Area") + 
    theme(axis.line=element_blank(),
          axis.text.x=element_blank(),
          axis.text.y=element_blank(),
          axis.ticks=element_blank(),
          axis.title.x=element_blank(),
          axis.title.y=element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.border = element_blank(),
          panel.background = element_blank(),
          legend.title=element_blank())