Mapping stuff
May, 2015
Refs:
# 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: ggplot2europe.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=falseeurope.limits## lon lat
## 1 60.64878 80.58823
## 2 24.08464 34.83469
## 3 -31.26192 39.45479
## 4 56.00000 74.00000plot(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_MoresbyNow 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: mapsmap <- 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=falseggmap(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 CR2Compute 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 7arrivals <- 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 7airportD <- 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 6airportA <- 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: TRUElibrary(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=falsem0 <- 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: survivaltxtVal <- 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=falseplot(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=falseggmap(pt_map)