# Maps with R, and polygon boundaries

December 21, 2011
By

(This article was first published on Freakonometrics - Tag - R-english, and kindly contributed to R-bloggers)

With R, it is extremely easy to draw maps. Let us start with something simple, like French regions. Baptiste mentioned on his blog that shapefiles can be downloaded from http://ign.fr/ website. Hence, if you extract the zip file, it is possible to get claims frequency per region (as done in the course ACT2040),

> library(maptools)
> library(maps)
> region<-tapply(baseFREQ[,"nbre"],
+ as.factor(baseFREQ[,"region"]),sum)/
+ tapply(baseFREQ[,"exposition"],
+ as.factor(baseFREQ[,"region"]),sum)
> depFREQ=rep(NA,nrow(departements))
> names(depFREQ)=as.character(
+ departements$CODE_REG) > for(nom in names(region)){ + depFREQ[names(depFREQ)==nom] = + region[nom]} > plot(departements,col=gray((depFREQ-.05)*20)) > legend(166963,6561753,legend=seq(1,0,by=-.1)/20+.05, + fill=gray(seq(1,0,by=-.1)),cex=1.25, bty="n") Another application is on earthquakes. It is possible to use shapefiles of tectonic plates contour, and to relate earthquakes to plates. Shapefiles can be found on http://www.colorado.edu/ (here). First, we can extract the shapes of the tectonic plates > plates = readShapePoly("plates.shp", + proj4string=CRS("+proj=longlat")) > PP=SpatialPolygons2PolySet(plates) Consider Montreal, > montreal=c(-73.600,45.500) Given that specific location, it is possible to use the following code to get the associated plate, > PLATE.loc=function(pt){ + K=NA + for(k in 1:17){ + c=point.in.polygon(pt[1], pt[2], + PP[PP$PID==k,c("X")],PP[PP$PID==k,c("Y")], + mode.checked=FALSE) + if(c>0){K=k} + } + return(K)} > abline(v=montreal[1],col="red") > abline(h=montreal[2],col="red") > PLATE.loc(montreal) [1] 1 and then to plot the associated tectonic plate very easily > PLATE=function(k0){ + library(maps) + map("world") + polygon(PP[PP$PID==k0,c("X")],PP[PP$PID==k0,c("Y")], + col="red") + for(k in (1:17)[-k0]){polygon(PP[PP$PID==k,c("X")],
+ PP[PP\$PID==k,c("Y")],col="light blue")}
> PLATE(PLATE.loc(montreal))

Those code were used in the paper written with Mathieu, and that will be presented on January 30th at the Geotop seminar.