Intro
The first thing we need to do to use GSV images for urban analytics is to sample points, prepare location data (i.e., coordinates and headings) and download the images. The literature uses various methods to sample points to download GSV images from. 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 learn how to do it in a different way at the end.
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)
library(progress)
library(nominatimlite)
ttm()Step 1. Select a study area.
Let’s find a small city that is suitable for the exercise. You can choose a city you like. I chose Chamblee.
cities_ga <- tigris::places('GA', year = 2024) %>%
filter(LSAD == 25)
metro_atl <- tigris::core_based_statistical_areas(year = 2024) %>%
filter(str_detect(NAME, 'Atlanta'))
cities_atl <- cities_ga[metro_atl,] %>%
mutate(temp_id = sample(1:n(), n()))tm_shape(metro_atl) + tm_borders(lwd = 3) +
tm_shape(cities_atl) + tm_polygons(col = 'temp_id', palette = 'inferno',
alpha = 0.5, lwd = 1.5) +
tm_shape(cities_atl) + tm_text('NAME') +
tm_view(set.view = c(-84.37, 33.75, 10))Step 2. Get OSM data and clean it.
2-1. Download OSM.
Download OSM data, convert it to an sfnetwork object, and clean it.
# Bounding Box for the city you chose.
my_bb <- nominatimlite::geo_lite_sf("Chamblee", points_only = F) %>%
st_bbox()
# Get OSM data.
osm_road <- opq(bbox = my_bb) %>%
add_osm_feature(key = "highway",
value = c("motorway", "trunk", "primary",
"secondary", "tertiary", "unclassified",
"residential")) %>%
osmdata_sf() %>%
osm_poly2line()
# Convert the OSM line to sfnetwork and clean it.
net <- osm_road$osm_lines %>%
select(osm_id, highway) %>%
sfnetworks::as_sfnetwork(directed = FALSE) %>%
activate("edges") %>%
filter(!edge_is_multiple()) %>% # remove duplicated edges
filter(!edge_is_loop()) %>% # remove loops
convert(., sfnetworks::to_spatial_subdivision) %>% # subdivide edges
convert(., sfnetworks::to_spatial_smooth) # delete pseudo nodes## Warning: to_spatial_subdivision assumes attributes are constant over geometries2-2. Extract edges.
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 > 100) %>%
# Add a unique ID for each edge
mutate(edge_id = seq(1,nrow(.)))Step 3. Get coordinates and headings.
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 illustration
test_edge_n <- 12
test_edge <- edges %>% filter(edge_id == test_edge_n) %>%
st_geometry()
# View it
test_edge %>% st_coordinates() %>%
as.data.frame() %>%
st_as_sf(coords = c("X", "Y"), crs = 4326) %>%
with(
tm_basemap("OpenStreetMap") +
tm_shape(.) + tm_dots() +
tm_shape(.[1,]) + tm_dots(col = "red") + # Start = red
tm_shape(.[nrow(.),]) + tm_dots(col = "yellow") # End = yellow
)
3-1. Sample points and get azimuth.
Extract start, mid, and end points of each segment. For those points, calculate the azimuth.
# -----------------------------------------------------------
# First intersection
# First two points from a line
start_p <- test_edge %>%
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 <- test_edge %>%
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_p3 <- test_edge %>%
.[[1]] %>%
st_line_sample(sample = c(0.45, 0.5, 0.55)) %>% # 0.45 and 0.55 are for azimuth calculation
st_cast("POINT") %>%
st_coordinates()
mid_p <- mid_p3[2,]
mid_azi <- atan2(mid_p3[2,"X"] - mid_p3[1, "X"],
mid_p3[2,"Y"] - mid_p3[1, "Y"])*180/pi
mid_azi2 <- ifelse(mid_azi < 180, mid_azi + 180, mid_azi - 180)3-2. Define a function.
Define a function that performs Step 3-1.
getAzimuth <- 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_p3 <- line %>%
st_line_sample(sample = c(0.45, 0.5, 0.55)) %>%
st_cast("POINT") %>%
st_coordinates()
mid_p <- mid_p3[2,]
mid_azi <- atan2(mid_p3[2,"X"] - mid_p3[1, "X"],
mid_p3[2,"Y"] - mid_p3[1, "Y"])*180/pi
mid_azi2 <- ifelse(mid_azi < 180, mid_azi + 180, mid_azi - 180)
# 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"], mid_azi2,
"end", end_p[2,"X"], end_p[2,"Y"], end_azi))
}3-3. Apply the function to edges.
Apply the function to all street segments.
edges_azi <- edges %>%
st_geometry() %>%
map_df(getAzimuth, .progress = T)
edges_azi <- edges_azi %>%
bind_cols(edges %>%
st_drop_geometry() %>% # same as `st_set_geometry(NULL)`
slice(rep(1:nrow(edges),each=4))) %>%
st_as_sf(coords = c("X", "Y"), crs = 4326, remove=FALSE) %>%
mutate(node_id = seq(1, nrow(.)))
st_write(edges_azi, paste0("C:/Users/benso/Documents/Georgia Tech/CP8883", "/", "edges_azi_chamblee.geojson"))## Writing layer `edges_azi_chamblee' to data source
## `C:/Users/benso/Documents/Georgia Tech/CP8883/edges_azi_chamblee.geojson' using driver `GeoJSON'
## Writing 5364 features with 9 fields and geometry type Point.edges_azi## Simple feature collection with 5364 features and 9 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -84.33123 ymin: 33.83633 xmax: -84.26031 ymax: 33.92455
## Geodetic CRS: WGS 84
## # A tibble: 5,364 × 10
## type X Y azi osm_id highway length edge_id
## * <chr> <dbl> <dbl> <dbl> <chr> <chr> <dbl> <int>
## 1 start -84.3 33.9 -162. 9164345 residential 110. 1
## 2 mid1 -84.3 33.9 69.4 9164345 residential 110. 1
## 3 mid2 -84.3 33.9 249. 9164345 residential 110. 1
## 4 end -84.3 33.9 205. 9164345 residential 110. 1
## 5 start -84.3 33.9 124. 9164886 residential 128. 2
## 6 mid1 -84.3 33.9 114. 9164886 residential 128. 2
## 7 mid2 -84.3 33.9 294. 9164886 residential 128. 2
## 8 end -84.3 33.9 284. 9164886 residential 128. 2
## 9 start -84.3 33.9 20.0 9165137 residential 118. 3
## 10 mid1 -84.3 33.9 10.2 9165137 residential 118. 3
## # ℹ 5,354 more rows
## # ℹ 2 more variables: geometry <POINT [°]>, node_id <int>Export the result
edges_azito a geojson file (e.g., edges_azi_{city_name}.geojson) usingst_write(). We will need this data later.
You may as well use the Chamblee data I created.
# Download images
edges_azi <- st_read("https://raw.githubusercontent.com/ujhwang/urban-analytics-2024/main/Lab/module_4/edges_azi_chamblee.geojson")Step 4. Get GSV images.
4-1. Define a function.
Define a function that formats request URL and download images.
getImage <- function(iterrow){
type = iterrow$type
location <- paste0(iterrow$Y %>% round(5), ",", iterrow$X %>% round(5))
heading <- iterrow$azi %>% round(1)
edge_id <- iterrow$edge_id
node_id <- iterrow$node_id
highway <- iterrow$highway
key <- Sys.getenv("google_api") # your Google API key
endpoint <- "https://maps.googleapis.com/maps/api/streetview"
request <- glue::glue("{endpoint}?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("GSV_Images", fname) # you may need to change this
download.file(request, fpath, mode = 'wb')
}4-2. Download the images.
Download the images using the custom function
getImage(). KEEP IN MIND THAT THIS WILL TAKE
A SIGNIFICANT AMOUNT OF TIME TO COMPLETE. FOR CHAMBLEE, IT TOOK ALMOST
30 MINUTES.ADDITIONALLY, WHEN RUNNING
THE FOLLOWING CODE CHUNK, MONITOR YOUR API USAGE IN THE GOOGLE DEVELOPER
CONSOLE TO AVOID UNEXPECTED CHARGES
# Create a progress bar to monitor the progress
pb <- txtProgressBar(min = 0, max = nrow(edges_azi), style = 3)## | | | 0%# Loop!
for (i in seq(1, nrow(edges_azi))){
getImage(edges_azi[i,])
setTxtProgressBar(pb, i)
}## | | | 1% | |= | 1% | |= | 2% | |== | 2% | |== | 3% | |== | 4% | |=== | 4% | |=== | 5% | |==== | 5% | |==== | 6% | |===== | 6% | |===== | 7% | |===== | 8% | |====== | 8% | |====== | 9% | |======= | 9% | |======= | 10% | |======= | 11% | |======== | 11% | |======== | 12% | |========= | 12% | |========= | 13% | |========= | 14% | |========== | 14% | |========== | 15% | |=========== | 15% | |=========== | 16% | |============ | 16% | |============ | 17% | |============ | 18% | |============= | 18% | |============= | 19% | |============== | 19% | |============== | 20% | |============== | 21% | |=============== | 21% | |=============== | 22% | |================ | 22% | |================ | 23% | |================ | 24% | |================= | 24% | |================= | 25% | |================== | 25% | |================== | 26% | |=================== | 26% | |=================== | 27% | |=================== | 28% | |==================== | 28% | |==================== | 29% | |===================== | 29% | |===================== | 30% | |===================== | 31% | |====================== | 31% | |====================== | 32% | |======================= | 32% | |======================= | 33% | |======================= | 34% | |======================== | 34% | |======================== | 35% | |========================= | 35% | |========================= | 36% | |========================== | 36% | |========================== | 37% | |========================== | 38% | |=========================== | 38% | |=========================== | 39% | |============================ | 39% | |============================ | 40% | |============================ | 41% | |============================= | 41% | |============================= | 42% | |============================== | 42% | |============================== | 43% | |============================== | 44% | |=============================== | 44% | |=============================== | 45% | |================================ | 45% | |================================ | 46% | |================================= | 46% | |================================= | 47% | |================================= | 48% | |================================== | 48% | |================================== | 49% | |=================================== | 49% | |=================================== | 50% | |=================================== | 51% | |==================================== | 51% | |==================================== | 52% | |===================================== | 52% | |===================================== | 53% | |===================================== | 54% | |====================================== | 54% | |====================================== | 55% | |======================================= | 55% | |======================================= | 56% | |======================================== | 56% | |======================================== | 57% | |======================================== | 58% | |========================================= | 58% | |========================================= | 59% | |========================================== | 59% | |========================================== | 60% | |========================================== | 61% | |=========================================== | 61% | |=========================================== | 62% | |============================================ | 62% | |============================================ | 63% | |============================================ | 64% | |============================================= | 64% | |============================================= | 65% | |============================================== | 65% | |============================================== | 66% | |=============================================== | 66% | |=============================================== | 67% | |=============================================== | 68% | |================================================ | 68% | |================================================ | 69% | |================================================= | 69% | |================================================= | 70% | |================================================= | 71% | |================================================== | 71% | |================================================== | 72% | |=================================================== | 72% | |=================================================== | 73% | |=================================================== | 74% | |==================================================== | 74% | |==================================================== | 75% | |===================================================== | 75% | |===================================================== | 76% | |====================================================== | 76% | |====================================================== | 77% | |====================================================== | 78% | |======================================================= | 78% | |======================================================= | 79% | |======================================================== | 79% | |======================================================== | 80% | |======================================================== | 81% | |========================================================= | 81% | |========================================================= | 82% | |========================================================== | 82% | |========================================================== | 83% | |========================================================== | 84% | |=========================================================== | 84% | |=========================================================== | 85% | |============================================================ | 85% | |============================================================ | 86% | |============================================================= | 86% | |============================================================= | 87% | |============================================================= | 88% | |============================================================== | 88% | |============================================================== | 89% | |=============================================================== | 89% | |=============================================================== | 90% | |=============================================================== | 91% | |================================================================ | 91% | |================================================================ | 92% | |================================================================= | 92% | |================================================================= | 93% | |================================================================= | 94% | |================================================================== | 94% | |================================================================== | 95% | |=================================================================== | 95% | |=================================================================== | 96% | |==================================================================== | 96% | |==================================================================== | 97% | |==================================================================== | 98% | |===================================================================== | 98% | |===================================================================== | 99% | |======================================================================| 99% | |======================================================================| 100%ZIP THE DOWNLOADED IMAGES AND NAME IT ‘gsv_images.zip’ FOR THE SUBSEQUENT ANALYSIS.
(Optional) Another way to sample points.
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 are twice as many as the points for actual location.
# Extract a curvy line
z <- edges %>% filter(edge_id == test_edge_n)
tm_basemap("OpenStreetMap") + 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(32616) %>% st_line_sample(sample = sample_actual) %>% st_cast("POINT")
sampled_azi <- z %>% st_transform(32616) %>% st_line_sample(sample = sample_azi) %>% st_cast("POINT")
# View them
tm_basemap("OpenStreetMap") +
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
)
}