Visualizing spatial data over a map (Putting it all together)

Pulling up maps in R is a fun distraction, but let’s put it together with some data.

qmplot sidenote

ggmap likes to build its plots in a function called qmplot which is the ggplot equivalent of qplot. Is that a bad thing? It appears that’s the way things are done with this package.

qplot is known as the “quick and dirty” way to quickly produce graphs, but it runs into issues when you try to do more complicated things.

example below of qplot vs ggplot (user kohske on stackoverflow)

qplot(x,y, geom=“line”) # I will use this ggplot(data.frame(x,y), aes(x,y)) + geom_line() # verbose

d <- data.frame(x, y)

qplot(x, y, data=d, geom=“line”) ggplot(d, aes(x,y)) + geom_line() # I will use this

qmap is a “wrapper” for ggmap and get_map: qmap(location = “X”, …)

Example with the crime dataset

The ggmap package has a ‘crime’ dataset built in; the dataset contains the crimes for Houston from Januaray to August (2010), geocoded with Google Maps.

head(crime)

##                      time     date hour premise            offense  beat
## 82729 2010-01-01 01:00:00 1/1/2010    0     18A             murder 15E30
## 82730 2010-01-01 01:00:00 1/1/2010    0     13R            robbery 13D10
## 82731 2010-01-01 01:00:00 1/1/2010    0     20R aggravated assault 16E20
## 82732 2010-01-01 01:00:00 1/1/2010    0     20R aggravated assault  2A30
## 82733 2010-01-01 01:00:00 1/1/2010    0     20A aggravated assault 14D20
## 82734 2010-01-01 01:00:00 1/1/2010    0     20R           burglary 18F60
##           block    street type suffix number   month    day
## 82729 9600-9699   marlive   ln      -      1 january friday
## 82730 4700-4799 telephone   rd      -      1 january friday
## 82731 5000-5099  wickview   ln      -      1 january friday
## 82732 1000-1099   ashland   st      -      1 january friday
## 82733 8300-8399    canyon           -      1 january friday
## 82734 9300-9399     rowan   ln      -      1 january friday
##                       location           address       lon      lat
## 82729    apartment parking lot   9650 marlive ln -95.43739 29.67790
## 82730 road / street / sidewalk 4750 telephone rd -95.29888 29.69171
## 82731        residence / house  5050 wickview ln -95.45586 29.59922
## 82732        residence / house   1050 ashland st -95.40334 29.79024
## 82733                apartment       8350 canyon -95.37791 29.67063
## 82734        residence / house     9350 rowan ln -95.54830 29.70223

#View(crime) #for full

We can see the data contains a lot of useful info: dates, types of crimes, locations by type of place, locations by street, locations by longitude/latitude.

First, let’s get an overview of the crimes on the map. Using qmplot, we put in longitude and latitude for the x, and y parameters, and specify the data as the crime dataset. This plots the crimes indiscriminantly.

# to see all crimes
qmplot(lon, lat, data = crime, maptype = "toner-lite", colour = I("red"),size = I(0.9),alpha=.3) +
  theme(legend.position="none") #to remove the .3 from the legend (pointless)

## Using zoom = 7...

## Map from URL : http://tile.stamen.com/toner-lite/7/28/49.png

## Map from URL : http://tile.stamen.com/toner-lite/7/29/49.png

## Map from URL : http://tile.stamen.com/toner-lite/7/30/49.png

## Map from URL : http://tile.stamen.com/toner-lite/7/31/49.png

## Map from URL : http://tile.stamen.com/toner-lite/7/28/50.png

## Map from URL : http://tile.stamen.com/toner-lite/7/29/50.png

## Map from URL : http://tile.stamen.com/toner-lite/7/30/50.png

## Map from URL : http://tile.stamen.com/toner-lite/7/31/50.png

## Map from URL : http://tile.stamen.com/toner-lite/7/28/51.png

## Map from URL : http://tile.stamen.com/toner-lite/7/29/51.png

## Map from URL : http://tile.stamen.com/toner-lite/7/30/51.png

## Map from URL : http://tile.stamen.com/toner-lite/7/31/51.png

## Map from URL : http://tile.stamen.com/toner-lite/7/28/52.png

## Map from URL : http://tile.stamen.com/toner-lite/7/29/52.png

## Map from URL : http://tile.stamen.com/toner-lite/7/30/52.png

## Map from URL : http://tile.stamen.com/toner-lite/7/31/52.png

## Map from URL : http://tile.stamen.com/toner-lite/7/28/53.png

## Map from URL : http://tile.stamen.com/toner-lite/7/29/53.png

## Map from URL : http://tile.stamen.com/toner-lite/7/30/53.png

## Map from URL : http://tile.stamen.com/toner-lite/7/31/53.png

## Map from URL : http://tile.stamen.com/toner-lite/7/28/54.png

## Map from URL : http://tile.stamen.com/toner-lite/7/29/54.png

## Map from URL : http://tile.stamen.com/toner-lite/7/30/54.png

## Map from URL : http://tile.stamen.com/toner-lite/7/31/54.png

## Warning: `panel.margin` is deprecated. Please use `panel.spacing` property
## instead

## Warning: Removed 5 rows containing missing values (geom_point).

Working with the data

Suppose we are interested in violent crimes in the downtown Houston area.

A slick method to isolate this could involve dplyr to create a ‘new pipe’ (similar to defining a new class method)

`%notin%` = function(setA, setB) !(setA %in% setB) #credit to kahle, checks that setA is not in setB

# reduce crime to violent crimes in downtown houston
violent_crimes = crime %>% filter(
    offense %notin% c("auto theft", "theft", "burglary"),
    -95.39681 <= lon & lon <= -95.34188,
    29.73631 <= lat & lat <=  29.78400
  ) %>% 
  mutate(
    offense = fct_drop(offense), #drops unused levels, does not drop NA levels that have values, will prevent errors
    offense = fct_relevel(offense,  #reorders factor levels by hand
                          c("robbery", "aggravated assault", "rape", "murder")
    )
  )

## Here is a traditional way to code it without pipes:
#
# qmplot(lon, lat, data = crime)
# 
# 
# # only violent crimes
# violent_crimes = subset(crime,
#   offense != "auto theft" &
#   offense != "theft" &
#   offense != "burglary"
# )
# 
# # rank 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
# )
# 
# theme_set(theme_bw())
# 
# qmplot(lon, lat, data = violent_crimes, colour = offense,
#   size = I(3.5), alpha = I(.6), legend = "topleft")
# 
# qmplot(lon, lat, data = violent_crimes, geom = c("point","density2d"))
# qmplot(lon, lat, data = violent_crimes) + facet_wrap(~ offense)
# qmplot(lon, lat, data = violent_crimes, extent = "panel") + facet_wrap(~ offense)
# qmplot(lon, lat, data = violent_crimes, extent = "panel", colour = offense, darken = .4) +
#   facet_wrap(~ month)

qmplot(lon, lat, data = violent_crimes, maptype = "toner-lite", colour = I("red"),size = I(0.9),alpha=.3) +
  theme(legend.position="none")

## Using zoom = 14...

## Map from URL : http://tile.stamen.com/toner-lite/14/3850/6770.png

## Map from URL : http://tile.stamen.com/toner-lite/14/3851/6770.png

## Map from URL : http://tile.stamen.com/toner-lite/14/3852/6770.png

## Map from URL : http://tile.stamen.com/toner-lite/14/3853/6770.png

## Map from URL : http://tile.stamen.com/toner-lite/14/3850/6771.png

## Map from URL : http://tile.stamen.com/toner-lite/14/3851/6771.png

## Map from URL : http://tile.stamen.com/toner-lite/14/3852/6771.png

## Map from URL : http://tile.stamen.com/toner-lite/14/3853/6771.png

## Map from URL : http://tile.stamen.com/toner-lite/14/3850/6772.png

## Map from URL : http://tile.stamen.com/toner-lite/14/3851/6772.png

## Map from URL : http://tile.stamen.com/toner-lite/14/3852/6772.png

## Map from URL : http://tile.stamen.com/toner-lite/14/3853/6772.png

## Map from URL : http://tile.stamen.com/toner-lite/14/3850/6773.png

## Map from URL : http://tile.stamen.com/toner-lite/14/3851/6773.png

## Map from URL : http://tile.stamen.com/toner-lite/14/3852/6773.png

## Map from URL : http://tile.stamen.com/toner-lite/14/3853/6773.png

## Warning: `panel.margin` is deprecated. Please use `panel.spacing` property
## instead

Different visualizations

In the previous example, the geom parameter is set by default to “point”. Let’s see it with the countour and bin geometries:

qmplot(x=lon, y=lat, data = violent_crimes, maptype = "toner-2011", geom = "density2d", colour = I("red"),size=I(1))

## Using zoom = 14...

## Map from URL : http://tile.stamen.com/toner-2011/14/3850/6770.png

## Map from URL : http://tile.stamen.com/toner-2011/14/3851/6770.png

## Map from URL : http://tile.stamen.com/toner-2011/14/3852/6770.png

## Map from URL : http://tile.stamen.com/toner-2011/14/3853/6770.png

## Map from URL : http://tile.stamen.com/toner-2011/14/3850/6771.png

## Map from URL : http://tile.stamen.com/toner-2011/14/3851/6771.png

## Map from URL : http://tile.stamen.com/toner-2011/14/3852/6771.png

## Map from URL : http://tile.stamen.com/toner-2011/14/3853/6771.png

## Map from URL : http://tile.stamen.com/toner-2011/14/3850/6772.png

## Map from URL : http://tile.stamen.com/toner-2011/14/3851/6772.png

## Map from URL : http://tile.stamen.com/toner-2011/14/3852/6772.png

## Map from URL : http://tile.stamen.com/toner-2011/14/3853/6772.png

## Map from URL : http://tile.stamen.com/toner-2011/14/3850/6773.png

## Map from URL : http://tile.stamen.com/toner-2011/14/3851/6773.png

## Map from URL : http://tile.stamen.com/toner-2011/14/3852/6773.png

## Map from URL : http://tile.stamen.com/toner-2011/14/3853/6773.png

## Warning: `panel.margin` is deprecated. Please use `panel.spacing` property
## instead

qmplot(x=lon, y=lat, data = violent_crimes, maptype = "toner-2011", geom = "bin2d")

## Using zoom = 14...

## Warning: `panel.margin` is deprecated. Please use `panel.spacing` property
## instead

Neither of these graphs are particularly great, so let’s see if we can improve on them. The countour is a bit hard to tell what areas have higher crime, and the bins are differentiated by hue, but it clashes a bit with the map itself. Trying to colour the bin puts a border around it, which only makes the hue harder to distinguish (optical illusion).

Heatmap

To make the countour map more useful, we can assign a gradient. Let’s look at the robberies:

robberies = violent_crimes %>% filter(offense=='robbery')

qmplot(lon, lat, data = violent_crimes, geom = "blank", 
       zoom = 15, maptype = "toner-background", legend = "bottomright"
) +
  stat_density_2d(aes(fill = ..level..), geom = "polygon", alpha = .35, colour = NA) +
  scale_fill_gradient2("Robbery\nHeatmap", low = "white", mid = "yellow", high = "red", midpoint = 650)

## 49 tiles needed, this may take a while (try a smaller zoom).

## Map from URL : http://tile.stamen.com/toner-background/15/7700/13541.png

## Map from URL : http://tile.stamen.com/toner-background/15/7701/13541.png

## Map from URL : http://tile.stamen.com/toner-background/15/7702/13541.png

## Map from URL : http://tile.stamen.com/toner-background/15/7703/13541.png

## Map from URL : http://tile.stamen.com/toner-background/15/7704/13541.png

## Map from URL : http://tile.stamen.com/toner-background/15/7705/13541.png

## Map from URL : http://tile.stamen.com/toner-background/15/7706/13541.png

## Map from URL : http://tile.stamen.com/toner-background/15/7700/13542.png

## Map from URL : http://tile.stamen.com/toner-background/15/7701/13542.png

## Map from URL : http://tile.stamen.com/toner-background/15/7702/13542.png

## Map from URL : http://tile.stamen.com/toner-background/15/7703/13542.png

## Map from URL : http://tile.stamen.com/toner-background/15/7704/13542.png

## Map from URL : http://tile.stamen.com/toner-background/15/7705/13542.png

## Map from URL : http://tile.stamen.com/toner-background/15/7706/13542.png

## Map from URL : http://tile.stamen.com/toner-background/15/7700/13543.png

## Map from URL : http://tile.stamen.com/toner-background/15/7701/13543.png

## Map from URL : http://tile.stamen.com/toner-background/15/7702/13543.png

## Map from URL : http://tile.stamen.com/toner-background/15/7703/13543.png

## Map from URL : http://tile.stamen.com/toner-background/15/7704/13543.png

## Map from URL : http://tile.stamen.com/toner-background/15/7705/13543.png

## Map from URL : http://tile.stamen.com/toner-background/15/7706/13543.png

## Map from URL : http://tile.stamen.com/toner-background/15/7700/13544.png

## Map from URL : http://tile.stamen.com/toner-background/15/7701/13544.png

## Map from URL : http://tile.stamen.com/toner-background/15/7702/13544.png

## Map from URL : http://tile.stamen.com/toner-background/15/7703/13544.png

## Map from URL : http://tile.stamen.com/toner-background/15/7704/13544.png

## Map from URL : http://tile.stamen.com/toner-background/15/7705/13544.png

## Map from URL : http://tile.stamen.com/toner-background/15/7706/13544.png

## Map from URL : http://tile.stamen.com/toner-background/15/7700/13545.png

## Map from URL : http://tile.stamen.com/toner-background/15/7701/13545.png

## Map from URL : http://tile.stamen.com/toner-background/15/7702/13545.png

## Map from URL : http://tile.stamen.com/toner-background/15/7703/13545.png

## Map from URL : http://tile.stamen.com/toner-background/15/7704/13545.png

## Map from URL : http://tile.stamen.com/toner-background/15/7705/13545.png

## Map from URL : http://tile.stamen.com/toner-background/15/7706/13545.png

## Map from URL : http://tile.stamen.com/toner-background/15/7700/13546.png

## Map from URL : http://tile.stamen.com/toner-background/15/7701/13546.png

## Map from URL : http://tile.stamen.com/toner-background/15/7702/13546.png

## Map from URL : http://tile.stamen.com/toner-background/15/7703/13546.png

## Map from URL : http://tile.stamen.com/toner-background/15/7704/13546.png

## Map from URL : http://tile.stamen.com/toner-background/15/7705/13546.png

## Map from URL : http://tile.stamen.com/toner-background/15/7706/13546.png

## Map from URL : http://tile.stamen.com/toner-background/15/7700/13547.png

## Map from URL : http://tile.stamen.com/toner-background/15/7701/13547.png

## Map from URL : http://tile.stamen.com/toner-background/15/7702/13547.png

## Map from URL : http://tile.stamen.com/toner-background/15/7703/13547.png

## Map from URL : http://tile.stamen.com/toner-background/15/7704/13547.png

## Map from URL : http://tile.stamen.com/toner-background/15/7705/13547.png

## Map from URL : http://tile.stamen.com/toner-background/15/7706/13547.png

## Warning: `panel.margin` is deprecated. Please use `panel.spacing` property
## instead

Faceting different categories, and not using qmplot

We can also look at different things, depending on what we’re interested in:

#faceting by crime
(
qmplot(lon, lat, data = violent_crimes, maptype = "toner-background", colour = offense) + 
  facet_wrap(~ offense)
)

## Using zoom = 14...

## Map from URL : http://tile.stamen.com/toner-background/14/3850/6770.png

## Map from URL : http://tile.stamen.com/toner-background/14/3851/6770.png

## Map from URL : http://tile.stamen.com/toner-background/14/3852/6770.png

## Map from URL : http://tile.stamen.com/toner-background/14/3853/6770.png

## Map from URL : http://tile.stamen.com/toner-background/14/3850/6771.png

## Map from URL : http://tile.stamen.com/toner-background/14/3851/6771.png

## Map from URL : http://tile.stamen.com/toner-background/14/3852/6771.png

## Map from URL : http://tile.stamen.com/toner-background/14/3853/6771.png

## Map from URL : http://tile.stamen.com/toner-background/14/3850/6772.png

## Map from URL : http://tile.stamen.com/toner-background/14/3851/6772.png

## Map from URL : http://tile.stamen.com/toner-background/14/3852/6772.png

## Map from URL : http://tile.stamen.com/toner-background/14/3853/6772.png

## Map from URL : http://tile.stamen.com/toner-background/14/3850/6773.png

## Map from URL : http://tile.stamen.com/toner-background/14/3851/6773.png

## Map from URL : http://tile.stamen.com/toner-background/14/3852/6773.png

## Map from URL : http://tile.stamen.com/toner-background/14/3853/6773.png

## Warning: `panel.margin` is deprecated. Please use `panel.spacing` property
## instead

#faceting by day

#First specify your map
houston = get_googlemap('houston texas', zoom=15, maptype='roadmap')

## Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=houston+texas&zoom=15&size=640x640&scale=2&maptype=roadmap&sensor=false

## Information from URL : http://maps.googleapis.com/maps/api/geocode/json?address=houston%20texas&sensor=false

ggmap(houston)

#then add the layers like we would in ggplot
HoustonMap = ggmap(houston, base_layer = ggplot(aes(x = lon, y = lat),
                                                 data = violent_crimes))
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) +
  theme(legend.position="none")

## Warning: Removed 467 rows containing non-finite values (stat_density2d).

The map we supplied here (roadmap with google) is probably a bit too colourful, but this is just to illustrate the capabilities.

The following was VERY experimental, and I couldn’t get it to work out properly… perhaps too many locations, but the idea is that we can include more information by adding an additional piece to the legend –> We let the offense be coded by colour, and the goal was to let the type of location be denoted by a particular symbol.

theme_set(theme_bw())
houmap = qmap("houston",zoom=13, colour='bw', legend ='topleft')

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

## Information from URL : http://maps.googleapis.com/maps/api/geocode/json?address=houston&sensor=false

## Warning: `panel.margin` is deprecated. Please use `panel.spacing` property
## instead

(houmap+
    geom_point(data=violent_crimes, aes(x=lon,y=lat,colour=offense,shape=location))
)

## Warning: The shape palette can deal with a maximum of 6 discrete values
## because more than 6 becomes difficult to discriminate; you have
## 60. Consider specifying shapes manually if you must have them.

## Warning: Removed 616 rows containing missing values (geom_point).

(houmap+
    geom_point(data=violent_crimes, aes(x=lon,y=lat,colour=offense,shape=location))+
    theme(legend.position="none")
)

## Warning: The shape palette can deal with a maximum of 6 discrete values
## because more than 6 becomes difficult to discriminate; you have
## 60. Consider specifying shapes manually if you must have them.

## Warning: Removed 616 rows containing missing values (geom_point).

We can also put the data into bins if we are interested in an overview (if we are more interested in zones).

(
houmap+
    stat_bin2d(
      data=violent_crimes,
      aes(x=lon,y=lat,colour=offense,fill=offense),
      size=0.6, bins= 90, alpha=.4
    )
)

#compare to closer, we can use less bins
houmap = qmap("houston",zoom=14, colour='bw', legend ='topleft')

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

## Information from URL : http://maps.googleapis.com/maps/api/geocode/json?address=houston&sensor=false

## Warning: `panel.margin` is deprecated. Please use `panel.spacing` property
## instead

(houmap+
    geom_point(data=violent_crimes, aes(x=lon,y=lat,colour=offense))
)

## Warning: Removed 11 rows containing missing values (geom_point).

(
  houmap+
    stat_bin2d(
      data=violent_crimes,
      aes(x=lon,y=lat,colour=offense,fill=offense),
      size=0.6, bins= 30, alpha=.4
    )
)

## Warning: Removed 11 rows containing non-finite values (stat_bin2d).

## Warning: Removed 8 rows containing missing values (geom_tile).

Depending on the level of detail we want, we can manually adjust the zoom, and the bin number.

Note, it is important that as you zoome in, you’ll want to reduce the number of bins, or you may as well just use the dot-plot since each crime will fit uniquely in its own bin.

4crime example

Alex(ks)’s Bushby and Jovic

March 20, 2018

packages