How to create multilayered maps in Stata and R from shapefiles

During the last decade several packages in both Stata and R have been developed in order to allow users to manage geographical data. In this post I show how to create maps by superimposing different shapefiles with our favorite statistical software. In this example a map of Lebanon is created. Districts (Cazas), roads and terrorist attacks are the three layers used in this mapping exercise.

Stata. To the best of my knowledge the most complete suite of commands for spatial data analysis in Stata is the one from Pisati (2001). Since Stata cannot directly deal with shapefiles, in this exercise the command shp2dta (Crow 2013) is used to convert .shp and .dbf files into Stata’s .dta. Subsequently the spmap command (Pisati 2008) is used to aggregate layers and plot the map.

* Title: Mapping Lebanon with Stata
* Author: Elio Amicarelli
* Software version: Stata/SE 12.0
* Date: 12 December 2013

clear all 
cd C:\Users\Utentew\Desktop\Lebanon\Stata 

* Install relevant packages 
ssc install spmap 
ssc install shp2dta 

* Conversion from .shp to .dta 
shp2dta using LBN_adm2, database(lebdb) coordinates(lebcd) genid(id)
/* "LBN_adm" is my administrative boundaries shapefile, it enters the shp2dta without the suffix. Two .dta files will be created from this shapefile: a database file with name lebdb and a coordinate file I call lebcd. Below the same procedure is applied to other shapes I want to add in my map. */ 
shp2dta using LBN_roads, database(lebroadb) coordinates(lebroacd) genid(id) 
shp2dta using attacks, database(lebterrdb) coordinates(lebterrcd) genid(id) 

* Put layers together on the same map 
use "lebterrcd.dta", clear 
generate _id = _n 
spmap using "lebcd.dta", id(_id) fcolor(white) /// 
point(x(_X) y(_Y) size(*1.2) fcolor(red) ocolor(black) osize(*0.5) /// 
 legenda(on) leglabel(Attack)) ///
line(data("lebroacd.dta")color(orange) legenda(on) leglabel(roads)) /// 
title("Terrorist attacks", size(*0.8)) /// 
subtitle("Lebanon (1970-2011)" " ", size(*0.8)) /// 
legend(region(lcolor(black) fcolor(white)) /// 
 title("Legend", size(*0.4) bexpand justification(center)) position (4))


Stata output - click to enlarge


R. The R library offers a wide choice of cartographic packages including maptools, rgdal, mapproj, Rmap and RArcInfo. For this example I decided to use rgdal (Bivand et al. 2013). Unlike Stata, R can directly interface with the shapefile format so we do not need to make any conversion.

# Title: Mapping Lebanon with R
# Author: Elio Amicarelli
# Software version: R 3.0.1
# Date: 12 December 2013

# Install relevant packages
library(rgdal)

# Import the shapefiles
counties.rg <- readOGR("C:/Users/Utentew/Desktop/Lebanon/R","LBN_adm2")
attacks.rg <- readOGR("C:/Users/Utentew/Desktop/Lebanon/R", "attacks")
roads.rg <- readOGR("C:/Users/Utentew/Desktop/Lebanon/R", "LBN_roads")

# Generate a simple map showing all three layers
plot(counties.rg, axes=TRUE, border="black", bg="white")
lines(roads.rg, col="orange", lwd=1.0)
points(attacks.rg, col="black", bg="red", pch=21, cex =1.0)
title(main="R - Terrorist attacks",cex.main = 1, font.main=1)
mtext("Lebanon (1970-2011)", cex= 0.80, font=1)
legend("bottomright",inset=c(0.1,0.1),c("Border","Road","Attack"),cex=0.6, lty=c(0,1,0),pch=c(22,26,16),col=c("black","orange","red"),pt.cex=1)


R output - click to enlarge

It is particularly satisfiyng to create maps with professional statistical software, isn't it? Obviously this is just a very simple example of the spatial data analysis capacity developed in Stata and R. Indeed I will come back on this subject, in the meanwhile enjoy maps!

Data:

Districs and roads shapefiles are from http://diva-gis.org/

Terrorist attacks data are from the Global Terrorism Database http://www.start.umd.edu/gtd/

References:

Crow K. 2013. SHP2DTA: Stata module to converts shape boundary files to Stata datasets, http://EconPapers.repec.org/RePEc:boc:bocode:s456718

Bivand R. et al. 2013. Package 'rgdal': Bindings for the Geospatial Data Abstracion Library, http://cran.r-project.org/web/packages/rgdal/rgdal.pdf 

Pisati M. 2001. sg162: Tools for spatial data analysis. Stata Technical Bulletin 60: 21{37. In Stata Technical Bulletin Reprints, vol. 10, 277{298. College Station, TX: Stata Press, http://www.stata.com/products/stb/journals/stb60.pdf

Pisati M. 2008. SPMAP: Stata module to visualize spatial data, http://EconPapers.repec.org/RePEc:boc:bocode:s456812

EA

No comments :

Post a Comment