Natural Earth – world maps – Excellent world shapefile maps with national boundaries, first order admin (provinces, states, etc.), disputed areas, coastlines, populated areas and more at various scale levels.
United Nations – UN-SALB (Second Administrative Level Boundaries) data set project. These maps are free but copyright rests with the UN. The maps are validated by the National Mapping Agencies (NMA) of each UN Member State. This website requires registration to download the maps.
DIVA-GIS – Free country level maps, including administrative boundaries, roads, railroads, altitude and land cover.
Geocommons – A huge database of user submitted maps. Search for any kind of map, then click on the “shapefile” button to the right of the map (if available) to download the shapefiles.
World Bank – Country shapefiles can be downloaded through the World Bank’s Mapping platform as follows: Click on ‘explore countries’ to select a country, click on ‘download’ in the top-right corner, and select ‘Indicators (shapefile)’. This option is not available for all countries.
Map Library – Country shapefile maps for Africa in the public domain.
MapCruzin – A great source for all kinds of shapefiles.
Open Street Map – A growing collection of shapefile map downloads by continent, region and country. These maps are community edited and are not always complete. Shapefiles by country are also available from VDS tech.
US Census Bureau – USA County and District Cartographic Boundary Files.

Drawing maps

1.1 rworldmap, World Map and countries

library(rworldmap) # start with the entire world newmap <- getMap(resolution = "low") plot(newmap, main = "World") # crop to the area desired (outside US) # (can use maps.google.com, right-click, drop lat/lon markers at corners) plot(newmap , xlim = c(-139.3, -58.8) # if you reverse these, the world gets flipped , ylim = c(13.5, 55.7) , asp = 1 # different aspect projections , main = "US from worldmap" )

1.2 ggmap, World Map and countries

library(ggplot2) map.world <- map_data(map = "world") # map = name of map provided by the maps package. # These include county, france, italy, nz, state, usa, world, world2. str(map.world) # how many regions length(unique(map.world$region)) # how many group polygons (some regions have multiple parts) length(unique(map.world$group)) p1 <- ggplot(map.world, aes(x = long, y = lat, group = group)) p1 <- p1 + geom_polygon() # fill areas p1 <- p1 + labs(title = "World, plain") #print(p1) p2 <- ggplot(map.world, aes(x = long, y = lat, group = group, colour = region)) p2 <- p2 + geom_polygon() # fill areas p2 <- p2 + theme(legend.position="none") # remove legend with fill colours p2 <- p2 + labs(title = "World, colour borders") #print(p2) p3 <- ggplot(map.world, aes(x = long, y = lat, group = group, fill = region)) p3 <- p3 + geom_polygon() # fill areas p3 <- p3 + theme(legend.position="none") # remove legend with fill colours p3 <- p3 + labs(title = "World, filled regions") #print(p3) p4 <- ggplot(map.world, aes(x = long, y = lat, group = group, colour = region)) p4 <- p4 + geom_path() # country outline, instead p4 <- p4 + theme(legend.position="none") # remove legend with fill colours p4 <- p4 + labs(title = "World, path outlines only") #print(p4) library(gridExtra) grid.arrange(p1, p2, p3, p4, ncol=2, main="ggmap examples")

1.3 ggmap, New Mexico and Albuquerque

library(ggmap) library(mapproj) map <- get_map( location = "New Mexico" # google search string , zoom = 7 # larger is closer , maptype = "hybrid" # map type ) p <- ggmap(map) p <- p + labs(title = "NM hybrid") print(p) # some options are cute, but not very informative map <- get_map( location = "Albuquerque, New Mexico" # google search string , zoom = 10 # larger is closer , maptype = "watercolor" # map type ) p <- ggmap(map) p <- p + labs(title = "Albuquerque watercolor") print(p)

1.2 Adding data to map underlay

1.2.1 Points

# identify some points around campus dat <- read.table(text = " location lat long MathStat 35.08396 -106.62410 Ducks 35.08507 -106.62238 SC1Class 35.08614 -106.62349 Biology 35.08243 -106.62296 CSEL 35.08317 -106.62414 ", header = TRUE) ## Sometimes the watercolor style can look nice. # get map layer map <- get_map( location = "University of New Mexico" # google search string , zoom = 16 # larger is closer , maptype = "watercolor" # map type ) # plot map p <- ggmap(map) p <- p + geom_point(data = dat, aes(x = long, y = lat, shape = location, colour = location), size = 7) p <- p + geom_text(data = dat, aes(x = long, y = lat, label = location), hjust = -0.2) # legend positioning, removing grid and axis labeling p <- p + theme( legend.position = "none" # remove legend , panel.grid.major = element_blank() , panel.grid.minor = element_blank() , axis.text = element_blank() , axis.title = element_blank() , axis.ticks = element_blank() ) p <- p + labs(title = "UNM SC1 locations") print(p) # Let s say I started in my office in Math & Stat, # then visited with the Ducks, # then taught the SC1 class, # then walked over to Biology, # then finished by picking up a book in the CSEL library. ## Satellite view with points plotted from get_googlemap() # the points need to be called "x" and "y" to get the google markers and path dat.pts <- data.frame(x = dat$long, y = dat$lat) # get map layer map <- get_googlemap( "University of New Mexico" # google search string , zoom = 16 # larger is closer , maptype = "satellite" # map type , markers = dat.pts # markers for map , path = dat.pts # path, in order of points , scale = 2 ) # plot map p <- ggmap(map , extent = "device" # remove white border around map , darken = 0.2 # darken map layer to help points stand out ) p <- p + geom_text(data = dat, aes(x = long, y = lat, label = location), hjust = -0.2, colour = "white", size = 6) # legend positioning, removing grid and axis labeling p <- p + theme( legend.position = c(0.05, 0.05) # put the legend inside the plot area , legend.justification = c(0, 0) , legend.background = element_rect(colour = F, fill = "white") , legend.key = element_rect(fill = F, colour = F) , panel.grid.major = element_blank() , panel.grid.minor = element_blank() , axis.text = element_blank() , axis.title = element_blank() , axis.ticks = element_blank() ) p <- p + labs(title = "UNM Walk around campus") print(p)

1.3 Incidence and density maps

str(crime) # Extract location of crimes in houston violent_crimes <- subset(crime, ((offense != "auto theft") & (offense != "theft") & (offense != "burglary"))) # rank violent crimes violent_crimes$offense <- factor(violent_crimes$offense , levels = c("robbery", "aggravated assault" , "rape", "murder")) # restrict to downtown violent_crimes <- subset(violent_crimes, ((-95.39681 <= lon) & (lon <= -95.34188) & (29.73631 <= lat) & (lat <= 29.784))) map <- get_map( location = "Houston TX" , zoom = 14 , maptype = "roadmap" , color = "bw" # make black & white so color is data ) p <- ggmap(map) p <- p + geom_point(data = violent_crimes , aes(x = lon, y = lat, size = offense, colour = offense)) # legend positioning, removing grid and axis labeling p <- p + theme( legend.position = c(0.0, 0.7) # put the legend inside the plot area , legend.justification = c(0, 0) , legend.background = element_rect(colour = F, fill = "white") , legend.key = element_rect(fill = F, colour = F) , panel.grid.major = element_blank() , panel.grid.minor = element_blank() , axis.text = element_blank() , axis.title = element_blank() , axis.ticks = element_blank() ) print(p) # 2D density plot p <- ggmap(map) overlay <- stat_density2d(data = violent_crimes , aes(x = lon, y = lat, fill = ..level.. , alpha = ..level..) , size = 2, bins = 4, geom = "polygon") p <- p + overlay p <- p + scale_fill_gradient("ViolentnnCrimennDensity") p <- p + scale_alpha(range = c(0.4, 0.75), guide = FALSE) p <- p + guides(fill = guide_colorbar(barwidth = 1.5, barheight = 10)) #p <- p + inset(grob = ggplotGrob(ggplot() + overlay + theme_inset()) # , xmin = -95.35836, xmax = Inf, ymin = -Inf, ymax = 29.75062) print(p) p <- p + facet_wrap( ~ day, nrow = 2) print(p)

1.4 Minard’s map, modern

library(ggplot2) library(plyr) troops <- read.table("http://stat405.had.co.nz/data/minard-troops.txt", header=T) cities <- read.table("http://stat405.had.co.nz/data/minard-cities.txt", header=T) russia <- map_data("world", region = "USSR") p <- ggplot(troops, aes(long, lat)) p <- p + geom_polygon(data = russia, aes(x = long, y = lat, group = group), fill = "white") p <- p + geom_path(aes(size = survivors, colour = direction, group = group), lineend = "round") p <- p + geom_text(data = cities, aes(label = city), size = 3) p <- p + scale_size(range = c(1, 6) , breaks = c(1, 2, 3) * 10^5 , labels = c(1, 2, 3) * 10^5) p <- p + scale_colour_manual(values = c("bisque2", "grey10")) p <- p + xlab(NULL) p <- p + ylab(NULL) p <- p + coord_equal(xlim = c(20, 40), ylim = c(50, 60)) print(p)

1.5 Choropleth maps

library(maps) library(ggplot2) library(plyr) # make fake choropleth data newmexico <- map("county", regions = "new mexico", plot = FALSE, fill = TRUE) newmexico <- fortify(newmexico) newmexico <- ddply(newmexico, "subregion", function(df){ mutate(df, fake = rnorm(1)) } # make standard ggplot map (without geom_map) p <- ggplot(newmexico, aes(x = long, y = lat, group = group, fill = fake)) p <- p + geom_polygon(colour = "white", size = 0.3) print(p) # Now, a fancier map using ggmap... library(ggmap) p <- qmap( New Mexico , zoom = 7, maptype = satellite , legend = topleft ) p <- p + geom_polygon(data = newmexico , aes(x = long, y = lat, group = group, fill = fake) , color = white , alpha = .75, size = .2) # Add some city names, by looking up their location cities <- c("Albuquerque NM", "Las Cruces NM", "Rio Rancho NM", "Santa Fe NM", "Roswell NM", "Farmington NM", "South Valley NM", "Clovis NM", "Hobbs NM", "Alamogordo NM", "Carlsbad NM", "Gallup NM", "Los Alamos NM") cities_locs <- geocode(cities) cities_locs$city <- cities p <- p + geom_text(data = cities_locs, aes(label = city) , color = yellow , size = 3) print(p)
Reference: