By Ken Steif

Replication of data science workflows is as critical in government as it is in any domain.

Much of Philadelphia’s vast open data resources are updated nightly or weekly, providing intelligence at high space/time resolutions.

Typically, my students develop algorithms that are less replicable, beginning their workflow by downloading a snapshot of data from an open data site in csv or shapefile form.

This semester, I have tried to give them tools to plug their code directly in to open data APIs – a timely feed of city administrative data. To do so, I’ve introduced R packages including ckanr, esri2sf, jsonLite and others.

I wanted to demonstrate a particular introductory workflow that I recently adopted for a series of data-driven neighborhood reports.

The goal was to create a very simple workflow in R, allowing a user to specify a PostGRES query; feed that query in to one of several Philadelphia open data APIs which are supported by Carto; wrangle that data to Simple Features and then create a series of maps and plots.

I demonstrate the workflow briefly below – one that can be reused for many of the Philadelphia Open Data APIs. The talented folks at the Office of Open Data and Digital Transformation have created helpful vignettes for each API. For example, here is the page for the crime data.

Begin by loading the requisite libraries, and specifying a ‘map theme’ that we can use to customize our data visualizations.

Setup library(lubridate) library(tidyverse) library(sf) library(jsonlite) library(esri2sf) mapTheme <- function(base_size = 12) { theme( text = element_text( color = "black"), plot.title = element_text(size = 14,colour = "black"), plot.subtitle=element_text(face="italic"), plot.caption=element_text(hjust=0), axis.ticks = element_blank(), panel.background = element_blank(),axis.title = element_blank(), axis.text = element_blank(), axis.title.x = element_blank(), axis.title.y = element_blank(), panel.grid.minor = element_blank(), panel.border = element_rect(colour = "black", fill=NA, size=2) ) } 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 library ( lubridate ) library ( tidyverse ) library ( sf ) library ( jsonlite ) library ( esri2sf ) mapTheme < - function ( base_size = 12 ) { theme ( text = element_text ( color = "black" ) , plot . title = element_text ( size = 14 , colour = "black" ) , plot . subtitle = element_text ( face = "italic" ) , plot . caption = element_text ( hjust = 0 ) , axis . ticks = element_blank ( ) , panel . background = element_blank ( ) , axis . title = element_blank ( ) , axis . text = element_blank ( ) , axis . title . x = element_blank ( ) , axis . title . y = element_blank ( ) , panel . grid . minor = element_blank ( ) , panel . border = element_rect ( colour = "black" , fill = NA , size = 2 ) ) }

The below query is written in PostGRES, not R. If copied and pasted in to a browser, it returns raw crime data in json format – like so.

Setup https://phl.carto.com/api/v2/sql?q=SELECT * FROM incidents_part1_part2 WHERE (dispatch_date_time >= current_date - 30) AND (dc_dist = '18') AND Text_General_Code in ('Aggravated Assault No Firearm', 'Aggravated Assault Firearm', 'Burglary Residential', 'Motor Vehicle Theft') 1 2 3 4 5 6 7 https : / / phl . carto . com / api / v2 / sql ? q = SELECT * FROM incidents_part1_part2 WHERE ( dispatch_date_time >= current_date - 30 ) AND ( dc_dist = '18' ) AND Text_General_Code in ( 'Aggravated Assault No Firearm' , 'Aggravated Assault Firearm' , 'Burglary Residential' , 'Motor Vehicle Theft' )

Note that we use the Police District (‘dc_dist’) here to query the data spatially. Grab the resulting url and paste it into the ‘fromJSON’ command as below. Note the browser has recoded the POSTGRES into a url that can be read in to R.

Setup crimeRecentJSON <- fromJSON("https://phl.carto.com/api/v2/sql?q=SELECT%20*%20FROM%20incidents_part1_part2%20WHERE%20(dispatch_date_time%20%3E=%20current_date%20-%2030)%20AND%20(dc_dist%20=%20%2718%27)%20AND%20Text_General_Code%20in%20(%27Aggravated%20Assault%20No%20Firearm%27,%20%27Aggravated%20Assault%20Firearm%27,%20%27Burglary%20Residential%27,%20%27Motor%20Vehicle%20Theft%27)") 1 crimeRecentJSON < - fromJSON ( "https://phl.carto.com/api/v2/sql?q=SELECT%20*%20FROM%20incidents_part1_part2%20WHERE%20(dispatch_date_time%20%3E=%20current_date%20-%2030)%20AND%20(dc_dist%20=%20%2718%27)%20AND%20Text_General_Code%20in%20(%27Aggravated%20Assault%20No%20Firearm%27,%20%27Aggravated%20Assault%20Firearm%27,%20%27Burglary%20Residential%27,%20%27Motor%20Vehicle%20Theft%27)" )

Next, the json can be converted into a Simple Features data frame. We create new X and Y coordinates fields; a ‘yearQuarter’ time period field; and recode the crime type field into a ‘Crime_Type’ which will read easier on a legend.

Setup crimeRecent <- as.data.frame(crimeRecentJSON$rows) %>% select(point_x,point_y,dispatch_date,text_general_code,dispatch_date_time) %>% na.omit %>% st_as_sf(coords = c("point_x", "point_y"), crs = 4326) %>% mutate( X = st_coordinates(.)[,1], Y = st_coordinates(.)[,2], yearQuarter = paste(year(dispatch_date),"-", quarter(dispatch_date),sep=''), year = year(dispatch_date), Crime_Type = text_general_code) 1 2 3 4 5 6 7 8 9 10 11 crimeRecent < - as . data . frame ( crimeRecentJSON $ rows ) % > % select ( point_x , point_y , dispatch_date , text_general_code , dispatch_date_time ) % > % na . omit % > % st_as_sf ( coords = c ( "point_x" , "point_y" ) , crs = 4326 ) % > % mutate ( X = st_coordinates ( . ) [ , 1 ] , Y = st_coordinates ( . ) [ , 2 ] , yearQuarter = paste ( year ( dispatch_date ) , "-" , quarter ( dispatch_date ) , sep = '' ) , year = year ( dispatch_date ) , Crime_Type = text_general_code )

Recall, we downloaded the data for the entire 18th Police District, but we actually need data at a smaller neighborhood geography. As such, the code below uses the esri2sf package to download the neighborhood boundary shapefile from OpenDataPhilly, selects the Spruce Hill neighborhood and converts to Simple Features. It does the same for the street centerlines file and uses the former shapefile to clip the latter.

Setup #Get the neighborhood boundary thisRCO <- esri2sf("http://gis.phila.gov/arcgis/rest/services/PhilaGov/RCO/MapServer/0") %>% filter(ORGANIZATION_NAME == "Spruce Hill Community Association") #Get the streets that are just in the neighborhood streets <- esri2sf("https://services.arcgis.com/fLeGjb7u4uXqeF9q/ArcGIS/rest/services/Street_Centerline/FeatureServer/0") streets <- st_intersection(thisRCO, streets) #Sometimes, non-polyline geometries are returned. This gets rid of them. streets <- streets[st_is(streets, 'LINESTRING'),] 1 2 3 4 5 6 7 8 9 10 11 #Get the neighborhood boundary thisRCO < - esri2sf ( "http://gis.phila.gov/arcgis/rest/services/PhilaGov/RCO/MapServer/0" ) % > % filter ( ORGANIZATION_NAME == "Spruce Hill Community Association" ) #Get the streets that are just in the neighborhood streets < - esri2sf ( "https://services.arcgis.com/fLeGjb7u4uXqeF9q/ArcGIS/rest/services/Street_Centerline/FeatureServer/0" ) streets < - st_intersection ( thisRCO , streets ) #Sometimes, non-polyline geometries are returned. This gets rid of them. streets < - streets [ st_is ( streets , 'LINESTRING' ) , ]

Using the resulting data frames, we create the below basemap.

Now all that’s left is to layer the crime data. First we perform a spatial selection by selecting all crimes that are inside the ‘thisRCO’ boundary. We then map `thisRCO.crimesRecent’. Note that the labels pull directly from the data frame ensuring that this routine is replicable no matter which neighborhood is selected.

Setup #Perform a spatial selection thisRCO.crimesRecent <- crimeRecent[thisRCO,] ggplot() + geom_sf(data=thisRCO) + geom_sf(data=streets) + geom_point(data=thisRCO.crimesRecent, aes(x=X,y=Y,colour=Crime_Type), size=4) + labs(title="Assaults, burglaries & motor vehicle thefts, Last 30 days", subtitle=paste(thisRCO$ORGANIZATION_NAME,"Source: Philadelphia Police Department", sep="; "), caption=paste(Sys.Date()-30, Sys.Date(), sep=" through ")) + annotate("text", label="S. 42nd St.", x=-75.2065, y=39.9528, angle=83, fontface=2, colour="dodgerblue4",size=3) + mapTheme() 1 2 3 4 5 6 7 8 9 10 11 12 13 14 #Perform a spatial selection thisRCO . crimesRecent < - crimeRecent [ thisRCO , ] ggplot ( ) + geom_sf ( data = thisRCO ) + geom_sf ( data = streets ) + geom_point ( data = thisRCO . crimesRecent , aes ( x = X , y = Y , colour = Crime_Type ) , size = 4 ) + labs ( title = "Assaults, burglaries & motor vehicle thefts, Last 30 days" , subtitle = paste ( thisRCO $ ORGANIZATION_NAME , "Source: Philadelphia Police Department" , sep = "; " ) , caption = paste ( Sys . Date ( ) - 30 , Sys . Date ( ) , sep = " through " ) ) + annotate ( "text" , label = "S. 42nd St." , x = - 75.2065 , y = 39.9528 , angle = 83 , fontface = 2 , colour = "dodgerblue4" , size = 3 ) + mapTheme ( )

If you make one simple change to your PostGRES query (‘current_date >= n’), you can download historical data as well. Here is the crime trend for the last 8 years.

That’s it. As I mentioned, this simple workflow can be adopted for a host of datasets including zoning, Licenses & Inspections, real estate transactions, 311, parking data and more. To check out these and more, head over to OpenDataPhilly and look for the “API Documentation” associated with a given dataset.

Ken Steif, PhD is the founder of Urban Spatial. He is also the director of the Master of Urban Spatial Analytics program at the University of Pennsylvania. You can follow him on Twitter @KenSteif.

March 3, 2018