Below is how to get map data for extraction and plotting in R, including how to add points. The R solution was developed by John Robinson for our lab.
One way is to just go to http://www.aquarius.ifm-geomar.de/ but not quite as nice or clean as finished products in R unless you're good with Illustrator. So you've got two good options.
Marine
Go to http://rimmer.ngdc.noaa.gov/mgg/coast/getcoast.html and click on "Use Java Map for Lat/Lon". Make your map area and accept it to return to the Coastline Extractor page. I chose "World Vector Shoreline" which is at 1:250,000 scale. Put the format in "Splus" and click SUBMIT to save as a .dat file. Note if you are covering a large area, it may say Warning: 500,000 point limit reached. No more coast extracted!! in which case you may need to adjust your resolution or bounding box.
You can then modify these R commands, written by John Robinson:
#Map of East Coast, with a box around GCE-LTER and an inset map of GCE
#Be sure to set the working directory to the place where you have the LoRes.dat files saved
setwd("/Users/jpw/Desktop");
#As the code presently stands the window needs to be resized immediately after plotting, not sure how to do this in code... drag it out to about 600 x 600 and it looks a little better!
#Reads in the two .dat files, one for the east coast and another for the GCE-LTER (note this is changed now to grab carib data)
#smallmaplines = read.table("LowResGCE.dat", header=FALSE);
bigmaplines = read.table("carib.dat", header=FALSE);
#Plotting the Coast
par(pty="s");
plot(bigmaplines, type="l", xlab="Degrees West Longitude",ylab="Degrees North Latitude", bty="l", asp=1);
#Vectors for lat and long of sites, enter long first (in vec)
#Input values for these vectors to get points for our other sample sites
vec = c(-76,-84,-85,-82.25,-79.75,-64.8,-67,-67);
tor = c(24.25,10,10,9.33,9.5,18.33,17.96,18.8);
points(vec,tor,pch=19, cex=1.3);
#This part of the code generates a smaller window where we want to plot our insert graph
#The dimensions of this window are set based on the range of values in the smallmaplines file
xrange = range(smallmaplines[,2], na.rm=TRUE);
yrange = range(smallmaplines[,1], na.rm=TRUE);
aspect = abs(diff(yrange))/abs(diff(xrange));
#Par makes a new plot region inside the old one with the dimensions specified
par(fig=c(0.99-0.4, 0.99, 0.15, 0.15+0.4 *aspect), mar=rep(0,4), new=TRUE);
plot.new();
#Sets the limits of the plot window created
plot.window(xlim=xrange, ylim=yrange);
plot(smallmaplines, type="l", lwd=0.5, tcl=0, ann=FALSE, col.axis= "white");
#Lat and longitude for the GCE sites and adding points for those sites
GCEvec = c(-81.419506, -81.304664, -81.215763, -81.364286, -81.341052, -81.280885, -81.335514, -81.275928);
GCEtor = c(31.540321, 31.542803, 31.534771, 31.456326, 31.430951, 31.382560, 31.348796, 31.481690);
points(GCEvec,GCEtor,pch=19, cex=1.3);
#Putting small white numbers on each of the GCE sites
names = c(1,2,3,4,5,6,9,10);
for(i in 1:8){
- text(GCEvec[i], (GCEtor[i]+0.004), names[i], col="white", cex=0.85)

Freshwater
Hey guys,
Had some success with making maps in R using ESRI shapefiles from the US Wild and Scenic Rivers website (http://www.rivers.gov/maps.html). I'm attaching an R script that uses the zip file from that website to plot a window on East Tennessee (see the figure). Thought this might be useful for future stuff.
John
#Maps based on the GIS map provided by the National Wild and Scenic Rivers System @ http://www.rivers.gov/maps.html
#Save the .zip file and set the working directories of the hydrom and states folders below
#Everything afte the FINALS_GIS is the construction as it stands in the .zip file as of 6/13
library(maptools);
#Read in sample site information
setwd("/Users/John/Desktop");
sites = read.table("Sites.txt", header=TRUE);
#These vectors give x and y lim values for the plot
#These NEED / DONT NEED to be a square
#Not sure, seems like yes to me, or at least to match the fig dimensions?
longs = c(-85.1, -83.1);
lats = c(34.75, 36.75);
#Note that these lines take a bit (a large bit) of time, reading in lots of data (for the whole U.S.)
#STATES
setwd("/Users/John/Desktop/FINALS_GIS/Data_WSR/NatlAtlas/statesp020");
states = readShapeLines("statesp020");
#HYDROM
setwd("/Users/John/Desktop/FINALS_GIS/Data_WSR/NatlAtlas/hydrom020");
rivs1 = readShapeLines("hydrogp020");
rivs2 = readShapeLines("hydrogl020");
#Want to remove all the bullshit whitespace from our map, R loves the margins
par(mar = c(0,0,0,0));
#This is the outline of all the US states if xlim and ylim are not set by the vectors above
#Otherwise, it's a nice plot with heavy lines for state boundaries
plot(states, xlim = longs, ylim = lats, lwd = 2, lty="dotted");
#This starts drawing lakes and rivers in with polygons and lines
#I think it looks nice with 0.5 weight.
plot(rivs1, add=T, lwd = 0.5);
plot(rivs2, add=T, lwd = 0.5);
#I'd recommend adding in place names with the locator(), sample sites should go in easy with GPS.
points(sites$Long, sites$Lat, pch = 19);
text(sites$Long, sites$Lat, sites$Pop, pos=sites$pos, cex=0.75);
lines(c(-84.96807, -84.96807+(2*5.39956803*(1/60))), c(36.39265, 36.39265), lty = "solid", lwd = 1.6);
text(x= -84.87724,y=36.39265,"20 km",pos=3);
text(x = -84.70934, y = 36.21498 , "TN", cex = 1.5);
text(x = -84.49190, y = 34.89201, "GA", cex = 1.5);
text(x = -83.53132, y = 35.19816, "NC", cex = 1.5);