loadLibs = function(libs="RgoogleMaps"){
  for (l in libs) 
    suppressMessages(suppressWarnings(require(l,character.only =TRUE, quietly=TRUE)))
}
if (!require(RgoogleMaps))  install.packages("RgoogleMaps", repos="http://R-Forge.R-project.org")

pckgsForExercises = c("geosphere", "ggmap", "RgoogleMaps", "png")
pckgs2Compile = c("sp", "rgdal", "leaflet", "loa", "latticeExtra")

pckgs2load = c(pckgsForExercises,pckgs2Compile)

for (pkg in pckgs2load) library(pkg,character.only =TRUE) 


OFFLINE=FALSE # use previously downloaded maps/data
set.seed(123)

if (!dir.exists("data")) dir.create("data")
dataDir = "H:/DropboxHWR/InternationalVisits/Pune2017/Teaching/MapsInR/ForStudents/data"

data(meuse)
LevCols = rev(heat.colors(5))[cut(meuse$lead,breaks=c(0,100,200,300,400,Inf),labels=FALSE)]

Spatial Data in R

The motivation to want to draw on map tiles directly within R is twofold. While the R environment boasts a long history of spatial analysis tools and libraries, the plots have traditionally been created without any spatial context that a map would provide.

The meuse data set is the “Hello world” equivalent of the sp package. I quote:

It gives locations and topsoil heavy metal concentrations, along with a number of soil and landscape variables at the observation locations, collected in a flood plain of the river Meuse, near the village of Stein (NL). Heavy metal concentrations are from composite samples of an area of approximately 15 m x 15m.

Note the “strange” coordinate system

head(meuse[,1:8])
##        x      y cadmium copper lead zinc  elev       dist
## 1 181072 333611    11.7     85  299 1022 7.909 0.00135803
## 2 181025 333558     8.6     81  277 1141 6.983 0.01222430
## 3 181165 333537     6.5     68  199  640 7.800 0.10302900
## 4 181298 333484     2.6     81  116  257 7.655 0.19009400
## 5 181307 333330     2.8     48  117  269 7.480 0.27709000
## 6 181390 333260     3.0     61  137  281 7.791 0.36406700
coordinates(meuse) <- c("x", "y")
plot(meuse, bg=LevCols, col="grey",main="Lead Distribution", pch=21, axes=TRUE)

There are times when such a plot will be the preferred choice: a clean and simple graph of the locations of interest with no clutter and no distractions. The analyst can focus on the pattern itself and the marker attributes.

However, a lot of the modern massive data sets gathered such as location information from mobile devices, surveys, crime data, vehicle tracks, demographic data, etc. require a map based spatial context for even the most basic data explorations. In those settings, the somewhat narrow or ``blind’’ exploration of spatial data on a blank background can be rather limiting and often leads to less insight than would be possible had the data been graphed on a map canvas. Anecdotally, we dare to speculate that the British phyisican John Snow might have been slower to identify contaminated drinking water as the cause of the cholera epidemic in 1854, had he not used a dot map to illustrate the cluster of cholera cases around one of the pumps in London.

Google/Bing/OSM maps and tiles

Although the Mercator projection significantly distorts scale and area (particularly near the poles), it has two important properties that outweigh the scale distortion:

  1. It’s a conformal projection, which means that it preserves the shape of relatively small objects. This is especially important when showing aerial imagery, because we want to avoid distorting the shape of buildings. Square buildings should appear square, not rectangular.

  2. It’s a cylindrical projection, which means that north and south are always straight up and down, and west and east are always straight left and right.

At the lowest level of detail (Level 1), the map is 512 x 512 pixels. At each successive level of detail, the map width and height grow by a factor of 2: Level 2 is 1024 x 1024 pixels, Level 3 is 2048 x 2048 pixels, Level 4 is 4096 x 4096 pixels, and so on.

For example, at level 3, the pixel coordinates range from (0, 0) to (2047, 2047), like this:

Pixel coordinates at level 3

Pixel coordinates at level 3

To optimize the performance of map retrieval and display, the rendered map is cut into tiles of \(256 x 256\) pixels each. As the number of pixels differs at each level of detail, so does the number of tiles:

map width = map height = \(2^{level}\) tiles

Each tile is given XY coordinates ranging from (0, 0) in the upper left to \((2^{level}-1, 2^{level}-1)\) in the lower right.

Quadtree structure

Quadtree structure

Example: The XY coordinates of the SIG institute at Symbiosis University at zoom level 13 are:

(SIU <- getGeoCode("GeoInformatics, Symbiosis University, Pune, India"))
##      lat      lon 
## 18.53338 73.83364
(SIU_XY= LatLon2XY(SIU[1],SIU[2],zoom=13))
## $Tile
##        X    Y
## lon 5776 3666
## 
## $Coords
##            x        y
## lon 32.14468 181.6492

fetch map tile and plot it:

try({
download.file(paste0("http://mt1.google.com/vt/lyrs=m&x=",SIU_XY$Tile[1,"X"],"&y=",SIU_XY$Tile[1,"Y"],"&z=13"),"mtSIU.png", mode="wb")
})

img = readPNG("mtSIU.png")
par(pty="s", mar=rep(0,4))
plot(0:256,0:256, type='n', axes = FALSE)
rasterImage(img, 0,0, 256,256)
points(SIU_XY$Coords[1,"x"], 256-SIU_XY$Coords[1,"y"],col=2,pch=20)

Stitching map tiles versus static maps

The RgoogleMaps package

The Google Static Maps API (https://developers.google.com/maps/documentation/static-maps/intro) lets you embed a Google Maps image on your webpage without requiring JavaScript or any dynamic page loading. The Google Static Map service creates your map based on URL parameters sent through a standard HTTP request and returns the map as an image you can display on your web page.

This package serves two purposes:

  1. Provide a comfortable R interface to query the Google server for static maps
  2. Use the map as a background image to overlay plots within R. This requires proper coordinate scaling.

The first step naturally will be to download a static map from the Google server. A simple example is shown in this Figure:

par(pty="s")
map1 <- GetMap(markers = paste0(
"&markers=color:blue|label:S|40.702147,-74.015794",
"&markers=color:green|label:G|40.711614,-74.012318",
"&markers=color:red|color:red|label:C|40.718217,-73.998284"),
center=NULL,destfile = "MyTile1.png",size=c(640,640),
maptype = "terrain")

PlotOnStaticMap(map1)

No data, just maps

The 2-step process from above (fetch a map -> plot) can be simplified:

par(pty="s")
mapBG1 = plotmap("Brandenburg Gate, Berlin", zoom = 15)

mapBG2 = plotmap("Brandenburg Gate, Berlin", zoom = 16, maptype="satellite")

Show traffic on bing maps

try({
par(pty="s")
mapBG3 = GetBingMap(center="Brandenburg Gate, Berlin", zoom=12, extraURL="&mapLayer=TrafficFlow", apiKey=apiKey, destfile="BerlinTraffic.png")
PlotOnStaticMap(mapBG3)
})

Customize google maps

no highways:

par(pty="s")

ManHatMap <- GetMap(center="Lower Manhattan", zoom=15,                  extraURL="&style=feature:road.highway|visibility:off")
  PlotOnStaticMap(ManHatMap)

Plots with spatial data

data(incidents)
col=as.numeric(incidents$Category)
par(pty="s")

mapSF_Z16 = plotmap(lat, lon, zoom = 16, col = col, pch=20, data = incidents)

zoom 15

par(pty="s")
mapSF_Z15 = plotmap(lat, lon, zoom = 15, col = col, pch=20, data = incidents)

zoom 13

par(pty="s")
mapSF_Z13 = plotmap(lat, lon, zoom = 13, col = col, pch=20, data = incidents, alpha = 0.7)

Working offline

It would be wasteful to have to fetch a new map from the map server for each new plot! Instead, we pass the map object to the next calls:

par(pty="s")
SundayCrimes = subset(incidents, DayOfWeek=="Sunday")
col=as.numeric(SundayCrimes$violent)
tmp=plotmap(lat, lon, mapSF_Z13, col = col, pch=20, data = SundayCrimes, alpha = 0.5)

meuse data with map background

  #RgoogleMaps
  crs = CRS("+init=epsg:28992") # set original projection
  data("meuse")
  coordinates(meuse) <- ~x+y
  proj4string(meuse) <- crs
  
  merc = CRS("+init=epsg:3857")
  WGS84 = CRS("+init=epsg:4326")
  #this funcion requires rgdal to be installed!
  meuse2 = spTransform(meuse, WGS84)
  meuseMap <- GetMap(rev(rowMeans(bbox(meuse2))), zoom=13, destfile="data/meuseMap.png")  
  
  meuseMap2 <- get_map(location=rowMeans(bbox(meuse2)), zoom=13)   # get Google map

LevCols = rev(heat.colors(5))[cut(meuse2$lead,breaks=c(0,100,200,300,400,Inf),labels=FALSE)]

meuseLatLon=as.data.frame(meuse2)
par(pty="s")
  

ggmap(meuseMap2) + 
  geom_point(data=meuseLatLon, aes(x,y,fill=lead), 
             color="grey70", size=2.5, shape=21)+
  scale_fill_gradientn(colours=rev(heat.colors(5)))

Density maps

# only violent crimes
data("crime")
violent_crimes <- subset(crime,
 offense != "auto theft" & offense != "theft" & offense != "burglary")

# order violent crimes
violent_crimes$offense <- factor(violent_crimes$offense,
 levels = c("robbery", "aggravated assault", "rape", "murder"))

# restrict to downtown
violent_crimes <- subset(violent_crimes,
 -95.39681 <= lon & lon <= -95.34188 &
 29.73631 <= lat & lat <= 29.78400)

houston <- get_map("houston", zoom = 14)
HoustonMap <- ggmap(houston)
HoustonMap +
stat_density2d(
aes(x = lon, y = lat, fill = ..level.., alpha = ..level..),
size = 2, bins = 4, data = violent_crimes,
geom = "polygon"
)
## Warning: Removed 11 rows containing non-finite values (stat_density2d).

Facets:

HoustonMap +
stat_density2d(aes(x = lon, y = lat, fill = ..level.., alpha = ..level..),
bins = 5, geom = "polygon",
data = violent_crimes) +
scale_fill_gradient(low = "black", high = "red") +
facet_wrap(~ day)
## Warning: Removed 11 rows containing non-finite values (stat_density2d).

Insets in RgoogleMaps

Crime Cluster

Conditional Plotting

data(incidents)

cols <- c("yellow", "darkred")
GoogleMap(~lat*lon|DayOfWeek,         #plot conditioned         
          data=incidents, map=NULL,
          groups=incidents$violent,                                                 
          col=cols, cex=0.1, alpha=0.3,                 
          panel=panel.loaPlot2,                         
          key.groups=list(main="Violent", 
                          cex=1, col=cols))    

loa_mapSF <- getMapArg()

Heat Maps:

v <- ifelse(incidents$violent==TRUE, "Violent", "Non-violent")
#weekend versus weekday:
WE = ifelse(incidents$DayOfWeek %in% c("Saturday", "Sunday"),"weekend", "weekday")

useOuterStrips(GoogleMap(~lat*lon|WE+v,         #plot conditioned         
          data=incidents, map=loa_mapSF,
          panel=panel.kernelDensity,                          #surface plot control                                                                
          col.regions=c("lightyellow", "darkred"),            #
          alpha.regions=0.5,                                  #                       
          col=1,                                              #surface calculation    
          n=200, at=c(0.5,1,2.5,5,10,25,50), 
          scales=list(draw=FALSE), xlab="", ylab=""))     

#useOuterStrips(trellis.last.object())

From static to interactive maps

Interactive Web Maps with the JavaScript ‘Leaflet’ Library

  m = leaflet() %>%  addTiles()
  m %>% leaflet::setView(SIU[2], SIU[1], zoom = 14) %>% addMarkers(SIU[2], SIU[1])

Ability to store map tiles offline

  m = leaflet::leaflet() %>% 
    addTiles( urlTemplate = "http:/localhost:8000/mapTiles/OSM/{z}_{x}_{y}.png")
  m = m %>% leaflet::setView(-73.99733, 40.73082 , zoom = 16)
  m = m %>% leaflet::addMarkers(-73.99733, 40.73082 )
  m

Posting map applications on the Web