Home / Tutorials

Map Projections in R

The Earth is a spheroid - it is round like a ball or sphere but flattened by slightly by the centrifugal force of rotation. The Earth has a circumference of around 24,900 miles, but is around 27 miles wider than it is tall.

The Earth exists in three-dimensions but, other than globes, most representations of the earth, such as maps, are two dimensional. A projection is a set of mathematical transformations used to represent the three-dimensional world in two dimensions.

Projections in R are handled by the RGDAL library and are usually specified using PROJ.4 strings that can be given as parameters to the CRS() function, and in turn used when either creating spatial objects from coordinates, or transforming the coordinates to another coordinate reference system.

The data for this example is in 2017-state-data.zip.

To get a proj4 string For R, you can search SpatialReference.org.

library(rgdal)

# Load the shapefile
states = readOGR(".", "2017-state-data", stringsAsFactors=F)

# Select only the contiguous 48 states
states = states[((states$ST != "AK") & (states$ST != "HI")),]

# Transform to Albers Equal Area Conic to make cartographically proper
usa_albers = CRS("+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=37.5 +lon_0=-96 +x_0=0 +y_0=0 
	+ellps=GRS80 +datum=NAD83 +units=m +no_defs")
states = spTransform(states, usa_albers)

# Plot a red-state/blue-state choropleth
colors = ifelse(states$WIN2012 == "Obama", "navy", "red")
plot(states, col = colors)

legend(x="bottomleft", legend=c("Obama","Romney"), 
	title="2012 US Presidential Election", 
	pch=15, col=c("navy","red"))

Example US Choropleth Using an Albers Equal Area Conic Projection

The following is a list of variables with the PROJ.4 strings for commonly used projections.

An extensive list of PROJ.4 strings is available at spatialreference.org.

# EPSG:4326 Unprojected WGS 84 latitude and longitude
# NEVER USED THIS IN A PRINT MAP
wgs84 = CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")

# EPSG: 3395 World Mercator
world_mercator = CRS("+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs")

# EPSG: 3857 Spherical Mercator
spherical_mercator = 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 +nadgrids=@null +wktext  +no_defs")

# EPSG: 54030 World Robinson
robinson = CRS("+proj=robin +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs")

# EPSG: 102008 North America Albers Equal Area Conic
usa_albers = CRS("+proj=aea +lat_1=20 +lat_2=60 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 
	+datum=NAD83 +units=m +no_defs")

# EPSG:3082 US/Texas-Centric Lambert Conformal
usa_lambert = CRS("+proj=lcc +lat_1=27.5 +lat_2=35 +lat_0=18 +lon_0=-100 
	+x_0=1500000 +y_0=5000000 +ellps=GRS80 
	+towgs84=0,0,0,0,0,0,0 +units=m +no_defs")

# EPSG: 32610 UTM zone 10 North (Washington state, WGS 84)
utm10n = CRS("+proj=utm +zone=10 +datum=WGS84 +units=m +no_defs")

# EPSG:3692 State Plane Coordinate System - NAD 83 NRS 2007 Washington South, US feet
washington = CRS("+proj=lcc +lat_1=47.33333333333334 +lat_2=45.83333333333334 
	+lat_0=45.33333333333334 +lon_0=-120.5 +x_0=500000.0001016001 +y_0=0 
	+ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=us-ft +no_defs")