Chapter 4 GIS Analysis Toolkit
This chapter requires the following packages:
library(sf)
library(readr)
library(readxl)
library(tidyr)
library(cancensus)
library(osmdata)
library(mapview)
library(dplyr)The previous chapter was dedicated to getting the baseline skills necessary to load spatial data into R and understand the different data structures. In this section we are going to go over the basic geoprocessing operations that are commonly encountered in applied research with spatial data. Geoprocessing is a commonly used term for describing a suite of tools and operations that manipulate an existing spatial dataset to create new or augmented versions of the original data.
4.1 Attribute Queries
One of the most basic spatial data operations is to create subsets based on the values of an attribute (i.e a column). Let’s retrieve Metro Vancouver data using the cancensus package that was introduced in the previous chapter, and examine the attributes associated with each dissemination area. Note that the ID for Metro Vancouver is 5915. Feel free to augment the code below to retrieve a different subset of Canadian census data.
find_census_vectors("Median total income",type = "total", dataset = "CA16",query_type = "semantic")#options(cancensus.api_key = "your_api_key") # uncomment and input your API Key here to enable retrieval of data from Statistics Canada## # A tibble: 9 x 4
## vector type label details
## <chr> <fct> <chr> <chr>
## 1 v_CA16_2~ Total Median total income in 2015 among recipients ($) Income; Individuals; Total - Income statistics in 2015 for the population aged 15 years and over in private house~
## 2 v_CA16_2~ Total Median total income of households in 2015 ($) Income; Households; Total - Income statistics in 2015 for private households by household size - Median total inc~
## 3 v_CA16_2~ Total Median total income of one-person households in 2015 ($) Income; Households; Total - Income statistics in 2015 for one-person private households - Median total income of ~
## 4 v_CA16_2~ Total Median total income of two-or-more-person households in 20~ Income; Households; Total - Income statistics in 2015 for two-or-more-person private households - Median total in~
## 5 v_CA16_2~ Total Median total income of economic families in 2015 ($) Income; Families; Total - Income statistics in 2015 for economic families in private households - Median total in~
## 6 v_CA16_2~ Total Median total income of couple economic families without ch~ Income; Families; Total - Income statistics in 2015 for couple economic families without children or other relati~
## 7 v_CA16_2~ Total Median total income of couple economic families with child~ Income; Families; Total - Income statistics in 2015 for couple economic families with children in private househo~
## 8 v_CA16_2~ Total Median total income of lone-parent economic families in 20~ Income; Families; Total - Income statistics in 2015 for lone-parent economic families in private households - Med~
## 9 v_CA16_2~ Total Median total income in 2015 for persons aged 15 years and ~ Income; Families; Total - Income statistics in 2015 for persons aged 15 years and over not in economic families i~mv_census <- get_census(dataset='CA16', #specify the census of interest
regions=list(CSD="5915"),#specify the region of interest
vectors=c("v_CA16_2397"), # specify the attribute data of interest
labels = "short",# here we specify short variabhle names, otherwise they get very long
level='DA', #specify the census geography
geo_format = "sf", #specify the geographic data format as output
quiet = TRUE #turn off download progress messages (set to FALSE to see messages)
)
mv_census## Simple feature collection with 3450 features and 13 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -123.4319 ymin: 49.00192 xmax: -122.4086 ymax: 49.57428
## Geodetic CRS: WGS 84
## First 10 features:
## Population Households GeoUID Type CD_UID Shape Area Dwellings CSD_UID CT_UID CMA_UID Region Name Area (sq km) v_CA16_2397 geometry
## 1 368 159 59150004 DA 5915 0.85230 172 5915055 9330133.01 59933 West Vancouver 0.85230 94592 MULTIPOLYGON (((-123.2421 4...
## 2 503 205 59150005 DA 5915 0.16598 221 5915055 9330133.01 59933 West Vancouver 0.16598 102997 MULTIPOLYGON (((-123.2789 4...
## 3 432 196 59150006 DA 5915 0.15593 200 5915055 9330133.01 59933 West Vancouver 0.15593 85504 MULTIPOLYGON (((-123.2773 4...
## 4 546 212 59150007 DA 5915 0.71085 221 5915055 9330133.01 59933 West Vancouver 0.71085 133803 MULTIPOLYGON (((-123.289 49...
## 5 596 224 59150008 DA 5915 0.65167 247 5915055 9330133.01 59933 West Vancouver 0.65167 167936 MULTIPOLYGON (((-123.2815 4...
## 6 399 136 59150009 DA 5915 0.56986 143 5915055 9330133.01 59933 West Vancouver 0.56986 146432 MULTIPOLYGON (((-123.2797 4...
## 7 831 298 59150010 DA 5915 0.80527 314 5915055 9330133.01 59933 West Vancouver 0.80527 144725 MULTIPOLYGON (((-123.2766 4...
## 8 512 173 59150012 DA 5915 0.44894 182 5915055 9330133.02 59933 West Vancouver 0.44894 129792 MULTIPOLYGON (((-123.2429 4...
## 9 360 118 59150013 DA 5915 0.27287 122 5915055 9330133.02 59933 West Vancouver 0.27287 163328 MULTIPOLYGON (((-123.2429 4...
## 10 1070 413 59150014 DA 5915 32.12439 455 5915055 9330133.02 59933 West Vancouver 32.12439 123136 MULTIPOLYGON (((-123.195 49...Notice how our census data has attributes that describe geographical regions we are interested in via the column Region Name. While there are several ways to subset rows of data (spatial or otherwise) based on the values of columns, the filter() function from the dplyr package integrates nicely with sf objects, is intuitive and easy to implement. In the filter function we can construct a boolean expression (a statement that evaluates to either TRUE or FALSE) to determine which rows (which with sf objects, correspond to a spatial feature) we want to keep in our new dataset. Example functions and operators include:
==(is equal to)>(greater than)<(less than)is.na()(is equal to NA[ no value])
We can combine these functions using operators such as
&(and, to combine)|(or)!(does not equal)
Let’s try an attribute query that will isolate (or filter) the dissemination areas that are within the City of Vancouver and has a median household income > $175,000.
cov_high_income <- mv_census %>%
filter(`Region Name`=="Vancouver" & v_CA16_2397 >175000)In the above code we are stating that we want to create a new spatial dataset called cov_wealthy. THis new object is based on the original dataset mv_census where, using the filter() function, we only keep the rows (corresponding to a polygon) where the Region Name column value is equal to the string "Vancouver" AND the column describing the median household income v_CA16_2397 is greater than $175,000. Let’s see if this worked by quickly mapping both the original data and the new data using mapview(). Click on the red polygons and check to see if the value for v_CA16_2397 is greater than 175000.
mapview(mv_census) + mapview(cov_high_income,col.regions = "red")
We can see that this attribute query has worked perfectly and that the layer coloured in red only contains dissemination areas within the City of Vancouver and where the median household income is greater than $175,000.
Now let’s do a slightly more complex attribute query using data retrieved from OpenStreetMaps using the osmdata package where we only want separated bicycle infrastructure data (E.g. lines representing the location of separated infrastructure). First let’s load data for the city
van_bb <- getbb("Vancouver BC")
streets <- van_bb %>%
opq()%>%
add_osm_feature(key = "highway") %>%
osmdata_sf()
streets_sf <- streets$osm_linesIn the osm data there are several attributes (columns) that could be used to designate separated bicycling infrastructure including the “highway” column and the “bicycle” column. Here we construct a query that will select lines based on the different values that could be plausibly be used to describe separated bicycling infrastructure. In plain engligh this looks like.
Select the features where the following is TRUE:
highwayequals “cycleway” ORhighwayequals “path” &bicycleequals “yes” ORhighwayequals “path” &bicycleequals “designated” ORhighwayequals “footway” &bicycleequals “yes” ORhighwayequals “footway” &bicycleequals “designated”
bicycle <- streets_sf %>%
filter(highway=="cycleway" | (highway=="path"&bicycle=="yes") | (highway=="path"&bicycle=="designated") | (highway=="footway"&bicycle=="yes") | (highway=="footway"&bicycle=="designated"))
mapview(streets_sf) + mapview(bicycle,color = "red") 
In the above map we can see our selection highlighted in red.
4.2 Spatial Queries
What happens if we wanted to isolate spatial data based on the location, but the data don’t have existing attributes that are associated with that location? For example, what if our census data did not have the Region Name column? One way we could solve this problem would be to find another spatial dataset that represents the location of interest, and perform what is called a spatial query. Spatial queries refer to creating new data or subsets of data based purely on the location of different features, in relation to another dataset. For example, if we wanted to select all the crashes within a certain distance of specific infrastructure types, we’d need to perform a spatial query. More formally, spatial queries are typically based on either topological spatial relationships or distance relationships. Topological relationships describe the relationship in terms of whether any two spatial data a and b intersects, touches, contains among others. While distance relationships are about how far away certain objects are from each other. Let’s get the flavour for a spatial query by going through a few examples.
First, let’s answer the first question I posed in the above paragraph. How would we isolate the Dissmimination Areas within the City of Vancouver, if our dataset did not have an existing attribute that describes which region the DA fell in? We will need another spatial dataset which represents the boundary of the City of Vancouver. Such a shapefile was included in the course data package in Chapter 1.
van_boundary <- read_sf("D:/Github/crash-course-gis-r/data/CoV_boundary.shp") # use the file location of where you have stored your data
mapview(mv_census) + mapview(van_boundary,col.regions = "red")We can see from our map above that we want to keep those DAs that are within the red polygon (the City of Vancouver boundary). We can use functions from sf to evaluate the topological spatial relationships between these two datasets (See https://r-spatial.github.io/sf/reference/geos_binary_pred.html for a full list of operations). We are going to use the function st_interscts() to evaluate which polygons from the mv_census dataset intersect with the spatial features of the van_boundary dataset. First let’s just run the function setting the argument “sparse” = FALSE and see what the output is. We are setting the arugment sparse to false for teaching purposes, but I encourage you to read (Chapter 4)[https://geocompr.robinlovelace.net/spatial-operations.html#introduction-1] of Geocomputation with R for a greater background on the technical details.
st_intersects(mv_census,van_boundary,sparse = FALSE)## Error in st_geos_binop("intersects", x, y, sparse = sparse, prepared = prepared, : st_crs(x) == st_crs(y) is not TRUEShoot! We get an error message indicating that our two datasets are in two different coordinate systems! This is an important lesson to learn. In order to perform any kind of spatial operation involving more than one dataset, we need to ensure they are both in the same coordinate system.
Let’s troubleshoot this. First let’s check the CRS information of each dataset and compare.
st_crs(mv_census)#check CRS## Coordinate Reference System:
## User input: 4326
## wkt:
## GEOGCS["WGS 84",
## DATUM["WGS_1984",
## SPHEROID["WGS 84",6378137,298.257223563,
## AUTHORITY["EPSG","7030"]],
## AUTHORITY["EPSG","6326"]],
## PRIMEM["Greenwich",0,
## AUTHORITY["EPSG","8901"]],
## UNIT["degree",0.0174532925199433,
## AUTHORITY["EPSG","9122"]],
## AXIS["Latitude",NORTH],
## AXIS["Longitude",EAST],
## AUTHORITY["EPSG","4326"]]st_crs(van_boundary)#check CRS## Coordinate Reference System:
## User input: NAD83 / UTM zone 10N
## wkt:
## PROJCRS["NAD83 / UTM zone 10N",
## BASEGEOGCRS["NAD83",
## DATUM["North American Datum 1983",
## ELLIPSOID["GRS 1980",6378137,298.257222101,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4269]],
## CONVERSION["UTM zone 10N",
## METHOD["Transverse Mercator",
## ID["EPSG",9807]],
## PARAMETER["Latitude of natural origin",0,
## ANGLEUNIT["Degree",0.0174532925199433],
## ID["EPSG",8801]],
## PARAMETER["Longitude of natural origin",-123,
## ANGLEUNIT["Degree",0.0174532925199433],
## ID["EPSG",8802]],
## PARAMETER["Scale factor at natural origin",0.9996,
## SCALEUNIT["unity",1],
## ID["EPSG",8805]],
## PARAMETER["False easting",500000,
## LENGTHUNIT["metre",1],
## ID["EPSG",8806]],
## PARAMETER["False northing",0,
## LENGTHUNIT["metre",1],
## ID["EPSG",8807]]],
## CS[Cartesian,2],
## AXIS["(E)",east,
## ORDER[1],
## LENGTHUNIT["metre",1]],
## AXIS["(N)",north,
## ORDER[2],
## LENGTHUNIT["metre",1]],
## ID["EPSG",26910]]Based on the CRS information above, we see that the mv_census dataset is “unprojected” and is in WGS1984, while van_boundary is in the projected coordinated system of NAD83/UTM Zone 10N. Since we need our data to be projected, let’s tranform the mv_census object to be in the same projection as van_boundary.
crs_to_use <- st_crs(van_boundary) # store crs information from the van_boundary dataset
crs_epsg <- crs_to_use$epsg # store the epsg code from the van_boundary dataset
mv_census_proj <- st_transform(mv_census,crs=crs_epsg)Now let’s check to see if they are both in the same coordinate system
st_crs(van_boundary)==st_crs(mv_census_proj)## [1] TRUEGreat! We are rolling now. Let’s try the spatial operation of st_intersects() again, and see what the output is
st_intersects_output <- st_intersects(mv_census_proj,van_boundary,sparse = FALSE)
head(st_intersects_output)## [,1]
## [1,] FALSE
## [2,] FALSE
## [3,] FALSE
## [4,] FALSE
## [5,] FALSE
## [6,] FALSEWhen we set the argument “sparse” = FALSE, what is returned is matrix where each row corresponds to a spatial feature within the mv_census_proj object and if the value is TRUE than it means that the feature intersects with the van_boundary geometry and if FALSE than it means it does not. We can combine this output with the the filter() function to return only the polygons that intersect with the van_boundary polygon.
census_intersects <- filter(mv_census_proj,st_intersects(x = mv_census_proj, y = van_boundary, sparse = FALSE))
mapview(van_boundary,col.region = "red") + mapview(census_intersects) Another important lesson to learn here! It looks like our spatial query returned the dissemination areas just on the borders of the van_boundary file. This is a common occurrence in GIS analysis where boundaries from different data sources don’t exactly line up. You can see that the borders of the red polygon representing the City of Vancouver boundaries don’t line up with borders of our dissemination areas, causing us to select polygons that intersect with the Vancouver boundary layer, but are technically in Richmond, Burnaby and UBC. This is why it is important to map out your geoprocessing to do a spot check on whether what you wanted to do, actually occurred.
Now let’s brainstorm on how we can fix this problem. Short of manually fixing the boundaries to ensure that they line up with the dissemination area boundaries, we need to get a bit creative. I’m going to introduce you to the concept of a centroid, something that can be quite handy in situations like these. A centroid in GIS speak, is the “average” coordinate within a polygon, that is a single point (x,y location) that is representative of the location of a polygon.
mv_centroids <- st_centroid(mv_census_proj)
mapview(van_boundary,col.regions="red") + mapview(mv_census_proj) + mapview(mv_centroids,col.regions="black") Make sure to explore this map by zooming in and clicking on the points and the polygons they fall on. See if you can find any polygons where the centroid doesn’t fall within the boundaries of its polygon. Zoom into areas around the border of the red polygon and see if any of the centroid fall within it.
When we are working with two polygon datasets whose boundaries don’t align in the way we want, often using the centroids to represent a polygon location can be a good solution to this problem. Instead of selecting disseminiation area polygons whose boundaries intersect with the Vancouver boundary, let’s select the dissmeination area polygons whose centroid intersects with the Vancouver boundary layer. We can do this by combining filter(), st_intersect() and st_centroid().
census_intersects_cntrd <- filter(mv_census_proj,st_intersects(x = st_centroid(mv_census_proj), y = van_boundary, sparse = FALSE))
mapview(van_boundary,col.region = "red") + mapview(census_intersects_cntrd) There we have it. Explore the output in the map to ensure that we have selected the correct polygons.
This section was designed to give you some intuition as to how to perform a spatial query. There are lots of different ways of conducting spatial queries to produce the spatial dataset you need. Check out a (full list of functions)[https://r-spatial.github.io/sf/reference/geos_binary_pred.html] to give you an idea of what other topological and distance relationships you can query using sf.
4.3 Attribute Joins
Often data are provided in tabular format that are associated with a certain geographic location such as census geographies (census tracts, dissemination areas etc.) or neighbourhoods, or other geographic identifiers. For example in your downloaded data there is an excel spread sheet called bc_scores_quintiles_EN.xlsx which is originally from Statistics Canada that provides a score for different dimensions deprivation for each dissemination area in the province of British Columbia. In the figure below the column that refers to a unique identifier for each dissemination area in the province is named “PRCDDA.”
Figure 4.1: Tabular data with unique geographic identifier
That same variable is available within the census data we have previously retrieved in the chapter. As can be seen in the figure below, within the census data retrieved using cancensus the dissemination area unique identifier is named “GeoUID.” With this knowledge in mind we can join the spread sheet data that lacks any geometry to our sf census data using this unique identifier.
Figure 4.2: Geographic data with unique geographic identifier
First we need to use the function read_xlsx() from the readxl package to load in the .xlsx excel spreadsheet data.
deprivation_da <- read_xlsx("D:/Github/crash-course-gis-r/data/bc_scores_quintiles_EN.xlsx",sheet = 1)Then we use the function left_join() from the dplyrpackage. The left_join() function takes two datasets and returns all the rows from x and only the rows from y where there is a matching value.
](images/left_join.png)
Figure 4.3: Excample of left join operation for two datasets.source
While there are many kind of ways to join tables together based on common attributes (see here), here we are interested in ensuring that we keep all of the geographic data and only the data on deprivation for the dissemination areas we are interested in.
In this function we make sure to input our geographic sf object first to indicate that is the dataset for which we wish to keep all the records (e.g. circle A in the above figure), and the a-spatial tabular data second to indicate that is the data we wish to join matching records. We can indicate that our unique identifier we use to match the records have different names in each dataset. The other alternative is to rename one of the variables to match the other, but this works also.
vancouver_da_dep <- left_join(census_intersects_cntrd,deprivation_da, by = c("GeoUID"="PRCDDA"))Let’s take a look at the data and then map it out to ensure it worked.
vancouver_da_dep## Simple feature collection with 993 features and 23 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 483690.9 ymin: 5449535 xmax: 498329.2 ymax: 5462381
## Projected CRS: NAD83 / UTM zone 10N
## First 10 features:
## Population Households GeoUID Type CD_UID Shape Area Dwellings CSD_UID CT_UID CMA_UID Region Name Area (sq km) v_CA16_2397 Province Dissemination area (DA) Population
## 1 632 254 59150307 DA 5915 0.29926 273 5915022 9330053.02 59933 Vancouver 0.29926 90965 British Columbia 669.4691
## 2 501 203 59150308 DA 5915 0.10956 223 5915022 9330053.02 59933 Vancouver 0.10956 93611 British Columbia 470.3942
## 3 745 269 59150309 DA 5915 0.11190 299 5915022 9330053.02 59933 Vancouver 0.11190 84736 British Columbia 728.7991
## 4 536 283 59150310 DA 5915 0.10940 310 5915022 9330053.02 59933 Vancouver 0.10940 45184 British Columbia 396.3275
## 5 532 181 59150311 DA 5915 0.08094 199 5915022 9330053.02 59933 Vancouver 0.08094 82517 British Columbia 623.7097
## 6 562 203 59150312 DA 5915 0.08706 211 5915022 9330053.02 59933 Vancouver 0.08706 74368 British Columbia 530.6063
## 7 1088 402 59150313 DA 5915 0.17284 418 5915022 9330052.02 59933 Vancouver 0.17284 61184 British Columbia 954.5501
## 8 556 170 59150328 DA 5915 0.23224 181 5915022 9330052.02 59933 Vancouver 0.23224 80896 British Columbia 558.0631
## 9 959 333 59150329 DA 5915 0.30854 356 5915022 9330052.02 59933 Vancouver 0.30854 66688 British Columbia 891.8479
## 10 506 174 59150330 DA 5915 0.07820 189 5915022 9330052.02 59933 Vancouver 0.07820 82176 British Columbia 460.2466
## Ethno-cultural composition Quintiles Ethno-cultural composition Scores Situational vulnerability Quintiles Situational vulnerability Scores Economic dependency Quintiles
## 1 3 -0.2744621 2 -0.4505756 1
## 2 3 -0.2350062 1 -0.7415757 2
## 3 4 0.5496419 1 -0.6846734 2
## 4 4 0.8662144 5 0.3722661 2
## 5 4 0.6849209 2 -0.4435429 2
## 6 5 1.0570569 3 -0.2977577 1
## 7 4 0.5790914 5 0.6069958 3
## 8 5 1.4650933 3 -0.2751484 1
## 9 5 2.4115451 4 0.0428343 4
## 10 5 0.9880264 3 -0.2410908 1
## Economic dependency Scores Residential instability Quintiles Residential instability Scores geometry
## 1 -0.9301352 4 0.1246316 MULTIPOLYGON (((498308 5459...
## 2 -0.6086799 2 -0.5748885 MULTIPOLYGON (((498301.6 54...
## 3 -0.6227231 3 -0.2030643 MULTIPOLYGON (((497939.5 54...
## 4 -0.4802208 5 1.9569124 MULTIPOLYGON (((498298.6 54...
## 5 -0.5483774 3 -0.4832458 MULTIPOLYGON (((498128.4 54...
## 6 -0.7775863 3 -0.1416693 MULTIPOLYGON (((498298.6 54...
## 7 -0.0186630 4 0.2009225 MULTIPOLYGON (((498259 5458...
## 8 -1.4517890 3 -0.2034135 MULTIPOLYGON (((498102 5458...
## 9 0.1864334 3 -0.3812451 MULTIPOLYGON (((497801.1 54...
## 10 -0.7402528 2 -0.6010522 MULTIPOLYGON (((497312 5458...mapview(vancouver_da_dep, zcol = "Residential instability Scores")4.4 Buffers
Buffers are a geometric operation which creates a polygon based on a specified distance around individual spatial features. Buffers can be constructed around points, line and polygons.
](images/GIS_Buffer.png)
Figure 4.4: Diagram showing a buffer (green) of a point, line and polygon feature (red).Source
Let’s illustrate with a few examples.
First let’s load some data. Recall the code from last chapter to convert the .csv data on crash events to a spatial object in R. Below we repeat that code.
icbc_crash <- read_csv("D:/Github/crash-course-gis-r/data/ICBC_cyclist_crashes.csv") %>%# use the file location of where you have stored your data
filter(!is.na(Latitude)) %>% # remove the recrods which have missing values
st_as_sf(.,coords = c("Longitude","Latitude"),crs = 4326) %>%
mutate(uniqid = row_number()) #assign unique identifier to each crash record##
## -- Column specification ------------------------------------------------------------------------------------------------------------------------------------------------------------------------
## cols(
## .default = col_character(),
## `Date Of Loss Year` = col_double(),
## `Crash Count` = col_double(),
## Latitude = col_double(),
## Longitude = col_double(),
## `Number of Records` = col_double()
## )
## i Use `spec()` for the full column specifications.Now on your own select the crashes that occur within the City of Vancouver and assign it to a new spatial object called icbc_crash_van.
HINT you can do either an attribute query or a spatial query and you will need to convert the layer to a projected coordinate system.
The output should look like this:
mapview(icbc_crash_van)Buffers can be created using the function st_buffer() from the sf package and inputting the sf object that you wish to create buffers around, and the size of the buffer in a unit of length (e.g. metres). Check the units of your map projection to ensure you know what numbers to input for distance. Here let’s create a 200m buffer around each crash event in our dataset.
st_crs(icbc_crash_van)$units # my map projection uses metres as a unit - ## [1] "m"crash_buffer_icbc <- st_buffer(icbc_crash_van,dist = 150)
crash_buffer_icbc## Simple feature collection with 3373 features and 26 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 484177.4 ymin: 5449829 xmax: 498479.2 ymax: 5462176
## Projected CRS: NAD83 / UTM zone 10N
## # A tibble: 3,373 x 27
## `Crash Breakdown~ `Date Of Loss Ye~ `Crash Severity` `Cyclist Flag` `Day Of Week` `Derived Crash Con~ `Intersection Cr~ `Month Of Year` `Motorcycle Fla~ `Municipality Nam~ `Parking Lot Fl~
## * <chr> <dbl> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr>
## 1 LOWER MAINLAND 2018 CASUALTY Yes SATURDAY SINGLE VEHICLE Yes MAY No VANCOUVER No
## 2 LOWER MAINLAND 2017 PROPERTY DAMAGE ~ Yes THURSDAY UNDETERMINED Yes MARCH No VANCOUVER No
## 3 LOWER MAINLAND 2015 CASUALTY Yes FRIDAY SINGLE VEHICLE Yes JUNE No VANCOUVER No
## 4 LOWER MAINLAND 2015 CASUALTY Yes WEDNESDAY SINGLE VEHICLE Yes JUNE No VANCOUVER No
## 5 LOWER MAINLAND 2016 CASUALTY Yes THURSDAY SINGLE VEHICLE Yes NOVEMBER No VANCOUVER No
## 6 LOWER MAINLAND 2019 PROPERTY DAMAGE ~ Yes WEDNESDAY SINGLE VEHICLE Yes JANUARY No VANCOUVER No
## 7 LOWER MAINLAND 2016 CASUALTY Yes FRIDAY SINGLE VEHICLE Yes AUGUST No VANCOUVER No
## 8 LOWER MAINLAND 2018 CASUALTY Yes FRIDAY SINGLE VEHICLE Yes MAY No VANCOUVER No
## 9 LOWER MAINLAND 2016 CASUALTY Yes WEDNESDAY SINGLE VEHICLE Yes AUGUST No VANCOUVER No
## 10 LOWER MAINLAND 2017 CASUALTY Yes TUESDAY SINGLE VEHICLE Yes NOVEMBER No VANCOUVER No
## # ... with 3,363 more rows, and 16 more variables: Pedestrian Flag <chr>, Region <chr>, Street Full Name (ifnull) <chr>, Time Category <chr>, Crash Breakdown <chr>,
## # Cross Street Full Name <chr>, Heavy Veh Flag <chr>, Mid Block Crash <chr>, Municipality Name <chr>, Municipality With Boundary <chr>, Road Location Description <chr>,
## # Street Full Name <chr>, Crash Count <dbl>, Number of Records <dbl>, geometry <POLYGON [m]>, uniqid <int>mapview(crash_buffer_icbc)Note how the buffer retains the attributes of each original crash event!
We can also create buffers for lines, and polygons. For example, lets create 20m buffers around our separated infrastructure data we created earlier in the chapter.
bicycle_buffer <- st_buffer(bicycle,dist=20)## Warning in st_buffer.sfc(st_geometry(x), dist, nQuadSegs, endCapStyle = endCapStyle, : st_buffer does not correctly buffer longitude/latitude data## dist is assumed to be in decimal degrees (arc_degrees).Shoot! Coordinate systems problems at it again! Remember that our data have to be projected in order to spatial operations on them. We know what to do from here though. Let’s convert it to be the same projection as the rest of our data NAD83/UTM Zone 10N.
bicycle_proj <- st_transform(bicycle,crs=st_crs(icbc_crash_van)$epsg)
bicycle_buffer <- st_buffer(bicycle_proj,dist=20)
mapview(bicycle_buffer)4.5 Union
A union of a single layer can be used to combine the geometry into a single unit using st_union(). This removes all of the attributes however.
](images/Union.jpg)
Figure 4.5: Diagram showing a union of two polygon layers.Source
Let’s union the buffer data from our crashes.
crash_buffer_union <- st_union(crash_buffer_icbc)
mapview(crash_buffer_union)4.6 Intersection
](images/intersection.png)
Figure 4.6: Diagram showing intersection of two polygon layers.Source
intersect_van_census_crash_buffer <- st_intersection(census_intersects_cntrd,crash_buffer_icbc)
intersect_van_census_crash_buffer## Simple feature collection with 13303 features and 39 fields
## Geometry type: GEOMETRY
## Dimension: XY
## Bounding box: xmin: 484205.8 ymin: 5449833 xmax: 498329.2 ymax: 5462176
## Projected CRS: NAD83 / UTM zone 10N
## First 10 features:
## Population Households GeoUID Type CD_UID Shape.Area Dwellings CSD_UID CT_UID CMA_UID Region.Name Area..sq.km. v_CA16_2397 Crash.Breakdown.2 Date.Of.Loss.Year Crash.Severity
## 912 1093 529 59153583 DA 5915 0.03968 607 5915022 9330059.06 59933 Vancouver 0.03968 33344 LOWER MAINLAND 2018 CASUALTY
## 990 4685 2596 59154001 DA 5915 0.28495 2960 5915022 9330059.13 59933 Vancouver 0.28495 72082 LOWER MAINLAND 2018 CASUALTY
## 989 1668 930 59154000 DA 5915 0.37695 1081 5915022 9330059.14 59933 Vancouver 0.37695 80811 LOWER MAINLAND 2017 PROPERTY DAMAGE ONLY
## 990.1 4685 2596 59154001 DA 5915 0.28495 2960 5915022 9330059.13 59933 Vancouver 0.28495 72082 LOWER MAINLAND 2017 PROPERTY DAMAGE ONLY
## 423 779 413 59150755 DA 5915 0.05163 454 5915022 9330057.01 59933 Vancouver 0.05163 17740 LOWER MAINLAND 2015 CASUALTY
## 424 203 76 59150756 DA 5915 0.02422 109 5915022 9330057.01 59933 Vancouver 0.02422 56960 LOWER MAINLAND 2015 CASUALTY
## 425 1123 631 59150757 DA 5915 0.07878 726 5915022 9330057.01 59933 Vancouver 0.07878 28832 LOWER MAINLAND 2015 CASUALTY
## 426 499 254 59150758 DA 5915 0.02731 261 5915022 9330057.01 59933 Vancouver 0.02731 28032 LOWER MAINLAND 2015 CASUALTY
## 772 558 262 59153191 DA 5915 0.05637 262 5915022 9330056.01 59933 Vancouver 0.05637 66389 LOWER MAINLAND 2015 CASUALTY
## 773 552 366 59153192 DA 5915 0.03022 375 5915022 9330056.01 59933 Vancouver 0.03022 25280 LOWER MAINLAND 2015 CASUALTY
## Cyclist.Flag Day.Of.Week Derived.Crash.Congifuration Intersection.Crash Month.Of.Year Motorcycle.Flag Municipality.Name..ifnull. Parking.Lot.Flag Pedestrian.Flag Region
## 912 Yes SATURDAY SINGLE VEHICLE Yes MAY No VANCOUVER No No LOWER MAINLAND
## 990 Yes SATURDAY SINGLE VEHICLE Yes MAY No VANCOUVER No No LOWER MAINLAND
## 989 Yes THURSDAY UNDETERMINED Yes MARCH No VANCOUVER No No LOWER MAINLAND
## 990.1 Yes THURSDAY UNDETERMINED Yes MARCH No VANCOUVER No No LOWER MAINLAND
## 423 Yes FRIDAY SINGLE VEHICLE Yes JUNE No VANCOUVER No No LOWER MAINLAND
## 424 Yes FRIDAY SINGLE VEHICLE Yes JUNE No VANCOUVER No No LOWER MAINLAND
## 425 Yes FRIDAY SINGLE VEHICLE Yes JUNE No VANCOUVER No No LOWER MAINLAND
## 426 Yes FRIDAY SINGLE VEHICLE Yes JUNE No VANCOUVER No No LOWER MAINLAND
## 772 Yes WEDNESDAY SINGLE VEHICLE Yes JUNE No VANCOUVER No No LOWER MAINLAND
## 773 Yes WEDNESDAY SINGLE VEHICLE Yes JUNE No VANCOUVER No No LOWER MAINLAND
## Street.Full.Name..ifnull. Time.Category Crash.Breakdown Cross.Street.Full.Name Heavy.Veh.Flag Mid.Block.Crash Municipality.Name Municipality.With.Boundary
## 912 EXPO BLVD 18:00-20:59 SATURDAY ABBOTT ST N N VANCOUVER VANCOUVER
## 990 EXPO BLVD 18:00-20:59 SATURDAY ABBOTT ST N N VANCOUVER VANCOUVER
## 989 EXPO BLVD 15:00-17:59 THURSDAY CARRALL ST N N VANCOUVER VANCOUVER
## 990.1 EXPO BLVD 15:00-17:59 THURSDAY CARRALL ST N N VANCOUVER VANCOUVER
## 423 E GEORGIA ST 09:00-11:59 FRIDAY GORE AVE N N VANCOUVER VANCOUVER
## 424 E GEORGIA ST 09:00-11:59 FRIDAY GORE AVE N N VANCOUVER VANCOUVER
## 425 E GEORGIA ST 09:00-11:59 FRIDAY GORE AVE N N VANCOUVER VANCOUVER
## 426 E GEORGIA ST 09:00-11:59 FRIDAY GORE AVE N N VANCOUVER VANCOUVER
## 772 ADANAC ST 15:00-17:59 WEDNESDAY COMMERCIAL DR N N VANCOUVER VANCOUVER
## 773 ADANAC ST 15:00-17:59 WEDNESDAY COMMERCIAL DR N N VANCOUVER VANCOUVER
## Road.Location.Description Street.Full.Name Crash.Count Number.of.Records uniqid geometry
## 912 ABBOTT ST & EXPO BLVD & PAT QUINN WAY EXPO BLVD 1 1 42 POLYGON ((492282.1 5458528,...
## 990 ABBOTT ST & EXPO BLVD & PAT QUINN WAY EXPO BLVD 1 1 42 POLYGON ((492282.1 5458528,...
## 989 CARRALL ST & EXPO BLVD EXPO BLVD 1 1 43 POLYGON ((492168.2 5458279,...
## 990.1 CARRALL ST & EXPO BLVD EXPO BLVD 1 1 43 POLYGON ((492167.4 5458271,...
## 423 E GEORGIA ST & GORE AVE E GEORGIA ST 1 1 45 POLYGON ((492947.4 5458322,...
## 424 E GEORGIA ST & GORE AVE E GEORGIA ST 1 1 45 POLYGON ((492947.4 5458322,...
## 425 E GEORGIA ST & GORE AVE E GEORGIA ST 1 1 45 POLYGON ((492947.9 5458274,...
## 426 E GEORGIA ST & GORE AVE E GEORGIA ST 1 1 45 POLYGON ((493012.9 5458522,...
## 772 ADANAC ST & COMMERCIAL DR ADANAC ST 1 1 46 POLYGON ((494757.3 5458221,...
## 773 ADANAC ST & COMMERCIAL DR ADANAC ST 1 1 46 POLYGON ((494903.4 5458315,...Now we can see that each polygon has information from both the census data and the crash data. We could now map the information from overlapping census data based on the location of each buffer. It takes the overlapping polygons and combines the attribute data, removing areas without overlapping data. Let’s plot income for illustration
plot(intersect_van_census_crash_buffer["v_CA16_2397"])
There are a lot of overlapping polygons in this dataset. Let’s simplify these data by averaging the census income data for each unique crash buffer. Recall that each icbc crash had a unique identifer we assigned it when we loaded it. The code below essentially is taking all of the different census data that overlapped with a crash buffer and averaging the income from that dissemination area. More formally, we can think about it as for each polygon that has the same uniqid, dissolve the geometry so that they all share the same extent and take the average value of v_CA16_2397 and assign it to the new geography. This will necessarily only leave 1 income value for each crash buffer, and “dissolve” the boundaries of buffers that were on the border of multiple census areas. This might make more sense by comparing the intersection plot with the dissolved/aggregated plot below.
average_income_by_crash_buffer <- intersect_van_census_crash_buffer %>%
group_by(uniqid) %>%
summarise(average_income = mean(v_CA16_2397,na.rm=TRUE))
plot(average_income_by_crash_buffer["average_income"])
4.7 Spatial Joins
A spatial join will likely be one of the most common types of overlay analyses you will employ in your analyses. In a spatial join the goal is to combine the attribute data from one geographic dataset (the join dataset) and append it to the geometry of another geographic dataset (the target dataset), based on the spatial relationship between the two datasets (figure 4.7. In a spatial join, the join dataset typically has 1 to 2 more dimensions that of the target features (recall that polygons are 2 dimensional, lines are 1 dimensional and points are 0 dimensional).
Usually the goal in a spatial join is to summarize or aggregate attribute information of either points that fall on/within the geometry of polygon or line features, or, to summarize the attribute information of lines that fall within polygons. In a transportation safety context, a good example would be to aggregate crash point data to the neighbourhood it occured within (point-in-polygon), or the number of crashes that occur on a road segment (points-on-lines). Another transportation example would be to aggregate the road network features by neighbourhoods (lines-in polygons).
For example in Figure 4.7 we have the crash point map where each crash has an ID. We also have a neighbourhood map with the names of neighbourhoods. We want to join the neighbourhood names to the point data, based on the co-location of the points and the neighbourhoods. We define the crash points dataset as the “target dataset.” The target dataset is the one where we want to preserve the existing geometry and attributes. We define the neighbourhood polygon dataset as the “join dataset.” THe join dataset contains the attribute information we want to join to the attribute table of the target dataset. In short, when we run a spatial join we preserve the geometry of the target dataset and add the attribute information of the join dataset.
Figure 4.7: Diagram showing spatial join of polygons to points
ta.
In R we can use the st_join() function from the sf package to accomplish a spatial join. In this function the arguments that are input are st_join(x,y,join) where x refers to a target dataset (which needs to be an sf object), y refers to the join dataset. The default spatial relationship to use to define which feature attributes from the join dataset will be added to to which features from the target dataset is based on an intersection (i.e. st_intersects()). This can be modified to join attributes based on a different spatial relationship such as st_contains() or st_is_within_distance() among others. See the reference page for the full list.
Let’s go a through a few examples based on different purposes.
4.7.0.1 Points-in-polygon
In this example our end goal is to map the number of ICBC recorded bicycle crashes by dissemination area boundaries for the city of Vancouver. We will use the datasets we loaded earlier in the chapter:
vancouver_da_dep: Dissemination Area census geographies that include information on deprivation levels (see Attribute Join section above)icbc_crash_van: ICBC bicycle crash data for the City of Vancouver (see Buffer section above)
To enable counting of crashes within each dissemination area we are going to use a spatial join within a series of steps that will also draw upon what we learned about attribute joins and some data summary functionality from dplyr.
The first step is to join the DA level information to each point based on whether the point intersects with that polygon. Essentially, we want to take all of the DA level information and join it to each crash event based on the colocation of the crash and census data.
icbc_crash_van_join_da <- st_join(icbc_crash_van,vancouver_da_dep)
head(icbc_crash_van_join_da)## Simple feature collection with 6 features and 49 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 492170.9 ymin: 5458315 xmax: 494873.9 ymax: 5458438
## Projected CRS: NAD83 / UTM zone 10N
## # A tibble: 6 x 50
## `Crash Breakdown~ `Date Of Loss Ye~ `Crash Severity` `Cyclist Flag` `Day Of Week` `Derived Crash Cong~ `Intersection Cr~ `Month Of Year` `Motorcycle Fla~ `Municipality Nam~ `Parking Lot Fl~
## <chr> <dbl> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr>
## 1 LOWER MAINLAND 2018 CASUALTY Yes SATURDAY SINGLE VEHICLE Yes MAY No VANCOUVER No
## 2 LOWER MAINLAND 2017 PROPERTY DAMAGE ~ Yes THURSDAY UNDETERMINED Yes MARCH No VANCOUVER No
## 3 LOWER MAINLAND 2015 CASUALTY Yes FRIDAY SINGLE VEHICLE Yes JUNE No VANCOUVER No
## 4 LOWER MAINLAND 2015 CASUALTY Yes WEDNESDAY SINGLE VEHICLE Yes JUNE No VANCOUVER No
## 5 LOWER MAINLAND 2016 CASUALTY Yes THURSDAY SINGLE VEHICLE Yes NOVEMBER No VANCOUVER No
## 6 LOWER MAINLAND 2019 PROPERTY DAMAGE ~ Yes WEDNESDAY SINGLE VEHICLE Yes JANUARY No VANCOUVER No
## # ... with 39 more variables: Pedestrian Flag <chr>, Region <chr>, Street Full Name (ifnull) <chr>, Time Category <chr>, Crash Breakdown <chr>, Cross Street Full Name <chr>,
## # Heavy Veh Flag <chr>, Mid Block Crash <chr>, Municipality Name <chr>, Municipality With Boundary <chr>, Road Location Description <chr>, Street Full Name <chr>, Crash Count <dbl>,
## # Number of Records <dbl>, geometry <POINT [m]>, uniqid <int>, Population <int>, Households <int>, GeoUID <chr>, Type <fct>, CD_UID <chr>, Shape Area <dbl>, Dwellings <int>, CSD_UID <chr>,
## # CT_UID <chr>, CMA_UID <chr>, Region Name <fct>, Area (sq km) <dbl>, v_CA16_2397 <dbl>, Province <chr>, Dissemination area (DA) Population <dbl>,
## # Ethno-cultural composition Quintiles <dbl>, Ethno-cultural composition Scores <dbl>, Situational vulnerability Quintiles <dbl>, Situational vulnerability Scores <dbl>,
## # Economic dependency Quintiles <dbl>, Economic dependency Scores <dbl>, Residential instability Quintiles <dbl>, Residential instability Scores <dbl>In the attribute table of our new object icbc_crash_van_join_da we can see that each point now has information for the dissemination area with which it is located (e.g. intersects with). A successful spatial join! But how do we go from here to map the number of crashes by dissemination area? The next step is to take our spatially joined data and use dplyr to group the data by the dissemination area unique identifier, add up the Crash Count variable for each row with the same unique id. Right now our data is still an sf point object. Now that the spatial join is done, we don’t need to keep the geometry associated with this feature so we can remove it using the st_drop_geometry() function. This will speed up the computation time for the functions performed on grouped data.
In the code chunk below we are grouping our data by the dissemination area unique id we just joined to our point data using the group_by() function and summing the Crash Count variable using the summarise() function.
no_icbc_crashes_by_da <- icbc_crash_van_join_da %>%
st_drop_geometry() %>%
group_by(GeoUID) %>%
summarise(total_crashes = sum(`Crash Count`))The output is a table with one row for each dissemination area (unique GeoUID)
head(no_icbc_crashes_by_da)## # A tibble: 6 x 2
## GeoUID total_crashes
## <chr> <dbl>
## 1 59150310 1
## 2 59150311 1
## 3 59150313 4
## 4 59150329 3
## 5 59150330 1
## 6 59150332 1Now we have counts of crashes by the unique identifier of the dissemination area it occured in. Now we can simply join this table back to the polygon data using a left_join().
vancouver_da_dep_crashes <- left_join(vancouver_da_dep,no_icbc_crashes_by_da,by="GeoUID")
head(vancouver_da_dep_crashes)## Simple feature collection with 6 features and 24 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 497754.1 ymin: 5458297 xmax: 498323.1 ymax: 5460083
## Projected CRS: NAD83 / UTM zone 10N
## Population Households GeoUID Type CD_UID Shape Area Dwellings CSD_UID CT_UID CMA_UID Region Name Area (sq km) v_CA16_2397 Province Dissemination area (DA) Population
## 1 632 254 59150307 DA 5915 0.29926 273 5915022 9330053.02 59933 Vancouver 0.29926 90965 British Columbia 669.4691
## 2 501 203 59150308 DA 5915 0.10956 223 5915022 9330053.02 59933 Vancouver 0.10956 93611 British Columbia 470.3942
## 3 745 269 59150309 DA 5915 0.11190 299 5915022 9330053.02 59933 Vancouver 0.11190 84736 British Columbia 728.7991
## 4 536 283 59150310 DA 5915 0.10940 310 5915022 9330053.02 59933 Vancouver 0.10940 45184 British Columbia 396.3275
## 5 532 181 59150311 DA 5915 0.08094 199 5915022 9330053.02 59933 Vancouver 0.08094 82517 British Columbia 623.7097
## 6 562 203 59150312 DA 5915 0.08706 211 5915022 9330053.02 59933 Vancouver 0.08706 74368 British Columbia 530.6063
## Ethno-cultural composition Quintiles Ethno-cultural composition Scores Situational vulnerability Quintiles Situational vulnerability Scores Economic dependency Quintiles
## 1 3 -0.2744621 2 -0.4505756 1
## 2 3 -0.2350062 1 -0.7415757 2
## 3 4 0.5496419 1 -0.6846734 2
## 4 4 0.8662144 5 0.3722661 2
## 5 4 0.6849209 2 -0.4435429 2
## 6 5 1.0570569 3 -0.2977577 1
## Economic dependency Scores Residential instability Quintiles Residential instability Scores total_crashes geometry
## 1 -0.9301352 4 0.1246316 NA MULTIPOLYGON (((498308 5459...
## 2 -0.6086799 2 -0.5748885 NA MULTIPOLYGON (((498301.6 54...
## 3 -0.6227231 3 -0.2030643 NA MULTIPOLYGON (((497939.5 54...
## 4 -0.4802208 5 1.9569124 1 MULTIPOLYGON (((498298.6 54...
## 5 -0.5483774 3 -0.4832458 1 MULTIPOLYGON (((498128.4 54...
## 6 -0.7775863 3 -0.1416693 NA MULTIPOLYGON (((498298.6 54...Our new column total_crashes has NA values for some polygons. This is because if a dissemination area did not have crashes occur within their boundaries they would not have been attribtued to a crash in our point dataset during the spatial join. The NA value actually then represents a 0. So let’s change the NA values to 0 using the replace_na() function from the tidyr package.
vancouver_da_dep_crashes <- vancouver_da_dep_crashes %>%
mutate(total_crashes = replace_na(total_crashes,0))Now let’s map the number of crashes from the new dataset vancouver_da_dep_crashes and overlay the original crash data icbc_crash_van
mapview(vancouver_da_dep_crashes,zcol = "total_crashes") + mapview(icbc_crash_van)4.7.0.2 Points-on-line
In a points-on-line spatial join we are applying the same concept as the previous example, but we instead want to join attribute data from line data to our point data. This type of operation is a bit trickier to implement because lines only 1-dimensions. So if a point does not fall exactly onto it’s geometry it is harder to join based on an intersection. There are a few ways around this as is the case with almost any kind of geoprocessing task. The simplest way would be to change spatial relationship from the default st_intersects to something like st_is_within_distance and specify a distance tolerance for considering a line to be attributed to a point. This is great if you are okay with a one-to-many spatial join where you create a new point object for every line segment within this distance. However, if we wanted to ensure a one-to-one join where you only join the attribute data from one line segment to one point - you will need to find a work around. One way to do this is to “snap” each point data to the nearest road segment, and then use the within st_join(x,y,join = st_is_within_distance, dist = tiny_value_like_1cm) to conduct the spatial join.
sf does not have a built in function for snapping points to nearest line segment but we can create one ourselves. Below is a function called st_snap_points() we define manually. It has 4 arguments, x for the point data, y for the line data, max_dist is the distance tolerance value, and namevar for specifying any uniqid from th original dataset. Feel free to ignore the details of the function below but make sure to run it before going through the rest of the section.
st_snap_points <- function(x, y, namevar, max_dist = 1000) {
# this evaluates the length of the data
if (inherits(x, "sf")) n = nrow(x)
if (inherits(x, "sfc")) n = length(x)
# this part:
# 1. loops through every piece of data (every point)
# 2. snaps a point to the nearest line geometries
# 3. calculates the distance from point to line geometries
# 4. retains only the shortest distances and generates a point at that intersection
out = do.call("c",
lapply(seq(n), function(i) {
nrst = st_nearest_points(st_geometry(x)[i], y)
nrst_len = st_length(nrst)
nrst_mn = which.min(nrst_len)
if (as.vector(nrst_len[nrst_mn]) > max_dist) return(st_geometry(x)[i])
return(st_cast(nrst[nrst_mn], "POINT")[2])
})
)
# this part converts the data to a dataframe and adds a named column of your choice
out_xy <- st_coordinates(out) %>% as.data.frame()
out_xy <- out_xy %>%
mutate({{namevar}} := x[[namevar]]) %>%
st_as_sf(coords=c("X","Y"), crs=st_crs(x), remove=FALSE)
return(out_xy)
}Now that we have defined this function let’s try out another scenario. Similar to the section above, our end goal is to map the number of bicycle crashes (icbc_van_crash) that occur on the separated bicycle infrastructure network in Vancouver (bicycle_proj). The steps do so are as follows:
- Snap the crash locations to the nearest line segment based on a reasonable distance threshold using the
st_snap_pointsfunction - Perform a spatial join where we join the attributes of the nearest line segment to each point
- Aggregate the attribute data of the new spatially joined data by a unique identifier of the line segment
- Join the aggregated data back to the line segment data based on the unique identifier
Snap Points
Set a max distance of 20m to snap points. Outside of 20m then we don’t consider them occurring on an infrastructure segment and the position of the event will stay the same. This function may take a bit of time to run.
icbc_snp <- st_snap_points(x = icbc_crash_van,#crash van
y = bicycle_proj, #network infrasructure data
max_dist = 20,
namevar = "uniqid") #this only returns the geometry - doesn't preserve attributes
icbc_snp <- left_join(icbc_snp,
icbc_crash_van %>%
st_drop_geometry(),
by = "uniqid")#return the attributes to the snapped locations
mapview(icbc_snp) + mapview(bicycle_proj)Spatial Join
Now for each point let’s count the number of crashes that occur on that segment based on our 20m threshold. To do so we will join the unique id of the nearest road segment to the point data using st_join() and the st_within_distance() function. Since we already have snapped the points to the nearest road segment there is only one road segment assigned to each point and we specify a tiny distance to consider a spatial join.
icbc_snp_bicycle_inf <- st_join(icbc_snp,
bicycle_proj["osm_id"], #select only the column "osm_id"
join = st_is_within_distance, dist = 0.001) #for each point assign the LIXID that it falls on by indicating a very small distance tolerance (here it is 1mm)
icbc_snp_bicycle_inf %>%
select(uniqid,`Crash Count`,osm_id) # look at the results## Simple feature collection with 3405 features and 3 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 484327.3 ymin: 5449979 xmax: 498329.2 ymax: 5462021
## Projected CRS: NAD83 / UTM zone 10N
## First 10 features:
## uniqid Crash Count osm_id geometry
## 1 42 1 <NA> POINT (492170.9 5458438)
## 2 43 1 368058222 POINT (492279.2 5458372)
## 3 45 1 <NA> POINT (492940 5458424)
## 4 46 1 368210443 POINT (494873.4 5458316)
## 5 47 1 368210443 POINT (494873.4 5458316)
## 6 48 1 368058222 POINT (492279.2 5458372)
## 7 51 1 748846437 POINT (490680.9 5458413)
## 8 52 1 <NA> POINT (492170.9 5458438)
## 9 56 1 827958314 POINT (491713.7 5458318)
## 10 57 1 540329767 POINT (490452.9 5458326)Note how only the points that were snapped to bicycle infrastructure locations have a uniqid assigned from the bicycle infrastructure data. The NA values indicate no lines were within a distance of 0.001 meters (i.e. equivalent to 20m since we snapped the points to the nearest segment based on a 20m threshold)
Spatial Aggregation
Now let’s count the number of points with same osm_id.
no_crashes_by_osm_id <- icbc_snp_bicycle_inf %>%
group_by(osm_id) %>%
summarise(n_crashes = sum(`Crash Count`)) %>%
st_drop_geometry() # drop the geometry from this - we only need the uniqid for the bicycle infrastrucutre data and the countsAttribute Join
Now we can join this data back to the bicycle infrastructure polyline data. Polylines without any crashes will be assigned an NA value because they don’t have a matching record in the join. We then simply replace the NA values with a 0 after the join using the replace_na() from the tidyr package.
bicycle_infra_crash_count <- left_join(bicycle_proj,no_crashes_by_osm_id , by = "osm_id") %>%
mutate(n_crashes = replace_na(n_crashes,0))Let’s map out the number of crashes by bicycle infrastructure and layer on the original point data and zoom into areas with collision events and see if the results make sense.
mapview(bicycle_infra_crash_count,zcol="n_crashes") + mapview(icbc_snp)4.7.0.3 Lines-in-polygon
Here let’s assume a scenario where the end goal is to calculate the length of separated bicycle infrastructure within each dissemination area in the City of Vancouver. Recall we created the following datasets from earlier in the chapter:
bicycle_proj: separated bicycle infrastructure polylines for the City of Vancouver in UTM Zone 10Nvancouver_da_dep: dissemination area boundaries for the City of Vancouver in UTM Zone 10N
The main problem is that as of now, the line segments for bicycle infrastructure do not line up with the boundary of disseminatin areas. This concept is illustrated in the left panel of Figure 4.7 where we have one polyline that crosses multiple boundaries. We need to cut the line segments where they intersect with the boundaries of dissemination areas such as in Figure 4.7 using st_intersection() so that we can then calculate the length of each individual segment using st_length().
Figure 4.8: Intersection of a single polyline (blue) by polygons (grey)
From there we can sum the length of each line segment by the dissemination area’s unique identfier (GeoUID) in a few different ways, including a spatial join. Let’s try this out so we can see some of the problems that can arise with this type of data and operation.
First, let’s create the “cut” version of our bicycle infrastructure layer using st_intersection(). I am going to do this with subsets of original datasets that only contain the unique identifier as attribute data to make our data exploration a little cleaner.
## Warning: attribute variables are assumed to be spatially constant throughout all geometriesWe can then calculate the length of each line segment using a function called st_length().
bicycle_cut <- bicycle_cut %>%
mutate(length_m = st_length(.))
bicycle_cut## Simple feature collection with 1473 features and 3 fields
## Geometry type: GEOMETRY
## Dimension: XY
## Bounding box: xmin: 483720 ymin: 5449700 xmax: 498317.2 ymax: 5462376
## Projected CRS: NAD83 / UTM zone 10N
## First 10 features:
## osm_id GeoUID geometry length_m
## 35260598 35260598 59150307 LINESTRING (498105.5 545987... 157.80005 [m]
## 35260599 35260599 59150307 LINESTRING (497911 5459475,... 32.86117 [m]
## 36691655 36691655 59150307 LINESTRING (497911 5459475,... 163.86383 [m]
## 41840546 41840546 59150307 LINESTRING (498180.4 545988... 76.01887 [m]
## 75573893 75573893 59150307 LINESTRING (498226.7 545987... 44.49847 [m]
## 210638528 210638528 59150307 LINESTRING (498190.6 545997... 84.65905 [m]
## 210638530 210638530 59150307 LINESTRING (498235.3 545998... 46.27267 [m]
## 336365962 336365962 59150307 LINESTRING (497942.1 545946... 56.43094 [m]
## 340103260 340103260 59150307 LINESTRING (497944.5 545952... 86.09529 [m]
## 380516622 380516622 59150307 LINESTRING (498317.2 545999... 84.11612 [m]Notice how the intersection also attributed the GeoUID from the Dissemination Area boundary. We could then use that information to summarise the length of polylines within each GeoUID.
But let’s take a look at our data.
mapview(bicycle_cut) + mapview(vancouver_da_sub)59150308
59150309
59150310
59150311
59150312
59150313
59150328
59150329
59150330
59150331
59150332
59150333
59150334
59150335
59150336
59150337
59150338
59150339
59150340
59150341
59150342
59150343
59150344
59150345
59150346
59150347
59150348
59150349
59150350
59150351
59150352
59150353
59150354
59150355
59150356
59150357
59150358
59150359
59150360
59150361
59150362
59150363
59150364
59150365
59150366
59150367
59150368
59150369
59150370
59150371
59150372
59150373
59150375
59150376
59150377
59150378
59150379
59150380
59150381
59150382
59150383
59150384
59150385
59150386
59150387
59150388
59150389
59150390
59150391
59150392
59150393
59150394
59150395
59150396
59150397
59150398
59150399
59150400
59150401
59150402
59150403
59150404
59150405
59150406
59150407
59150408
59150409
59150410
59150412
59150414
59150415
59150416
59150417
59150419
59150420
59150421
59150422
59150423
59150424
59150425
59150426
59150427
59150428
59150429
59150430
59150431
59150432
59150433
59150434
59150435
59150436
59150437
59150438
59150439
59150440
59150441
59150442
59150443
59150444
59150445
59150446
59150447
59150448
59150449
59150450
59150451
59150452
59150453
59150454
59150455
59150457
59150459
59150460
59150461
59150462
59150463
59150465
59150466
59150468
59150469
59150470
59150471
59150472
59150473
59150474
59150475
59150476
59150477
59150478
59150479
59150480
59150481
59150482
59150483
59150484
59150485
59150486
59150487
59150488
59150489
59150490
59150491
59150492
59150493
59150494
59150495
59150496
59150497
59150498
59150499
59150500
59150503
59150504
59150505
59150506
59150507
59150508
59150509
59150510
59150511
59150512
59150513
59150514
59150515
59150516
59150517
59150518
59150519
59150520
59150521
59150522
59150523
59150524
59150525
59150526
59150527
59150528
59150529
59150530
59150531
59150532
59150533
59150534
59150535
59150536
59150537
59150538
59150539
59150540
59150541
59150542
59150543
59150544
59150545
59150546
59150547
59150548
59150549
59150550
59150551
59150552
59150553
59150554
59150555
59150556
59150557
59150558
59150559
59150560
59150561
59150562
59150563
59150564
59150565
59150566
59150567
59150568
59150569
59150570
59150571
59150572
59150573
59150575
59150577
59150578
59150579
59150580
59150581
59150582
59150583
59150584
59150585
59150586
59150587
59150588
59150589
59150590
59150591
59150592
59150593
59150594
59150595
59150596
59150597
59150598
59150599
59150600
59150601
59150602
59150603
59150604
59150605
59150606
59150607
59150608
59150609
59150610
59150611
59150612
59150613
59150614
59150615
59150616
59150617
59150618
59150619
59150620
59150621
59150622
59150623
59150624
59150625
59150626
59150627
59150628
59150629
59150630
59150631
59150632
59150633
59150634
59150635
59150636
59150637
59150638
59150639
59150640
59150641
59150642
59150643
59150644
59150645
59150646
59150647
59150648
59150649
59150650
59150651
59150652
59150653
59150654
59150655
59150656
59150657
59150658
59150659
59150660
59150661
59150662
59150663
59150664
59150665
59150666
59150667
59150668
59150669
59150670
59150671
59150672
59150673
59150674
59150675
59150676
59150677
59150678
59150679
59150680
59150681
59150682
59150683
59150684
59150685
59150686
59150687
59150688
59150689
59150690
59150691
59150692
59150693
59150694
59150695
59150696
59150697
59150698
59150699
59150700
59150701
59150702
59150703
59150704
59150705
59150706
59150707
59150708
59150709
59150710
59150711
59150712
59150713
59150714
59150715
59150716
59150717
59150718
59150719
59150720
59150721
59150722
59150723
59150724
59150725
59150726
59150727
59150728
59150729
59150730
59150731
59150732
59150733
59150734
59150735
59150736
59150737
59150738
59150739
59150740
59150741
59150742
59150743
59150744
59150745
59150746
59150747
59150748
59150749
59150750
59150751
59150752
59150753
59150754
59150755
59150756
59150757
59150758
59150759
59150760
59150761
59150762
59150763
59150777
59150779
59150780
59150782
59150784
59150785
59150786
59150787
59150788
59150790
59150791
59150792
59150793
59150795
59150816
59150817
59150819
59150820
59150821
59150822
59150823
59150824
59150827
59150828
59150831
59150832
59150837
59150840
59150841
59150845
59150846
59150850
59150859
59150861
59150863
59150864
59150865
59150866
59150871
59150872
59150873
59150874
59150875
59150876
59150877
59150878
59150879
59150880
59150881
59150882
59150883
59150884
59150885
59150886
59150887
59150888
59150889
59150890
59150891
59150892
59150893
59150894
59150895
59150896
59150897
59150898
59150899
59150900
59150901
59150902
59150903
59150904
59150905
59150906
59150907
59150908
59150909
59150910
59150911
59150912
59150913
59150914
59150915
59150916
59150917
59150918
59150919
59150920
59150922
59150923
59150924
59150925
59150926
59150927
59150928
59150929
59150930
59150931
59150932
59150933
59150934
59150935
59150953
59150954
59150955
59150956
59150957
59150958
59150959
59150960
59150961
59150962
59150963
59150964
59150965
59150966
59150967
59150968
59150969
59150970
59150971
59150972
59150973
59150974
59150975
59150976
59150977
59150978
59150979
59150980
59150981
59150982
59150983
59150984
59150985
59150986
59150987
59151198
59151199
59151200
59151201
59151202
59151203
59151204
59151205
59151206
59151207
59151208
59151209
59151211
59151213
59151214
59151215
59151216
59151217
59151219
59151220
59151221
59151222
59151223
59151224
59151225
59151226
59151227
59151228
59151229
59151230
59151231
59151234
59151236
59151351
59151352
59151353
59151354
59151355
59151356
59151357
59151358
59151359
59151360
59151361
59151362
59151363
59151364
59151365
59151366
59151367
59151368
59151369
59151370
59151371
59151372
59151373
59151407
59151408
59151409
59151410
59151411
59151412
59151413
59151414
59151415
59151416
59151417
59151418
59151419
59151420
59151423
59151424
59151425
59151426
59151427
59151428
59151429
59151430
59151432
59151433
59151434
59151435
59151436
59151437
59151438
59151439
59151440
59151441
59151442
59151443
59151444
59151448
59151449
59151450
59151451
59151452
59151453
59151454
59151455
59151456
59151457
59151458
59151459
59151460
59151461
59151463
59151464
59151465
59151467
59151469
59151470
59151471
59151472
59151473
59151474
59151475
59151476
59151477
59151478
59151479
59151480
59151481
59151482
59151483
59151484
59151485
59151486
59152941
59152942
59152943
59152944
59152945
59152946
59153075
59153076
59153077
59153078
59153079
59153080
59153081
59153082
59153083
59153084
59153085
59153086
59153087
59153088
59153089
59153090
59153091
59153092
59153093
59153094
59153095
59153096
59153097
59153098
59153123
59153124
59153132
59153133
59153134
59153135
59153136
59153137
59153138
59153139
59153140
59153141
59153142
59153144
59153145
59153146
59153151
59153152
59153153
59153154
59153155
59153156
59153157
59153158
59153159
59153160
59153161
59153162
59153163
59153164
59153165
59153166
59153167
59153168
59153171
59153172
59153181
59153182
59153183
59153184
59153185
59153186
59153187
59153188
59153189
59153190
59153191
59153192
59153193
59153194
59153195
59153196
59153197
59153198
59153199
59153200
59153201
59153202
59153203
59153204
59153205
59153206
59153207
59153208
59153209
59153210
59153211
59153212
59153213
59153214
59153215
59153216
59153217
59153218
59153219
59153220
59153221
59153222
59153223
59153224
59153225
59153226
59153227
59153228
59153229
59153230
59153231
59153232
59153246
59153249
59153250
59153251
59153252
59153253
59153254
59153255
59153256
59153257
59153258
59153259
59153261
59153262
59153263
59153264
59153265
59153266
59153267
59153268
59153269
59153270
59153271
59153272
59153273
59153274
59153275
59153276
59153277
59153278
59153279
59153280
59153281
59153282
59153283
59153284
59153285
59153286
59153287
59153288
59153289
59153290
59153291
59153292
59153293
59153294
59153295
59153296
59153297
59153298
59153299
59153300
59153301
59153302
59153303
59153304
59153305
59153306
59153307
59153308
59153370
59153371
59153372
59153373
59153394
59153400
59153464
59153466
59153467
59153469
59153470
59153471
59153472
59153473
59153474
59153476
59153477
59153479
59153480
59153481
59153482
59153483
59153490
59153491
59153492
59153498
59153499
59153500
59153539
59153540
59153552
59153553
59153562
59153578
59153579
59153580
59153581
59153582
59153583
59153584
59153585
59153586
59153587
59153588
59153589
59153590
59153591
59153668
59153669
59153670
59153671
59153672
59153673
59153674
59153675
59153676
59153677
59153678
59153679
59153680
59153681
59153682
59153683
59153684
59153685
59153686
59153687
59153688
59153689
59153690
59153691
59153692
59153799
59153800
59153801
59153805
59153806
59153807
59153809
59153810
59153812
59153814
59153815
59153816
59153818
59153821
59153825
59153826
59153831
59153832
59153833
59153834
59153835
59153836
59153838
59153839
59153843
59153845
59153846
59153847
59153849
59153851
59153852
59153955
59153956
59153986
59153987
59153988
59153989
59153990
59153991
59153994
59153995
59153996
59153997
59154000
59154001
59154002
59154026
59154027
Note how some bicycle infrastructure line segments follow the border of dissemination areas. One of the problems with trying to summarise transportation data within census polygons is that often the boundaries of census areas follow roads! So technically a given road belongs to more than one dissemination area. The intersection we performed above only assigned the attribute data from the neighbourhoods that intersected with a given road, a one-to-one spatial join. One way to deal with this would be to use the st_within_distance within st_join and create a polyline for each dissemination area that is within a certain distance of a line segment. This is called a one-to-many spatial join that will create multiple line segments for each dissemination where that spatial relationship is true. Let’s try it. If a line segment is within 5 metres of a dissemination area, join the attribute data from that dissemination area to the line segment. I chose 5 m to allow for some inaccuracy in the way the line data were digitized (e.g. a line segment may suppose to follow a boundary exactly but is actually a few metres off to one side etc.)
cut_join_da <- st_join(bicycle_cut, vancouver_da_sub,join = st_is_within_distance,dist=5)
cut_join_da## Simple feature collection with 2585 features and 4 fields
## Geometry type: GEOMETRY
## Dimension: XY
## Bounding box: xmin: 483720 ymin: 5449700 xmax: 498317.2 ymax: 5462376
## Projected CRS: NAD83 / UTM zone 10N
## First 10 features:
## osm_id GeoUID.x length_m GeoUID.y geometry
## 35260598 35260598 59150307 157.80005 [m] 59150307 LINESTRING (498105.5 545987...
## 35260598.1 35260598 59150307 157.80005 [m] 59153215 LINESTRING (498105.5 545987...
## 35260599 35260599 59150307 32.86117 [m] 59150307 LINESTRING (497911 5459475,...
## 36691655 36691655 59150307 163.86383 [m] 59150307 LINESTRING (497911 5459475,...
## 41840546 41840546 59150307 76.01887 [m] 59150307 LINESTRING (498180.4 545988...
## 75573893 75573893 59150307 44.49847 [m] 59150307 LINESTRING (498226.7 545987...
## 210638528 210638528 59150307 84.65905 [m] 59150307 LINESTRING (498190.6 545997...
## 210638530 210638530 59150307 46.27267 [m] 59150307 LINESTRING (498235.3 545998...
## 336365962 336365962 59150307 56.43094 [m] 59150307 LINESTRING (497942.1 545946...
## 340103260 340103260 59150307 86.09529 [m] 59150307 LINESTRING (497944.5 545952...Note how the first two rows of our data have the same osm_id but different values for GeoUID.y. What this means is that for this particular line segment there are multiple dissemination areas that fullfill the spatial join criteria (is within 5 metres). Now we let’s summarise the length_m column by the GeoUID.y column.
length_by_da <- cut_join_da %>%
st_drop_geometry() %>%
group_by(GeoUID.y) %>%
summarise(total_length_km = sum(as.numeric(length_m))/1000)
head(length_by_da)## # A tibble: 6 x 2
## GeoUID.y total_length_km
## <chr> <dbl>
## 1 59150307 1.34
## 2 59150309 0.218
## 3 59150312 0.0560
## 4 59150313 0.0939
## 5 59150328 0.171
## 6 59150329 1.30Now let’s join this back to our dissemination area boundaries, replace the NAs with 0s and map the results!
vancouver_da_sub <- left_join(vancouver_da_sub,length_by_da , by = c("GeoUID"="GeoUID.y")) %>%
mutate(total_length_km = replace_na(total_length_km,0))
vancouver_da_sub## Simple feature collection with 993 features and 2 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 483690.9 ymin: 5449535 xmax: 498329.2 ymax: 5462381
## Projected CRS: NAD83 / UTM zone 10N
## First 10 features:
## GeoUID total_length_km geometry
## 1 59150307 1.33568195 MULTIPOLYGON (((498308 5459...
## 2 59150308 0.00000000 MULTIPOLYGON (((498301.6 54...
## 3 59150309 0.21754280 MULTIPOLYGON (((497939.5 54...
## 4 59150310 0.00000000 MULTIPOLYGON (((498298.6 54...
## 5 59150311 0.00000000 MULTIPOLYGON (((498128.4 54...
## 6 59150312 0.05595795 MULTIPOLYGON (((498298.6 54...
## 7 59150313 0.09386538 MULTIPOLYGON (((498259 5458...
## 8 59150328 0.17094843 MULTIPOLYGON (((498102 5458...
## 9 59150329 1.29789720 MULTIPOLYGON (((497801.1 54...
## 10 59150330 0.00000000 MULTIPOLYGON (((497312 5458...mapview(vancouver_da_sub,zcol = "total_length_km") + mapview(bicycle_sub)103703105
104017171
104017200
105430161
106147360
106147365
10716909
111473385
111476393
111476395
113244334
113244335
113244336
113244337
113244340
113244341
113244347
113244349
113530128
115548889
115548894
115548899
115551776
115937368
115937372
115937376
115937386
115939809
115939816
115939820
115939827
115939839
115939846
115939852
115939861
115939863
115939882
115939885
115939886
115939888
115939900
115939901
115939905
115939919
115939920
115939924
115939925
116061590
116061596
116061597
116061606
116061621
116061622
116061624
116061646
116061665
116268241
116420526
116527975
116527979
116527981
116527982
116528908
116530434
116530438
116785439
116785444
116785447
116785452
116785453
116785456
116785460
116785461
116785462
116785471
116785474
116785477
116785480
116785491
116785494
116785508
116895611
116896363
116896364
116896366
116896369
116896370
116896395
116896401
116896405
116901106
116901117
116901143
116901159
116901161
116901169
116901174
116901175
116901189
116901196
117025317
117025318
117025319
117025320
117025321
117158480
117158494
117158495
117158496
117158498
117158499
117158505
117158510
117158520
117158522
117158524
117158525
117158527
117158530
117158536
117158538
117158539
117158544
117240591
117240605
117240617
117240620
117240628
117240640
117253401
117253402
117432867
117434878
117438891
117438892
117438894
117438897
117441217
117441219
118001597
118001607
118001611
118001613
118512109
118755262
118755291
118755292
119017978
119082683
119287310
119416940
119497568
119721430
119721431
119721433
119945024
120072705
120076089
120254690
120254701
120254703
120254730
120254732
120629449
121142936
121142937
121303930
121464609
121969172
124221070
125101809
125101810
125101811
125101812
125961205
125961208
125961217
125961234
125961243
126121346
126121443
128922681
128922682
136798884
136798888
136798893
136954730
137188671
142867204
149964049
154604519
158096937
159120756
159120764
159602117
159602122
159733979
159733985
167699458
167699460
167699474
167699475
167699477
167699479
167699482
167699483
171318256
171318665
171319122
172879289
174530914
174530915
185723076
190670675
209720964
209720965
209720966
209720968
210638528
210638530
220771277
222406667
222406668
223882644
22748108
22797667
228267241
22977244
230767157
230767158
230785269
232195215
23253924
23445307
241309375
24262171
248512403
255690779
255691343
256979886
258407325
258500297
258500303
259483127
26144993
26148108
279674920
286005317
286005464
287980354
28843819
28843820
28885937
28885938
28885939
28885940
28885941
28885943
28990478
28990716
290605242
290605250
290605256
290605266
290605273
290605277
29159329
29159394
29159633
292247972
292256507
294771945
295319449
295320813
295556516
295758726
295758727
295946666
295946668
297474257
298944326
302330395
302351830
302351831
302351832
304330884
304625748
30515621
305660601
305660602
305660608
305660609
305661211
307519154
307519157
307830374
307830375
307832359
30877356
309672081
314897632
319794427
320898812
320898825
322170298
322170299
327680822
334630810
336022374
336022942
33614302
336143339
336145997
336146171
336147286
336147662
336147663
336148158
336149034
336150160
336150161
336150534
336151145
336151146
336151186
336151212
336152904
336153031
336154719
336154738
336365962
336367020
336751642
336751643
336789306
336803184
336847509
336850008
336852774
336856341
336858408
336869485
336920236
336920237
337094804
337096482
337096894
337096927
337124793
337145339
337145340
337145341
337146058
337685824
337685826
337685830
337690233
337781607
337785030
338061608
338061609
338096912
338097900
338098189
338172095
338222476
338224261
338224263
338249439
338249440
338252554
338273163
338420820
338543609
338543614
338545883
338550030
339060621
339071119
340040619
340043117
340103260
340206883
340383378
340383379
340499518
340499519
340504714
341777090
341777093
346365689
35260598
35260599
35261341
354927072
357020124
357020125
358443047
358444420
358445049
358999768
359146644
35916647
35916651
359442309
359442311
359442313
359442315
359442688
359530859
360078666
360208600
360209022
360209024
360209386
360209388
360217193
360622582
361350681
361493594
361493595
361493845
361494010
361494011
361502535
361503814
361503815
361503816
361503817
361503818
361503819
361503820
361503875
361503876
361503877
361503878
361503880
361504326
363048942
363049007
363049009
363061588
363061589
363106002
363395764
363395766
363403138
363403139
363404507
363404581
363404582
363575804
363649494
363653528
363653529
363686272
363686274
363686510
36541268
36544891
36661000
36685463
36685489
36686826
36687246
36691654
36691655
36693635
36693642
36696019
367243448
36747459
36747746
36747750
36748091
36749981
367508899
367517778
367521930
367521931
367521963
36760149
36799387
368058222
368061172
368071150
368072528
368072529
36809173
368210443
368210444
368242970
36826535
36836280
36837997
36837999
368542358
368542359
368662929
368662932
368662933
368704135
36870938
36870942
368848682
36916130
369188715
369448967
369448968
36971913
36974363
36980456
36981919
36981920
36982309
370001880
37002218
37002219
37004623
37008964
37010042
37010044
370250108
37060337
370625583
370625586
37183383
37183384
37183385
37183386
37393070
374638951
374640036
374640037
374640515
375740755
376902798
376902800
376913078
376913080
379035685
379036259
379070932
379238007
379238008
379240422
379240423
379241368
379241369
379241371
379241372
379241373
379241375
379841289
379841290
380516622
381179592
381179593
383369305
383400043
383507077
383922158
383922162
384612221
384612223
385937656
389622143
389626533
389626534
389802479
389802484
389803812
389808100
389809705
389810449
389810451
389810452
389810456
391304585
391304586
391307535
391307536
391307537
391307538
391307540
391307542
391307544
393094031
39413845
398451942
398452650
398624648
39865209
39865218
39865222
39865229
39865233
40114662
402957423
40866546
40866552
408750120
408750121
408750122
408752084
40879853
409010861
41840546
41998485
41998486
41998487
41998492
41998493
421205814
421205815
42240571
422996754
424687805
429215457
429215459
431472888
436187403
437892174
437892532
437892533
437892534
437892535
437892546
437892939
437893146
437893864
437903187
437903188
437903189
437903190
437903191
438021295
438021975
438021976
439048348
43933078
44032491
447113476
448422875
450120196
45388564
45405212
457361746
457361747
457361748
457958562
460000575
460000576
460000578
460000580
460000584
465693508
4684064
468598004
4704933
47405380
475646248
475657872
475660411
476267898
476639340
476896738
476903048
477219539
477446962
487839273
489891205
489891207
489891211
489891218
489891221
489891223
489998650
491470422
492566528
49527152
49549360
49549363
49551583
49551868
49607820
49617609
496385508
49670225
497243697
497243698
497243699
497248382
497248383
49742836
497473841
497473843
49828813
49828819
49829204
49829212
499514386
499771914
500434368
500479271
500479272
500479273
50097995
50163490
503132613
503135267
503135268
503135269
503144467
503153729
503153731
503837937
503837938
504283287
504287837
504310510
505104675
505106308
506878610
510923680
510923681
511621373
511798517
51320843
51320845
514641651
514641652
514641653
514641654
514641655
514641656
51715442
51715448
51744538
51744557
51744559
51744593
518945617
521085315
521085319
521085322
521085329
521085332
521085336
521085339
521085342
521085343
521085344
52250426
52250429
526619907
527856883
530058743
53559237
53559260
538938313
540329767
5405080
541458751
54146362
550586926
550590503
552074822
552483034
563593926
567741896
56864257
58177324
58177326
58177332
58473342
58473343
586130245
588231507
589067903
592187684
592423618
596647511
606976115
606976118
606976119
607510884
607511914
607511915
608820557
608820558
610027438
610027439
610027440
61512385
615537606
615537607
62075466
62075467
621252878
621252879
62125826
62125827
621945344
621952387
621954166
622456090
626052402
631770722
634085474
636369371
639854987
639854988
639854989
641420431
641420447
643633773
644523271
653681127
658893551
661602531
673226421
676472319
676472320
677269115
677269116
678304531
681296736
681755674
681758647
682033606
684905424
685714750
685783904
685783905
685997602
690838917
691136520
691136522
691136523
693432045
693596154
694771136
696267599
697811523
699819938
699935828
699983092
699983093
700454853
700779788
700779789
700779790
700779793
700779794
70954668
70954671
70954672
70954674
70954677
70954679
71493588
71493594
729225829
729774097
729774099
729774100
730871282
731947170
731971184
731987072
732475903
73901009
73902061
740852282
74112853
74112862
74263740
74265305
74267956
74267959
74267962
74267967
74267970
74267972
74267974
74267975
746220721
74819440
74837530
74837539
74837543
74837551
74837554
74837556
748822399
748827656
748846436
748846437
750371203
750371206
750371207
75138132
75138177
75138182
75138196
751431203
751432982
751432987
75337223
75573893
75576268
75576270
75578008
75581675
75581702
757048226
75775688
75778322
75778327
75778328
75778330
75840914
76146118
76146120
76146127
76146129
76150423
76159666
76320583
76320622
76320625
76320652
76456081
76456085
76457094
76461710
764646059
768687624
771410371
771410372
77451356
774858742
774858743
776371399
776371400
776823054
782465259
782546366
785505589
785505590
789694693
789707415
789707421
789707423
789707425
790555076
790555077
791974722
791974723
796357054
797670684
798169853
798169856
803443589
803443590
803446768
803446770
804075658
805818821
80632623
80632625
806884833
809208845
809414002
809414003
809414004
809414010
809414011
809414012
809414016
809420567
809420571
81141119
81141120
81141123
815082906
815082916
81607344
816658462
816658463
816658464
816658620
816658621
816730795
817154138
818882898
818892753
818909352
820667373
821139852
821162576
82249874
824881442
825548178
825553878
825553879
825553880
825606442
825606443
825606444
825630957
825630958
826263586
826606988
827789218
827955833
827955834
827958010
827958011
827958314
827958315
827958316
827958317
827958445
827959007
827959008
828014626
828239103
828239104
828239117
828239120
828243449
828265744
828266366
828266367
828572601
828572603
82872016
830183824
830931208
831017791
831602139
831602141
831604045
831631121
831867177
831867178
83428213
83428214
83428215
83467624
835077990
835077991
835077992
835077993
835077994
835077995
835077996
835077997
835077998
835077999
83580863
83580871
83581201
836641455
836678217
838362439
839209124
839209128
839209130
839209131
839519397
839522323
842888496
845543871
845756374
845768532
845771505
846228524
846279016
84672142
847008541
850028422
850333887
85240509
85240514
85240520
85240523
85240524
85244180
853378466
85535825
856851806
856851807
856862174
856862175
856862176
856862177
856862178
859101445
859102478
859102479
859304114
859311424
859311426
859311429
859311431
859311432
859311433
860572471
860575835
862842076
863070919
863811844
863811845
863811846
863811847
863819776
863819777
863819779
863819780
863819781
863819782
863819786
863819787
863819788
863819791
863819792
863819795
863819796
863822654
863847084
863847086
864003835
865426194
867056949
867130223
867130224
872769239
872769240
872769241
872769242
872769246
872769247
872769248
872769249
874455936
874455937
874461788
874461789
874461790
874461791
874465022
874465023
874465024
874465025
874465026
874858897
876351626
877408899
877408901
877408903
877750117
878019455
878019456
87914022
87914036
880711654
880719215
880871611
880871612
88285291
884248889
885271854
88808788
888511057
89118741
891357219
89353357
89353365
89353372
89353378
89353415
89353421
89367597
899345983
900490294
900896674
905690152
905955372
905955373
905955376
905955377
907496777
907950878
908209929
911681217
913907337
913907338
914040718
914040720
915132409
915297457
915973911
915973912
916027710
916054075
916054076
916112197
916833991
916833992
916833994
917969890
917969896
917969897
917969898
917969899
917969900
918813454
919972793
919972794
919972795
919972796
919972798
919974826
919974827
919974830
92147661
921561009
921561010
921564449
92159343
92159348
92159349
921818859
921818861
921818862
92187712
921892891
921896328
922784577
922784578
922784579
922784580
922784581
922784582
922784583
922784584
922784586
922784588
922784589
922784591
922968299
924643737
924643738
925032669
925032670
926432727
926432728
926916045
926954993
926954994
926963103
927883927
927980391
928464732
928473654
929831446
930466265
931230698
931912468
931912469
931912470
931912471
935168401
935193928
935193929
935193930
935193931
939119936
939119937
940225673
944777727
944981743
95132019
96024722
96024723
96024729
96024734
96168096
96168097
96168098
96168100
96168102
96921154
96993481
97278473
97337324
97337373
98020497
99144095
We could also map the density instead:
vancouver_da_sub <- vancouver_da_sub %>%
mutate(total_area_km2 = as.numeric(st_area(.)/1000/1000),
total_bike_dens = total_length_km/ total_area_km2
)
mapview(vancouver_da_sub,zcol = "total_bike_dens") + mapview(bicycle_sub)103703105
104017171
104017200
105430161
106147360
106147365
10716909
111473385
111476393
111476395
113244334
113244335
113244336
113244337
113244340
113244341
113244347
113244349
113530128
115548889
115548894
115548899
115551776
115937368
115937372
115937376
115937386
115939809
115939816
115939820
115939827
115939839
115939846
115939852
115939861
115939863
115939882
115939885
115939886
115939888
115939900
115939901
115939905
115939919
115939920
115939924
115939925
116061590
116061596
116061597
116061606
116061621
116061622
116061624
116061646
116061665
116268241
116420526
116527975
116527979
116527981
116527982
116528908
116530434
116530438
116785439
116785444
116785447
116785452
116785453
116785456
116785460
116785461
116785462
116785471
116785474
116785477
116785480
116785491
116785494
116785508
116895611
116896363
116896364
116896366
116896369
116896370
116896395
116896401
116896405
116901106
116901117
116901143
116901159
116901161
116901169
116901174
116901175
116901189
116901196
117025317
117025318
117025319
117025320
117025321
117158480
117158494
117158495
117158496
117158498
117158499
117158505
117158510
117158520
117158522
117158524
117158525
117158527
117158530
117158536
117158538
117158539
117158544
117240591
117240605
117240617
117240620
117240628
117240640
117253401
117253402
117432867
117434878
117438891
117438892
117438894
117438897
117441217
117441219
118001597
118001607
118001611
118001613
118512109
118755262
118755291
118755292
119017978
119082683
119287310
119416940
119497568
119721430
119721431
119721433
119945024
120072705
120076089
120254690
120254701
120254703
120254730
120254732
120629449
121142936
121142937
121303930
121464609
121969172
124221070
125101809
125101810
125101811
125101812
125961205
125961208
125961217
125961234
125961243
126121346
126121443
128922681
128922682
136798884
136798888
136798893
136954730
137188671
142867204
149964049
154604519
158096937
159120756
159120764
159602117
159602122
159733979
159733985
167699458
167699460
167699474
167699475
167699477
167699479
167699482
167699483
171318256
171318665
171319122
172879289
174530914
174530915
185723076
190670675
209720964
209720965
209720966
209720968
210638528
210638530
220771277
222406667
222406668
223882644
22748108
22797667
228267241
22977244
230767157
230767158
230785269
232195215
23253924
23445307
241309375
24262171
248512403
255690779
255691343
256979886
258407325
258500297
258500303
259483127
26144993
26148108
279674920
286005317
286005464
287980354
28843819
28843820
28885937
28885938
28885939
28885940
28885941
28885943
28990478
28990716
290605242
290605250
290605256
290605266
290605273
290605277
29159329
29159394
29159633
292247972
292256507
294771945
295319449
295320813
295556516
295758726
295758727
295946666
295946668
297474257
298944326
302330395
302351830
302351831
302351832
304330884
304625748
30515621
305660601
305660602
305660608
305660609
305661211
307519154
307519157
307830374
307830375
307832359
30877356
309672081
314897632
319794427
320898812
320898825
322170298
322170299
327680822
334630810
336022374
336022942
33614302
336143339
336145997
336146171
336147286
336147662
336147663
336148158
336149034
336150160
336150161
336150534
336151145
336151146
336151186
336151212
336152904
336153031
336154719
336154738
336365962
336367020
336751642
336751643
336789306
336803184
336847509
336850008
336852774
336856341
336858408
336869485
336920236
336920237
337094804
337096482
337096894
337096927
337124793
337145339
337145340
337145341
337146058
337685824
337685826
337685830
337690233
337781607
337785030
338061608
338061609
338096912
338097900
338098189
338172095
338222476
338224261
338224263
338249439
338249440
338252554
338273163
338420820
338543609
338543614
338545883
338550030
339060621
339071119
340040619
340043117
340103260
340206883
340383378
340383379
340499518
340499519
340504714
341777090
341777093
346365689
35260598
35260599
35261341
354927072
357020124
357020125
358443047
358444420
358445049
358999768
359146644
35916647
35916651
359442309
359442311
359442313
359442315
359442688
359530859
360078666
360208600
360209022
360209024
360209386
360209388
360217193
360622582
361350681
361493594
361493595
361493845
361494010
361494011
361502535
361503814
361503815
361503816
361503817
361503818
361503819
361503820
361503875
361503876
361503877
361503878
361503880
361504326
363048942
363049007
363049009
363061588
363061589
363106002
363395764
363395766
363403138
363403139
363404507
363404581
363404582
363575804
363649494
363653528
363653529
363686272
363686274
363686510
36541268
36544891
36661000
36685463
36685489
36686826
36687246
36691654
36691655
36693635
36693642
36696019
367243448
36747459
36747746
36747750
36748091
36749981
367508899
367517778
367521930
367521931
367521963
36760149
36799387
368058222
368061172
368071150
368072528
368072529
36809173
368210443
368210444
368242970
36826535
36836280
36837997
36837999
368542358
368542359
368662929
368662932
368662933
368704135
36870938
36870942
368848682
36916130
369188715
369448967
369448968
36971913
36974363
36980456
36981919
36981920
36982309
370001880
37002218
37002219
37004623
37008964
37010042
37010044
370250108
37060337
370625583
370625586
37183383
37183384
37183385
37183386
37393070
374638951
374640036
374640037
374640515
375740755
376902798
376902800
376913078
376913080
379035685
379036259
379070932
379238007
379238008
379240422
379240423
379241368
379241369
379241371
379241372
379241373
379241375
379841289
379841290
380516622
381179592
381179593
383369305
383400043
383507077
383922158
383922162
384612221
384612223
385937656
389622143
389626533
389626534
389802479
389802484
389803812
389808100
389809705
389810449
389810451
389810452
389810456
391304585
391304586
391307535
391307536
391307537
391307538
391307540
391307542
391307544
393094031
39413845
398451942
398452650
398624648
39865209
39865218
39865222
39865229
39865233
40114662
402957423
40866546
40866552
408750120
408750121
408750122
408752084
40879853
409010861
41840546
41998485
41998486
41998487
41998492
41998493
421205814
421205815
42240571
422996754
424687805
429215457
429215459
431472888
436187403
437892174
437892532
437892533
437892534
437892535
437892546
437892939
437893146
437893864
437903187
437903188
437903189
437903190
437903191
438021295
438021975
438021976
439048348
43933078
44032491
447113476
448422875
450120196
45388564
45405212
457361746
457361747
457361748
457958562
460000575
460000576
460000578
460000580
460000584
465693508
4684064
468598004
4704933
47405380
475646248
475657872
475660411
476267898
476639340
476896738
476903048
477219539
477446962
487839273
489891205
489891207
489891211
489891218
489891221
489891223
489998650
491470422
492566528
49527152
49549360
49549363
49551583
49551868
49607820
49617609
496385508
49670225
497243697
497243698
497243699
497248382
497248383
49742836
497473841
497473843
49828813
49828819
49829204
49829212
499514386
499771914
500434368
500479271
500479272
500479273
50097995
50163490
503132613
503135267
503135268
503135269
503144467
503153729
503153731
503837937
503837938
504283287
504287837
504310510
505104675
505106308
506878610
510923680
510923681
511621373
511798517
51320843
51320845
514641651
514641652
514641653
514641654
514641655
514641656
51715442
51715448
51744538
51744557
51744559
51744593
518945617
521085315
521085319
521085322
521085329
521085332
521085336
521085339
521085342
521085343
521085344
52250426
52250429
526619907
527856883
530058743
53559237
53559260
538938313
540329767
5405080
541458751
54146362
550586926
550590503
552074822
552483034
563593926
567741896
56864257
58177324
58177326
58177332
58473342
58473343
586130245
588231507
589067903
592187684
592423618
596647511
606976115
606976118
606976119
607510884
607511914
607511915
608820557
608820558
610027438
610027439
610027440
61512385
615537606
615537607
62075466
62075467
621252878
621252879
62125826
62125827
621945344
621952387
621954166
622456090
626052402
631770722
634085474
636369371
639854987
639854988
639854989
641420431
641420447
643633773
644523271
653681127
658893551
661602531
673226421
676472319
676472320
677269115
677269116
678304531
681296736
681755674
681758647
682033606
684905424
685714750
685783904
685783905
685997602
690838917
691136520
691136522
691136523
693432045
693596154
694771136
696267599
697811523
699819938
699935828
699983092
699983093
700454853
700779788
700779789
700779790
700779793
700779794
70954668
70954671
70954672
70954674
70954677
70954679
71493588
71493594
729225829
729774097
729774099
729774100
730871282
731947170
731971184
731987072
732475903
73901009
73902061
740852282
74112853
74112862
74263740
74265305
74267956
74267959
74267962
74267967
74267970
74267972
74267974
74267975
746220721
74819440
74837530
74837539
74837543
74837551
74837554
74837556
748822399
748827656
748846436
748846437
750371203
750371206
750371207
75138132
75138177
75138182
75138196
751431203
751432982
751432987
75337223
75573893
75576268
75576270
75578008
75581675
75581702
757048226
75775688
75778322
75778327
75778328
75778330
75840914
76146118
76146120
76146127
76146129
76150423
76159666
76320583
76320622
76320625
76320652
76456081
76456085
76457094
76461710
764646059
768687624
771410371
771410372
77451356
774858742
774858743
776371399
776371400
776823054
782465259
782546366
785505589
785505590
789694693
789707415
789707421
789707423
789707425
790555076
790555077
791974722
791974723
796357054
797670684
798169853
798169856
803443589
803443590
803446768
803446770
804075658
805818821
80632623
80632625
806884833
809208845
809414002
809414003
809414004
809414010
809414011
809414012
809414016
809420567
809420571
81141119
81141120
81141123
815082906
815082916
81607344
816658462
816658463
816658464
816658620
816658621
816730795
817154138
818882898
818892753
818909352
820667373
821139852
821162576
82249874
824881442
825548178
825553878
825553879
825553880
825606442
825606443
825606444
825630957
825630958
826263586
826606988
827789218
827955833
827955834
827958010
827958011
827958314
827958315
827958316
827958317
827958445
827959007
827959008
828014626
828239103
828239104
828239117
828239120
828243449
828265744
828266366
828266367
828572601
828572603
82872016
830183824
830931208
831017791
831602139
831602141
831604045
831631121
831867177
831867178
83428213
83428214
83428215
83467624
835077990
835077991
835077992
835077993
835077994
835077995
835077996
835077997
835077998
835077999
83580863
83580871
83581201
836641455
836678217
838362439
839209124
839209128
839209130
839209131
839519397
839522323
842888496
845543871
845756374
845768532
845771505
846228524
846279016
84672142
847008541
850028422
850333887
85240509
85240514
85240520
85240523
85240524
85244180
853378466
85535825
856851806
856851807
856862174
856862175
856862176
856862177
856862178
859101445
859102478
859102479
859304114
859311424
859311426
859311429
859311431
859311432
859311433
860572471
860575835
862842076
863070919
863811844
863811845
863811846
863811847
863819776
863819777
863819779
863819780
863819781
863819782
863819786
863819787
863819788
863819791
863819792
863819795
863819796
863822654
863847084
863847086
864003835
865426194
867056949
867130223
867130224
872769239
872769240
872769241
872769242
872769246
872769247
872769248
872769249
874455936
874455937
874461788
874461789
874461790
874461791
874465022
874465023
874465024
874465025
874465026
874858897
876351626
877408899
877408901
877408903
877750117
878019455
878019456
87914022
87914036
880711654
880719215
880871611
880871612
88285291
884248889
885271854
88808788
888511057
89118741
891357219
89353357
89353365
89353372
89353378
89353415
89353421
89367597
899345983
900490294
900896674
905690152
905955372
905955373
905955376
905955377
907496777
907950878
908209929
911681217
913907337
913907338
914040718
914040720
915132409
915297457
915973911
915973912
916027710
916054075
916054076
916112197
916833991
916833992
916833994
917969890
917969896
917969897
917969898
917969899
917969900
918813454
919972793
919972794
919972795
919972796
919972798
919974826
919974827
919974830
92147661
921561009
921561010
921564449
92159343
92159348
92159349
921818859
921818861
921818862
92187712
921892891
921896328
922784577
922784578
922784579
922784580
922784581
922784582
922784583
922784584
922784586
922784588
922784589
922784591
922968299
924643737
924643738
925032669
925032670
926432727
926432728
926916045
926954993
926954994
926963103
927883927
927980391
928464732
928473654
929831446
930466265
931230698
931912468
931912469
931912470
931912471
935168401
935193928
935193929
935193930
935193931
939119936
939119937
940225673
944777727
944981743
95132019
96024722
96024723
96024729
96024734
96168096
96168097
96168098
96168100
96168102
96921154
96993481
97278473
97337324
97337373
98020497
99144095