Intro

Welcome! Here, I’ll demonstrate a quick way to link the British Columbia (BC) Community Health Service Area (CHSA) shapefile, with any given point data (E.g. imagine you have health events by longitude / latidude). We’ll do this through a spatial join!

I have this published on RPubs as well, you can check that out here: https://rpubs.com/jdsimkin04

For RPubs sake, here’s the github repo: https://github.com/jdsimkin04/bc_chsa_spatial_join

What is the BC CHSA shapefile?

The CHSA is the smallest health administrative area in BC. It was created to facilitate community level health analyses. The shapefile is openly available through the BC Gov data catalogue. It also contains attribute data - contextual information related to that CHSA. These include but are not limited too:

  • Corresponding health administrative areas (LHA, HSDA, HA)
  • Census population
  • Area km2
  • Level of urbanization

Spatial join rural-urban class to point data

Download the CHSA map file

We can do this a few ways… you could go to the BC Gov data catalogue and download the shapefile. Or you could use the download.file and unzip functions shown below.

Libraries

install these libraries if you don’t have them

library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──
## ✔ ggplot2 3.3.6     ✔ purrr   0.3.4
## ✔ tibble  3.1.7     ✔ dplyr   1.0.9
## ✔ tidyr   1.2.0     ✔ stringr 1.4.0
## ✔ readr   2.1.2     ✔ forcats 0.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
library(sf)
## Linking to GEOS 3.9.1, GDAL 3.4.2, PROJ 8.2.1; sf_use_s2() is TRUE
library(tmap)
library(tmaptools)
library(kableExtra)
## 
## Attaching package: 'kableExtra'
## The following object is masked from 'package:dplyr':
## 
##     group_rows
library(rmapshaper)
## Registered S3 method overwritten by 'geojsonlint':
##   method         from 
##   print.location dplyr

Download and unzip

Download the file from the website using the download.file and unzip function to unzip into a folder called map

url <- "https://catalogue.data.gov.bc.ca/dataset/68f2f577-28a7-46b4-bca9-7e9770f2f357/resource/f89f99b0-ca68-41e2-afc4-63fdc0edb666/download/chsa_2018.zip"
download.file(url, "chsa_2018.zip")
unzip("chsa_2018.zip", exdir = paste0(getwd(), "/map"))

Visualize

read the map files using the st_read function

chsa <-
st_read("map") %>% 
  st_make_valid() %>% 
  st_transform(., crs = 4326) %>% 
  simplify_shape(.)

visualize the chsa map

tmap_mode("view")
## tmap mode set to interactive viewing
tm_shape(chsa) +
  tm_polygons("HSDA_Name",
              id = "CHSA_Name",
              palette = "PRGn",  
              legend.show = FALSE)

see the data elements of the chsa map, the CHSA_UR_Cl is the rural urban classification

chsa %>% 
  st_drop_geometry() %>% 
  str()
## 'data.frame':    218 obs. of  20 variables:
##  $ CHSA_CD   : chr  "1110" "1120" "1130" "1140" ...
##  $ CHSA_Name : chr  "Fernie" "Cranbrook" "Kimberley" "Windermere" ...
##  $ CHSA_Title: chr  "1110 Fernie" "1120 Cranbrook" "1130 Kimberley" "1140 Windermere" ...
##  $ LHA_CD    : chr  "111" "112" "113" "114" ...
##  $ LHA_Name  : chr  "Fernie" "Cranbrook" "Kimberley" "Windermere" ...
##  $ LHA_Title : chr  "111 Fernie" "112 Cranbrook" "113 Kimberley" "114 Windermere" ...
##  $ LHA_CD1997: chr  "001" "002" "003" "004" ...
##  $ HSDA_CD   : chr  "11" "11" "11" "11" ...
##  $ HSDA_ID   : chr  "11 EK" "11 EK" "11 EK" "11 EK" ...
##  $ HSDA_Name : chr  "East Kootenay" "East Kootenay" "East Kootenay" "East Kootenay" ...
##  $ HSDA_Title: chr  "11 East Kootenay" "11 East Kootenay" "11 East Kootenay" "11 East Kootenay" ...
##  $ HA_CD     : chr  "1" "1" "1" "1" ...
##  $ HA_ID     : chr  "1 IHA" "1 IHA" "1 IHA" "1 IHA" ...
##  $ HA_Name   : chr  "Interior" "Interior" "Interior" "Interior" ...
##  $ HA_Title  : chr  "1 Interior" "1 Interior" "1 Interior" "1 Interior" ...
##  $ CHSA_UR_Cl: chr  "5 Rural Hub" "4 Small Urban" "6 Rural" "6 Rural" ...
##  $ CHSA_Pop16: int  15531 26248 9178 9487 12634 6851 3209 25523 13710 4501 ...
##  $ CHSA_Area : num  8044 4474 4346 10981 3794 ...
##  $ Latitude  : num  49.4 49.6 49.8 50.5 49.2 ...
##  $ Longitude : num  -115 -116 -116 -116 -117 ...

Explore CHSA_UR_Cl

Let’s look a bit closer at the CHSA_UR_Cl variable

chsa %>% 
  st_drop_geometry() %>% 
  count(CHSA_UR_Cl) %>% 
  mutate(proportion = round(n/sum(n)*100,1))
##       CHSA_UR_Cl  n proportion
## 1 1 Metropolitan 60       27.5
## 2  2 Large Urban 23       10.6
## 3 3 Medium Urban 19        8.7
## 4  4 Small Urban 18        8.3
## 5    5 Rural Hub 25       11.5
## 6        6 Rural 57       26.1
## 7       7 Remote 16        7.3
tmap_mode("view")
## tmap mode set to interactive viewing
tm_shape(chsa) +
  tm_polygons("CHSA_UR_Cl",
              id = "CHSA_Name",
              palette = "PRGn")

Sample point data

Here’s a sample of point data that we want to assign CHSA rural urban class. I got these randomly from google maps.

spatial_point_data <-
tibble(
  case_id = c(1:10),
  long = c(-122.79245404269491,
                -122.50750051771003,
                -121.47061705091572,
                -120.78489105416672,
                -119.46708958311152,
                -122.12044919737846,
                -123.9735411710232,
                -127.48038830551505,
                -122.66126844286325,
                -131.11218885196038),
  lat = c(
    49.1152475836367,  
49.117737238190344,
49.38845365756505, 
50.09960249743014, 
49.87203980873458, 
52.13143107578074, 
49.130988612031395,
50.69492877221613, 
58.80790367189857, 
57.88986797544227)
) %>% 
  st_as_sf(., coords = c("long", "lat"), crs = 4326, agr = "constant") %>% 
    st_transform(., crs = 4326)

Visualize

Viz the point data

tmap_mode("view")
## tmap mode set to interactive viewing
tm_shape(spatial_point_data) +
  tm_dots(size = 0.1,
          col = "black")

Viz teh point data on the CHSAs

tmap_mode("view")
## tmap mode set to interactive viewing
tm_shape(chsa) +
  tm_polygons("HSDA_Name",
              id = "CHSA_Name",
              palette = "PRGn",  
              legend.show = FALSE) +
    tm_shape(spatial_point_data) +
  tm_dots(size = 0.2,
          col = "black")

Spatial join rural urban class to the point data

link cohort sf and chsa map

int <-
  st_intersects(spatial_point_data, chsa, prepared = TRUE) #Wherever there are intersection between the two layers (the CHSA Map and the study data), we can link them! and store them in "int"

linked_data <-
spatial_point_data %>%
  mutate(
    CHSA_Name = as.character(chsa$CHSA_Name[as.numeric(as.character(int))]),
    CHSA_UR_Cl = as.character(chsa$CHSA_UR_Cl[as.numeric(as.character(int))]),
    LHA_Name = as.character(chsa$LHA_Name[as.numeric(as.character(int))]),
    HSDA_Name = as.character(chsa$HSDA_Name[as.numeric(as.character(int))]),
    HA_Name = as.character(chsa$HA_Name[as.numeric(as.character(int))])
    ) %>%
  st_drop_geometry() #drop the geometry info from the file so we can analyze the data like a norma lfile

Final linked dataset

… and Voila! See the spatially joined point data with CHSA attributes.

linked_data %>% 
  kable() %>% 
  kable_styling()
case_id CHSA_Name CHSA_UR_Cl LHA_Name HSDA_Name HA_Name
1 Panorama 1 Metropolitan Surrey Fraser South Fraser
2 North Langley Township 6 Rural Langley Fraser South Fraser
3 Hope 5 Rural Hub Hope Fraser East Fraser
4 Merritt 6 Rural Merritt Thompson Cariboo Shuswap Interior
5 Downtown Kelowna 2 Large Urban Central Okanagan Okanagan Interior
6 Williams Lake/East Cariboo 6 Rural Cariboo/Chilcotin Thompson Cariboo Shuswap Interior
7 Nanaimo West/Rural 6 Rural Greater Nanaimo Central Vancouver Island Vancouver Island
8 Port Hardy/Port Alice 5 Rural Hub Vancouver Island North North Vancouver Island Vancouver Island
9 Fort Nelson Population Centre 5 Rural Hub Fort Nelson Northeast Northern
10 Telegraph Creek 7 Remote Telegraph Creek Northwest Northern