Visualizing Baltimore 2: Vacant Property and Some More Crime

One of the key predictors in my model for this crime project I'm working on is vacant houses and lots. I'll speak to some findings about the relationship between levels of the different types of crime and vacant property in a later post. But I wanted to put some of these images up now, before I'm done with that, after a conversation at a party tonight.

This plot is of 15,928 vacant buildings and 17,169 vacant lots (according to the datasets here) across the city of Baltimore:

Here are visualizations of the 2-dimensional kernel density estimates for both of them. A density estimate essentially gives values at every point on a plane that communicate how close that point is to how many observations of the variable or point process you care about. So the more red, the more vacant properties are clustered together in that area.

And here are kernel density visualizations for homicide and aggravated assault:

And if you're interested, all the data is from here and here's the code:

## gis libraries
library(spBayes)
library(MBA)
library(geoR)
library(fields)
library(sp)
library(maptools)
library(rgdal)
library(classInt)
library(lattice)
library(xtable)
library(spatstat)
library(splancs)

## Other packages
library(ggplot2)
library(foreign)
library(stringr)
library(lubridate)
library(plyr)
library(xtable)
library(scales)
library(RColorBrewer)
library(grid)
library(ggmap)
library(gridExtra)
library(ggmcmc)

setwd('/home/rmealey/Dropbox/school/gisClass/FinalProject')
options(digits=10)

Save <- function(projName){
	savehistory(paste(projName,'.Rhistory',sep=''))
	save.image(paste(projName,'.RData',sep=''))
}

sv <- function() Save('FinalProject')

########################################################################
# City Boundary Shape File
city_df <- read.dbf('Baltcity_20Line/baltcity_line.dbf')
city_shp <- readOGR(dsn='Baltcity_20Line', layer='baltcity_line')
origProj <- city_shp@proj4string ## Store original projection
#city_shp = spTransform(city_shp,CRS("+proj=longlat +datum=WGS84"))
city_pl_df <- fortify(city_shp, region='LABEL')
cityLineCoords <- data.frame(city_shp@lines[[1]]@Lines[[1]]@coords)
cityLinePoly <- Polygon(cityLineCoords)
cityLinePolys <- Polygons(list(cityLinePoly), ID='cityline')
cityLineSpPoly <- SpatialPolygons(list(cityLinePolys),proj4string=origProj)

cityLineCoords <- cityLineCoords[,c(2,1)])

########################################################################
## Neighborhood Shape Files read in v1
nbhds_df <- read.dbf('Neighborhood_202010/nhood_2010.dbf')
nbhds_shp <- readOGR(dsn='Neighborhood_202010', layer='nhood_2010')
origProj <- nbhds_shp@proj4string ## Store original projection
#nbhds_shp = spTransform(nbhds_shp,CRS("+proj=longlat +datum=WGS84"))
nbhds_pl_df <- fortify(nbhds_shp, region='LABEL')

########################################################################
## Utility Functions

## Read lat/lng coords function
str2LatLong <- function(in_df){
	latlng <- str_replace(str_replace(in_df$Location.1,'\$$',''),')','')
	latlng <- str_split(latlng,', ')
	latlng_df <- ldply(latlng[in_df$Location.1 != ''])
	out_df <- in_df
	out_df$lat <- as.numeric(latlng_df[,1])
	out_df$long <- as.numeric(latlng_df[,2])
	return(out_df)
}

## convert projection function
convProj <- function(in_df,in_proj,out_proj){
	latlong <- in_df[,c('long','lat')]
	latlong_spdf <- SpatialPoints(latlong, 
	proj4string=in_proj)
	latlong_spdf <-  spTransform(latlong_spdf,out_proj)
	latlong_spdf_coords <- coordinates(latlong_spdf)
	out_df <- in_df
	out_df$long <- latlong_spdf_coords[,1]
	out_df$lat <- latlong_spdf_coords[,2]
	return(out_df)
}

########################################################################
## Preprocess Crime Data
crimeData <- read.csv('OpenDataSets/BPD_Part_1_Victim_Based_Crime_Data.csv')
crimeData_NoCoords <- crimeData[crimeData$Location.1 == '',]
crimeData <- crimeData[crimeData$Location.1 != '',]

## Get and convert projection
crimeData <- str2LatLong(crimeData)
## Incidents already in correct proj
crimeData_ProjOrig <- crimeData[crimeData$lat>100,]
crimeData <- crimeData[crimeData$lat<100,]
inProj <- CRS("+proj=longlat +datum=WGS84")
outProj <- origProj
crimeData <- convProj(crimeData, inProj, outProj)

## Parse Dates
crimeData$crimeDate2 <- parse_date_time(
	crimeData$crimeDate,
	orders='%m/%d/%Y'
)

## Get Burglary Incidents
burg_df <- subset(crimeData, description=='BURGLARY')
## Hold Out 2012 Incidents
burg_df_ho <- subset(burg_df, year(crimeDate2) == '2012') 
burg_df <- subset(burg_df, year(crimeDate2) != '2012') 
ggplot(data=burg_df, aes(x=long,y=lat)) + geom_point() + coord_equal()

## Get Street Robbery Incidents
robbStr_df <- subset(crimeData, description=="ROBBERY - STREET")
## Hold Out 2012 Incidents
robbStr_df_ho <- subset(robbStr_df, year(crimeDate2) == '2012') 
robbStr_df <- subset(robbStr_df, year(crimeDate2) != '2012') 
ggplot(data=robbStr_df, aes(x=long,y=lat)) + geom_point() + coord_equal()

## Homicide
homic_df <- subset(crimeData, description=='HOMICIDE')
## Hold Out 2012 Incidents
homic_df_ho <- subset(homic_df, year(crimeDate2) == '2012') 
homic_df <- subset(homic_df, year(crimeDate2) != '2012') 
ggplot(data=homic_df, aes(x=long,y=lat)) + geom_point() + coord_equal()

## Aggravated Assault
aggrAslt_df <- subset(crimeData, description=='AGG. ASSAULT')
## Hold Out 2012 Incidents
aggrAslt_df_ho <- subset(aggrAslt_df, year(crimeDate2) == '2012') 
aggrAslt_df <- subset(aggrAslt_df, year(crimeDate2) != '2012') 
ggplot(data=aggrAslt_df, aes(x=long,y=lat)) + geom_point() + coord_equal()


########################################################################
# Religous Building Locations
relig_df <- read.csv('geocoded/Religious_Buildings_gc.csv')
## Remove na rows
relig_df <- relig_df[complete.cases(relig_df),]
inProj <- CRS("+proj=longlat +datum=WGS84")
outProj <- origProj
relig_df <- convProj(relig_df, inProj, outProj)

########################################################################
# Police Station Locations
police_df <- read.csv('geocoded/Police_Stations_gc.csv')
inProj <- CRS("+proj=longlat +datum=WGS84")
outProj <- origProj
police_df <- convProj(police_df, inProj, outProj)

########################################################################
# Hospitals Locations
hospitals_df <- read.csv('geocoded/Hospitals.csv')
inProj <- CRS("+proj=longlat +datum=WGS84")
outProj <- origProj
hospitals_df <- convProj(hospitals_df, inProj, outProj)

########################################################################
# CCTV Locations
cams_df <- read.csv('OpenDataSets/CCTV_Locations.csv')
cams_df <- str2LatLong(cams_df)
inProj <- CRS("+proj=longlat +datum=WGS84")
outProj <- origProj
cams_df <- convProj(cams_df, inProj, outProj)
cams_df$type <- "CCTV Camera"

########################################################################
# Vacant Buildings
vacantBuildings_df <- read.csv('OpenDataSets/Vacant_Buildings.csv')
vacantBuildings_df <- str2LatLong(vacantBuildings_df)
inProj <- CRS("+proj=longlat +datum=WGS84")
outProj <- origProj
vacantBuildings_df <- convProj(vacantBuildings_df, inProj, outProj)
vacantBuildings_df$type <- 'Vacant Building'

########################################################################
# Vacant Lots
vacantLots_df <- read.csv('OpenDataSets/Vacant_Lots.csv')
vacantLots_df <- str2LatLong(vacantLots_df)
inProj <- CRS("+proj=longlat +datum=WGS84")
outProj <- origProj
vacantLots_df <- convProj(vacantLots_df, inProj, outProj)
vacantLots_df$type <- 'Vacant Lot'

########################################################################
## Get kernel density estimates
kde2dRange <- c(apply(burg_df[,c('long','lat')], 2, range))

getKde <- function(in_df, N=400, Lims=kde2dRange){
	pts <- as.matrix(in_df[,c('long','lat')])
	dens <- kde2d(pts[,1],pts[,2], n=N, lims=Lims)
	dens_df <- data.frame(expand.grid(dens$x, dens$y), z = c(dens$z))
	colnames(dens_df) <- c('x','y','z')
	return(dens_df)
}

plotKde2d <- function(in_df){
	fillCols <- rev(brewer.pal(11,'Spectral'))
	return(
		ggplot() + 
		geom_tile(data = in_df, aes(x=x, y=y, fill=z, group=1)) + 
		scale_fill_gradientn(colours=fillCols) + 
		theme_bw() +
		coord_equal()
	)
}

saveKde2Plot <- function(plotDf, plotName, plotTitle,titlCol='white'){
	# https://github.com/wch/ggplot2/wiki/New-theme-system
	new_theme_empty <- theme_bw()
	new_theme_empty$line <- element_blank()
	new_theme_empty$rect <- element_blank()
	new_theme_empty$strip.text <- element_blank()
	new_theme_empty$axis.text <- element_blank()
	new_theme_empty$plot.title <- element_blank()
	new_theme_empty$axis.title <- element_blank()
	new_theme_empty$legend.position <- 'none'
	new_theme_empty$plot.margin <- structure(c(0, 0, -1, -1), unit = "lines", valid.unit = 3L, class = "unit")

	nbhds_pl_df2 <- nbhds_pl_df[,c('long','lat','group')]
	colnames(nbhds_pl_df2) <- c('x','y','group')
	plotKde2d(plotDf) + 
		geom_path(data=nbhds_pl_df2,aes(x=x,y=y,
		group=group),color='black',alpha=0.4) +
		new_theme_empty + 
		annotate("text", x = 1405000, y = 568000,
		label=plotTitle, size=8, color=titlCol)
		ggsave(paste('img/',plotName,'.png', sep=''))
}

## Get all simple gaussian 2d kernel density estimates
burgDens <- getKde(burg_df) 				## Burglary, 7
robbStrDens <- getKde(robbStr_df) 			## Street Robbery, 7
homicDens <- getKde(homic_df) 				## Homicide, 7
aggrAsltDens <- getKde(aggrAslt_df) 			## Aggr Assault, 7
hospitalsDens <- getKde(hospitals_df)		## Hospitals
policeDens <- getKde(police_df)				## Police Stations
religDens <- getKde(relig_df)				## Religous Buildings
camsDens <- getKde(cams_df)					## Cameras, 1
vacBldDens <- getKde(vacantBuildings_df)	## Vacant Buildings, 5
vacLotsDens <- getKde(vacantLots_df)		## Vacant Lots, 6

## plot densities
saveKde2Plot(burgDens, 'BurglaryKde2d', 'Burglary\n Density')
saveKde2Plot(robbStrDens, 'StreetRobberyKde2d', 'Street\n Robbery\n Density')
saveKde2Plot(homicDens, 'HomicideKde2d', 'Homicide\n Density')
saveKde2Plot(aggrAsltDens, 'aggrAsltKde2d', 'Aggravated\n Assault\n Density')
saveKde2Plot(hospitalsDens, 'HospitalKde2d', 'Hospital\n Location\n Density')
saveKde2Plot(policeDens, 'PoliceKde2d', 'Police\n Station\n Density')
saveKde2Plot(religDens, 'ReligiousKde2d', 'Religous\n Building\n Density')
saveKde2Plot(camsDens, 'CCTVCamsKde2d', 'CCTV\n Cameras\n Density')
saveKde2Plot(vacBldDens, 'VacBldgKde2d', 'Vacant\n Building\n Density')
saveKde2Plot(vacLotsDens, 'vacLotsKde2d', 'Vacant\n Lot\n Density')