Install and Load Packages
The following code chunk installs the required R packages and loads them onto RStudio environment. sf, an R package specially designed to handle geospatial data in simple feature objects.
show
packages = c('tidyverse', 'sf', 'tmap', 'lubridate', 'clock',
'sftime', 'rmarkdown')
for (p in packages){
if(!require(p, character.only = T)){
install.packages(p)
}
library(p,character.only = T)
}
Importing wkt Data
Well-known text (WKT) is a human readable representation for spatial objects like points, lines, or enclosed areas on a map.
In the code chunk below, read_sf() of sf
package is used to parse School.csv into R as an sf data.frame.
show
schools <- read_sf('data/Schools.csv',
options= 'GEOM_POSSIBLE_NAMES=location')
buildings <- read_sf('data/Buildings.csv',
options= 'GEOM_POSSIBLE_NAMES=location')
apartments <- read_sf('data/Apartments.csv',
options= 'GEOM_POSSIBLE_NAMES=location')
employers <- read_sf('data/Employers.csv',
options= 'GEOM_POSSIBLE_NAMES=location')
restaurants <- read_sf('data/Restaurants.csv',
options= 'GEOM_POSSIBLE_NAMES=location')
pubs <- read_sf('data/Pubs.csv',
options= 'GEOM_POSSIBLE_NAMES=location')
Data Overview
After importing the data file into R, it is important for us to review the data object.
show
print(schools)
Simple feature collection with 4 features and 4 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: -4701.463 ymin: 1607.984 xmax: -376.7505 ymax: 6556.032
CRS: NA
# A tibble: 4 × 5
schoolId monthlyCost maxEnrollment location
<chr> <chr> <chr> <POINT>
1 0 12.81244502 242 (-376.7505 1607.984)
2 450 91.14351385 418 (-2597.448 3194.155)
3 900 38.00537955 394 (-2539.158 6556.032)
4 1350 73.19785215 384 (-4701.463 5141.763)
# … with 1 more variable: buildingId <chr>show
print(buildings)
Simple feature collection with 1042 features and 4 fields
Geometry type: POLYGON
Dimension: XY
Bounding box: xmin: -4762.191 ymin: -30.08359 xmax: 2650 ymax: 7850.037
CRS: NA
# A tibble: 1,042 × 5
buildingId location buildingType maxOccupancy
<chr> <POLYGON> <chr> <chr>
1 1 ((350.0639 4595.666, 390.0633… Commercial ""
2 2 ((-1926.973 2725.611, -1948.1… Residental "12"
3 3 ((685.6846 1552.131, 645.9985… Commercial ""
4 4 ((-976.7845 4542.382, -1053.2… Commercial ""
5 5 ((1259.306 3572.727, 1299.255… Residental "2"
6 6 ((478.8969 1082.484, 473.6596… Commercial ""
7 7 ((-1920.823 615.7447, -1960.8… Residental ""
8 8 ((-3302.657 5394.354, -3301.5… Commercial ""
9 9 ((-600.5789 4429.228, -495.95… Commercial ""
10 10 ((-68.75908 5379.924, -28.782… Residental "5"
# … with 1,032 more rows, and 1 more variable: units <chr>Plotting Building Footprint Map: Tmap
The code chunk below plots the building polygon features by using
tm_polygon().
show
tmap_mode('view')
tm_shape(buildings) +
tm_polygons(col = 'grey60',
size = 1,
border.col = 'black',
border.lwd = 1)
show
tmap_mode('plot')
Plotting A Composite Map
show
tm_shape(buildings)+
tm_polygons(col = "grey60",
size = 1,
border.col = "black",
border.lwd = 1) +
tm_shape(employers) +
tm_dots(col = '#003366', size = 0.5)
The Task: Plot a composite map by combining buildings, apartments, employers, pubs, restaurants, and schools.
show
tm_shape(buildings) +
tm_polygons(col = "grey60",
size = 1,
border.col = "black",
border.lwd = 1) +
tm_shape(employers) +
tm_dots(col = "#003366", size = 0.5) +
tm_shape(apartments) +
tm_dots(col = '#daa520', size = 0.3, alpha= 0.7 ) +
tm_shape(restaurants) +
tm_dots(col = '#ff7f50', size = 0.5) +
tm_shape(pubs) +
tm_dots(col = '#00ff00', size = 0.5) +
tm_shape(schools) +
tm_dots(col = '#20b2aa', size = 0.5)
show
logs <- read_sf('data/ParticipantStatusLogs1.csv',
options= 'GEOM_POSSIBLE_NAMES=currentLocation')
show
logs_selected <- logs %>%
mutate(Timestamp = date_time_parse(timestamp,zone = '',
format = '%Y-%m-%dT%H:%M:%S')) %>%
mutate(day = get_day(Timestamp)) %>%
filter(currentMode == 'Transport')
show
write_rds(logs_selected, 'data/rds/logs_selected.rds')
show
logs_selected <- read_rds('data/rds/logs_selected.rds')
show
hex <- st_make_grid(buildings,
cellsize=100,
square=FALSE) %>%
st_sf() %>%
rowid_to_column('hex_id')
plot(hex)
show
points_in_hex <- st_join(logs_selected,
hex,
join= st_within) %>%
st_set_geometry(NULL) %>%
count(name = 'pointCount', hex_id)
head(points_in_hex)
# A tibble: 6 × 2
hex_id pointCount
<int> <int>
1 169 35
2 212 56
3 225 21
4 226 94
5 227 22
6 228 45show
tm_shape(hex_combined %>%
filter(pointCount > 0))+
tm_fill("pointCount",
n = 8,
style = "quantile") +
tm_borders(alpha = 0.1)
show
logs_path <- logs_selected %>%
group_by(participantId, day) %>%
summarize(m = mean(Timestamp),
do_union=FALSE) %>%
st_cast("LINESTRING")
show
logs_path_selected <- logs_path %>%
filter(participantId == 0)
tm_shape(buildings)+
tm_polygons(col = "grey60",
size = 1,
border.col = "black",
border.lwd = 1) +
tm_shape(logs_path_selected) +
tm_lines(col = '#ffa500')
