-
Notifications
You must be signed in to change notification settings - Fork 2
/
for_graphics.R
369 lines (283 loc) · 18.3 KB
/
for_graphics.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
library(rgdal)
library(rgeos)
library(tidyverse)
library(tigris)
library(ggplot2)
library(geosphere)
library(tidycensus)
source("../scripts/keys.R")
library(viridis)
library(ggmap)
library(DT)
library(leaflet)
library(readr)
library(lubridate)
library(dplyr)
library(stringr)
library(sp)
library(knitr)
source("keys.R")
# Brings in a census API key from the keys.R script (which is hidden)
# Sign up here: http://api.census.gov/data/key_signup.html
census_api_key(key)
# Alternatively: census_api_key("your_key_goes_here")
# Pulling the shapefile for New Jersey's Hudson county census tracts
nj_h <- tracts("NJ", county="Hudson", cb=F)
# Figuring out the latitude and longitudes for the center of each census tract
nj_h_centroids <- SpatialPointsDataFrame(gCentroid(nj_h, byid=TRUE),
nj_h@data, match.ID=FALSE)
nj_h_centroids <- as.data.frame(nj_h_centroids)
nj_h_centroids <- select(nj_h_centroids, GEOID, x, y)
# Setting up some blank columns
nj_h_centroids$distance_proj1 <- 0
nj_h_centroids$distance_proj2 <- 0
# Determining the distance in miles between the project site and the center of each census tract in Hudson county
for (i in 1:nrow(nj_h_centroids)) {
nj_h_centroids$distance_proj1[i] <- distm (c(-74.063644, 40.734330), c(nj_h_centroids$x[i], nj_h_centroids$y[i]), fun = distHaversine)[,1]/ 1609
nj_h_centroids$distance_proj2[i] <- distm (c(-74.03577, 40.72), c(nj_h_centroids$x[i], nj_h_centroids$y[i]), fun = distHaversine)[,1]/ 1609
}
# 1 Journal Square (2017)
# Census tracts included by planners in this project
# GEOIDs: 34017001900, 34017004600, 34017005300, 34017006600, 34017006700, 34017007100
# Census tracts immediately adjacent to the project
# GEIODs: 19, 20 18, 17, 9.02, 12.02, 71
# Getting the unemployment figures and total labor workforce for all tracts in Hudson county
hudson_unemployment_2013 <- get_acs(geography="tract", endyear=2013, variables= c("B23025_005E", "B23025_002E"), county = "Hudson", state="NJ")
# Why 2009-2013? That's the ACS 5-year data that developers used to justify their project to officials, according to documents
# Discarding the margin of error
hudson_unemployment_2013$moe <- NULL
# Tidying up the structure of the new data frame
hudson_unemployment_2013 <- spread(hudson_unemployment_2013,variable, estimate)
# Calculating the unemployment rate
hudson_unemployment_2013$per_un <- round(hudson_unemployment_2013[,4]/hudson_unemployment_2013[,3]*100,2)
# Filtering out the blank census tracts
hudson_unemployment_2013 <- filter(hudson_unemployment_2013, GEOID!="34017006900" & GEOID!="34017980100")
# Array of group-designated census tracts
proj1 <- c("34017001900", "34017004600", "34017005300", "34017006600", "34017006700", "34017007100")
# Filtering out the tracts data specific to the project
hudson_unemployment_2013_sm <- filter(hudson_unemployment_2013, GEOID %in% proj1)
# Cleaning up the data frame
colnames(hudson_unemployment_2013_sm) <- c("GEOID", "name", "total", "unemployed", "unemp_rate")
hudson_unemployment_2013_sm$un_rate <- hudson_unemployment_2013_sm$unemp_rate$B23025_005
hudson_unemployment_2013_sm$unemp_rate <- NULL
# Adding a label column
hudson_unemployment_2013_sm$radius <- "gerrymandered"
# Array of census tracts immediately adjacent to the project
proj1_half <- c("34017001900", "34017002000", "34017001800", "34017001700", "34017000902", "34017001202", "34017007100")
# Filtering out the tracts data specific to the project
hudson_unemployment_2013_sm_half <- filter(hudson_unemployment_2013, GEOID %in% proj1_half)
# Cleaning up the data frame
colnames(hudson_unemployment_2013_sm_half) <- c("GEOID", "name", "total", "unemployed", "unemp_rate")
hudson_unemployment_2013_sm_half$un_rate <- hudson_unemployment_2013_sm_half$unemp_rate$B23025_005
hudson_unemployment_2013_sm_half$unemp_rate <- NULL
# Adding a label column
chudson_unemployment_2013_sm_half$radius <- "adjacent tracts"
# Combining the two census tracts data frames
hudson_unemployment_sm <- rbind(hudson_unemployment_2013_sm, hudson_unemployment_2013_sm_half)
# Cleaning up the census tract names
hudson_unemployment_sm$name <- gsub(",.*", "", hudson_unemployment_sm$name)
# More cleaning
hudson_unemployment_2013$un_rate <- hudson_unemployment_2013$per_un$B23025_005
hudson_unemployment_2013$per_un <- NULL
# Prepping the shape file to be joined with data
# (This process is necessary when visualizing with ggplot)
nj_hf <- fortify(nj_h, region="GEOID")
# Joining census data of project to shape file
nj_hud <- left_join(nj_hf, hudson_unemployment_sm, by=c("id"="GEOID"))
# Filter out the census tracts with no data
nj_hud <- filter(nj_hud, !is.na(radius))
# Joining census data of all tracts in county to shape file
nj_all <- left_join(nj_hf, hudson_unemployment_2013, by=c("id"="GEOID"))
# Slicing up the unemployment data (for all tracts) into buckets
nj_all$nj_mid <- cut(nj_all$un_rate, breaks = seq(0, 30, by=5))
# colors <- colorRampPalette(c("white", "red"))(length(levels(nj_all$nj_mid)))
# Slicing up the unemployment data (for project tracts) into buckets
nj_hud$nj_mid <- cut(nj_hud$un_rate, breaks = seq(0, 30, by=5))
#colors <- colorRampPalette(c("white", "red"))(length(levels(nj_hud$nj_mid)))
# Mapping project-only census tract data
nj_map <- ggplot()
nj_map <- nj_map + geom_polygon(data=nj_hud, aes(x=long, y=lat, group=group, fill=factor(nj_mid)), color="black", size=.5)
nj_map <- nj_map + facet_wrap(~radius)
nj_map <- nj_map + coord_map()
nj_map <- nj_map + scale_fill_manual(drop=FALSE, values=c("blue", "yellow", "orange", "red", "green", "purple"), na.value="#EEEEEE", name="Unemployment rate")
nj_map <- nj_map + theme_nothing(legend=TRUE)
nj_map <- nj_map + labs(x=NULL, y=NULL, title="1 Journal Square (2017)")
nj_map <- nj_map + theme(panel.grid.major = element_line(colour = NA))
nj_map <- nj_map + theme(text = element_text(size=15))
nj_map <- nj_map + theme(plot.title=element_text(face="bold", hjust=.4))
nj_map <- nj_map + theme(plot.subtitle=element_text(face="italic", size=9, margin=margin(l=20)))
nj_map <- nj_map + theme(plot.caption=element_text(size=12, margin=margin(t=12), color="#7a7d7e", hjust=0))
nj_map <- nj_map + theme(legend.key.size = unit(1, "cm"))
nj_map <- nj_map + annotate("segment", x = -74.063644, xend = -74.035, y = 40.734330, yend = 40.734330, colour = "tomato", size=.5)
nj_map <- nj_map + annotate("point", x = -74.063644, y = 40.734330, colour = "lightblue", size = 1)
nj_map <- nj_map + annotate("text", x = -74.01, y = 40.734330, label = "1 Journal Square", size=3, colour="gray30")
print(nj_map)
# Saving as a pdf
ggsave("map1a.pdf", nj_map, device="pdf")
# Mapping all census tract borders in the county and layering on the project-only census tract data
nj_map <- ggplot()
nj_map <- nj_map + geom_polygon(data=nj_hf, aes(x=long, y=lat, group=group), fill=NA, color="black", size=.1)
nj_map <- nj_map + geom_polygon(data=nj_hud, aes(x=long, y=lat, group=group, fill=factor(nj_mid)), color="black", size=.5)
nj_map <- nj_map + facet_wrap(~radius)
nj_map <- nj_map + coord_map()
nj_map <- nj_map + scale_fill_manual(drop=FALSE, values=c("blue", "yellow", "orange", "red", "green", "purple"), na.value="#EEEEEE", name="Unemployment rate")
nj_map <- nj_map + theme_nothing(legend=TRUE)
nj_map <- nj_map + labs(x=NULL, y=NULL, title="1 Journal Square (2017)")
nj_map <- nj_map + theme(panel.grid.major = element_line(colour = NA))
nj_map <- nj_map + theme(text = element_text(size=15))
nj_map <- nj_map + theme(plot.title=element_text(face="bold", hjust=.4))
nj_map <- nj_map + theme(plot.subtitle=element_text(face="italic", size=9, margin=margin(l=20)))
nj_map <- nj_map + theme(plot.caption=element_text(size=12, margin=margin(t=12), color="#7a7d7e", hjust=0))
nj_map <- nj_map + theme(legend.key.size = unit(1, "cm"))
nj_map <- nj_map + annotate("segment", x = -74.063644, xend = -74.035, y = 40.734330, yend = 40.734330, colour = "tomato", size=.5)
nj_map <- nj_map + annotate("point", x = -74.063644, y = 40.734330, colour = "lightblue", size = 1)
nj_map <- nj_map + annotate("text", x = -74.01, y = 40.734330, label = "1 Journal Square", size=3, colour="gray30")
print(nj_map)
# Saving output as a pdf
ggsave("map1b.pdf", nj_map, device="pdf")
# Mapping all census tract data in the county and layering on the project-only census tract data
nj_map <- ggplot()
nj_map <- nj_map + geom_polygon(data=nj_all, aes(x=long, y=lat, group=group, fill=factor(nj_mid)), color="gray", size=.3)
nj_map <- nj_map + geom_polygon(data=nj_hud, aes(x=long, y=lat, group=group, fill=factor(nj_mid)), color="black", size=.6)
nj_map <- nj_map + facet_wrap(~radius)
nj_map <- nj_map + coord_map()
nj_map <- nj_map + scale_fill_manual(drop=FALSE, values=c("blue", "yellow", "orange", "red", "green", "purple"), na.value="#EEEEEE", name="Unemployment rate")
nj_map <- nj_map + theme_nothing(legend=TRUE)
nj_map <- nj_map + labs(x=NULL, y=NULL, title="1 Journal Square (2017)")
nj_map <- nj_map + theme(panel.grid.major = element_line(colour = NA))
nj_map <- nj_map + theme(text = element_text(size=15))
nj_map <- nj_map + theme(plot.title=element_text(face="bold", hjust=.4))
nj_map <- nj_map + theme(plot.subtitle=element_text(face="italic", size=9, margin=margin(l=20)))
nj_map <- nj_map + theme(plot.caption=element_text(size=12, margin=margin(t=12), color="#7a7d7e", hjust=0))
nj_map <- nj_map + theme(legend.key.size = unit(1, "cm"))
nj_map <- nj_map + annotate("segment", x = -74.063644, xend = -74.035, y = 40.734330, yend = 40.734330, colour = "tomato", size=.5)
nj_map <- nj_map + annotate("point", x = -74.063644, y = 40.734330, colour = "lightblue", size = 1)
nj_map <- nj_map + annotate("text", x = -74.01, y = 40.734330, label = "1 Journal Square", size=3, colour="gray30")
print(nj_map)
ggsave("map1c.pdf", nj_map, device="pdf")
# Summarizing the table
# Unemployment rate for the census tracts picked by the project group
# versus unemployment rate for census tracts only adjacent to the project
sm_table <- hudson_unemployment_sm %>%
group_by(radius) %>%
summarize(average_unemployment=round(mean(un_rate, na.rm=T),2), median_unemployment=round(median(un_rate, na.rm=T),2)) %>%
arrange(average_unemployment)
kable(sm_table)
## OK, SAME AS ABOVE BUT NOW FOR A NEW PROJECT SITE
# 65 Bay Street (2015)
# Census tracts included by planners in this project
# GEOIDs: 34017004400, 34017004500, 34017004600, 34017005200, 34017005300, 34017005500, 34017005600, 34017005801, 34017006000, 34017006100, 34017006200, 34017006400, 34017006500, 34017006700, 34017006800, 34017007600
# Census tracts immediately adjacent to the project
# GEOIDS: 77, 78, 70, 64, 75, 74, 76
# Getting the unemployment figures and total labor workforce for all tracts in Hudson county
hudson_unemployment_2012 <- get_acs(geography="tract", endyear=2012, variables= c("B23025_005E", "B23025_002E"), county = "Hudson", state="NJ")
# Why 2008-2012? That's the ACS 5-year data that developers used to justify their project to officials, according to documents
# Discarding the margin of error
hudson_unemployment_2012$moe <- NULL
# Tidying up the structure of the new data frame
hudson_unemployment_2012 <- spread(hudson_unemployment_2012,variable, estimate )
# Calculating the unemployment rate
hudson_unemployment_2012$per_un <- round(hudson_unemployment_2012[,4]/hudson_unemployment_2012[,3]*100,2)
# Filtering out the blank census tracts
hudson_unemployment_2012 <- filter(hudson_unemployment_2012, GEOID!="34017006900" & GEOID!="34017980100")
# Array of group-designated census tracts
proj1 <- c("34017004400", "34017004500", "34017004600", "34017005200", "34017005300", "34017005500", "34017005600", "34017005801", "34017006000", "34017006100", "34017006200", "34017006400", "34017006500", "34017006700", "34017006800", "34017007600")
# Filtering out the tracts data specific to the project
hudson_unemployment_2012_sm <- filter(hudson_unemployment_2012, GEOID %in% proj1)
# Cleaning up the data frame
colnames(hudson_unemployment_2012_sm) <- c("GEOID", "name", "total", "unemployed", "unemp_rate")
hudson_unemployment_2012_sm$un_rate <- hudson_unemployment_2012_sm$unemp_rate$B23025_005
hudson_unemployment_2012_sm$unemp_rate <- NULL
# Adding a label column
hudson_unemployment_2012_sm$radius <- "gerrymandered"
# Array of census tracts immediately adjacent to the project
proj1_half <- c("34017007700", "34017007800", "34017007000", "34017006400", "34017007500", "34017007400", "34017007600")
# Filtering out the tracts data specific to the project
hudson_unemployment_2012_sm_half <- filter(hudson_unemployment_2012, GEOID %in% proj1_half)
# Cleaning up the data frame
colnames(hudson_unemployment_2012_sm_half) <- c("GEOID", "name", "total", "unemployed", "unemp_rate")
hudson_unemployment_2012_sm_half$un_rate <- hudson_unemployment_2012_sm_half$unemp_rate$B23025_005
hudson_unemployment_2012_sm_half$unemp_rate <- NULL
# Adding a label column
hudson_unemployment_2012_sm_half$radius <- "adjacent tracts"
# Combining the two census tracts data frames
hudson_unemployment_sm <- rbind(hudson_unemployment_2012_sm, hudson_unemployment_2012_sm_half)
# Cleaning up the census tract names
hudson_unemployment_sm$name <- gsub(",.*", "", hudson_unemployment_sm$name)
# Prepping the shape file to be joined with data
# (This process is necessary when visualizing with ggplot)
nj_hf <- fortify(nj_h, region="GEOID")
# Joining census data of project to shape file
nj_hud <- left_join(nj_hf, hudson_unemployment_sm, by=c("id"="GEOID"))
# Filter out the census tracts with no data
nj_hud <- filter(nj_hud, !is.na(radius))
# Slicing up the unemployment data (for project tracts) into buckets
nj_hud$nj_mid <- cut(nj_hud$un_rate, breaks = seq(0, 30, by=5))
#colors <- colorRampPalette(c("white", "red"))(length(levels(nj_hud$nj_mid)))
# Mapping project-only census tract data
nj_map <- ggplot()
nj_map <- nj_map + geom_polygon(data=nj_hud, aes(x=long, y=lat, group=group, fill=factor(nj_mid)), color="black", size=.5)
nj_map <- nj_map + facet_wrap(~radius)
nj_map <- nj_map + coord_map()
nj_map <- nj_map + scale_fill_manual(drop=FALSE, values=c("blue", "yellow", "orange", "red", "green", "purple"), na.value="#EEEEEE", name="Unemployment rate")
nj_map <- nj_map + theme_nothing(legend=TRUE)
nj_map <- nj_map + labs(x=NULL, y=NULL, title="65 Bay Street (2015)")
nj_map <- nj_map + theme(panel.grid.major = element_line(colour = NA))
nj_map <- nj_map + theme(text = element_text(size=15))
nj_map <- nj_map + theme(plot.title=element_text(face="bold", hjust=.4))
nj_map <- nj_map + theme(plot.subtitle=element_text(face="italic", size=9, margin=margin(l=20)))
nj_map <- nj_map + theme(plot.caption=element_text(size=12, margin=margin(t=12), color="#7a7d7e", hjust=0))
nj_map <- nj_map + theme(legend.key.size = unit(1, "cm"))
nj_map <- nj_map + annotate("segment", x = -74.03577, xend = -74.005, y = 40.72, yend = 40.72, colour = "tomato", size=.5)
nj_map <- nj_map + annotate("point", x = -74.03577, y = 40.72, colour = "lightblue", size = 1)
nj_map <- nj_map + annotate("text", x = -73.98217, y = 40.72, label = "65 Bay Street", size=3, colour="gray30")
print(nj_map)
ggsave( "map2a.pdf", nj_map, device="pdf")
# Mapping all census tract borders in the county and layering on the project-only census tract data
nj_map <- ggplot()
nj_map <- nj_map + geom_polygon(data=nj_hf, aes(x=long, y=lat, group=group), fill=NA, color="black", size=.1)
nj_map <- nj_map + geom_polygon(data=nj_hud, aes(x=long, y=lat, group=group, fill=factor(nj_mid)), color="black", size=.5)
nj_map <- nj_map + facet_wrap(~radius)
nj_map <- nj_map + coord_map()
nj_map <- nj_map + scale_fill_manual(drop=FALSE, values=c("blue", "yellow", "orange", "red", "green", "purple"), na.value="#EEEEEE", name="Unemployment rate")
nj_map <- nj_map + theme_nothing(legend=TRUE)
nj_map <- nj_map + labs(x=NULL, y=NULL, title="65 Bay Street (2015)")
nj_map <- nj_map + theme(panel.grid.major = element_line(colour = NA))
nj_map <- nj_map + theme(text = element_text(size=15))
nj_map <- nj_map + theme(plot.title=element_text(face="bold", hjust=.4))
nj_map <- nj_map + theme(plot.subtitle=element_text(face="italic", size=9, margin=margin(l=20)))
nj_map <- nj_map + theme(plot.caption=element_text(size=12, margin=margin(t=12), color="#7a7d7e", hjust=0))
nj_map <- nj_map + theme(legend.key.size = unit(1, "cm"))
nj_map <- nj_map + annotate("segment", x = -74.03577, xend = -74.005, y = 40.72, yend = 40.72, colour = "tomato", size=.5)
nj_map <- nj_map + annotate("point", x = -74.03577, y = 40.72, colour = "lightblue", size = 1)
nj_map <- nj_map + annotate("text", x = -73.98217, y = 40.72, label = "65 Bay Street", size=3, colour="gray30")
print(nj_map)
ggsave("map2b.pdf", nj_map, device="pdf")
# Mapping all census tract data in the county and layering on the project-only census tract data
nj_map <- ggplot()
nj_map <- nj_map + geom_polygon(data=nj_all, aes(x=long, y=lat, group=group, fill=factor(nj_mid)), color="gray", size=.3)
nj_map <- nj_map + geom_polygon(data=nj_hud, aes(x=long, y=lat, group=group, fill=factor(nj_mid)), color="black", size=.6)
nj_map <- nj_map + facet_wrap(~radius)
nj_map <- nj_map + coord_map()
nj_map <- nj_map + scale_fill_manual(drop=FALSE, values=c("blue", "yellow", "orange", "red", "green", "purple"), na.value="#EEEEEE", name="Unemployment rate")
nj_map <- nj_map + theme_nothing(legend=TRUE)
nj_map <- nj_map + labs(x=NULL, y=NULL, title="1 Journal Square (2017)")
nj_map <- nj_map + theme(panel.grid.major = element_line(colour = NA))
nj_map <- nj_map + theme(text = element_text(size=15))
nj_map <- nj_map + theme(plot.title=element_text(face="bold", hjust=.4))
nj_map <- nj_map + theme(plot.subtitle=element_text(face="italic", size=9, margin=margin(l=20)))
nj_map <- nj_map + theme(plot.caption=element_text(size=12, margin=margin(t=12), color="#7a7d7e", hjust=0))
nj_map <- nj_map + theme(legend.key.size = unit(1, "cm"))
nj_map <- nj_map + annotate("segment", x = -74.063644, xend = -74.035, y = 40.734330, yend = 40.734330, colour = "tomato", size=.5)
nj_map <- nj_map + annotate("point", x = -74.063644, y = 40.734330, colour = "lightblue", size = 1)
nj_map <- nj_map + annotate("text", x = -74.01, y = 40.734330, label = "1 Journal Square", size=3, colour="gray30")
print(nj_map)
ggsave("map2c.pdf", nj_map, device="pdf")
# Summarizing the table
# Unemployment rate for the census tracts picked by the project group
# versus unemployment rate for census tracts only adjacent to the project
sm_table <- hudson_unemployment_sm %>%
group_by(radius) %>%
summarize(average_unemployment=round(mean(un_rate, na.rm=T),2), median_unemployment=round(median(un_rate, na.rm=T),2)) %>%
arrange(average_unemployment)
kable(sm_table)