Maps with `R`

Dhafer Malouche

http://dhafermalouche.net

Outline

Importing the map of Tunisia into R: shp files, raster, ggplot2, googleVis, Leaflet
Statistics and maps
Choropleth maps, plot, ggplot2, spplot, Leaflet

Importing the map of Tunisia into `R`

Using `shp` files

Available shp files on the web:

https://github.com/Groupe-ElementR/FreeCarto/tree/master/SCRIPTS/data

download .prj, .sbn, sbx, .shp, .shp.xml, shx files and put them in the same directory.

Importing `shp` into `R`

> library(maptools)
> library(sp)
> library(shapefiles)
> fdc <- readShapePoly("/Users/dhafermalouche/Documents/Teaching/AdvancedR/geom/Tunisie_snuts4.shp")

Importing `shp` into `R`

> plot(fdc)

Importing `shp` into `R`

Using `raster` package

> library(raster)
> fdc0<- getData(name="GADM", country="TUN", level=0)

Using `raster` package

> plot(fdc0)

Using `raster` package

A level 1 (By gouvernorat) map

> fdc1<- getData(name="GADM", country="TUN", level=1)

A level 1 (By gouvernorat) map

> plot(fdc1)

A level 1 (By gouvernorat) map

A level 2 (By delegation) map

> fdc2<- getData(name="GADM", country="TUN", level=2)

A level 2 (By delegation) map

> plot(fdc2)

A level 2 (By delegation) map

Import data on Gouvernorats into `R`

> library(googleVis)
> library(readr)
> library(rgdal)
> op <- options(gvis.plot.tag="chart")
> tn <- read_delim("/Users/dhafermalouche/Documents/Teaching/AdvancedR/tunisie_data_del_2011.csv", 
+                    ";", escape_double = FALSE, trim_ws = TRUE)
> colnames(tn)
 [1] "del"       "del_nom"   "gou"       "gou_nom"   "reg"      
 [6] "reg_nom"   "zon"       "zon_nom"   "cap_gou"   "POPTO2004"
[11] "POPCO2004" "NBLOG2004" "IMMIG9404" "EMIG9404"  "POPTO2010"
[16] "SUP2010"   "IDRVA2011"

Transforming the map in `html` code (only gouvernorat)

> g1 <- gvisGeoChart(tn, locationvar = "gou",hovervar = "gou_nom",
+                           options=list(region="TN", 
+                                 displayMode="regions",
+                                 resolution="provinces",
+                                 width=800, height=600))
> print(g1, file="g1.html")
> cat(g1$html$chart, file = "g1.html")

Transforming the map in `html` code (only gouvernorat)

With `leaflet` package

> library(leaflet)
> mm<-leaflet(data = fdc2) %>%
+     addProviderTiles("CartoDB.Positron")%>%
+   addPolygons(fillColor = "#444444", 
+               weight = 1, 
+               smoothFactor = 0.5,
+               opacity = 1.0, fillOpacity = 0.5,
+               highlightOptions = highlightOptions(color = "white",
+                                                   weight = 2,
+                                                   bringToFront = TRUE))

With `leaflet` package

Adding pop-ups to the map

> mypopup=fdc2$NAME_2
> mm<-leaflet(data = fdc2) %>%
+     addProviderTiles("CartoDB.Positron")%>%
+   addPolygons(fillColor = "red", 
+               color="black",
+               weight = 1, 
+               smoothFactor = 0.5,
+               popup = mypopup,
+               opacity = 0.5, fillOpacity = 0.5,
+               highlightOptions = highlightOptions(color = "white",
+                                                   weight = 2,
+                                                   bringToFront = TRUE))

Adding pop-ups to the map

With `ggplot2`

> library(ggplot2)
> tun_df <- fortify(fdc,region = "id")
> colnames(tun_df)
[1] "long"  "lat"   "order" "hole"  "piece" "id"    "group"

With `ggplot2`

>  p<-ggplot(tun_df, aes(x=long, y=lat, group=group)) +
+     geom_polygon(aes(fill="pink"),color = "grey50")+
+     labs(x="",y="")+ theme_bw()+
+     coord_fixed()
> p<-p+theme(axis.line=element_blank(),
+      axis.text.x=element_blank(),
+        axis.text.y=element_blank(),
+        axis.ticks=element_blank(),
+        axis.title.x=element_blank(),
+        axis.title.y=element_blank(),
+      panel.grid.major = element_blank(),
+      panel.grid.minor = element_blank(),
+      panel.border = element_blank(),
+      panel.background = element_blank())
>  p<-p+ theme(legend.position="none")

bubbles on the map!

Representing populations by delegation

Importing the data

> data_del <- read.csv("tunisie_data_del_2011.csv", header = TRUE, sep = ";", 
+    dec = ",", encoding = "latin1")

Extract the coordinates of the centers of the delegations in Tunisia from the fdc object of the slide #3

> pt <- cbind(fdc@data[, "id"], as.data.frame(coordinates(fdc)))
> dim(pt)
[1] 263   3

bubbles on the map!

Merge pt (containing the coordinates) and the data containing the variable

> colnames(pt) <- c("id", "x", "y")
> i=match(pt[, "id"], data_del[, "del"])
> pt <- data.frame(pt, data_del[i, ])
> pt$var <- pt$POPTO2010

bubbles on the map!

We determine the maximum extension of the base map, the bbox function gives the max and min coordinates of the base map

> x1 <- bbox(fdc)[1]
> y1 <- bbox(fdc)[2]
> x2 <- bbox(fdc)[3]
> y2 <- bbox(fdc)[4]

bubbles on the map!

The maximum area of the map \(M\)

> M <- (x2 - x1) * (y2 - y1)

Sum of the variable \(s=\sum_i x_i\)

> sc <- sum(pt$var, na.rm = TRUE)
> sc 
[1] 10489116

bubbles on the map!

Computing the rays of the circles

\[r_i=\sqrt{\displaystyle\frac{kx_iM}{s\,\pi}}\] Hence the area of one circle is \(a_i=\pi r_i^2= \pi \displaystyle\frac{kx_iM}{s\,\pi}=\displaystyle\frac{kM}{s}x_i\)

Then \(\sum_i a_i= kM\), we choose \(k=.2\) (20% of the total surface of the map)

> k <- 0.2  
> pt$size <- sqrt((pt$var * k * M/sc)/pi)

bubbles on the map!

We draw then the maps and add the bubbles

> plot(fdc, border = "white", col = "grey")
> symbols(pt[, c("x", "y")], circles = pt$size, add = TRUE, bg = "red", inches = FALSE)

bubbles on the map!

Title of the legend

> LegTitle <- "Popuplation"

Size of the circles according to the quantiles

> rLeg <- quantile(pt$size, c(1, 0.9, 0.25, 0), type = 1, na.rm = TRUE);rLeg
> rVal <- quantile(pt$var, c(1, 0.9, 0.25, 0), type = 1, na.rm = TRUE);rVal
> l <- data.frame(x = x1, y = y1)
> head(l)

bubbles on the map!

The coordinate of the legend

> xinit <- l$x + rLeg[1]
> xinit
> ypos <- l$y + rLeg
> ypos

Adding the circles of the legend

> symbols(x = rep(xinit, 4), y = ypos, circles = rLeg, 
+         add = TRUE, bg = "red", 
+      inches = FALSE)
> 
> text(x = rep(xinit, 4) + rLeg[1] * 1.2, y = (l$y + (2 * rLeg)), rVal, cex = 0.3, 
+      srt = 0, adj = 0)

bubbles on the map!

Adding segments indicating the correspondence size of the circle and the population

> for (i in 1:4) {
+      segments(xinit, (l$y + (2 * rLeg[i])), xinit + rLeg[1] * 1.1, (l$y + (2 * 
+          rLeg[i]))) }
>  text(x = xinit - rLeg[1], y = (l$y + (2 * rLeg[1])), LegTitle, adj = c(0, 0), 
+      cex = 0.7)

bubbles on the map!

Title

> title(main = "Population, 2010",  cex.sub = 0.7)

Scale

> xscale <- x2
> yscale <- y1
> sizescale <- 50000
> labelscale <- "50km"
> SpatialPolygonsRescale(layout.scale.bar(), offset = c(xscale, yscale), scale = sizescale, 
+    fill = c("black"), plot.grid = F)
> text(xscale + sizescale/2, yscale, paste(labelscale, "\n\n", sep = ""), cex = 0.7)

Arrow for the orientation

> xarrow <- x1
> yarrow <- y2 - (y2 - y1)/10
> SpatialPolygonsRescale(layout.north.arrow(2), offset = c(xarrow, yarrow), scale = 50000, 
+      plot.grid = F)

Choropleth maps with `plot`

Matching the names of the delegations between data sets (map data and the one containing the variable)

> i=match(fdc@data[, "id"], data_del[, "del"])
> fdc@data <- data.frame(fdc@data, data_del[i,])
> colnames(fdc@data)
 [1] "id"        "id_snuts3" "id_snuts2" "id_snuts1" "id_snuts0"
 [6] "del"       "del_nom"   "gou"       "gou_nom"   "reg"      
[11] "reg_nom"   "zon"       "zon_nom"   "cap_gou"   "POPTO2004"
[16] "POPCO2004" "NBLOG2004" "IMMIG9404" "EMIG9404"  "POPTO2010"
[21] "SUP2010"   "IDRVA2011"
> summary(fdc@data$IDRVA2011)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.0000  0.1885  0.3120  0.3179  0.4220  1.0000

We will represent the Regional Development Index (IDR)

> summary(fdc@data$IDRVA2011)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.0000  0.1885  0.3120  0.3179  0.4220  1.0000

Choropleth maps with `plot`

Creating clusters from the variable (we choose here 8 clusters)

> nbclass <- 8
> var<-fdc@data$IDRVA2011
> library(classInt)
> distr <- classIntervals(var, nbclass, style = "quantile")$brks
> library(RColorBrewer)
> colours <- brewer.pal(nbclass, "YlOrRd")
> colMap <- colours[(findInterval(fdc@data$IDRVA2011, distr, all.inside = TRUE))]

Choropleth maps with `plot`

drawing the map

> plot(fdc, col = colMap, border = "black", lwd = 1)

Choropleth maps with `plot`

Adding legend

> legend(x = "topright", 
+         legend = leglabs(round(distr, 2),
+                          over = "plus de", under = "moins de"), 
+         fill = colours, bty = "n", pt.cex = 1,
+        cex = 0.7, title = "indice 0-1")
> title(main = "Regional Developpment Index", 
+     cex.sub = 0.7)

Choropleth maps with `plot`

Case study: Comparing Nahdha voters between 2011 and 2014

Importing data into `R`

Import the v csv file

> library(readr)
> data_delegation_tunisia <- read_csv("data_delegation_tunisia.csv")

We need only the variables del (delegation codes), Nahdha_2011 and Nahdha_2014. Let’s extract them!

> dt=data_delegation_tunisia[,c("del","Nahdha_2011","Nahdha_2014")]

Merging the data with `fdc` (see slide 4)

> i=match(fdc@data$id,dt$del)
> fdc@data$Nahdha_2014=dt$Nahdha_2014[i]
> fdc@data$Nahdha_2011=dt$Nahdha_2011[i]

Preparing data

Imputing missing values (to discussed)

> library(e1071)

Attaching package: 'e1071'
The following object is masked from 'package:raster':

    interpolate
> x=matrix(fdc@data$Nahdha_2014,ncol=1)
> x=impute(x)
> fdc@data$Nahdha_2014=x
> y=matrix(fdc@data$Nahdha_2011,ncol=1)
> y=impute(y)
> fdc@data$Nahdha_2011=y

Drawing two maps: 2011 vs 2014

Make it ready for ggplot2: data for polygons

> library(ggplot2)
> tun_df <- fortify(fdc,region = "id")
> dim(tun_df)
[1] 3490    7

On “fortify” data

> library(reshape2)
> dt=fdc@data[,c(1,6,7)]
> colnames(dt)=c("id","2014","2011")
> dt=melt(dt,id.vars = "id")
> colnames(dt)=c("id","Year","Votes")
> tun_df_all=merge(tun_df,dt,by="id")
> colnames(tun_df_all)
[1] "id"    "long"  "lat"   "order" "hole"  "piece" "group" "Year"  "Votes"

Drawing two maps: 2011 vs 2014

drawing the map

> library(RColorBrewer)
>  p<-ggplot(tun_df_all, aes(x=long, y=lat, group=group)) +
+    geom_polygon(aes(fill=Votes),color = "grey50")+
+    scale_fill_gradientn(colours=brewer.pal(5,"OrRd"),name="Votes")+
+    labs(x="",y="")+ theme_bw()+facet_grid(~Year)+
+    coord_fixed()

customizing the map

>  p<-p+theme(axis.line=element_blank(),
+        axis.text.x=element_blank(),
+        axis.text.y=element_blank(),
+        axis.ticks=element_blank(),
+        axis.title.x=element_blank(),
+        axis.title.y=element_blank(),panel.grid.major = element_blank(),
+      panel.grid.minor = element_blank(),
+      panel.border = element_blank(),
+      panel.background = element_blank())
> p<-p+ theme(legend.position="right")
> p

Drawing two maps: 2011 vs 2014

Displaying the `Difference` variable

Method 1: using spplot package from lattice package (All delegations)

> library(lattice)
> fdc@data$Difference=fdc@data$Nahdha_2014-fdc@data$Nahdha_2011
> spplot(fdc,"Difference",
+        main="Difference of number of voters\n 2014 and 2011",sub="")

Displaying the `Difference` variable

We will display the delegations where Nahdha party gained voters (positive)

> i=which(fdc@data$Difference>0)
> spplot(fdc[i,],"Difference",main="Difference of number of voters\n 2014 and 2011",sub="")

Using `Leaflet` package

We consider fdc2 object import in the slide 13

> i=match(fdc2$HASC_2,data_delegation_tunisia$HASC_2)
> fdc2a=cbind.Spatial(fdc2,data_delegation_tunisia[i,])

We will display the rate of the unemployment in 2004 (according to INS)

> xv=fdc2a@data$taux.de.chomage.en.2004
> summary(xv)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
   4.50   10.60   14.50   15.34   19.10   39.00      15

Using `Leaflet` package

Palette of colors

> library(leaflet)
> MyPaletteColor <- colorQuantile("YlGn", NULL, 7)

Pop ups

> my_popup <- paste0("<strong>",fdc2a@data$delegation,"</strong>",xv,"%")

The code leaflet

> mm<-leaflet(data = fdc2a) %>%
+     addProviderTiles("CartoDB.Positron") %>%
+     addPolygons(fillColor = ~MyPaletteColor(xv), 
+                 fillOpacity = 0.8, 
+                 color = "#BDBDC3", 
+                 weight = 1, 
+                 popup = my_popup)
> mm

Maps with R

Outline

Importing the map of Tunisia into R

Using shp files

Importing shp into R

Importing shp into R

Importing shp into R

Using raster package

Using raster package

Using raster package

A level 1 (By gouvernorat) map

A level 1 (By gouvernorat) map

A level 1 (By gouvernorat) map

A level 2 (By delegation) map

A level 2 (By delegation) map

A level 2 (By delegation) map

With googleVis

Import data on Gouvernorats into R

Transforming the map in html code (only gouvernorat)

Transforming the map in html code (only gouvernorat)

With leaflet package

With leaflet package

Adding pop-ups to the map

Adding pop-ups to the map

With ggplot2

With ggplot2

With ggplot2

Statistics and maps

bubbles on the map!

bubbles on the map!

bubbles on the map!

bubbles on the map!

bubbles on the map!

bubbles on the map!

bubbles on the map!

bubbles on the map!

bubbles on the map!

bubbles on the map!

bubbles on the map!

bubbles on the map!

bubbles on the map!

bubbles on the map!

Choropleth maps

Choropleth maps with plot

Choropleth maps with plot

Choropleth maps with plot

Choropleth maps with plot

Choropleth maps with plot

Choropleth maps with plot

Case study: Comparing Nahdha voters between 2011 and 2014

Importing data into R

Merging the data with fdc (see slide 4)

Preparing data

Drawing two maps: 2011 vs 2014

Drawing two maps: 2011 vs 2014

Drawing two maps: 2011 vs 2014

Displaying the Difference variable

Displaying the Difference variable

Using Leaflet package

Using Leaflet package

Unemployment rate in 2004

Maps with `R`

Importing the map of Tunisia into `R`

Using `shp` files

Importing `shp` into `R`

Importing `shp` into `R`

Importing `shp` into `R`

Using `raster` package

Using `raster` package

Using `raster` package

With `googleVis`

Import data on Gouvernorats into `R`

Transforming the map in `html` code (only gouvernorat)

Transforming the map in `html` code (only gouvernorat)

With `leaflet` package

With `leaflet` package

With `ggplot2`

With `ggplot2`

With `ggplot2`

Choropleth maps with `plot`

Choropleth maps with `plot`

Choropleth maps with `plot`

Choropleth maps with `plot`

Choropleth maps with `plot`

Choropleth maps with `plot`

Importing data into `R`

Merging the data with `fdc` (see slide 4)

Displaying the `Difference` variable

Displaying the `Difference` variable

Using `Leaflet` package

Using `Leaflet` package