Skip to content

How can I plot a shapefile with geom_map from ggplot2

4 messages · Manuel Spínola, Roman Luštrik, O'Hanlon, Simon J

Original
O'Hanlon, Simon J
# 
Dear Manuel,
I had the same problem as you.

I googled: "plotting a shapefile with ggplot2"

First hit:

https://github.com/hadley/ggplot2/wiki/plotting-polygon-shapefiles

Using this guide it only take a couple of very minor adjustments to plot using geom_map() rather than geom_polygon(). My example below plots country borders of West African countries using the shapefiles provided by the Natural Earth website.

Note that on the tutorial webpage they use gpclibPermit() to use gpclib for polygon union operations in the fortify function. You should omit this and the fortify method will now select rgeos to union polygons by default (works from ggplot2 0.9.1 onwards - possibly 0.9 as well). The example below uses both geom_map() for outlines and geom_map() for filled polygons.


#Load relevant libraries
library(ggplot2)
library(maptools)
library(rgeos)


#	Download and unzip the Natural Earth Country Polygon data
oldwd <- getwd()
tmp <- tempdir()
setwd(tmp)
url <- "http://www.naturalearthdata.com/http//www.naturalearthdata.com/download/50m/cultural/50m-admin-0-countries.zip"
dest <- file.path(tmp,"tmp.zip")
download.file(url,dest)	#File is 1.3Mb
unzip(dest)

#	Read in the data to a SpatialPolygonsDataFrame and subset to West African countries
wld <- readShapePoly("ne_50m_admin_0_countries")
wa <- wld[c(21,22,43,44,81,82,84,144,151,158,160,190,195,213),]

#	Convert to data.frame for plotting with ggplot2
wa at data$id <- rownames(wa at data) # Explicitly identifies attribute rows by the .dbf offset
wa.df <- as.data.frame(wa)
gpclibPermit() # required for fortify method if rgeos is unavailable
wa.fort <- fortify(wa, region="id")
wa.gg <- join(wa.fort, wa.df,by="id")

#	Create the plot
p <- ggplot(wa.gg)+
	geom_map(map=wa.gg,aes(long,lat,map_id=id))+
	geom_path(colour="white",aes(long,lat,group=group),size=1.2)+
	coord_equal()+
	scale_x_continuous(name=expression(paste("Longitude (",degree,")")),limits=c(-18,4),expand=c(0,0))+
	scale_y_continuous(name=expression(paste("Latitude (",degree,")")),limits=c(4,17),expand=c(0,0))
print(p)
unlink(tmp,recursive=T)
setwd(oldwd)

Simon


--------------------------------
Simon O'Hanlon, BSc MSc
Department of Infectious Disease Epidemiology
Imperial College London
St. Mary's Hospital
London
W2 1PG