Maps with R, and polygon boundaries
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
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 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")} + 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.
R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.