Administrative regions map of a country with ggmap and ggplot2

I don't know the spatial level of administrative areas you require, but here's two waysto read in shapefile data and .RData formats from Global Administrative Areas (gadm.org), and converting them into data frames for use in ggplot2. Also, in order to replicate the U.S. map, you will need to plot the administrative area names located at the polygon centroids.

library(ggplot2)library(rgdal)

Method 1. SpatialPolygonDataFrames stored as .RData format

# Data from the Global Administrative Areas# 1) Read in administrative area level 2 dataload("/Users/jmuirhead/Downloads/PAK_adm2.RData")pakistan.adm2.spdf <- get("gadm")

Method2. Shapefile format read in with rgdal::readOGR

pakistan.adm2.spdf <- readOGR("/Users/jmuirhead/Downloads/PAK_adm", "PAK_adm2",  verbose = TRUE, stringsAsFactors = FALSE)

Creating a data.frame from the spatialPolygonDataframes and merging with a data.frame containing the information on unemployment, for example.

pakistan.adm2.df <- fortify(pakistan.adm2.spdf, region = "NAME_2")# Sample dataframe of unemployment infounemployment.df <- data.frame(id= unique(pakistan.adm2.df[,'id']),  unemployment = runif(n = length(unique(pakistan.adm2.df[,'id'])), min = 0, max = 25))pakistan.adm2.df <- merge(pakistan.adm2.df, unemployment.df, by.y = 'id', all.x = TRUE)

Extracting names and centoids of administrative areas for plotting

# Get centroids of spatialPolygonDataFrame and convert to dataframe# for use in plotting  area names. pakistan.adm2.centroids.df <- data.frame(long = coordinates(pakistan.adm2.spdf)[, 1],    lat = coordinates(pakistan.adm2.spdf)[, 2]) # Get names and id numbers corresponding to administrative areaspakistan.adm2.centroids.df[, 'ID_2'] <- pakistan.adm2.spdf@data[,'ID_2']pakistan.adm2.centroids.df[, 'NAME_2'] <- pakistan.adm2.spdf@data[,'NAME_2']

Create ggplot with labels for administrative areas

p <- ggplot(pakistan.adm2.df, aes(x = long, y = lat, group = group)) + geom_polygon(aes(fill = cut(unemployment,5))) +geom_text(data = pakistan.adm2.centroids.df, aes(label = NAME_2, x = long, y = lat, group = NAME_2), size = 3) + labs(x=" ", y=" ") + theme_bw() + scale_fill_brewer('Unemployment Rate (Jan 2011)', palette  = 'PuRd') + coord_map() + theme(panel.grid.minor=element_blank(), panel.grid.major=element_blank()) + theme(axis.ticks = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank()) + theme(panel.border = element_blank())print(p)

Have you tried http://www.diva-gis.org/gdata? That should give you Shapefiles with state/district boundaries for Pakistan.

The maptools package has functions that let you read and convert Shapefiles to a dataframe that can be used to lay out the boundary polygons.