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.
OpenEdition suggests that you cite this post as follows:
Arthur Charpentier (December 21, 2011). Maps with R, and polygon boundaries. Freakonometrics. Retrieved October 5, 2024 from https://doi.org/10.58079/oujo
baseFREQ=read.table(
“http://freakonometrics.free.fr/baseFREQ.csv”,
sep=”;”,header=TRUE,encoding=”latin1″)
(based on http://freakonometrics.hypotheses.org/978)
but the data does not include used above columns like baseFREQ$region … (instead:
[1] “numeropol” “debut_pol” “fin_pol” “freq_paiement”
[5] “langue” “type_prof” “alimentation” “type_territoire”
[9] “utilisation” “presence_alarme” “marque_voiture” “sexe”
[13] “cout1” “cout2” “cout3” “cout4”
[17] “nbsin” “exposition” “cout” “age”
[21] “duree_permis” “annee_vehicule”),
so the example does not look being replicable…
Sorry if I went (too) quickly on that one. The database used can be, e.g., the following,
(actually, it is an extraction of the one here to get that graph, and that is used in a forthcoming book on R for actuarial science… but I’ll post more on that soon).
in the first example, the baseFREQ object is not defined in that code. It would be very helpful for those (me included) who’re not very well versed on R mapping. Thannk you in advance!