Choropleths are a staple of geospatial visualization. The idea of this example is to introduced those basics in stroke context.

The example loads the data, estimates annual stroke incidence based on data from the NEMISIS study, performs a join, calculates some geospatial measures, filters postcodes based on distance to the hospital and displays the results.

0. Package loading

library(tidyverse)
library(sf)
library(units)
library(tmaptools)

1. Loading census and boundary data

postcodeboundariesAUS <- sf::read_sf(
  here::here("ABSData", 
            "Boundaries", 
            "POA_2016_AUST.shp"))

basicDemographicsVIC <- readr::read_csv(
  here::here("ABSData", 
            "2016 Census GCP Postal Areas for VIC", 
            "2016Census_G01_VIC_POA.csv"))

2. Combine demographics and spatial data

## Join the demographics and shape tables, retaining victoria only
## use postcode boundaries as the reference data frame so that coordinate
## reference system is retained.
basicDemographicsVIC <- right_join(postcodeboundariesAUS, basicDemographicsVIC, 
                                  by=c("POA_CODE" = "POA_CODE_2016"))

3. Geocode hospital location

## Location of hopsital providing acute stroke services
## address: 246 Clayton Rd, Clayton VIC, 3168
MMCLocation <- tmaptools::geocode_OSM("Monash Medical Centre, Clayton, Victoria, Australia", as.sf=TRUE)
MMCLocation
#> Simple feature collection with 1 feature and 7 fields
#> geometry type:  POINT
#> dimension:      XY
#> bbox:           xmin: 145.1207 ymin: -37.92093 xmax: 145.1207 ymax: -37.92093
#> epsg (SRID):    4326
#> proj4string:    +proj=longlat +datum=WGS84 +no_defs
#>                                                 query       lat      lon
#> 1 Monash Medical Centre, Clayton, Victoria, Australia -37.92093 145.1207
#>     lat_min   lat_max  lon_min  lon_max                   geometry
#> 1 -37.92098 -37.92088 145.1207 145.1208 POINT (145.1207 -37.92093)

4. Compute per-postcode stroke incidence

basicDemographicsVIC <- mutate(basicDemographicsVIC, 
                               Age_0_24_yr_P = Age_0_4_yr_P + Age_5_14_yr_P + 
                                 Age_15_19_yr_P + Age_20_24_yr_P)
basicDemographicsVIC <- mutate(basicDemographicsVIC, stroke_count_estimate = ( 
  Age_0_24_yr_P  * 5 +  
    Age_25_34_yr_P * 30 +   
    Age_35_44_yr_P * 44 +  
    Age_45_54_yr_P * 111 +  
    Age_55_64_yr_P * 299 +  
    Age_65_74_yr_P * 747 + 
    Age_75_84_yr_P * 1928 +  
    Age_85ov_P     * 3976) / 100000)

5. Compute distance to hospital from postcode

## Add some geospatial measures to the data frame
## Reproduce the existing area column as a demo.
basicDemographicsVIC <- mutate(basicDemographicsVIC, 
                               PostcodeArea=units::set_units(st_area(geometry), km^2))

## Distance to MMC
basicDemographicsVIC <- sf::st_transform( basicDemographicsVIC, crs = sf::st_crs( MMCLocation ) )
basicDemographicsVIC <- mutate(basicDemographicsVIC, 
                               DistanceToMMC=units::set_units(st_distance(geometry,MMCLocation)[,1], km))

## comments re auto-great-circle, straight line assumption, 
## distance to polygon or polygon-centroid ...
plot(basicDemographicsVIC["DistanceToMMC"])

6. Discard remote postcodes

## Make a small dataset for MMC surrounds.
basicDemographicsMMC <- filter(basicDemographicsVIC, DistanceToMMC < set_units(20, km))

7. Interactive display of the result

library(tmap)
tmap_mode("view")
#> tmap mode set to interactive viewing

MMCLocation <- mutate(MMCLocation, ID="Monash Medical Centre")
basicDemographicsMMC <- mutate(basicDemographicsMMC, Over65 = Age_65_74_yr_P + Age_75_84_yr_P + Age_85ov_P)
tm_shape(basicDemographicsMMC, name="Annual stroke counts") + 
  tm_polygons("stroke_count_estimate", id="POA_NAME", popup.vars=c("Cases"="stroke_count_estimate"), alpha=0.6) + 
  tm_shape(MMCLocation) + tm_markers() + 
  tm_basemap("OpenStreetMap")