R: Using map() and plot() with the same projection

Question

I am trying to plot a shapefile (SpatialPolygonDataFrame in R) of 180 census tracts on top of an outline of the US using map(). The .prj file for the census tracts lists the projection as follows:

 PROJCS["USA_Contiguous_Albers_Equal_Area_Conic",GEOGCS["GCS_North_American_1983", 
 DATUM["D_North_American_1983", 
 SPHEROID["GRS_1980",6378137.0,298.257222101]],
 PRIMEM["Greenwich",0.0], 
 UNIT["Degree",0.0174532925199433]], 
 PROJECTION["Albers"], 
 PARAMETER["False_Easting",0.0], 
 PARAMETER["False_Northing",0.0], 
 PARAMETER["Central_Meridian",-96.0],
 PARAMETER["Standard_Parallel_1",29.5], 
 PARAMETER["Standard_Parallel_2",45.5],
 PARAMETER["Latitude_Of_Origin",37.5],
 UNIT["Meter",1.0]]

And when I call proj4string() on the tracts, I get the projection in this format.

"+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 +ellps=GRS80 +towgs84=0,0,0"

I tried setting the projection within map() to match this .prj file, but when I go to plot my census tracts on top, they don't show up anywhere.

 map("usa", project="albers", par=c(29.5,45.5))
 plot(tracts, add=TRUE)

I feel I'm misunderstanding either how to add the projection properly, or how to over-plot polygons on top of map().

Answer 1

I think your best bet is to get the US outline from another source. I don't think the map function really handles projections quite the same way as the sp functions.

For example, your USA map with albers above - look at par()$usr or show the coordinates with axis(1);axis(2) - those are small numbers around 0 to plus or minus 1. Your projection string for your data has units of m - metres. maps has not scaled the USA to metres.

I suggest you get the raster package and use getData to get the USA outline from the GADM database. Instructions are there. Then you can project it to your data coordinates using spTransform from the rgdal and sp packages.