Traffic Maps and Rasters
Josh Pepper
10/29/2016
1 Load libraries
You’ll have to install these libraries if you don’t already have them
#stephanie traffic code
#load libraries
library(leaflet)
library(rgdal)
library(RColorBrewer)
library(raster)
library(gstat)
library(maptools)2 QGIS Multipart to Single Parts
I used the QGIS Geometry Tool to turn Multi parts into single parts. This converted Multipoints into individual points which makes it possible to import into R. 
3 Load traffic data
#load traffic data
trafficCounts <- readOGR("aadt.geojson", "OGRGeoJSON", disambiguateFIDs=T)## OGR data source with driver: GeoJSON
## Source: "aadt.geojson", layer: "OGRGeoJSON"
## with 13616 features
## It has 14 fieldssummary(trafficCounts)## Object of class SpatialPointsDataFrame
## Coordinates:
## min max
## coords.x1 -124.26486 -114.29893
## coords.x2 32.54431 42.00547
## Is projected: FALSE
## proj4string :
## [+proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs]
## Number of points: 13616
## Data attributes:
## Id DYNSEGPM PMSTATUS District
## Min. :0 ALA 880 R0 : 6 Valid:13616 Min. : 0.000
## 1st Qu.:0 LA 5 R44.924 : 6 1st Qu.: 4.000
## Median :0 KER 178 R1.959: 4 Median : 6.000
## Mean :0 LAS 36 R29.394: 4 Mean : 6.164
## 3rd Qu.:0 MNO 395 120.49: 4 3rd Qu.: 8.000
## Max. :0 SAC 160 35.045: 4 Max. :12.000
## (Other) :13588
## Route County Postmile
## 101 : 1008 LA :1556 0 : 564
## 5 : 849 SD : 912 R0 : 178
## 1 : 560 SBD : 712 L0 : 32
## 99 : 508 KER : 540 M0 : 10
## 10 : 312 ORA : 529 0.23 : 8
## 80 : 311 RIV : 504 1.25 : 8
## (Other):10068 (Other):8863 (Other):12816
## Descriptn Back_pk_h Back_pk_m
## JCT. RTE. 5 : 72 Min. : 0 Min. : 0
## BREAK IN ROUTE : 40 1st Qu.: 590 1st Qu.: 6000
## JCT. RTE. 101 : 40 Median : 2200 Median : 24100
## JCT. RTE. 99 : 28 Mean : 4979 Mean : 63786
## NEVADA STATE LINE : 27 3rd Qu.: 8100 3rd Qu.: 98000
## PLACER/NEVADA COUNTY LINE: 25 Max. :89000 Max. :665000
## (Other) :13384
## Back_AADT Ahead_pk_h Ahead_pk_m Ahead_AADT
## Min. : 0 Min. : 0 Min. : 0 Min. : 0
## 1st Qu.: 5000 1st Qu.: 590 1st Qu.: 6000 1st Qu.: 5000
## Median : 22000 Median : 2250 Median : 24200 Median : 22000
## Mean : 60477 Mean : 4955 Mean : 63742 Mean : 60425
## 3rd Qu.: 93000 3rd Qu.: 8100 3rd Qu.: 98000 3rd Qu.: 93000
## Max. :602600 Max. :31000 Max. :410000 Max. :377000
## 4 Make map
#make color scale for Ahead_AADT
spectral <- brewer.pal(11, "Spectral") %>% rev()
AADTColor <- colorNumeric(spectral, trafficCounts@data$Ahead_AADT)
#make map
leaflet() %>% addProviderTiles("CartoDB.Positron") %>%
addCircleMarkers(data=trafficCounts,
stroke=FALSE, radius=5, color = ~AADTColor(Ahead_AADT), fillOpacity=0.3,
popup=~paste("Description: ",Descriptn,
"<br>Route: ", Route,
"<br>Back AADT: ", Back_AADT,
"<br>Ahead AADT: ", Ahead_AADT,
sep="")) %>%
addLegend(title="AADT", pal = AADTColor, values = trafficCounts@data$Ahead_AADT, position="bottomleft")5 Make raster layer
#assign latitude and longitude base on coordinates
trafficCounts$latitude <- trafficCounts@coords[,2]
trafficCounts$longitude <- trafficCounts@coords[,1]
#make raster
sSp <- as(SpatialPoints(trafficCounts@data[,c("longitude","latitude")]), "ppp") # convert points to pp class
Dens <- density(sSp, adjust = 0.1) # create density object
r <- raster(Dens)
crs(r) <- sp::CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
pal.act <- colorNumeric(c("#ffffd400","#fed98e","#fe9929","#d95f0e","#993404"), values(r), na.color = "transparent", alpha=TRUE)
pal.bin <- colorBin(c("#ffffd400","#fed98e","#fe9929","#d95f0e","#993404"), values(r), na.color = "transparent", alpha=TRUE)
leaflet() %>% addProviderTiles("CartoDB.Positron") %>%
addRasterImage(r, opacity=0.7, colors=pal.act) %>%
addLegend(title="Traffic 'Exposure'", pal=pal.bin, values = values(r), position="bottomleft")Note that because you are “spreading” out the value of the AADT, the value is not a direct measure of the AADT. Rather, it is a value that is a proxy for how close they are to traffic. The units are somewhat meaningless. If you increase the adjust coeffecient in the above code, the exposure values will become more concentrated over less space. To illustrate this point, I’ll include a map below using an adjust coeffecient of 0.05. Notice the change in the raster, but also the corresponding doubling of the legend values.
6 Raster map of bay area only
#download box around bay area
bayBox <- readOGR("bay.geojson", "OGRGeoJSON", disambiguateFIDs=T)## OGR data source with driver: GeoJSON
## Source: "bay.geojson", layer: "OGRGeoJSON"
## with 1 features
## It has 0 fieldscrs(bayBox) <- crs(trafficCounts@proj4string)
#subset for bay area only
bay <- trafficCounts[bayBox,]
#make raster
sSp <- as(SpatialPoints(bay@data[,c("longitude","latitude")]), "ppp") # convert points to pp class
Dens <- density(sSp, adjust = 0.1) # create density object
r <- raster(Dens)
crs(r) <- sp::CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
#extract value at one point
point1 <- SpatialPoints(data.frame(-122.298, 37.867))
point2 <- SpatialPoints(data.frame(-122.168, 37.806))
#assign mean at point
value1 <- extract(r, point1)
value2 <- extract(r, point2)
#make color palette
pal.act <- colorNumeric(c("#ffffd400","#fed98e","#fe9929","#d95f0e","#993404"), values(r), na.color = "transparent", alpha=TRUE)
pal.bin <- colorBin(c("#ffffd400","#fed98e","#fe9929","#d95f0e","#993404"), values(r), na.color = "transparent", alpha=TRUE)
#make leaflet map
leaflet() %>% addProviderTiles("CartoDB.Positron") %>% addRasterImage(r, opacity=0.7, colors=pal.act) %>%
addPolygons(data=bayBox, fill=FALSE, color="grey") %>%
addMarkers(data=point1, popup=paste("Exposure:", round(value1))) %>%
addMarkers(data=point2, popup=paste("Exposure:", round(value2))) %>%
addLegend(title="Traffic 'Exposure'", pal=pal.bin, values = values(r), position="bottomleft")