[This article was first published on Obscure Analytics » Rstats, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
A few quick plots of West Baltimore neighborhoods, first Sandtown-Winchester:
and Harlem Park:
These aren’t very polished, I’ll put up better versions.
Here’s the code for those that want it:
## 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)
########################################################################
# 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
# Source:
## 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')
names(nbhds_shp@polygons) <- nbhds_shp@data$LABEL
## Neighborhood Shape Files read in v2 (from spatstat docs)
#nbhds_shp <- readShapePoly('Neighborhood_202010/nhood_2010.shp')
#nbhds_sp <- as(nbhds_shp, "SpatialPolygons")
#nbhds_owin <- as(nbhds_sp, "owin")
#centroids <- coordinates(nbhds_shp)
hoodNames <- 'Mount Vernon'
ggplot(data=nbhds_pl_df[nbhds_pl_df$id==hoodNames,],
aes(x=long, y=lat, group=group)) +
geom_path() +
ggtitle(hoodNames) +
coord_equal()
## plot actual city shape using empty nbhd boundaries
city_plot <- bound_plot +
geom_polygon(data=nbhds_pl_df, fill='white',color='white')
ggsave('img/emptyCity.png')
## plot nbhd boundaries
nbhds_plot <- bound_plot +
geom_polygon(data=nbhds_pl_df,color='gray',fill='white')
ggsave('img/nbhds.png')
########################################################################
# Parcel Shape Polygon Data
# Source:
parcel_df <- read.dbf('Parcel_Shp/parcel.dbf')
parcel_shp <- readOGR(dsn='Parcel_Shp', layer='parcel')
parcel_df <- data.frame(parcel_df, coordinates(parcel_shp))
parcel_mtrx <- as.matrix(coordinates(parcel_shp))
########################################################################
# Vacant Buildings
# Source:
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'
vacBld_mtrx <- as.matrix(vacantBuildings_df[,c('long','lat')])
########################################################################
# Vacant Lots
# Source:
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'
vacantLots_mtrx <- as.matrix(vacantLots_df[,c('long','lat')])
########################################################################
## Plot by neighborhood
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)
crime_mtrx <- as.matrix(crimeData[,c('long','lat')])
## 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()
########################################################################
# Vacant Lots
SandtownWinchester_df <- fortify(nbhds_shp@polygons[['Sandtown-Winchester']])
sw_mtr <- as.matrix(SandtownWinchester_df[,1:2])
sw_props <- data.frame(pip(parcel_mtrx, sw_mtr))
sw_vacB <- data.frame(pip(vacBld_mtrx, sw_mtr))
sw_vacL <- data.frame(pip(vacBld_mtrx, sw_mtr))
sw_crime <- data.frame(pip(crime_mtrx, sw_mtr))
sw_crime <- crimeData[rownames(sw_crime),]
sw_crime_2012 <- subset(sw_crime, year(crimeDate2)==2012)
colnames(sw_props) <- c('long','lat')
colnames(sw_vacB) <- c('long','lat')
colnames(sw_vacL) <- c('long','lat')
# 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 <- 'bottom'
new_theme_empty$plot.margin <- structure(c(0, 0, -1, -1), unit = "lines", valid.unit = 3L, class = "unit")
crimeCols <- brewer.pal(12,'Paired')
crimeTypes <- list('RAPE'=c(crimeCols[1],crimeCols[2]),
'ARSON'=c(crimeCols[1],crimeCols[2]),
'COMMON ASSAULT'=c(crimeCols[3],crimeCols[4]),
'AGG. ASSAULT'=c(crimeCols[3],crimeCols[4]),
'SHOOTING'=c(crimeCols[5],crimeCols[6]),
'HOMICIDE'=c(crimeCols[5],crimeCols[6]),
'ROBBERY - STREET'=c(crimeCols[7],crimeCols[8]),
'ROBBERY - CARJACKING'=c(crimeCols[7],crimeCols[8]),
'ROBBERY - RESIDENCE'=c(crimeCols[7],crimeCols[8]),
'ROBBERY - COMMERCIAL'=c(crimeCols[7],crimeCols[8]),
'BURGLARY'=c(crimeCols[9],crimeCols[10]),
'LARCENY'=c(crimeCols[9],crimeCols[10]),
'AUTO THEFT'=c(crimeCols[11],crimeCols[12]),
'LARCENY FROM AUTO'=c(crimeCols[11],crimeCols[12]))
crimeCols <- as.data.frame(t(data.frame(crimeTypes)))
col_cols <- crimeCols[,2]
names(col_cols) <- names(crimeTypes)
ggplot(data = SandtownWinchester_df) +
geom_polygon(aes(x=long, y=lat, group=group),
color='black', fill='white') +
geom_point(data = sw_props, aes(x=long, y=lat),
shape = 0, color = 'gray') +
geom_point(data = sw_vacB, aes(x=long, y=lat),
shape = 4, color = 'red') +
geom_point(data = sw_crime_2012,
aes(x=long, y=lat, color=description), shape = 'o',size=2) +
scale_color_manual(values = col_cols) +
coord_equal() +
annotate("text", x = 1415200, y = 598300,
label="Sandtown-\nWinchester\nVacant Properties\nand Crime",
size=6, color="black") +
new_theme_empty +
guides(color=guide_legend("",nrow=5)) +
ggsave('img/SandtownWinchesterVacantsandCrime.png')
########################################################################
# Vacant Lots
HarlemPark_df <- fortify(nbhds_shp@polygons[['Harlem Park']])
hp_mtr <- as.matrix(HarlemPark_df[,1:2])
hp_props <- data.frame(pip(parcel_mtrx, hp_mtr))
hp_vacB <- data.frame(pip(vacBld_mtrx, hp_mtr))
hp_vacL <- data.frame(pip(vacBld_mtrx, hp_mtr))
hp_crime <- data.frame(pip(crime_mtrx, hp_mtr))
hp_crime <- crimeData[rownames(hp_crime),]
hp_crime_2012 <- subset(hp_crime, year(crimeDate2)==2012)
colnames(hp_props) <- c('long','lat')
colnames(hp_vacB) <- c('long','lat')
colnames(hp_vacL) <- c('long','lat')
# 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 <- 'bottom'
new_theme_empty$plot.margin <- structure(c(0, 0, -1, -1), unit = "lines", valid.unit = 3L, class = "unit")
crimeCols <- brewer.pal(12,'Paired')
crimeTypes <- list('RAPE'=c(crimeCols[1],crimeCols[2]),
'ARSON'=c(crimeCols[1],crimeCols[2]),
'COMMON ASSAULT'=c(crimeCols[3],crimeCols[4]),
'AGG. ASSAULT'=c(crimeCols[3],crimeCols[4]),
'SHOOTING'=c(crimeCols[5],crimeCols[6]),
'HOMICIDE'=c(crimeCols[5],crimeCols[6]),
'ROBBERY - STREET'=c(crimeCols[7],crimeCols[8]),
'ROBBERY - CARJACKING'=c(crimeCols[7],crimeCols[8]),
'ROBBERY - RESIDENCE'=c(crimeCols[7],crimeCols[8]),
'ROBBERY - COMMERCIAL'=c(crimeCols[7],crimeCols[8]),
'BURGLARY'=c(crimeCols[9],crimeCols[10]),
'LARCENY'=c(crimeCols[9],crimeCols[10]),
'AUTO THEFT'=c(crimeCols[11],crimeCols[12]),
'LARCENY FROM AUTO'=c(crimeCols[11],crimeCols[12]))
crimeCols <- as.data.frame(t(data.frame(crimeTypes)))
col_cols <- crimeCols[,2]
names(col_cols) <- names(crimeTypes)
hpplot <- ggplot(data = HarlemPark_df) +
geom_polygon(aes(x=long, y=lat, group=group),
color='black', fill='white') +
geom_point(data = hp_props, aes(x=long, y=lat),
shape = 0, color = 'gray',size=5) +
geom_point(data = hp_vacB, aes(x=long, y=lat),
shape = 4, color = 'red',size=5) +
geom_point(data = hp_crime_2012,
aes(x=long, y=lat, color=description), shape = 'o',size=3) +
scale_color_manual(values = col_cols) +
coord_equal() +
annotate("text", x = 1416400, y = 594500,
label="Harlem Park\nVacant Properties\nand Crime",
size=4, color="black") +
new_theme_empty +
guides(color=guide_legend("",nrow=5))
ggsave('img/HarlemParkVacantsandCrime.png',width=11, height=8.5) |
To leave a comment for the author, please follow the link and comment on their blog: Obscure Analytics » Rstats.
R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.