Mapping stuff

Refs:

• Maps in R: Introduction - Drawing the map of Europe

• Maps in R: Plotting data points on a map

• Maps in R: choropleth maps

# download from http://cran.r-project.org/web/packages/rworldmap/index.html
# and high def map from http://cran.r-project.org/web/packages/rworldxtra/index.html
library(rworldmap)

## Loading required package: sp
## ### Welcome to rworldmap ###
## For a short introduction type :   vignette('rworldmap')

newmap <- getMap(resolution = "low")
plot(newmap)

Defining map limits

library(ggmap)

## Loading required package: ggplot2

europe.limits <- geocode(c(
  "CapeFligely,RudolfIsland,Franz Josef Land,Russia",
  "Gavdos,Greece",
  "Faja Grande,Azores",
  "SevernyIsland,Novaya Zemlya,Russia")
)

## Information from URL : http://maps.googleapis.com/maps/api/geocode/json?address=CapeFligely,RudolfIsland,Franz+Josef+Land,Russia&sensor=false
## Information from URL : http://maps.googleapis.com/maps/api/geocode/json?address=Gavdos,Greece&sensor=false
## Information from URL : http://maps.googleapis.com/maps/api/geocode/json?address=Faja+Grande,Azores&sensor=false
## Information from URL : http://maps.googleapis.com/maps/api/geocode/json?address=SevernyIsland,Novaya+Zemlya,Russia&sensor=false

europe.limits

##         lon      lat
## 1  60.64878 80.58823
## 2  24.08464 34.83469
## 3 -31.26192 39.45479
## 4  56.00000 74.00000

plot(newmap,
  xlim = range(europe.limits$lon),
  ylim = range(europe.limits$lat),
  asp = 1.0
)

Adding information on map

Get some data, in this case a dataset including airport coordinates:

# from http://openflights.org/data.html
airports <- read.csv("airports.dat", header = FALSE)
colnames(airports) <- c("ID", "name", "city", "country", "IATA_FAA", "ICAO", "lat", "lon", "altitude", "timezone", "DST")
head(airports)

##   ID                       name         city          country IATA_FAA
## 1  1                     Goroka       Goroka Papua New Guinea      GKA
## 2  2                     Madang       Madang Papua New Guinea      MAG
## 3  3                Mount Hagen  Mount Hagen Papua New Guinea      HGU
## 4  4                     Nadzab       Nadzab Papua New Guinea      LAE
## 5  5 Port Moresby Jacksons Intl Port Moresby Papua New Guinea      POM
## 6  6                 Wewak Intl        Wewak Papua New Guinea      WWK
##   ICAO       lat      lon altitude timezone DST                   NA
## 1 AYGA -6.081689 145.3919     5282       10   U Pacific/Port_Moresby
## 2 AYMD -5.207083 145.7887       20       10   U Pacific/Port_Moresby
## 3 AYMH -5.826789 144.2959     5388       10   U Pacific/Port_Moresby
## 4 AYNZ -6.569828 146.7262      239       10   U Pacific/Port_Moresby
## 5 AYPY -9.443383 147.2200      146       10   U Pacific/Port_Moresby
## 6 AYWK -3.583828 143.6692       19       10   U Pacific/Port_Moresby

Now place them in the map

plot(newmap, xlim = range(europe.limits$lon), ylim = range(europe.limits$lat), asp = 1.0)
points(airports$lon, airports$lat, col = "red", cex = .25)

Using ggmap

Get a map from Google Maps (there’s also other servers)

library(ggmap)
library(mapproj)

## Loading required package: maps

map <- get_map(location = 'Europe', zoom = 4)

## Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=Europe&zoom=4&size=640x640&scale=2&maptype=terrain&language=en-EN&sensor=false
## Information from URL : http://maps.googleapis.com/maps/api/geocode/json?address=Europe&sensor=false

ggmap(map)

Get arrival and departure info to attach to the airpot data:

# Also from http://openflights.org/data.html
routes <- read.csv("routes.dat", header=F)
colnames(routes) <- c("airline", "airlineID", "sourceAirport", "sourceAirportID", "destinationAirport", "destinationAirportID", "codeshare", "stops", "equipment")
head(routes)

##   airline airlineID sourceAirport sourceAirportID destinationAirport
## 1      2B       410           AER            2965                KZN
## 2      2B       410           ASF            2966                KZN
## 3      2B       410           ASF            2966                MRV
## 4      2B       410           CEK            2968                KZN
## 5      2B       410           CEK            2968                OVB
## 6      2B       410           DME            4029                KZN
##   destinationAirportID codeshare stops equipment
## 1                 2990               0       CR2
## 2                 2990               0       CR2
## 3                 2962               0       CR2
## 4                 2990               0       CR2
## 5                 4078               0       CR2
## 6                 2990               0       CR2

Compute how many arrivals and departures from each airport (the joining key is ID)

library(plyr)
departures <- ddply(routes, .(sourceAirportID), "nrow")
names(departures)[2] <- "flights"
head(departures)

##   sourceAirportID flights
## 1             \\N     239
## 2               1       5
## 3              10       2
## 4             100      45
## 5            1001       4
## 6            1004       7

arrivals <- ddply(routes, .(destinationAirportID), "nrow")
names(arrivals)[2] <- "flights"
head(arrivals)

##   destinationAirportID flights
## 1                  \\N     242
## 2                    1       5
## 3                   10       2
## 4                  100      44
## 5                 1001       4
## 6                 1004       7

airportD <- merge(airports, departures, by.x = "ID", by.y = "sourceAirportID")
head(airportD)

##   ID                       name         city          country IATA_FAA
## 1  1                     Goroka       Goroka Papua New Guinea      GKA
## 2  2                     Madang       Madang Papua New Guinea      MAG
## 3  3                Mount Hagen  Mount Hagen Papua New Guinea      HGU
## 4  4                     Nadzab       Nadzab Papua New Guinea      LAE
## 5  5 Port Moresby Jacksons Intl Port Moresby Papua New Guinea      POM
## 6  6                 Wewak Intl        Wewak Papua New Guinea      WWK
##   ICAO       lat      lon altitude timezone DST                   NA
## 1 AYGA -6.081689 145.3919     5282       10   U Pacific/Port_Moresby
## 2 AYMD -5.207083 145.7887       20       10   U Pacific/Port_Moresby
## 3 AYMH -5.826789 144.2959     5388       10   U Pacific/Port_Moresby
## 4 AYNZ -6.569828 146.7262      239       10   U Pacific/Port_Moresby
## 5 AYPY -9.443383 147.2200      146       10   U Pacific/Port_Moresby
## 6 AYWK -3.583828 143.6692       19       10   U Pacific/Port_Moresby
##   flights
## 1       5
## 2       8
## 3      10
## 4      11
## 5      52
## 6       6

airportA <- merge(airports, arrivals,   by.x = "ID", by.y = "destinationAirportID")

A ggmap is an object of class ggplot, so typical ggplot operations work with it

mapPoints <- ggmap(map) +
             geom_point(aes(x = lon, y = lat, size = sqrt(flights)), data = airportD, alpha = .5) +
             scale_size_area(breaks = sqrt(c(1, 5, 10, 50, 100, 500)), 
                             labels = c(1, 5, 10, 50, 100, 500), name = "departing routes")
mapPoints

## Warning in loop_apply(n, do.ply): Removed 2727 rows containing missing
## values (geom_point).

Using shapefiles with maptools

library(maptools) # read shapefiles (.shp)

## Checking rgeos availability: TRUE

library(ggplot2)
library(ggmap)

# from http://ec.europa.eu/eurostat/web/gisco/geodata/reference-data/administrative-units-statistical-units
# read administrative boundaries (change folder appropriately)
eurMap <- readShapePoly(fn="NUTS_2010_60M_SH/data/NUTS_RG_60M_2010")
plot(eurMap)

# convert shp data to data frame for ggplot
eurMapDf <- fortify(eurMap, region='NUTS_ID')

# read downloaded data (change folder appropriately)
eurEdu <- read.csv("educ_thexp_1_Data.csv", stringsAsFactors = F)
eurEdu$Value <- as.double(eurEdu$Value) #format as numeric

## Warning: NAs introduced by coercion

# merge map and data
eurEduMapDf <- merge(eurMapDf, eurEdu, by.x="id", by.y="GEO")
eurEduMapDf <- eurEduMapDf[order(eurEduMapDf$order),]

eurEduMapDf <- subset(eurEduMapDf, 
                      long > min(europe.limits$lon) & long < max(europe.limits$lon) & 
                      lat  > min(europe.limits$lat) & lat  < max(europe.limits$lat))

# ggplot mapping
# data layer
m0 <- ggplot(data=eurEduMapDf)
# empty map (only borders)
m1 <- m0 + geom_path(aes(x=long, y=lat, group=group), color='gray') + coord_equal()

# fill with education expenditure data
m2 <- m1 + geom_polygon(aes(x=long, y=lat, group=group, fill=Value))

# inverse order (to have visible borders)
m0 <- ggplot(data=eurEduMapDf)
m1 <- m0 + geom_polygon(aes(x=long, y=lat, group=group, fill=Value)) + coord_equal()
m2 <- m1 + geom_path(aes(x=long, y=lat, group=group), color='black')
m2

# over a GoogleMap (not working if not correctly projected)
map <- get_map(location = 'Europe', zoom=4)

## Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=Europe&zoom=4&size=640x640&scale=2&maptype=terrain&language=en-EN&sensor=false

m0 <- ggmap(map)
m1 <- m0 + geom_polygon(aes(x=long, y=lat, group=group, fill=Value), data=eurEduMapDf, alpha=.9)
m2 <- m1 + geom_path(aes(x=long, y=lat, group=group), data=eurEduMapDf, color='black')

# add text
library(doBy) # from http://cran.r-project.org/web/packages/doBy/index.html

## Loading required package: survival

txtVal <- summaryBy(long + lat + Value ~ id, data=eurEduMapDf, FUN=mean, keep.names=T)
m3 <- m2 + geom_text(aes(x=long, y=lat, label=Value), data=txtVal, col="yellow", cex=3)
m3

## Warning in loop_apply(n, do.ply): Removed 1654 rows containing missing
## values (geom_path).

## Warning in loop_apply(n, do.ply): Removed 6 rows containing missing values
## (geom_text).

Maps and functions for Portugal

PT_limits <- geocode(c(    # from http://pt.wikipedia.org/wiki/Pontos_extremos_de_Portugal
  "Cristoval,Portugal",                     # mais setentrional
  "Paradela,Portugal",                      # mais oriental
  "Cabo da Roca, Sintra,Portugal",          # mais ocidental
  "Cabo de Santa Maria,Algarve,Portugal")   # mais meridional
)

## Information from URL : http://maps.googleapis.com/maps/api/geocode/json?address=Cristoval,Portugal&sensor=false
## Information from URL : http://maps.googleapis.com/maps/api/geocode/json?address=Paradela,Portugal&sensor=false
## Information from URL : http://maps.googleapis.com/maps/api/geocode/json?address=Cabo+da+Roca,+Sintra,Portugal&sensor=false
## Information from URL : http://maps.googleapis.com/maps/api/geocode/json?address=Cabo+de+Santa+Maria,Algarve,Portugal&sensor=false

plot(eurMap, xlim = range(PT_limits$lon), ylim = range(PT_limits$lat), asp = 1.0)

# pt_map <- get_map(location = 'Portugal', zoom = 7, language = "PT-pt")
# this is a better fit:
pt_map <- get_map(location = c(lon=mean(PT_limits$lon), lat= mean(PT_limits$lat)), zoom = 7, language = "PT-pt", maptype="satellite")

## Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=39.643885,-8.489004&zoom=7&size=640x640&scale=2&maptype=satellite&language=PT-pt&sensor=false

ggmap(pt_map)

# http://www.dgterritorio.pt/cartografia_e_geodesia/cartografia/carta_administrativa_oficial_de_portugal__caop_/caop_em_vigor/

freguesias <- read.csv("freguesias.csv", header = TRUE)
head(freguesias)

##   FREGUESIA_COD  NUTSI_DSG NUTSII_DSG NUTSIII_DSG DISTRITO_ILHA_DSG
## 1        010103 CONTINENTE     CENTRO BAIXO VOUGA            AVEIRO
## 2        010109 CONTINENTE     CENTRO BAIXO VOUGA            AVEIRO
## 3        010112 CONTINENTE     CENTRO BAIXO VOUGA            AVEIRO
## 4        010119 CONTINENTE     CENTRO BAIXO VOUGA            AVEIRO
## 5        010121 CONTINENTE     CENTRO BAIXO VOUGA            AVEIRO
## 6        010122 CONTINENTE     CENTRO BAIXO VOUGA            AVEIRO
##   MUNICIPIO_DSG                                   FREGUESIA_DSG
## 1        ÁGUEDA                                  AGUADA DE CIMA
## 2        ÁGUEDA                                     FERMENTELOS
## 3        ÁGUEDA                              MACINHATA DO VOUGA
## 4        ÁGUEDA                                VALONGO DO VOUGA
## 5        ÁGUEDA       UNIÃO DAS FREGUESIAS DE ÁGUEDA E BORRALHA
## 6        ÁGUEDA UNIÃO DAS FREGUESIAS DE BARRÔ E AGUADA DE BAIXO
##   AREA_2014.Ha.
## 1       2839.31
## 2        858.20
## 3       3195.44
## 4       4320.11
## 5       3602.93
## 6       1019.01

municipios <- read.csv("municipios.csv", header = TRUE)
head(municipios)

##   MUNICIPIO_COD.NUTSI_DSG.NUTSII_DSG.NUTSIII_DSG.DISTRITO_ILHA_DSG.MUNICIPIO_DSG.AREA_2014.Ha..ALTITUDE_MAX.M..ALTITUDE_MIN.M.
## 1                                                              0101;CONTINENTE;CENTRO;BAIXO VOUGA;AVEIRO;ÁGUEDA;33527.44;758;4
## 2                                                   0102;CONTINENTE;CENTRO;BAIXO VOUGA;AVEIRO;ALBERGARIA-A-VELHA;15882.5;425;0
## 3                                                             0103;CONTINENTE;CENTRO;BAIXO VOUGA;AVEIRO;ANADIA;21663.48;525;13
## 4                                                     0104;CONTINENTE;NORTE;ENTRE DOURO E VOUGA;AVEIRO;AROUCA;32910.52;1220;50
## 5                                                               0105;CONTINENTE;CENTRO;BAIXO VOUGA;AVEIRO;AVEIRO;19757.57;78;0
## 6                                                         0106;CONTINENTE;NORTE;TÂMEGA;AVEIRO;CASTELO DE PAIVA;11500.53;690;25

Functions for mapping dstrict and municipe information of Portugal

The maps from Portugal were from Global Administrative Areas datafiles. Zip file with shapefiles from Portugal

Ref: http://statisticaconr.blogspot.it/2009/12/choropleth-map-in-r-coloriamo-le-mappe.html (In Italian)

library(maptools) # read shapefiles (.shp)

# district_data is a data frame with two columns, district and value
pt_distritos <- function(district_data, breaks, palette = "hot", 
                         show_label=TRUE, 
                         # the next args define how to handle non-available datapoints
                         show.NA=FALSE, val.NA=-1, col.NA="black", label.NA="NA") {
  
  # read shapefile file for districts
  library(maptools)
  ptMap <- readShapePoly(fn="PRT_adm1")
  
  init <- ifelse(show.NA,val.NA,min(municipe_data[,2]))

  df <- data.frame(districts = ptMap$NAME_1,
                   value     = rep(init,length(ptMap$NAME_1)))

  # populate df with district_data (the districts might be in different sorting order)
  for(i in 1:nrow(df)) {
    query <- as.character(district_data[,1])==as.character(df[i,1])
    if (any(query))
      df[i,2] = district_data[which(query),2]
  }
  
  if (any(df==-1))
    breaks = c(0,breaks)
  
  library(Hmisc)
  lev <- cut2(df$value, cuts=breaks)
  n.cols <- length(breaks)+1

  if(show_label && length(breaks)>3) { # does not work with small number of breaks (?)
    label_names <- levels(lev)
    if (show.NA && any(df==-1))
      label_names[1] <- "NA"
  } else
    label_names <- ""
  
  if (palette=="hot") {
    cols <- heat.colors(n.cols)
    if (show.NA && any(df==-1))
      cols <- c(col.NA,cols)
  }
  else if (palette == "grey")
    cols <- grey(seq(0, 1, length = n.cols))
  else if (palette == "rainbow") {
    cols <- rainbow(n.cols)
    if (show.NA && any(df==-1))
      cols <- c(col.NA,cols)
  }
 
  ptMap$value <- as.factor( as.character(as.numeric(lev)) )
  
  library(sp)
  spplot(ptMap, "value", col.regions=cols, 
         xlim = range(-10,-6), ylim = range(36.9,42.2), asp = 1.0,
         colorkey = list(labels = list(labels = label_names, width = 1, cex = 1),
                         space = "bottom"))
}

# Use eg:
library(maptools)
ptMap <- readShapePoly(fn="PRT_adm1")
district_data <- data.frame(district = ptMap$NAME_1, 
                            value    = rpois(20,10)) # make some fake data
district_data <- district_data[1:19,] # remove one district, for testing

pt_distritos(district_data, breaks=c(3,6,11,15), palette="rainbow", show.NA=TRUE)

## Loading required package: grid
## Loading required package: lattice
## Loading required package: Formula
## 
## Attaching package: 'Hmisc'
## 
## The following object is masked from 'package:maptools':
## 
##     label
## 
## The following objects are masked from 'package:plyr':
## 
##     is.discrete, summarize
## 
## The following objects are masked from 'package:base':
## 
##     format.pval, round.POSIXt, trunc.POSIXt, units

And the same for the municipal areas:

# district_data is a data frame with two columns, municipe and value
pt_municipios <- function(municipe_data, breaks, palette = "hot", 
                          show_label=TRUE, 
                          # the next args define how to handle non-available datapoints
                          show.NA=FALSE, val.NA=-1, col.NA="black", label.NA="NA") {
  
  # read shapefile file for districts
  library(maptools)
  ptMap <- readShapePoly(fn="PRT_adm2")
  
  init <- ifelse(show.NA,val.NA,min(municipe_data[,2]))
  
  df <- data.frame(municipios = ptMap$NAME_2,
                   value      = rep(init,length(ptMap$NAME_2)))

  # populate df with municipe_Data (the municipes might be in different sorting order)
  for(i in 1:nrow(df)) {
    query <- as.character(municipe_data[,1])==as.character(df[i,1])
    if (any(query))
      df[i,2] = municipe_data[which(query),2]
  }
  
  if (any(df==-1))
    breaks = c(0,breaks)
  
  library(Hmisc)
  lev <- cut2(df$value, cuts=breaks)
  n.cols <- length(breaks)+1

  if(show_label && length(breaks)>3) { # does not work with small number of breaks (?)
    label_names <- levels(lev)
    if (show.NA && any(df==-1))
      label_names[1] <- "NA"
  } else
    label_names <- ""
  
  if (palette=="hot") {
    cols <- heat.colors(n.cols)
    if (show.NA && any(df==-1))
      cols <- c(col.NA,cols)
  }
  else if (palette == "grey")
    cols <- grey(seq(0, 1, length = n.cols))
  else if (palette == "rainbow") {
    cols <- rainbow(n.cols)
    if (show.NA && any(df==-1))
      cols <- c(col.NA,cols)
  }
  
  ptMap$value <- as.factor( as.character(as.numeric(lev)) )
  
  library(sp)
  spplot(ptMap, "value", col.regions=cols, 
         xlim = range(-10,-6), ylim = range(36.9,42.2), asp = 1.0,
         colorkey = list(labels = list(labels = label_names, width = 1, cex = 1),
                         space = "bottom"))
}

# Use eg:
municipe_data <- read.csv("passivos.csv", header = TRUE)
head(municipe_data)

##           Municipio Passivo
## 1 Arcos de Valdevez    0.00
## 2           Caminha 2281.97
## 3           Melgaço 3257.48
## 4            Monção  400.00
## 5  Paredes de Coura 3081.12
## 6    Ponte da Barca 1291.75

pt_municipios(municipe_data, breaks=c(1e2,1e3,5e3,1e4,2e4), palette="hot", show.NA=TRUE)

Also check https://procomun.wordpress.com/2012/02/18/maps_with_r_1/ and https://procomun.wordpress.com/2012/02/20/maps_with_r_2/.