Maps with R, and polygon boundaries

[This article was first published on Freakonometrics - Tag - R-english, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
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).

https://i2.wp.com/freakonometrics.blog.free.fr/public/perso4/plate-tekto.gif?w=578

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.

To leave a comment for the author, please follow the link and comment on their blog: Freakonometrics - Tag - R-english.

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.

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)