In a previous post we looked at how to use D3 TopoJSON files with R and make some very D3-esque maps. I mentioned that one thing missing was moving Alaska & Hawaii a bit closer to the continental United States and this post shows you how to do that.
The D3 folks have it easy. They just use the built in modified Albers composite projection. We R folk have to roll up our sleeves a bit (but not much) to do the same. Thankfully, we can do most of the work with the elide (“ih lied”) function from the
We’ll start with some package imports:
library(maptools) library(mapproj) library(rgeos) library(rgdal) library(RColorBrewer) library(gthemes) library(ggplot2)
I’m using a GeoJSON file that I made from the 2013 US Census shapefile. I prefer GeoJSON mostly due to it being single file and the easy conversion to TopoJSON if I ever need to use the same map in a D3 context (I work with information security data most of the time, so I rarely have to use maps at all for the day job). I simplified the polygons a bit (passing
-simplify 0.01 to
ogr2ogr) to reduce processing time.
We read in that file and then transform the projection to Albers equal area and join the polygon ids to the shapefile data frame:
# https://www.census.gov/geo/maps-data/data/cbf/cbf_counties.html # read U.S. counties moderately-simplified GeoJSON file us <- readOGR(dsn="data/us.geojson", layer="OGRGeoJSON") # convert it to Albers equal area us_aea <- spTransform(us, CRS("+proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs")) us_aea@data$id <- rownames(us_aea@data)
Now, to move Alaska & Hawaii, we have to:
- extract them from the main shapefile data frame
- perform rotation, scaling and transposing on them
- ensure they have the right projection set
- merge them back into the original spatial data frame
elide function has parameters for all the direct shape munging, so we’ll just do that for both states. I took a peek at the D3 source code for the Albers projection to get a feel for the parameters. You can tweak those if you want them in other positions or feel the urge to change the Alaska rotation angle.
# extract, then rotate, shrink & move alaska (and reset projection) # need to use state IDs via # https://www.census.gov/geo/reference/ansi_statetables.html alaska <- us_aea[us_aea$STATEFP=="02",] alaska <- elide(alaska, rotate=-50) alaska <- elide(alaska, scale=max(apply(bbox(alaska), 1, diff)) / 2.3) alaska <- elide(alaska, shift=c(-2100000, -2500000)) proj4string(alaska) <- proj4string(us_aea) # extract, then rotate & shift hawaii hawaii <- us_aea[us_aea$STATEFP=="15",] hawaii <- elide(hawaii, rotate=-35) hawaii <- elide(hawaii, shift=c(5400000, -1400000)) proj4string(hawaii) <- proj4string(us_aea) # remove old states and put new ones back in; note the different order # we're also removing puerto rico in this example but you can move it # between texas and florida via similar methods to the ones we just used us_aea <- us_aea[!us_aea$STATEFP %in% c("02", "15", "72"),] us_aea <- rbind(us_aea, alaska, hawaii)
Rather than just show the resultant plain county map, we’ll add some data to it. The first example uses US drought data (from November 11th, 2014). Drought conditions are pretty severe in some states, but we’ll just show areas that have any type of drought (levels D0-D4). The color ramp shows the % of drought coverage in each county (you’ll need a browser that can display SVGs to see the maps):
# get some data to show...perhaps drought data via: # http://droughtmonitor.unl.edu/MapsAndData/GISData.aspx droughts <- read.csv("data/dm_export_county_20141111.csv") droughts$id <- sprintf("%05d", as.numeric(as.character(droughts$FIPS))) droughts$total <- with(droughts, (D0+D1+D2+D3+D4)/5) # get ready for ggplotting it... this takes a cpl seconds map <- fortify(us_aea, region="GEOID") # plot it gg <- ggplot() gg <- gg + geom_map(data=map, map=map, aes(x=long, y=lat, map_id=id, group=group), fill="#ffffff", color="#0e0e0e", size=0.15) gg <- gg + geom_map(data=droughts, map=map, aes(map_id=id, fill=total), color="#0e0e0e", size=0.15) gg <- gg + scale_fill_gradientn(colours=c("#ffffff", brewer.pal(n=9, name="YlOrRd")), na.value="#ffffff", name="% of County") gg <- gg + labs(title="U.S. Areas of Drought (all levels, % county coverage)") gg <- gg + coord_equal() gg <- gg + theme_map() gg <- gg + theme(legend.position="right") gg <- gg + theme(plot.title=element_text(size=16)) gg
While that shows Alaska & Hawaii in D3-Albers-style, it would be more convincing if we actually used the FIPS county codes on Alaska & Hawaii, so we’ll riff off the previous post and make a county land-mass area choropleth (since we have the land mass area data available in the GeoJSON file):
gg <- ggplot() gg <- gg + geom_map(data=map, map=map, aes(x=long, y=lat, map_id=id, group=group), fill="white", color="white", size=0.15) gg <- gg + geom_map(data=us_aea@data, map=map, aes(map_id=GEOID, fill=log(ALAND)), color="white", size=0.15) gg <- gg + scale_fill_gradientn(colours=c(brewer.pal(n=9, name="YlGn")), na.value="#ffffff", name="County Land\nMass Area (log)") gg <- gg + labs(title="U.S. County Area Choropleth (log scale)") gg <- gg + coord_equal() gg <- gg + theme_map() gg <- gg + theme(legend.position="right") gg <- gg + theme(plot.title=element_text(size=16)) gg
Now, you have one less reason to be envious of the D3 cool kids!
The code & shapefiles are available on github.
hi bob thanks for this! is there a smart way to do the same thing when the coloring itself requires two layers? (such as the
maskcommand used at http://stackoverflow.com/a/26052226/1759499)
Pingback: Adding Annotations to a State Map with Alaska and Hawaii | Bertplot
Pingback: Making Static & Interactive Maps With ggvis (+ using ggvis maps w/shiny) | rud.is
How do I save this back to a geojson file?
This should do it:
Also, that was a really great idea. It could make using this much easier. I’ll try to optimize the polygons and also see if the
choroplethrfolks want a version of it.
Pingback: Easier Composite U.S. Choropleths with albersusa | rud.is
Hi, i keep getting the error message: Error in unit(x, default.units) : ‘x’ and ‘units’ must have length > 0
after the code: gg <- gg + geommap(data=droughts, map=map, aes(mapid=id, fill=total), color=”#0e0e0e”, size=0.15)
I cant place what is causing the error message.