Plotting US maps with selected locations (longitude and latitude of locations are known) in R?

Question

I've been struggling with the following R plotting of the US map and wonder if any of you have suggestions on how the plot should be done.

Here is what I have tried:

US <- map("state",fill=TRUE, plot=FALSE)
US.names <- US$names
US.IDs <- sapply(strsplit(US.names,":"),function(x) x[1])
US_poly_sp <- map2SpatialPolygons(US,IDs=US.IDs,proj4string=CRS("+proj=longlat + datum=wgs84"))
plot(US_poly_sp,col="white",axes=TRUE)

With the above code, I can add points if I need to. For example, suppose I have a location with lat=35 and lon=-110, then I can simply use the next line since the latitude and longitude of the spatial polygon I created from the object US contains information on the location:

points(lon,lat)

How can I get a plot of the map that is curved, as generated from the following line and also add points which are given in terms of their latitude and longitude?

US <- map("state",proj="albers", par=c(30,40), plot=TRUE)

I know it has something to do with the projection of maps, but the moment I projected the map, the corresponding range of the latitude and longitude (US$range) changes.

Answer 1

The best resource is the ggmap and ggplot2 packages in R. Here is a short paper describing details. http://journal.r-project.org/archive/2013-1/kahle-wickham.pdf

There are also materials from a workshop on visualizing climate data at this site: http://streaming.stat.iastate.edu/~dicook/NCAR/

Answer 2

Instead of points use SpatialPoints. Then you can reproject layers as you wish using spTransform from the rgdal package.

library(maps)
library(maptools)
library(sp)
library(rgdal)
US <- map("state",fill=TRUE, plot=FALSE)
US.names <- US$names
US.IDs <- sapply(strsplit(US.names,":"),function(x) x[1])
US_poly_sp <- map2SpatialPolygons(US,IDs=US.IDs,proj4string=CRS("+proj=longlat + datum=wgs84"))

#point
point_sp <- SpatialPoints(data.frame(x=-110,y=35), proj4string=CRS(proj4string(US_poly_sp)))

# new CRS:
albers.crs <- 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")

# reproject
US_poly_sp.albers <- spTransform(US_poly_sp, CRS=albers.crs)
point_sp.albers <- spTransform(point_sp, CRS=albers.crs)

plot(US_poly_sp.albers,col="white",axes=TRUE)
points(point_sp.albers)