1 Libraries

library(knitr)
library(sf)
library(dplyr)
library(ggplot2)
library(bcdata)
library(ggspatial)
library(bcmaps)
library(here)
library(mapview)
library(units)

2 Tidy

2.1 Select and filter

starwars %>%
  select(name, height, hair_color, homeworld) %>%
  filter(homeworld == "Tatooine")
## # A tibble: 10 × 4
##    name               height hair_color  homeworld
##    <chr>               <int> <chr>       <chr>    
##  1 Luke Skywalker        172 blond       Tatooine 
##  2 C-3PO                 167 <NA>        Tatooine 
##  3 Darth Vader           202 none        Tatooine 
##  4 Owen Lars             178 brown, grey Tatooine 
##  5 Beru Whitesun Lars    165 brown       Tatooine 
##  6 R5-D4                  97 <NA>        Tatooine 
##  7 Biggs Darklighter     183 black       Tatooine 
##  8 Anakin Skywalker      188 blond       Tatooine 
##  9 Shmi Skywalker        163 black       Tatooine 
## 10 Cliegg Lars           183 brown       Tatooine

2.2 Group and summarize

# Grouping and Summarizing
starwars %>%
  tidyr::drop_na() %>% 
  filter(species == "Human") %>%
  group_by(homeworld) %>%
  summarise(mean_mass = mean(mass, na.rm = TRUE)) %>%
  arrange(desc(mean_mass))
## # A tibble: 11 × 2
##    homeworld    mean_mass
##    <chr>            <dbl>
##  1 Tatooine          96  
##  2 Haruun Kal        84  
##  3 Serenno           80  
##  4 Bespin            79  
##  5 Concord Dawn      79  
##  6 Socorro           79  
##  7 Corellia          78.5
##  8 Kamino            78.2
##  9 Stewjon           77  
## 10 Naboo             60  
## 11 Alderaan          49

3 Spatial

3.1 bcdata and points

airports <- bcdc_query_geodata("bc-airports") %>% collect()
## Warning: It is advised to use the permanent id ('76b1b7a3-2112-4444-857a-afccf7b20da8') rather than the name of the record ('bc-airports') to guard against future name changes.
airports
## Simple feature collection with 455 features and 41 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 406543.7 ymin: 367957.6 xmax: 1796645 ymax: 1689146
## Projected CRS: NAD83 / BC Albers
## # A tibble: 455 × 42
##    id       CUSTODIAN_ORG_DESCRI…¹ BUSINESS_CATEGORY_CL…² BUSINESS_CATEGORY_DE…³
##  * <chr>    <chr>                  <chr>                  <chr>                 
##  1 WHSE_IM… "Ministry of Forest, … airTransportation      Air Transportation    
##  2 WHSE_IM… "Ministry of Forest, … airTransportation      Air Transportation    
##  3 WHSE_IM… "Ministry of Forest, … airTransportation      Air Transportation    
##  4 WHSE_IM… "Ministry of Forest, … airTransportation      Air Transportation    
##  5 WHSE_IM… "Ministry of Forest, … airTransportation      Air Transportation    
##  6 WHSE_IM… "Ministry of Forest, … airTransportation      Air Transportation    
##  7 WHSE_IM… "Ministry of Forest, … airTransportation      Air Transportation    
##  8 WHSE_IM… "Ministry of Forest, … airTransportation      Air Transportation    
##  9 WHSE_IM… "Ministry of Forest, … airTransportation      Air Transportation    
## 10 WHSE_IM… "Ministry of Forest, … airTransportation      Air Transportation    
## # ℹ 445 more rows
## # ℹ abbreviated names: ¹​CUSTODIAN_ORG_DESCRIPTION, ²​BUSINESS_CATEGORY_CLASS,
## #   ³​BUSINESS_CATEGORY_DESCRIPTION
## # ℹ 38 more variables: OCCUPANT_TYPE_DESCRIPTION <chr>, SOURCE_DATA_ID <chr>,
## #   SUPPLIED_SOURCE_ID_IND <chr>, AIRPORT_NAME <chr>, DESCRIPTION <chr>,
## #   PHYSICAL_ADDRESS <chr>, ALIAS_ADDRESS <chr>, STREET_ADDRESS <chr>,
## #   POSTAL_CODE <chr>, LOCALITY <chr>, CONTACT_PHONE <chr>, …
mapview(airports, zcol = "NUMBER_OF_RUNWAYS")
st_geometry(airports)
## Geometry set for 455 features 
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 406543.7 ymin: 367957.6 xmax: 1796645 ymax: 1689146
## Projected CRS: NAD83 / BC Albers
## First 5 geometries:
## POINT (833323.9 1054950)
## POINT (1193727 381604.1)
## POINT (1194902 382257.7)
## POINT (1193719 382179.3)
## POINT (1198292 383563.6)
st_bbox(airports)
##      xmin      ymin      xmax      ymax 
##  406543.7  367957.6 1796645.0 1689145.9
st_crs(airports)
## Coordinate Reference System:
##   User input: NAD83 / BC Albers 
##   wkt:
## PROJCRS["NAD83 / BC Albers",
##     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["British Columbia Albers",
##         METHOD["Albers Equal Area",
##             ID["EPSG",9822]],
##         PARAMETER["Latitude of false origin",45,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8821]],
##         PARAMETER["Longitude of false origin",-126,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8822]],
##         PARAMETER["Latitude of 1st standard parallel",50,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8823]],
##         PARAMETER["Latitude of 2nd standard parallel",58.5,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8824]],
##         PARAMETER["Easting at false origin",1000000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8826]],
##         PARAMETER["Northing at false origin",0,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8827]]],
##     CS[Cartesian,2],
##         AXIS["(E)",east,
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["(N)",north,
##             ORDER[2],
##             LENGTHUNIT["metre",1]],
##     USAGE[
##         SCOPE["Province-wide spatial data management."],
##         AREA["Canada - British Columbia."],
##         BBOX[48.25,-139.04,60.01,-114.08]],
##     ID["EPSG",3005]]

3.2 Make a point, and transform

albers_centre <- st_sfc(st_point(c(-126, 54)), crs = 4326) %>% 
  st_transform(3005)
albers_centre
## Geometry set for 1 feature 
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 1e+06 ymin: 999610.6 xmax: 1e+06 ymax: 999610.6
## Projected CRS: NAD83 / BC Albers
## POINT (1e+06 999610.6)

3.3 Filter spaial

bc_cities_filtered <- bc_cities() %>% 
  filter(NAME == "Prince George")
## bc_cities was updated on 2024-01-13
bc_cities_filtered
## Simple feature collection with 1 feature and 10 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 1215604 ymin: 994614.2 xmax: 1215604 ymax: 994614.2
## Projected CRS: NAD83 / BC Albers
## # A tibble: 1 × 11
##   id       FCODE BCMJ_TAG NAME  CITY_TYPE LONG_TYPE POP_2000 POP_SOURCE OBJECTID
## * <chr>    <chr>    <int> <chr> <chr>     <chr>        <int> <chr>         <int>
## 1 WHSE_BA… AR05…       24 Prin… C         CITY         81326 BCSTATS          24
## # ℹ 2 more variables: SE_ANNO_CAD_DATA <chr>, geometry <POINT [m]>

3.4 Spatial plot

ggplot() + 
  geom_sf(data = bc_bound()) +
  geom_sf(data = albers_centre, colour = "blue", size = 5) +
  geom_sf(data = bc_cities_filtered) + 
  geom_sf_text(data = bc_cities_filtered, aes(label = NAME), 
               nudge_x = 100000, nudge_y = 30000) +
  theme_void()

3.5 Polygons

elec_bc <- bcdc_query_geodata("provincial-electoral-districts-electoral-boundaries-redistribution-2015") %>% collect() 

ggplot() + geom_sf(data = elec_bc, aes(fill = FEATURE_AREA_SQM ))

3.6 Spatial Join

cities_by_nr <- bc_cities() %>% 
  st_join(nr_districts(), join = st_intersects) %>%
  group_by(REGION_ORG_UNIT_NAME) %>%
  summarise(nr_pop_2000 = sum(POP_2000))
ggplot(cities_by_nr) +
  geom_col(aes(x = REGION_ORG_UNIT_NAME, 
               y = nr_pop_2000/1000000)) +
  coord_flip() +
  theme_minimal() + 
  labs(x = "", y = "Pop (in Millions)")

3.7 Buffers

bc_cities_buffer <- st_buffer(bc_cities(), dist = 20000)
bc_cities_buffer
## Simple feature collection with 166 features and 10 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 543199.8 ymin: 355965.9 xmax: 1818382 ymax: 1663181
## Projected CRS: NAD83 / BC Albers
## # A tibble: 166 × 11
##    id      FCODE BCMJ_TAG NAME  CITY_TYPE LONG_TYPE POP_2000 POP_SOURCE OBJECTID
##  * <chr>   <chr>    <int> <chr> <chr>     <chr>        <int> <chr>         <int>
##  1 WHSE_B… AR08…        1 Uclu… DM        DISTRICT…     1824 BCSTATS           1
##  2 WHSE_B… AR08…        2 Tofi… DM        DISTRICT…     1540 BCSTATS           2
##  3 WHSE_B… AR32…        3 Burn… VL        VILLAGE       1892 BCSTATS           3
##  4 WHSE_B… AR08…        4 Lill… DM        DISTRICT…     2987 BCSTATS           4
##  5 WHSE_B… AR32…        5 Lytt… VL        VILLAGE        318 BCSTATS           5
##  6 WHSE_B… AR05…        6 Merr… C         CITY          8075 BCSTATS           6
##  7 WHSE_B… AR08…        7 Inve… DM        DISTRICT…     2952 BCSTATS           7
##  8 WHSE_B… AR32…        8 New … U         VILLAGE …      800 OTHER             8
##  9 WHSE_B… AR32…        9 Ashc… VL        VILLAGE       1995 BCSTATS           9
## 10 WHSE_B… AR32…       10 Cach… VL        VILLAGE       1114 BCSTATS          10
## # ℹ 156 more rows
## # ℹ 2 more variables: SE_ANNO_CAD_DATA <chr>, geometry <POLYGON [m]>
ggplot() +
  geom_sf(data = bc_bound()) +
  geom_sf(data = bc_cities()) +
  geom_sf(data = bc_cities_buffer, fill = "green", alpha = 0.5)

4 Some Vector Formats

glimpse(bc_cities())
## Rows: 166
## Columns: 11
## $ id               <chr> "WHSE_BASEMAPPING.BC_MAJOR_CITIES_POINTS_500M.1", "WH…
## $ FCODE            <chr> "AR08750000", "AR08750000", "AR32700000", "AR08750000…
## $ BCMJ_TAG         <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16…
## $ NAME             <chr> "Ucluelet", "Tofino", "Burns Lake", "Lillooet", "Lytt…
## $ CITY_TYPE        <chr> "DM", "DM", "VL", "DM", "VL", "C", "DM", "U", "VL", "…
## $ LONG_TYPE        <chr> "DISTRICT MUNICIPALITY", "DISTRICT MUNICIPALITY", "VI…
## $ POP_2000         <int> 1824, 1540, 1892, 2987, 318, 8075, 2952, 800, 1995, 1…
## $ POP_SOURCE       <chr> "BCSTATS", "BCSTATS", "BCSTATS", "BCSTATS", "BCSTATS"…
## $ OBJECTID         <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16…
## $ SE_ANNO_CAD_DATA <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ geometry         <POINT [m]> POINT (1036351 441199.7), POINT (1008983 455516…

4.1 KML

bc_cities() %>% st_write("bc_cities.kml", append = F)
## Writing layer `bc_cities' to data source `bc_cities.kml' using driver `KML'
## Writing 166 features with 10 fields and geometry type Point.
glimpse(st_read("bc_cities.kml"))
## Reading layer `bc_cities' from data source `C:\Users\bevin\Desktop\bc_cities.kml' using driver `KML'
## Simple feature collection with 166 features and 2 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -133.7266 ymin: 48.375 xmax: -114.8762 ymax: 59.55726
## Geodetic CRS:  WGS 84
## Rows: 166
## Columns: 3
## $ Name        <chr> "Ucluelet", "Tofino", "Burns Lake", "Lillooet", "Lytton", …
## $ Description <chr> "", "", "", "", "", "", "", "", "", "", "", "", "", "", ""…
## $ geometry    <POINT [°]> POINT (-125.504 48.98594), POINT (-125.8771 49.11591…

4.2 SHP

bc_cities() %>% st_write("bc_cities.shp", append = F)
## Warning in abbreviate_shapefile_names(obj): Field names abbreviated for ESRI
## Shapefile driver
## Deleting layer `bc_cities' using driver `ESRI Shapefile'
## Writing layer `bc_cities' to data source `bc_cities.shp' using driver `ESRI Shapefile'
## Writing 166 features with 10 fields and geometry type Point.
glimpse(st_read("bc_cities.shp")) 
## Reading layer `bc_cities' from data source `C:\Users\bevin\Desktop\bc_cities.shp' using driver `ESRI Shapefile'
## Simple feature collection with 166 features and 10 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 563199.8 ymin: 375965.9 xmax: 1798382 ymax: 1643181
## Projected CRS: NAD83 / BC Albers
## Rows: 166
## Columns: 11
## $ id       <chr> "WHSE_BASEMAPPING.BC_MAJOR_CITIES_POINTS_500M.1", "WHSE_BASEM…
## $ FCODE    <chr> "AR08750000", "AR08750000", "AR32700000", "AR08750000", "AR32…
## $ BCMJ_TA  <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18…
## $ NAME     <chr> "Ucluelet", "Tofino", "Burns Lake", "Lillooet", "Lytton", "Me…
## $ CITY_TY  <chr> "DM", "DM", "VL", "DM", "VL", "C", "DM", "U", "VL", "VL", "DM…
## $ LONG_TY  <chr> "DISTRICT MUNICIPALITY", "DISTRICT MUNICIPALITY", "VILLAGE", …
## $ POP_200  <int> 1824, 1540, 1892, 2987, 318, 8075, 2952, 800, 1995, 1114, 416…
## $ POP_SOU  <chr> "BCSTATS", "BCSTATS", "BCSTATS", "BCSTATS", "BCSTATS", "BCSTA…
## $ OBJECTI  <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18…
## $ SE_ANNO  <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ geometry <POINT [m]> POINT (1036351 441199.7), POINT (1008983 455516.9), POI…

4.3 GPKG

bc_cities() %>% st_write("bc_cities.gpkg", append = F)
## Deleting layer `bc_cities' using driver `GPKG'
## Writing layer `bc_cities' to data source `bc_cities.gpkg' using driver `GPKG'
## Writing 166 features with 10 fields and geometry type Point.
glimpse(st_read("bc_cities.gpkg")) 
## Reading layer `bc_cities' from data source `C:\Users\bevin\Desktop\bc_cities.gpkg' using driver `GPKG'
## Simple feature collection with 166 features and 10 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 563199.8 ymin: 375965.9 xmax: 1798382 ymax: 1643181
## Projected CRS: NAD83 / BC Albers
## Rows: 166
## Columns: 11
## $ id               <chr> "WHSE_BASEMAPPING.BC_MAJOR_CITIES_POINTS_500M.1", "WH…
## $ FCODE            <chr> "AR08750000", "AR08750000", "AR32700000", "AR08750000…
## $ BCMJ_TAG         <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16…
## $ NAME             <chr> "Ucluelet", "Tofino", "Burns Lake", "Lillooet", "Lytt…
## $ CITY_TYPE        <chr> "DM", "DM", "VL", "DM", "VL", "C", "DM", "U", "VL", "…
## $ LONG_TYPE        <chr> "DISTRICT MUNICIPALITY", "DISTRICT MUNICIPALITY", "VI…
## $ POP_2000         <int> 1824, 1540, 1892, 2987, 318, 8075, 2952, 800, 1995, 1…
## $ POP_SOURCE       <chr> "BCSTATS", "BCSTATS", "BCSTATS", "BCSTATS", "BCSTATS"…
## $ OBJECTID         <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16…
## $ SE_ANNO_CAD_DATA <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ geom             <POINT [m]> POINT (1036351 441199.7), POINT (1008983 455516…

4.4 SQLITE

bc_cities() %>% st_write("bc_cities.sqlite", append = F)
## Deleting layer `bc_cities' using driver `SQLite'
## Writing layer `bc_cities' to data source `bc_cities.sqlite' using driver `SQLite'
## Writing 166 features with 10 fields and geometry type Point.
glimpse(st_read("bc_cities.sqlite")) 
## Reading layer `bc_cities' from data source 
##   `C:\Users\bevin\Desktop\bc_cities.sqlite' using driver `SQLite'
## Simple feature collection with 166 features and 10 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 563199.8 ymin: 375965.9 xmax: 1798382 ymax: 1643181
## Projected CRS: NAD83 / BC Albers
## Rows: 166
## Columns: 11
## $ id               <chr> "WHSE_BASEMAPPING.BC_MAJOR_CITIES_POINTS_500M.1", "WH…
## $ fcode            <chr> "AR08750000", "AR08750000", "AR32700000", "AR08750000…
## $ bcmj_tag         <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16…
## $ name             <chr> "Ucluelet", "Tofino", "Burns Lake", "Lillooet", "Lytt…
## $ city_type        <chr> "DM", "DM", "VL", "DM", "VL", "C", "DM", "U", "VL", "…
## $ long_type        <chr> "DISTRICT MUNICIPALITY", "DISTRICT MUNICIPALITY", "VI…
## $ pop_2000         <int> 1824, 1540, 1892, 2987, 318, 8075, 2952, 800, 1995, 1…
## $ pop_source       <chr> "BCSTATS", "BCSTATS", "BCSTATS", "BCSTATS", "BCSTATS"…
## $ objectid         <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16…
## $ se_anno_cad_data <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ GEOMETRY         <POINT [m]> POINT (1036351 441199.7), POINT (1008983 455516…

4.5 Benchmark

microbenchmark::microbenchmark(
  st_read("bc_cities.kml", quiet = T), 
  st_read("bc_cities.shp", quiet = T), 
  st_read("bc_cities.gpkg", quiet = T), 
  st_read("bc_cities.sqlite", quiet = T))
## Unit: milliseconds
##                                    expr     min       lq      mean   median
##     st_read("bc_cities.kml", quiet = T) 18.7737 19.21415 19.354246 19.33640
##     st_read("bc_cities.shp", quiet = T) 10.3774 10.64790 10.860667 10.76665
##    st_read("bc_cities.gpkg", quiet = T) 11.3659 11.75765 11.967699 11.90985
##  st_read("bc_cities.sqlite", quiet = T)  8.8311  9.18075  9.512164  9.36680
##        uq     max neval  cld
##  19.48795 19.9480   100 a   
##  11.08690 12.5766   100  b  
##  12.09145 12.9420   100   c 
##   9.54240 15.4946   100    d

5 Querying

5.1 Example

glimpse(
  st_read("bc_cities.gpkg", 
          query = "SELECT * FROM bc_cities WHERE NAME == 'Prince George'"))
## Reading query `SELECT * FROM bc_cities WHERE NAME == 'Prince George''
## from data source `C:\Users\bevin\Desktop\bc_cities.gpkg' using driver `GPKG'
## Simple feature collection with 1 feature and 10 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 1215604 ymin: 994614.2 xmax: 1215604 ymax: 994614.2
## Projected CRS: NAD83 / BC Albers
## Rows: 1
## Columns: 11
## $ id               <chr> "WHSE_BASEMAPPING.BC_MAJOR_CITIES_POINTS_500M.24"
## $ FCODE            <chr> "AR05500000"
## $ BCMJ_TAG         <int> 24
## $ NAME             <chr> "Prince George"
## $ CITY_TYPE        <chr> "C"
## $ LONG_TYPE        <chr> "CITY"
## $ POP_2000         <int> 81326
## $ POP_SOURCE       <chr> "BCSTATS"
## $ OBJECTID         <int> 24
## $ SE_ANNO_CAD_DATA <chr> NA
## $ geom             <POINT [m]> POINT (1215604 994614.2)
glimpse(
  st_read("bc_cities.sqlite", 
          query = "SELECT * FROM bc_cities WHERE city_type == 'DM'"))
## Reading query `SELECT * FROM bc_cities WHERE city_type == 'DM''
## from data source `C:\Users\bevin\Desktop\bc_cities.sqlite' using driver `SQLite'
## Simple feature collection with 54 features and 10 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 732355 ymin: 375965.9 xmax: 1798382 ymax: 1252807
## Projected CRS: NAD83 / BC Albers
## Rows: 54
## Columns: 11
## $ id               <chr> "WHSE_BASEMAPPING.BC_MAJOR_CITIES_POINTS_500M.1", "WH…
## $ fcode            <chr> "AR08750000", "AR08750000", "AR08750000", "AR08750000…
## $ bcmj_tag         <int> 1, 2, 4, 7, 11, 12, 15, 16, 21, 22, 23, 25, 28, 32, 3…
## $ name             <chr> "Ucluelet", "Tofino", "Lillooet", "Invermere", "Sparw…
## $ city_type        <chr> "DM", "DM", "DM", "DM", "DM", "DM", "DM", "DM", "DM",…
## $ long_type        <chr> "DISTRICT MUNICIPALITY", "DISTRICT MUNICIPALITY", "DI…
## $ pop_2000         <int> 1824, 1540, 2987, 2952, 4167, 4206, 700, 6275, 11533,…
## $ pop_source       <chr> "BCSTATS", "BCSTATS", "BCSTATS", "BCSTATS", "BCSTATS"…
## $ objectid         <int> 1, 2, 4, 7, 11, 12, 15, 16, 21, 22, 23, 25, 28, 32, 3…
## $ se_anno_cad_data <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ GEOMETRY         <POINT [m]> POINT (1036351 441199.7), POINT (1008983 455516…

5.2 Benchmark

microbenchmark::microbenchmark(
  st_read("bc_cities.shp", quiet = T) %>% filter(NAME == "Prince George"), 
  st_read("bc_cities.gpkg", quiet = T, 
          query = "SELECT * FROM bc_cities WHERE NAME == 'Prince George'"), 
  st_read("bc_cities.sqlite", quiet = T, 
          query = "SELECT * FROM bc_cities WHERE NAME == 'Prince George'"))
## Unit: milliseconds
##                                                                                                     expr
##                                  st_read("bc_cities.shp", quiet = T) %>% filter(NAME == "Prince George")
##    st_read("bc_cities.gpkg", quiet = T, query = "SELECT * FROM bc_cities WHERE NAME == 'Prince George'")
##  st_read("bc_cities.sqlite", quiet = T, query = "SELECT * FROM bc_cities WHERE NAME == 'Prince George'")
##      min       lq      mean   median       uq     max neval cld
##  12.2382 12.50815 12.907519 12.66560 13.36905 14.0081   100 a  
##   8.9347  9.34875  9.663555  9.55475  9.87500 15.6543   100  b 
##   6.7460  7.00035  7.280490  7.29095  7.50400  7.9120   100   c

6 Lines

railways <- bcdc_get_data("railway-track-line",
                          resource = "bf30d34e-1f6b-4034-a35c-1cf7c9707ae7")
## This object has 10055 records and requires 2 paginated requests to complete.
## Retrieving data
## Parsing data
railways
## Simple feature collection with 10055 features and 47 fields
## Geometry type: LINESTRING
## Dimension:     XY
## Bounding box:  xmin: 484759.6 ymin: 382665.3 xmax: 1816901 ymax: 1699824
## Projected CRS: NAD83 / BC Albers
## # A tibble: 10,055 × 48
##    id    RAILWAY_TRACK_ID NID   TRACK_SEGMENT_ID TRACK_NAME TRACK_CLASSIFICATION
##  * <chr>            <int> <chr> <chr>            <chr>      <chr>               
##  1 WHSE…                1 46a2… c96278ed49bd49b… None       Yard                
##  2 WHSE…                2 b94a… a4fbee15d38c402… None       Siding              
##  3 WHSE…                3 d9ac… 1efb8684e99c4a0… None       Main                
##  4 WHSE…                4 dbb4… 9d3e669022b343e… Kamloops   Yard                
##  5 WHSE…                5 2c48… ec79608c92434b6… None       Spur                
##  6 WHSE…                6 8557… 7558c69d02ab487… None       Main                
##  7 WHSE…                7 6a73… 73333f4fdb9140b… Canada Li… Main                
##  8 WHSE…                8 73ce… 34179415bc56498… None       Siding              
##  9 WHSE…                9 6723… c38ef77344a14aa… Prince Ge… Yard                
## 10 WHSE…               10 4471… 73803bdaca7a48c… Prince Ge… Yard                
## # ℹ 10,045 more rows
## # ℹ 42 more variables: REGULATOR <chr>, TRANSPORT_TYPE <chr>, USE_TYPE <chr>,
## #   GAUGE <chr>, NUMBER_OF_TRACKS <int>, ELECTRIFICATION <chr>, STATUS <chr>,
## #   DESIGN_SPEED_FREIGHT <int>, DESIGN_SPEED_PASSENGER <int>, SOURCE_ID <chr>,
## #   OPERATOR_ENGLISH_NAME <chr>, OPERATOR_SUBDIV_PORTION_START <dbl>,
## #   OPERATOR_SUBDIV_PORTION_END <dbl>, OWNER_NAME <chr>,
## #   TRACK_USER1_ENGLISH_NAME <chr>, TRACK_USER2_ENGLISH_NAME <chr>, …
railways %>%
  ggplot() + 
  geom_sf(aes(colour = USE_TYPE))

railways <- railways %>% 
  select(TRACK_NAME, TRACK_CLASSIFICATION, USE_TYPE)
railways
## Simple feature collection with 10055 features and 3 fields
## Geometry type: LINESTRING
## Dimension:     XY
## Bounding box:  xmin: 484759.6 ymin: 382665.3 xmax: 1816901 ymax: 1699824
## Projected CRS: NAD83 / BC Albers
## # A tibble: 10,055 × 4
##    TRACK_NAME          TRACK_CLASSIFICATION USE_TYPE                    geometry
##    <chr>               <chr>                <chr>               <LINESTRING [m]>
##  1 None                Yard                 Freight    (1208684 482248.3, 12087…
##  2 None                Siding               Freight a… (1326619 542880.3, 13266…
##  3 None                Main                 Freight a… (1417685 644318.3, 14174…
##  4 Kamloops            Yard                 Freight    (1398520 650507.3, 13985…
##  5 None                Spur                 Freight    (1151069 462783.5, 11510…
##  6 None                Main                 Freight a… (1504164 902711.4, 15040…
##  7 Canada Line         Main                 Passenger  (1210137 469036, 1210165…
##  8 None                Siding               Freight    (1240785 914620.1, 12407…
##  9 Prince George South Yard                 Freight    (1213653 990118.6, 12136…
## 10 Prince George South Yard                 Freight    (1213251 991000.4, 12132…
## # ℹ 10,045 more rows
railways %>% st_length() %>% head()
## Units: [m]
## [1]    37.01437   275.19940  7390.58200    86.90221    55.29159 12178.29057
railways <- railways %>% mutate(track_length = st_length(.))
long_railways <- railways %>% filter(as.numeric(track_length) > 5000)
long_railways
## Simple feature collection with 442 features and 4 fields
## Geometry type: LINESTRING
## Dimension:     XY
## Bounding box:  xmin: 484759.6 ymin: 382949.2 xmax: 1816692 ymax: 1699824
## Projected CRS: NAD83 / BC Albers
## # A tibble: 442 × 5
##    TRACK_NAME TRACK_CLASSIFICATION USE_TYPE                             geometry
##  * <chr>      <chr>                <chr>                        <LINESTRING [m]>
##  1 None       Main                 Freight and Tourist (1417685 644318.3, 14174…
##  2 None       Main                 Freight and Passen… (1504164 902711.4, 15040…
##  3 None       Main                 Freight and Tourist (1662986 756995, 1662937…
##  4 None       Main                 Freight and Tourist (1566502 712251.9, 15665…
##  5 None       Main                 Freight and Passen… (1403899 673631.1, 14038…
##  6 None       Main                 Freight and Passen… (1146903 466782.8, 11469…
##  7 None       Main                 Freight             (1085898 1073262, 108576…
##  8 None       Main                 Freight and Tourist (1578861 723160.5, 15788…
##  9 None       Main                 Freight and Passen… (730156.7 1025246, 73018…
## 10 None       Main                 Freight and Passen… (1356126 960868.2, 13560…
## # ℹ 432 more rows
## # ℹ 1 more variable: track_length [m]
library(units)
long_railways <- filter(railways, track_length > as_units(5000, "m"))
ggplot() +
  geom_sf(data = long_railways, aes(colour = as.numeric(track_length)))

7 Unions (via summarize)

nr_district <- bcdc_get_data('natural-resource-nr-district')
## You are accessing a simplified view of the data - see the catalogue record for details.
nr_region <- nr_district %>% 
  group_by(REGION_ORG_UNIT_NAME) %>% 
  summarise() # << defaults to union

ggplot() +
  geom_sf(data = nr_region, colour = "white", 
          aes(fill = REGION_ORG_UNIT_NAME)) 

nr_region_multi <- nr_district %>% 
  group_by(REGION_ORG_UNIT_NAME) %>% 
  summarise(do_union = FALSE) # <<--

nr_region_multi # MULTIPOLYGON! 
## Simple feature collection with 8 features and 1 field
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 273366.8 ymin: 359771.8 xmax: 1870587 ymax: 1735721
## Projected CRS: NAD83 / BC Albers
## # A tibble: 8 × 2
##   REGION_ORG_UNIT_NAME                                                  geometry
##   <chr>                                                       <MULTIPOLYGON [m]>
## 1 Cariboo Natural Resource Region           (((1027839 883453.2, 1023202 878748…
## 2 Kootenay-Boundary Natural Resource Region (((1490768 821176.4, 1490447 820935…
## 3 Northeast Natural Resource Region         (((1235405 1136771, 1235426 1136708…
## 4 Omineca Natural Resource Region           (((863577.5 1272237, 863508.1 12721…
## 5 Skeena Natural Resource Region            (((1064630 1052166, 1064398 1052332…
## 6 South Coast Natural Resource Region       (((1091120 523587.2, 1094036 519727…
## 7 Thompson-Okanagan Natural Resource Region (((1162743 654314.6, 1162555 653829…
## 8 West Coast Natural Resource Region        (((862927.5 562402.1, 866696 565781…
ggplot() +
  geom_sf(data = nr_region_multi, colour = "white", 
          aes(fill = REGION_ORG_UNIT_NAME)) 

nr_region %>% 
  mutate(region_area = st_area(geometry))
## Simple feature collection with 8 features and 2 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 273366.8 ymin: 359771.8 xmax: 1870587 ymax: 1735721
## Projected CRS: NAD83 / BC Albers
## # A tibble: 8 × 3
##   REGION_ORG_UNIT_NAME                                      geometry region_area
## * <chr>                                                <POLYGON [m]>       [m^2]
## 1 Cariboo Natural Resource Region          ((1381787 807659.5, 1381…     8.25e10
## 2 Kootenay-Boundary Natural Resource Regi… ((1708109 725049.3, 1708…     8.23e10
## 3 Northeast Natural Resource Region        ((1390137 1036196, 13900…     1.75e11
## 4 Omineca Natural Resource Region          ((1235426 1136708, 12353…     1.58e11
## 5 Skeena Natural Resource Region           ((973837.6 1168731, 9739…     2.63e11
## 6 South Coast Natural Resource Region      ((1270064 574122.1, 1270…     4.62e10
## 7 Thompson-Okanagan Natural Resource Regi… ((1398207 458006.1, 1389…     7.51e10
## 8 West Coast Natural Resource Region       ((1095742 517629, 109657…     1.57e11
library(bcdata)

nr_district <- bcdc_get_data('natural-resource-nr-district')
## You are accessing a simplified view of the data - see the catalogue record for details.
lines <- bcdc_get_data('bc-transmission-lines')

fires_2017 <- bcdc_query_geodata('bc-wildfire-fire-perimeters-historical') %>%
  filter(FIRE_YEAR == 2017) %>% 
  collect()

big_fires <- fires_2017 %>%
  filter(FIRE_NUMBER %in% c('C10784', 'C50647'))

bc <- bcdc_query_geodata('7-5m-provinces-and-states-the-atlas-of-canada-base-maps-for-bc') %>% 
  filter(ENGLISH_NAME == 'British Columbia') %>% 
  collect()

m <- ggplot() +
  geom_sf(data = bc, fill = "grey80") +
  geom_sf(data = nr_district, fill = "purple", alpha = 0.5) +
  geom_sf(data = big_fires, fill = "orange", alpha = 0.5) +
  geom_sf(data = lines, colour = "yellow")

m

8 Summary and Conclusion

  • Learned data manipulation and geospatial visualization.
  • Utilized libraries like sf, ggplot2, and mapview.
  • Explored bcdata and performed advanced spatial operations.