This is a brief demonstration of common data manipulation and mapping techniques using spatial analysis tools in R. The goal here is to go from knowing nothing about shapefiles, to being able to create meaningful and attractive maps. For this exercise, I use crime data from Portland, Oregon that was provided by the National Institute of Justice for their crime prediction competition. The data and R scripts for this exercise are available at its Github repository.
Loading Spatial Objects from Shape Files
Spatial data usually comes in the form of shapefiles, which can be loaded using the readOGR function from the rgdal package.
library(rgdal)
port <- readOGR(dsn = "data", layer = "Police_Districts_Portland")
crime <- readOGR(dsn = "data", layer = "NIJ_Nov2016_Crime")
The port object contains information describing the geographic area of the 60 policing districts in Portland. The crime object contains data relating to every crime reported in Portland during the month of November 2016.
Exploring Spatial Objects
If you’re used to regular R dataframes, you’ll find that spatial objects (represented by S4 objects) are a little weird. But this allows a lot of flexibility. The first thing to know is that each spatial object belongs to one of many different “classes”. The most elementary type of spatial data are Points, which denote single point locations. Lines are a set of points connected by straight line segments that can be used to denote things like rivers or state boundaries. Polygons are shapes that denote area, and a Grid is a collection of cells organized into a regular grid. Running the class() command shows that the port object is a Spatial Polygons Data Frame and the crime object is a Spatial Points Data Frame. We can plot spatial objects easily using base R plot functions.
par(mfrow=c(1,2))
plot(port); title(main = list("Portland Police Districts (Polygon Object)", cex=0.8))
plot(crime); title(main = list("November 2016 Crime Reports (Points Object)", cex=0.8))
The map on the left is made up of 60 polygons, each representing a policing district within Portland. The Polygon class is especially suited to describing regions or zones. The plot on the right shows a scattering of points that each describe the location of a particular crime reported in November 2016. The Points class is particularly suited to describing exact locations.
Structure of Spatial Objects
The summary() command is a useful way of exploring a spatial object.
summary(crime)
Object of class SpatialPointsDataFrame
Coordinates:
min max
coords.x1 7608809 7717689
coords.x2 651489 727287
Is projected: TRUE
proj4string :
[+proj=lcc +lat_1=44.33333333333334 +lat_2=46 +lat_0=43.66666666666666
+lon_0=-120.5 +x_0=2500000 +y_0=0 +ellps=GRS80 +units=ft +no_defs]
Number of points: 16834
Data attributes:
CATEGORY occ_date census_tra
BURGLARY : 95 2016/11/04: 645 Min. : 0
MOTOR VEHICLE THEFT: 270 2016/11/03: 625 1st Qu.: 2000
OTHER :13889 2016/11/05: 607 Median : 4700
STREET CRIMES : 2580 2016/11/08: 607 Mean : 6609
2016/11/06: 602 3rd Qu.: 8201
2016/11/01: 599 Max. :980000
(Other) :13149
The first line of the summary() output confirms that the crime object is a Spatial Points Data Frame. A spatial object is made up of a number of different “slots” that, together, completely describe the geographical coordinates of the shape relevant to its class (polygons in this case). The exact number and types of slots differ according to the class of the object.
The next few lines in the output above summarize the contents of the Coordinates slot, which holds the geographic coordinates of each crime. The displayed values are known as the “bounding box,” which gives the four points (or “corners”) that denote the spatial extent of the data. The spatial extent is basically the corners of the smallest rectangle that contains all the crime points. These values can be accessed with the bbox() command.
bbox(crime)
min max
coords.x1 7608809 7717689
coords.x2 651489 727287
The next lines describe the projection characteristics of the crime object. A projection is an arbitrary coordinate reference system that provides a basis upon which to measure distances between objects in spatial data. All geographic objects are based on some form of coordinate system, and there are many different types of projections. The proj4string slot contains a string describing the projection information of the crime object’s coordinates. We’ll see below that some packages require us to change the projection system to one of the more common systems. The next line states that there were 16,834 recorded crime incidents in November 2016.
In addition to the geographic information, spatial objects also have additional data attributes that are contained in the data slot, which is actually a traditional R data frame. The output above shows that the data Data Frame has three columns - one describing the Category of each crime, one stating the date on which the crime occurred and one stating the census tract of each crime. Each row in the dataframe corresponds to a particular crime. The number of rows is therefore equal to the total number of crimes. Each slot can be accessed on its own via the @ operator, as follows.
head(crime@data)
CATEGORY occ_date census_tra
1 STREET CRIMES 2016/11/01 5200
2 STREET CRIMES 2016/11/01 7202
3 STREET CRIMES 2016/11/01 10600
4 STREET CRIMES 2016/11/01 1101
5 STREET CRIMES 2016/11/01 2000
6 STREET CRIMES 2016/11/01 1201
Manipulating spatial data
It turns out that some of the crimes contained in crime occur outside the city of Portland defined by the port polygons. We can clip the crimes that occur outside the city using simple base R functions. Plotting the clipped and non-clipped data illustrates the point.
crime_clp <- crime[port, ] # "Clip" crime to stay in same area
par(mfrow=c(1,2))
plot(port); points(crime); title(main = list("Original Crime Data", cex=0.8))
plot(port); points(crime_clp); title(main = list("Clipped Crime Data", cex=0.8))
When working with spatial data, we often want to make “cross polygon” comparisons. The aggregate function is very useful for this because it allows us to perform operations “by” polygons. The following code counts the number of crimes contained in crime by the polygons defined in port.
crime_agg <- aggregate(x=crime["CATEGORY"],by=port,FUN=length) # Aggregate by district
summary(crime_agg)
Object of class SpatialPolygonsDataFrame
Coordinates:
min max
x 7604004.6 7701431.1
y 651315.6 733815.4
Is projected: TRUE
proj4string :
[+proj=lcc +lat_1=44.33333333333334 +lat_2=46 +lat_0=43.66666666666666
+lon_0=-120.5 +x_0=2500000 +y_0=0 +ellps=GRS80 +units=ft +no_defs]
Data attributes:
CATEGORY
Min. :109.0
1st Qu.:194.2
Median :280.0
Mean :279.8
3rd Qu.:345.8
Max. :512.0
The aggregate function identifies which port polygon (police district) each crime is committed in and groups them accordingly. Here we used the CATEGORY variable from the crime@data dataframe, although any variable could have been used here as we are merely counting how many points exist. The function then counts the number of crime points in each district, using the length function.
The summary(crime_agg) command shows that the new spatial object crime_agg is a polygon dataframe which has the same projection properties and “bounding box” as the port object. This is because the aggregate function simply counted the number of crimes whose x, y coordinates were in each of the polygons described by the port object. crime_agg only has one data attribute, however, which corresponds to the crime counts within each police district. Because there are only 60 police districts, we can actually just print the CATEGORY data column, which lists the crime count within each of the 60 districts. Each entry represents the number of crimes committed in the row corresponding to each polygon.
print(crime_agg$CATEGORY)
[1] 262 176 295 327 210 228 253 186 285 143 140 123 164 345 353 239 177 279 288 401 444
[22] 307 370 400 292 211 143 439 276 395 427 144 487 480 182 148 120 174 220 109 285 265
[43] 168 268 197 215 320 248 238 512 284 451 316 348 391 390 337 281 292 337
Plotting spatial data with Leaflet
There are a number of very useful packages for plotting spatial data such as ggmap and tmap, which I will explore in a later post. Here, we will use Leaflet. It isn’t the simplest package, but it offers a lot of power, is apparently the leading open-source mapping library and, in my opinion, looks cool. The first thing to do is to transform the projection of the coordinates, because some of the packages we will use require spatial data to be in a certain projection.
crime_agg <- spTransform(crime_agg, CRS("+init=epsg:4326")) # Reproject coordinates
We’re going to create a choropleth map. The colorBin function of the Leaflet package divides the values into a specified number of bins (5 in this case) and maps these values to colors following a palette, in this case a palette of Reds. You can find information about the other palettes by typing ?RColorBrewer.
library(leaflet)
qpal <- colorBin("Reds", crime_agg$CATEGORY, bins=5)
Leaflet’s syntax is very similar to the popular ggplot syntax in the sense that you add layers to the map in steps. The first is to create a leaflet object from the crime_agg spatial object. We then use the addPolygons layer to overlay the police district outlines onto the map. We use the pallette mapping created above in the fillColor option. Finally, we add a legend.
leaflet(crime_agg) %>%
addPolygons(stroke = TRUE,opacity = 1,fillOpacity = 0.5, smoothFactor = 0.5,
color="black",fillColor = ~qpal(CATEGORY),weight = 1) %>%
addLegend(values=~CATEGORY,pal=qpal,title="Crime Count")
Creating a spatial grid
When working with spatial data, it’s often useful to work within our own evenly-spaced grid. The raster package is a simple way to make grids.
library(raster)
e <- extent(bbox(port)) # define boundaries of object
r <- raster(e) # create raster object
dim(r) <- c(40, 40) # specify number of cells
projection(r) <- CRS(proj4string(port)) # give it the same projection as port
g <- as(r, 'SpatialPolygonsDataFrame') # convert into polygon
The first step in the code above is to define an extent object, which basically defines the “boundaries” of our raster grid. Here, I’ve used the boundaries of the port object accessed by the bbox command. The next line creates the raster object with our given extent. We then specify how many cells we’d like the raster to have. Here we specify the dimensions of the raster, which splits the extent into a 40 x 40 square grid. Finally, I give the raster cell the same projection orientation as port. The resulting r is a raster object, which we then convert into a spatial polygon object g. Inspecting g, we see that it has the same min and max x,y coordinates as port because we defined this in the extent above. There are no data attributes because we haven’t defined any.
summary(g)
Object of class SpatialPolygonsDataFrame
Coordinates:
min max
x 7604004.6 7701431.1
y 651315.6 733815.4
Is projected: TRUE
proj4string :
[+proj=lcc +lat_1=44.33333333333334 +lat_2=46 +lat_0=43.66666666666666
+lon_0=-120.5 +x_0=2500000 +y_0=0 +ellps=GRS80 +units=ft +no_defs]
Data attributes:
layer
Min. : NA
1st Qu.: NA
Median : NA
Mean :NaN
3rd Qu.: NA
Max. : NA
NA's :1600
We now clip the new grid to match Portland and perform the same aggregation procedure to count the number of crimes within our newly defined (square) polygons.
p <- g[port,]
crime_agg1 <- aggregate(x=crime["CATEGORY"],by=p,FUN=length)
Something to watch out for when using a fine grid, which are now a lot smaller than the initial police districts, is that some polygons may not have any crimes committed within them. The length function therefore takes the length of an empty vector in those cases, which returns a NA value. To avoid ugly NA parts on the map, we replace them with with zeros.
crime_agg1$CATEGORY[is.na(crime_agg1$CATEGORY)] <- 0
Note that the above is simple base R code for dealing with dataframes. That’s because the object crime_agg1$CATEGORY is simply an R dataframe, which happens to live inside a spatial object.
crime_agg1 <- spTransform(crime_agg1, CRS("+init=epsg:4326")) # reproject
qpal <- colorBin("Reds", crime_agg1$CATEGORY, bins=5) # define color bins
leaflet(crime_agg1) %>%
addPolygons(stroke = TRUE,opacity = 1,fillOpacity = 0.5, smoothFactor = 0.5,
color="black",fillColor = ~qpal(CATEGORY),weight = 0.5) %>%
addLegend(values=~CATEGORY,pal=qpal,title="Crime Count")
The above looks pretty close, but because each of the grids are square, its only a rough approximation of the portland area. Ideally, we’d like our grid to exactly overlay the Portland map. This means trimming some of the outer squares. Using the rgeos package, we first “melt” the 60 Portland polygons into a single portland polygon. We then use the gIntersection function to find the polygon that is enclosed by both the portland polygons and the grid p we created above.
library(rgeos)
# Melt Portland districts into one polygon
portland = gUnaryUnion(port, port$dummy)
# Take the intersection of the grid and the Portland map
map <- gIntersection(p,portland,byid = TRUE,drop_lower_td = TRUE)
With this new polygon spatial object, we can then perform the same calculations.
crime_agg2 <- aggregate(x=crime["CATEGORY"],by=map,FUN=length)
crime_agg2$CATEGORY[is.na(crime_agg2$CATEGORY)] <- 0
crime_agg2 <- spTransform(crime_agg2, CRS("+init=epsg:4326"))
qpal <- colorBin("Reds", crime_agg2$CATEGORY, bins=5)
plot3 <- leaflet(crime_agg2) %>%
addPolygons(stroke = TRUE,opacity = 1,fillOpacity = 0.5, smoothFactor = 0.5,
color="black",fillColor = ~qpal(CATEGORY),weight = 0.5) %>%
addLegend(values=~CATEGORY,pal=qpal,title="Crime Count")
Looking pretty good. Now we can add an open source version of a Portland map to see how our map looks on the actual city. We do this by appending the addTiles() function to our original plot
plot4 <- plot3 %>% addTiles()
Notice too that Leaflet plots are fully interactive. This means you can set the zoom using the setView() function as follows. We first find the center of portland using the gCentroid function, and then set the zoom. Default is 10, so this makes it slightly more zoomed in.
cent <- gCentroid(crime_agg2) # Find center of map
plot5 <- plot4 %>% setView(zoom = 11,lng=cent@coords[[1]], lat=cent@coords[[2]])
Saving Shape Files
The rgdal package that we used to read in spatial data also has a useful writeOGR function for saving shapefiles.
writeOGR(crime_agg2, dsn='output_data', layer='Portland', driver="ESRI Shapefile",
overwrite_layer=T)