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) > departements<-readShapeSpatial("DEPARTEMENT.SHP") > 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 onhttp://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")} + map("world",add=TRUE)} > 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.