Major Assignment 2
Ilsu
2022-11-13
Section 1. Choose your Census Tracts.
Walkable: census tract with GEOID 13121001002
Reasons: This census tract covers the main part of Georgia Tech campus.
Georgia Tech feels quite walkable because of 1) a large proportion of
road segments with wide sidewalk and/or dedicated bike lanes that makes
walking on the sidewalk safe from cars, 2) a lot of stop signs that
prevent cars from speeding, 3) many car-free pedestrian pathways that
enable great connectivity among various spots on campus by walking or
biking, and 4) many people walking and/or biking.
Unwalkable: census tract with GEOID
13121011802
Reasons: This census tract is adjacent to the south-western side of
Georgia Tech campus. A railroad divides this area and Georgia Tech.
Major roads in the census tract have narrow, poorly maintained sidewalk
with no or minimal buffer between sidewalk and roadway. Residential
roads in the area are more poorly maintained and some of them do not
have sidewalk. Making things worse, there are many empty lots and
ill-maintained houses that would make feel walking unsafe.
Section 2. OSM, GSV, and computer vision.
Fill out the template to complete the script.
library(tidyverse)
library(tidycensus)
library(osmdata)
library(sfnetworks)
library(units)
library(sf)
library(tidygraph)
library(tmap)
library(here)Step 1. Get OSM data and clean it.
The getbb() function, which we used in the class
material to download OSM data, isn’t suitable for downloading just two
Census Tracts. We will instead use an alternative method.
- Using tidycensus package, download the Census Tract polygon for Fulton and DeKalb counties.
- Extract two Census Tracts, each of which will be your most walkable and least walkable Census Tracts.
- Using their bounding boxes, get OSM data.
- Convert them into sfnetworks data and clean it.
# TASK ////////////////////////////////////////////////////////////////////////
# 1. Set up your api key here
census.api.key <- read_rds(here("data", "census_api_key.rds"))
census_api_key(census.api.key)## To install your API key for use in future sessions, run this function with `install = TRUE`.# //TASK //////////////////////////////////////////////////////////////////////
# =========== NO MODIFICATION ZONE STARTS HERE ===============================
# Download Census Tract polygon for Fulton and DeKalb
tract <- get_acs("tract",
variables = c('tot_pop' = 'B01001_001'),
year = 2020,
state = "GA",
county = c("Fulton", "DeKalb"),
geometry = TRUE)## Getting data from the 2016-2020 5-year ACS## Downloading feature geometry from the Census website. To cache shapefiles for use in future sessions, set `options(tigris_use_cache = TRUE)`.##
|
| | 0%
|
| | 1%
|
|= | 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%
|
|===================== | 30%
|
|====================== | 31%
|
|======================= | 33%
|
|======================== | 34%
|
|========================= | 36%
|
|========================== | 37%
|
|=========================== | 39%
|
|============================ | 40%
|
|============================= | 41%
|
|============================== | 42%
|
|=============================== | 44%
|
|================================ | 45%
|
|================================= | 47%
|
|================================== | 48%
|
|================================== | 49%
|
|=================================== | 50%
|
|==================================== | 51%
|
|==================================== | 52%
|
|===================================== | 53%
|
|====================================== | 54%
|
|======================================= | 56%
|
|======================================== | 57%
|
|========================================= | 58%
|
|========================================== | 59%
|
|=========================================== | 61%
|
|============================================ | 62%
|
|============================================= | 64%
|
|============================================== | 65%
|
|=============================================== | 67%
|
|================================================ | 68%
|
|================================================= | 70%
|
|================================================== | 71%
|
|=================================================== | 73%
|
|==================================================== | 74%
|
|===================================================== | 75%
|
|====================================================== | 76%
|
|======================================================= | 78%
|
|======================================================= | 79%
|
|========================================================= | 81%
|
|========================================================= | 82%
|
|========================================================== | 83%
|
|=========================================================== | 84%
|
|=========================================================== | 85%
|
|============================================================= | 87%
|
|============================================================= | 88%
|
|=============================================================== | 90%
|
|=============================================================== | 91%
|
|================================================================ | 91%
|
|================================================================= | 92%
|
|================================================================= | 93%
|
|=================================================================== | 95%
|
|=================================================================== | 96%
|
|===================================================================== | 98%
|
|===================================================================== | 99%
|
|======================================================================| 100%# =========== NO MODIFY ZONE ENDS HERE ========================================
# TASK ////////////////////////////////////////////////////////////////////////
# The purpose of this TASK is to create one bounding box for walkable Census Tract and another bounding box for unwalkable Census Tract.
# As long as you generate what's needed for the subsequent codes, you are good. The numbered list of tasks below is to provide some hints.
# 1. Write the GEOID of walkable & unwalkable Census Tracts. e.g., tr1_ID <- c("13121001205", "13121001206")
# 2. Extract the selected Census Tracts using tr1_ID & tr2_ID
# 3. Create their bounding boxes using st_bbox(), and
# 4. assign them to tract_1_bb and tract_1_bb, respectively.
# For the walkable Census Tract(s)
# 1.
tr1_ID <- "13121001002"
# 2~4
tract_1_bb <- tract %>%
filter(GEOID == tr1_ID) %>%
st_bbox()
# For the unwalkable Census Tract(s)
# 1.
tr2_ID <- "13121011802"
# 2~4
tract_2_bb <- tract %>%
filter(GEOID == tr2_ID) %>%
st_bbox()
# //TASK //////////////////////////////////////////////////////////////////////
# =========== NO MODIFICATION ZONE STARTS HERE ===============================
# Get OSM data for the two bounding box
osm_1 <- opq(bbox = tract_1_bb) %>%
add_osm_feature(key = 'highway',
value = c("motorway", "trunk", "primary",
"secondary", "tertiary", "unclassified",
"residential")) %>%
osmdata_sf() %>%
osm_poly2line()
osm_2 <- opq(bbox = tract_2_bb) %>%
add_osm_feature(key = 'highway',
value = c("motorway", "trunk", "primary",
"secondary", "tertiary", "unclassified",
"residential")) %>%
osmdata_sf() %>%
osm_poly2line()
# =========== NO MODIFY ZONE ENDS HERE ========================================
# TASK ////////////////////////////////////////////////////////////////////////
# 1. Convert osm_1 and osm_2 to sfnetworks objects (set directed = FALSE)
# 2. Clean the network by (1) deleting parallel lines and loops, (2) create missing nodes, and (3) remove pseudo nodes,
# 3. Add a new column named length using edge_length() function.
net1 <- osm_1$osm_lines %>%
sfnetworks::as_sfnetwork(directed = FALSE) %>%
activate("edges") %>%
filter(!edge_is_multiple()) %>%
filter(!edge_is_loop()) %>%
convert(., sfnetworks::to_spatial_subdivision) %>%
convert(., sfnetworks::to_spatial_smooth) %>%
mutate(length = edge_length()) ## Warning: to_spatial_subdivision assumes attributes are constant over geometriesnet2 <- osm_2$osm_lines %>%
sfnetworks::as_sfnetwork(directed = FALSE) %>%
activate("edges") %>%
filter(!edge_is_multiple()) %>%
filter(!edge_is_loop()) %>%
convert(., sfnetworks::to_spatial_subdivision) %>%
convert(., sfnetworks::to_spatial_smooth) %>%
mutate(length = edge_length()) ## Warning: to_spatial_subdivision assumes attributes are constant over geometries# //TASK //////////////////////////////////////////////////////////////////////
# =========== NO MODIFICATION ZONE STARTS HERE ===============================
# OSM for the walkable part
edges_1 <- net1 %>%
# Extract 'edges'
st_as_sf("edges") %>%
# Drop redundant columns
select(osm_id, highway, length) %>%
# Drop segments that are too short (100m)
mutate(length = as.vector(length)) %>%
filter(length > 50) %>%
# Add a unique ID for each edge
mutate(edge_id = seq(1,nrow(.)),
is_walkable = "walkable")
# OSM for the unwalkable part
edges_2 <- net2 %>%
# Extract 'edges'
st_as_sf("edges") %>%
# Drop redundant columns
select(osm_id, highway, length) %>%
# Drop segments that are too short (100m)
mutate(length = as.vector(length)) %>%
filter(length > 50) %>%
# Add a unique ID for each edge
mutate(edge_id = seq(1,nrow(.)),
is_walkable = "unwalkable")
# Merge the two
edges <- bind_rows(edges_1, edges_2)
# =========== NO MODIFY ZONE ENDS HERE ========================================Step 2. Define a function that performs Step 3.
get_azi <- function(line){
# This function takes one edge (i.e., a street segment) as an input and
# outputs a data frame with four points (start, mid1, mid2, and end) and their azimuth.
# TASK ////////////////////////////////////////////////////////////////////////
# 1. From `line` object, extract the coordinates using st_coordinates() and extract the first two rows.
# 2. Use atan2() function to calculate the azimuth in degree.
# Make sure to adjust the value such that 0 is north, 90 is east, 180 is south, and 270 is west.
# 1.
start_p <- line %>%
st_coordinates() %>%
.[1:2, 1:2]
# 2
start_azi <- atan2(start_p[2,"X"] - start_p[1, "X"],
start_p[2,"Y"] - start_p[1, "Y"])*180/pi
# //TASK //////////////////////////////////////////////////////////////////////
# TASK ////////////////////////////////////////////////////////////////////////
# Repeat what you did above, but for last two rows (instead of the first two rows).
# Remember to flip the azimuth so that the camera would be looking at the street that's being measured
end_p <- line %>%
st_coordinates() %>%
.[(nrow(.)-1):nrow(.),1:2]
end_azi <- atan2(end_p[2,"X"] - end_p[1, "X"],
end_p[2,"Y"] - end_p[1, "Y"])*180/pi
end_azi <- if (end_azi < 180) {end_azi + 180} else {end_azi - 180}
# //TASK //////////////////////////////////////////////////////////////////////
# TASK ////////////////////////////////////////////////////////////////////////
# 1. From `line` object, use st_line_sample() function to generate points at 0.45 and 0.55 locations. These two points will be used to calculate the azimuth.
# 2. Use st_case() function to convert 'MULTIPOINT' object to 'POINT' object.
# 3. Extract coordinates using st_coordinates().
# 4. Use atan2() functino to Calculate azimuth.
# 5. Use st_line_sample() again to generate a point at 0.5 location and get its coordinates. This point will be the location at which GSV image will be downloaded.
mid_p <- line %>%
st_line_sample(sample = c(0.45, 0.55)) %>%
st_cast("POINT") %>%
st_coordinates()
mid_azi <- atan2(mid_p[2,"X"] - mid_p[1, "X"],
mid_p[2,"Y"] - mid_p[1, "Y"])*180/pi
mid_p <- line %>%
st_line_sample(sample = 0.5) %>%
st_coordinates() %>%
.[1,1:2]
# //TASK //////////////////////////////////////////////////////////////////////
# =========== NO MODIFICATION ZONE STARTS HERE ===============================
return(tribble(
~type, ~X, ~Y, ~azi,
"start", start_p[1,"X"], start_p[1,"Y"], start_azi,
"mid1", mid_p["X"], mid_p["Y"], mid_azi,
"mid2", mid_p["X"], mid_p["Y"], ifelse(mid_azi < 180, mid_azi + 180, mid_azi - 180),
"end", end_p[2,"X"], end_p[2,"Y"], end_azi))
# =========== NO MODIFY ZONE ENDS HERE ========================================
}Step 3. Apply the function to all street segments
We can apply get_azi() function to the edges object. We
finally append edges object to make use of the columns in
edges object (e.g., is_walkable column). When
you are finished with this code chunk, you will be ready to download GSV
images.
# TASK ////////////////////////////////////////////////////////////////////////
# Apply get_azi() function to all edges.
# Remember that you need to pass edges object to st_geometry()
# before you apply get_azi()
endp_azi <- edges %>%
st_geometry() %>%
map_df(get_azi)
# //TASK //////////////////////////////////////////////////////////////////////
# =========== NO MODIFICATION ZONE STARTS HERE ===============================
endp <- endp_azi %>%
bind_cols(edges %>%
st_drop_geometry() %>%
slice(rep(1:nrow(edges),each=4))) %>%
st_as_sf(coords = c("X", "Y"), crs = 4326, remove=FALSE) %>%
mutate(node_id = seq(1, nrow(.)))
# =========== NO MODIFY ZONE ENDS HERE ========================================Step 4. Define a function that formats request URL and download images.
get_image <- function(iterrow){
# This function takes one row of endp and downloads GSV image using the information from endp.
# TASK ////////////////////////////////////////////////////////////////////////
# Finish this function definition.
# 1. Extract required information from the row of endp, including
# type (i.e., start, mid1, mid2, end), location, heading, edge_id, node_id, source (i.e., outdoor vs. default) and key.
# 2. Format the full URL and store it in furl. Refer to this page: https://developers.google.com/maps/documentation/streetview/request-streetview
# 3. Format the full path (including the file name) of the image being downloaded and store it in fpath
type <- iterrow$type
location <- paste0(iterrow$Y %>% round(4), ",", iterrow$X %>% round(4))
heading <- iterrow$azi %>% round(1)
edge_id <- iterrow$edge_id
node_id <- iterrow$node_id
highway <- iterrow$highway
key <- read_rds(here("data", "google_api_key.rds"))[1]
furl <- glue::glue("https://maps.googleapis.com/maps/api/streetview?size=640x640&location={location}&heading={heading}&fov=90&pitch=0&key={key}")
fname <- glue::glue("GSV-nid_{node_id}-eid_{edge_id}-type_{type}-Location_{location}-heading_{heading}.jpg") # Don't change this code for fname
fpath <- here("data", "downloaded_image_assignment_2", fname)
# //TASK //////////////////////////////////////////////////////////////////////
# =========== NO MODIFICATION ZONE STARTS HERE ===============================
# Download images
if (!file.exists(fpath)){
download.file(furl, fpath, mode = 'wb')
}
# =========== NO MODIFY ZONE ENDS HERE ========================================
}Step 5. Download GSV images
Before you download GSV images, make sure
the row number of endp is not too large! The row number of
endp will be the number of GSV images you will be
downloading. Before you download images, always double-check your Google
Cloud Console’s Billing tab to make sure that you will not go above the
free credit of $200 each month. The price is $7 per 1000 images.
# =========== NO MODIFICATION ZONE STARTS HERE ===============================
# Loop!
for (i in seq(1,nrow(endp))){
get_image(endp[i,])
}
# =========== NO MODIFY ZONE ENDS HERE ========================================Step 6. Apply computer vision
Now, you need to upload the images you downloaded to Google Drive. You should upload the images to the same folder that we used in class - the ‘demo_images’ folder in the root directory of Google Drive. Then, use Google Colab to apply a semantic segmentation model called Pyramid Scene Parsing Network.
Step 7. Merging the processed data back to R
Once all of the images are processed and saved in your Google Drive as a CSV file, download the CSV file and merge it back to edges.
# Read the downloaded CSV file from Google Drive
pspnet <- read.csv(here("data", "downloaded_image_assignment_2", "seg_output.csv"))
# =========== NO MODIFICATION ZONE STARTS HERE ===============================
# Join the pspnet object back to endp object using node_id as the join key.
pspnet_nodes <- endp %>% inner_join(pspnet, by="node_id") %>%
select(type, X, Y, node_id, building, sky, tree, road, sidewalk) %>%
mutate(across(c(building, sky, tree, road, sidewalk), function(x) x/(640*640)))
# =========== NO MODIFY ZONE ENDS HERE ========================================Section 3. Summarise and analyze the results.
At the beginning of this assignment, you defined one Census Tract as walkable and the other as unwalkable. The key to the following analysis is the comparison between walkable/unwalkable Census Tracts.
Analysis 1 - Create map(s) to visualize the spatial distribution of the streetscape.
You need to create maps of the proportion of building, sky, tree, road, and sidewalk for walkable and unwalkable areas. In total, you will have 10 maps (two Census Tracts times five categories).
Below the maps, provide a brief description of your findings from the maps.
# TASK ////////////////////////////////////////////////////////////////////////
# Create map(s) to visualize the `pspnet_nodes` objects.
# As long as you can deliver the message clearly, you can use any format/package you want.
node.id.walkable <- endp %>%
filter(is_walkable == "walkable") %>% pull(node_id)
node.id.unwalkable <- endp %>%
filter(is_walkable == "unwalkable") %>% pull(node_id)
tmap_mode("view") %>% suppressMessages()
tmap_arrange(tm_shape(pspnet_nodes %>% filter(node_id %in% node.id.walkable)) +
tm_dots(col = "building", style = "fixed", breaks = seq(0, 0.6, 0.1)),
tm_shape(pspnet_nodes %>% filter(node_id %in% node.id.unwalkable)) +
tm_dots(col = "building", style = "fixed", breaks = seq(0, 0.6, 0.1)),
tm_shape(pspnet_nodes %>% filter(node_id %in% node.id.walkable)) +
tm_dots(col = "sky", style = "fixed", breaks = seq(0, 0.6, 0.1)),
tm_shape(pspnet_nodes %>% filter(node_id %in% node.id.unwalkable)) +
tm_dots(col = "sky", style = "fixed", breaks = seq(0, 0.6, 0.1)),
tm_shape(pspnet_nodes %>% filter(node_id %in% node.id.walkable)) +
tm_dots(col = "tree", style = "fixed", breaks = seq(0, 0.6, 0.1)),
tm_shape(pspnet_nodes %>% filter(node_id %in% node.id.unwalkable)) +
tm_dots(col = "tree", style = "fixed", breaks = seq(0, 0.6, 0.1)),
tm_shape(pspnet_nodes %>% filter(node_id %in% node.id.walkable)) +
tm_dots(col = "road", style = "fixed", breaks = seq(0, 0.6, 0.1)),
tm_shape(pspnet_nodes %>% filter(node_id %in% node.id.unwalkable)) +
tm_dots(col = "road", style = "fixed", breaks = seq(0, 0.6, 0.1)),
tm_shape(pspnet_nodes %>% filter(node_id %in% node.id.walkable)) +
tm_dots(col = "sidewalk", style = "fixed", breaks = seq(0, 0.6, 0.1)),
tm_shape(pspnet_nodes %>% filter(node_id %in% node.id.unwalkable)) +
tm_dots(col = "sidewalk", style = "fixed", breaks = seq(0, 0.6, 0.1)),
ncol = 2, nrow = 5, sync = F)# //TASK //////////////////////////////////////////////////////////////////////Findings from Analysis 1
The maps on the left are for the walkable tract that I chose and those on the right are for unwalkable tract. Colors in a map are designated using quintiles break points.
- The shares of building pixels are really low in both walkable and
unwalkable tracts without no big difference between two. However, it is
possible to see that the values are slightly lower for the unwalkable
tract probably because of a low building density and a large proportion
of buildings being low-rise compared to the walkable tract.
- The difference in the shares of sky pixels between two tracts is not
prominent.
- It seems that the shares of tree are slightly higher in the
unwalkable tract than those of the walkable tract. This is somewhat
expected considering a larger proportion of the unwalkable census tract
is residential roads.
- The shares of road pixels look higher in the walkable area because
roads are wider in general in the walkable tract (Georgia Tech) than the
unwalkable tract with most of its roads are residential.
- The shares of side walk are really low in both walkable and unwalkable tracts without no big difference between two.
Note that it is possible that my observations would vary (slightly, I hope) with the break points used in maps, which makes it a good idea to check mean values and box plots.
Analysis 2 - Compare the means.
You need to calculate the mean of the proportion of building, sky, tree, road, and sidewalk for walkable and unwalkable areas. For example, you need to calculate the mean of building category for each of walkable and unwalkable Census Tracts. Then, you need to calculate the mean of sky category for each of walkable and unwalkable Census Tracts. In total, you will have 10 mean values. After the calculation, provide a brief description of your findings.
# TASK ////////////////////////////////////////////////////////////////////////
# Perform the calculation as described above.
# As long as you can deliver the message clearly, you can use any format/package you want.
pspnet_nodes %>%
st_drop_geometry() %>%
mutate(is_walkable = case_when(node_id %in% node.id.walkable ~ "Walkable", TRUE ~ "Unwalkable")) %>%
group_by(is_walkable) %>%
summarise(
number_of_nodes = n(),
mean_building = mean(building),
mean_sky = mean(sky),
mean_tree = mean(tree),
mean_road = mean(road),
mean_sidewalk = mean(sidewalk)
) ## # A tibble: 2 × 7
## is_walkable number_of_nodes mean_building mean_sky mean_tree mean_road mean_…¹
## <chr> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Unwalkable 1152 0.0530 0.317 0.139 0.320 0.0345
## 2 Walkable 716 0.105 0.292 0.117 0.349 0.0413
## # … with abbreviated variable name ¹​mean_sidewalk# //TASK //////////////////////////////////////////////////////////////////////Findings from Analysis 2
The most prominent difference between the two tracts are in the mean share of building pixels. The value of walkable tract (10.5%) is almost double that of unwalkable tract (5.3%). The walkable tract has slightly higher mean shares of road and sidewalk while the unwalkable tract has slightly higher mean shares for sky and tree. I would expect the differences would be more prominent if Google had street-view images for many of pedestrian/bicycler pathways inside Georgia Tech campus (in the walkable tract), which is not the case as shown in the maps presented for Analysis 1.
Analysis 3 - Draw boxplot
You need to calculate the mean of the proportion of building, sky, tree, road, and sidewalk for walkable and unwalkable areas. For example, you need to calculate the mean of building category for each of walkable and unwalkable Census Tracts. Then, you need to calculate the mean of sky category for each of walkable and unwalkable Census Tracts. In total, you will have 10 mean values. After the calculation, provide a brief description of your findings.
# TASK ////////////////////////////////////////////////////////////////////////
# Create boxplot(s) using geom_boxplot() function from ggplot2 package.
# You may find the code from mini-assignment 4 useful here.
building <- ggplot(pspnet_nodes %>%
st_drop_geometry() %>%
mutate(is_walkable = case_when(node_id %in% node.id.walkable ~ "Walkable", TRUE ~ "Unwalkable"))) +
geom_boxplot(mapping = aes(x = is_walkable, y = building))
sky <- ggplot(pspnet_nodes %>%
st_drop_geometry() %>%
mutate(is_walkable = case_when(node_id %in% node.id.walkable ~ "Walkable", TRUE ~ "Unwalkable"))) +
geom_boxplot(mapping = aes(x = is_walkable, y = sky))
tree <- ggplot(pspnet_nodes %>%
st_drop_geometry() %>%
mutate(is_walkable = case_when(node_id %in% node.id.walkable ~ "Walkable", TRUE ~ "Unwalkable"))) +
geom_boxplot(mapping = aes(x = is_walkable, y = tree))
road <- ggplot(pspnet_nodes %>%
st_drop_geometry() %>%
mutate(is_walkable = case_when(node_id %in% node.id.walkable ~ "Walkable", TRUE ~ "Unwalkable"))) +
geom_boxplot(mapping = aes(x = is_walkable, y = road))
sidewalk <- ggplot(pspnet_nodes %>%
st_drop_geometry() %>%
mutate(is_walkable = case_when(node_id %in% node.id.walkable ~ "Walkable", TRUE ~ "Unwalkable"))) +
geom_boxplot(mapping = aes(x = is_walkable, y = sidewalk))
ggpubr::ggarrange(building, sky, tree, road, sidewalk,
labels = c("Building", "Sky", "Tree", "Road", "Sidewalk"),
ncol = 5, nrow = 1, label.x = 0.25, label.y = 1)
# //TASK //////////////////////////////////////////////////////////////////////Findings from Analysis 3
Similar to what I observed in Analysis 2, the most distinguished difference between two tracts is from the share of building pixels. The median of walkable tract is slightly larger than the third quartile of unwalkable tract. For all the others, the walkable tract has slightly higher median values for road and sidewalk and slightly lower meadina values for tree and sky, which aligns with what I observed in Analysis 2 when comparing mean values.