Roadmap
The first thing we need to do to use GSV images for urban analytics is to prepare download points. The literature uses various methods to sample GSV images. For example, some studies sampled four images per street segments while some others downloaded four images (i.e., panorama) every 20 meters. We will be downloading four images per street segment, but you will also be introduced to how to do it in different ways.
NOTE: The method for sampling GSV images in this document is a simplified version.
library(tidyverse)
library(osmdata)
library(sfnetworks)
library(units)
library(sf)
library(tidygraph)
library(tmap)
library(here)Step 1. Get OSM data and clean it.
Download OSM data, convert it to sfnetworks object, and clean it. Detailed descriptions can be found in the previous week’s material (link).
# Bounding Box for Atlanta.
bb_atl <- getbb("Atlanta,GA")
# Get OSM data.
osm_road <- opq(bbox = bb_atl) %>%
add_osm_feature(key = 'highway',
value = c("motorway", "trunk", "primary",
"secondary", "tertiary", "unclassified",
"residential")) %>%
osmdata_sf() %>%
osm_poly2line()
# Convert the OSM line to sfnetworks and clean it.
net <- sfnetworks::as_sfnetwork(osm_road$osm_lines ,directed = FALSE) %>%
activate("edges") %>%
filter(!edge_is_multiple()) %>%
filter(!edge_is_loop()) %>%
convert(., sfnetworks::to_spatial_subdivision) %>%
convert(., sfnetworks::to_spatial_smooth) %>%
mutate(legnth = edge_length())## Warning: to_spatial_subdivision assumes attributes are constant over geometriesStep 2. Further cleaning specific to GSV
Extract “edges” from the cleaned network, add length column. Then delete segments that are too short (< 100m). Finally, add a unique ID for each edge.
edges <- net %>%
# Extract 'edges'
st_as_sf("edges") %>%
# Drop redundant columns
select(osm_id, highway) %>%
# Add length column
mutate(length = st_length(.) %>% unclass()) %>%
# Drop segments that are too short (100m)
filter(length > 50) %>%
# Add a unique ID for each edge
mutate(edge_id = seq(1,nrow(.)))
# nodes <- net %>%
# st_as_sf("nodes") %>%
# mutate(node_id = seq(1,nrow(.)))Step 3. Extract start, mid, and end points of each segment. For those points, calculate the azimuth.
We will download four images for street segment, one at the start of
a segment, two in the middle of a segment (each looking opposite to the
other), and one at the end of a segment. The code below, we will extract
a street segment (i.e., e below) and demonstrate how to
calculate azimuth for the start point, end point, and mid point.
# Select an edge for demo
test_edge = 23126
e <- edges %>% filter(edge_id == test_edge)
# View it
tmap_mode('view')## tmap mode set to interactive viewinge %>% st_coordinates() %>%
as.data.frame() %>%
st_as_sf(coords = c("X", "Y"), crs = 4326) %>%
with(
tm_shape(.) + tm_dots() +
tm_shape(.[1,]) + tm_dots(col = "red") + # Start = red
tm_shape(.[nrow(.),]) + tm_dots(col = "yellow") # End = yellow
)
# -----------------------------------------------------------
# First intersection
# First two points from a line
start_p <- e %>%
st_coordinates() %>%
.[1:2,1:2]
# Calculate the azimuth of the line connecting the two points
start_azi <- atan2(start_p[2,"X"] - start_p[1, "X"],
start_p[2,"Y"] - start_p[1, "Y"])*180/pi # 180/pi because trigonometry in R takes radians
# -----------------------------------------------------------
# The other intersection
# Last two points from a line
end_p <- e %>%
st_coordinates() %>%
.[(nrow(.)-1):nrow(.),1:2]
# Calculate the azimuth of the line connecting the two points
end_azi <- atan2(end_p[2,"X"] - end_p[1, "X"],
end_p[2,"Y"] - end_p[1, "Y"])*180/pi
# Flip the azimuth so that the camera would be looking back
end_azi <- if (end_azi < 180) {end_azi + 180} else {end_azi - 180}
# ----------------------------------------------------------
# mid point
mid_p <- e %>%
st_geometry() %>%
.[[1]] %>%
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 <- e %>%
st_geometry() %>%
.[[1]] %>%
st_line_sample(sample = 0.5) %>%
st_coordinates() %>%
.[1,1:2]Step 4. Define a function that performs Step 3.
get_azi <- function(line){
# end point 1 ----------------------------------------------
start_p <- line %>%
st_coordinates() %>%
.[1:2,1:2]
start_azi <- atan2(start_p[2,"X"] - start_p[1, "X"],
start_p[2,"Y"] - start_p[1, "Y"])*180/pi
# end point 2 ----------------------------------------------
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}
# mid point 1 ---------------------------------------------
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]
# return in data frame ------------------------------------
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))
}Step 5. Apply the function to all street segments
endp_azi <- edges %>%
st_geometry() %>%
map_df(get_azi)
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(.)))
endp## Simple feature collection with 161536 features and 9 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -84.61196 ymin: 33.56219 xmax: -84.10788 ymax: 34.06792
## Geodetic CRS: WGS 84
## # A tibble: 161,536 × 10
## type X Y azi osm_id highway length edge_id
## * <chr> <dbl> <dbl> <dbl> <chr> <chr> <dbl> <int>
## 1 start -84.3 33.6 177. 9106163 residential 109. 1
## 2 mid1 -84.3 33.6 177. 9106163 residential 109. 1
## 3 mid2 -84.3 33.6 357. 9106163 residential 109. 1
## 4 end -84.3 33.6 2.85 9106163 residential 109. 1
## 5 start -84.3 33.6 -178. 9106163 residential 110. 2
## 6 mid1 -84.3 33.6 -178. 9106163 residential 110. 2
## 7 mid2 -84.3 33.6 1.52 9106163 residential 110. 2
## 8 end -84.3 33.6 1.52 9106163 residential 110. 2
## 9 start -84.3 33.6 -178. 9106166 tertiary 78.9 3
## 10 mid1 -84.3 33.6 -178. 9106166 tertiary 78.9 3
## # … with 161,526 more rows, and 2 more variables: geometry <POINT [°]>,
## # node_id <int>
## # ℹ Use `print(n = ...)` to see more rows, and `colnames()` to see all variable names# st_write(endp, "C:/Users/bonwo/Downloads/endp.geojson", delete_dsn = T)
# st_write(edges, "C:/Users/bonwo/Downloads/edges.geojson", delete_dsn = T)Step 6. Define a function that formats request URL and download images.
get_image <- function(iterrow){
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 <- Sys.getenv("google_api")
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}-highway_{highway}.jpg")
fpath <- here("Lab", "module_3", "downloaded_image", fname)
download.file(furl, fpath, mode = 'wb')
}To reduce the number of images you’d be downloading, I pre-filtered the points so that it includes only Georgia Tech & Atlantic Station. Download endp_gt.geojson from here.
# Download images
edges_gt <- st_read(here("Lab", "module_3", "endp_atlantic_station.geojson"))
# Loop!
tic()
for (i in seq(1,nrow(edges_gt))){
get_image(edges_gt[i,])
}
toc()Another way to sample GSV images
Previously, we sampled four points per street segment. However, other studies have used different methods of sampling GSV images. One of the most widely used sampling method is to download one panoramic image at some fixed distance interval, such as 20 meters. This can be easily achieved by a small change to the code shown above.
Here, two sets of points are generated, (1) for actual location of points and (2) for the azimuth calculation. Notice that the points for azimuth calculation (i.e., 2) are twice as many as the points for actual location (i.e., 1).
# Extract a curvy line
z <- edges %>% filter(edge_id == test_edge)
tm_shape(z) + tm_lines()# calculate what proportion equals 40 in the given segment
prop <- 40/z$length
# Vector for actual points
sample_actual <- c(seq(0, 1, by=prop),1)
# Vector for two points before/after the actual points for azimuth calculation
sample_azi <- c(sample_actual-0.02, sample_actual+0.02) %>% sort()
# Sample points
sampled_actual <- z %>% st_transform(26967) %>% st_line_sample(sample = sample_actual) %>% st_cast("POINT")
sampled_azi <- z %>% st_transform(26967) %>% st_line_sample(sample = sample_azi) %>% st_cast("POINT")
# View them
tm_shape(sampled_actual) + tm_dots(col = "red", size=0.07) +
tm_shape(sampled_azi) + tm_dots(col = "blue", alpha = 0.5)# Calculate azimuth
point_azi <- sampled_actual %>%
st_sf() %>%
mutate(azi = NA)
for (i in 1:nrow(point_azi)){
j <- i*2 - 1
point_azi[i, "azi"] <- sampled_azi %>%
st_cast("POINT") %>%
.[j:(j+1)] %>%
st_coordinates() %>%
as.data.frame() %>%
with(
.[2,] - .[1,]
) %>%
with(
atan2(.[["X"]], .[["Y"]])*180/pi
)
}Merging the processed data back to R
# Download the output from the computer vision models
pspnet <- read.csv(here("Lab", "Module_3", "seg_output.csv"))
# Join them back to the template used to download images
pspnet_nodes <- endp %>% inner_join(pspnet, by="node_id")
# Map!
tm_shape(pspnet_nodes %>%
mutate(pct_tree = tree/(640*640))) +
tm_dots(col = "pct_tree", style="quantile")