Retail Location Analysis — Canada Edition

Integrating Statistics Canada and CensusMapper APIs for Demographic and Market Analysis

Author

Prepared for Canadian Analytics Team

Published

October 24, 2025

Flag of Canada

Introduction

In the US based model, the demographic attributes are based on radii around geocoded points. I didn’t see this functionality on the Canadian sites, albeit I know some of the commercial providers (STDB and ESRI) allow both radius and drive time searches (drive time is the gold standard). But I was experimenting with free resources.

Radius / buffer queries vs fixed geography

“When I pull data for the Canadian resources, will it be able to pull the data like my US-based code does on a radius around a geocoded point?”

Here’s what I found and how it compares:

StatCan’s API / Census Profile WDS is designed to retrieve data for predefined geographic units (e.g., provinces, census divisions, census subdivisions, census tracts, dissemination areas) as defined in the Canadian census geography structure. Statistics Canada +1

It does not natively support retrieving data for an arbitrary radius (e.g., 1-mile or 5-mile buffer) around a geocoded point in a single API call in the same way the U.S. system might handle via e.g., ACS with custom geographies or spatial queries.

What this means for the workflow (which uses radius buffers around geocoded sites) is:

You would likely need to geocode your point, then select the smallest geography units available (e.g., DAs or tracts) that cover your buffer, and then retrieve the variables for those units.

Optionally, you could aggregate or weight those units (by area overlap, population proportion, etc.) to approximate the buffer value. This is a common technique when the API does not directly support arbitrary buffers.

Because the US-based code already works with buffers and “Euclidean” or drive-time representations in your U.S. workflows, you may need to build logic to (a) fetch shapefiles (boundaries) from StatCan for the geographies, (b) perform the spatial intersection/aggregation in R/Quarto, and (c) then join the retrieved demographic variables.

Additionally, StatCan provides boundary shapefiles and spatial products: e.g., “Spatial information products” for the 2021 census geographies (boundary files, road network, etc.). Statistics Canada +1

Because the pipeline now focuses on 1- and 5-mile Euclidean buffers (for U.S. census data) and you probably want the same for Canada, you will need to treat the Canadian side as fixed-geography data + buffer-aggregation procedure, rather than an API built for radius queries. That is of course unless you can find a service that allows radii.

This guide explains how to replicate U.S. Census style demographic retrieval for Canadian markets using Statistics Canada and CensusMapper. It is written for analysts adapting buffer based demographic and market analytics to Canadian geography. Canada’s Census of Population is collected every five years. The most recent is 2021, so most variables reference that year.

Variable Equivalency Table

U.S. Metric U.S. Census Code Canadian Equivalent StatCan Table or Vector Notes
Total Population B01003_001 Total population Census 2021 Profile, Population 2021 Available at DA, CT, CMA and other levels
Median Household Income B19013_001 Median total income of households (2020 $) 98-10-0057-01 Closest conceptual match
Median Home Value B25077_001 Median value of dwellings, owner estimated 98-10-0257-01 Proxy for owner occupied value
Median Gross Rent B25046_001 Median monthly rent of tenants 98-10-0256-01 Parallel to U.S. gross rent

You can retrieve these measures with the Statistics Canada Web Data Service or with the CensusMapper API through the cancensus R package.

Accessing Canadian Data

1. Statistics Canada Web Data Service (WDS)

The Web Data Service is a REST API for StatCan tables.

Example:

https://api.statcan.gc.ca/census-recensement/profile/sdmx/rest/data/STC_CP,DF_PR/A5.2021A000235+2021A000224.1..1

Refer to the Statistics Canada Developer Guide for usage syntax and table IDs.

2. CensusMapper API via cancensus

The CensusMapper API provides simplified access to 2021 Census data at multiple geographic levels and integrates seamlessly with R.

install.packages("cancensus")
library(cancensus)
options(cancensus.api_key = Sys.getenv("CANCENSUS_API_KEY"))

# Example query: change regions and vectors for your study area
get_census(dataset = "CA21",
           regions  = list(CSD = "3520005"),
           vectors  = c("v_CA21_1",      # population
                        "v_CA21_2397",  # median household income
                        "v_CA21_4836",  # median dwelling value
                        "v_CA21_4882"), # median rent
           level    = "DA",
           geo_format = "sf")

Geographic Considerations

The APIs do not accept radius queries around a coordinate. Work with fixed census geographies instead, most often Census Tracts for cities and Dissemination Areas for finer detail. To analyze a radius, you will buffer your site and intersect that buffer with the census units, then compute weighted results.

Buffer and Overlay Functions

library(sf)
library(cancensus)

# Build a buffer and return intersecting census units with geometry
get_canadian_demographics <- function(lat, lon, radius_km = 1,
                                      vectors = c("v_CA21_1","v_CA21_2397","v_CA21_4836","v_CA21_4882"),
                                      region = list(CMA = "59933")) {
  site <- st_sfc(st_point(c(lon, lat)), crs = 4326)
  buf  <- st_transform(st_buffer(st_transform(site, 3347), radius_km * 1000), 4326)

  data <- get_census(dataset = "CA21",
                     regions = region,
                     vectors = vectors,
                     level = "DA",
                     geo_format = "sf")

  st_intersection(data, buf)
}

Aggregating to Buffer Level Metrics

Use area or population weights to summarize values inside the buffer. The example below uses population as a weight since it is a strong proxy for demand.

aggregate_to_radius <- function(df, pop_vector = "v_CA21_1") {
  stopifnot(pop_vector %in% names(df))
  df |>
    mutate(weight = !!rlang::sym(pop_vector) / sum(!!rlang::sym(pop_vector), na.rm = TRUE)) |>
    summarise(across(where(is.numeric), ~ weighted.mean(.x, weight, na.rm = TRUE)))
}

Illustrative output layout:

Metric Value
Total population 7,421
Median household income (2020 $) $86,700
Median dwelling value $648,000
Median rent $1,320

Leaflet Visualization

Interactive maps help verify overlay logic and communicate trade areas.

library(leaflet)

visualize_buffer <- function(lat, lon, radius_km = 1, region = list(CMA = "59933")) {
  site <- st_sfc(st_point(c(lon, lat)), crs = 4326)
  buf  <- st_transform(st_buffer(st_transform(site, 3347), radius_km * 1000), 4326)

  data <- get_census(dataset = "CA21",
                     regions = region,
                     vectors = c("v_CA21_1"),
                     level = "DA",
                     geo_format = "sf")

  leaflet() %>%
    addTiles() %>%
    addPolygons(data = data, color = "#999999", weight = 1, fillOpacity = 0.4, group = "DAs") %>%
    addPolygons(data = buf,  color = "#E31837", weight = 3, fillOpacity = 0.2, group = "Buffer") %>%
    addMarkers(lng = lon, lat = lat, popup = "Target Site") %>%
    addLegend(position = "bottomright",
              colors = c("#E31837", "#999999"),
              labels = c("Analysis buffer", "Dissemination areas"),
              opacity = 0.7,
              title = "Map legend") %>%
    addLayersControl(overlayGroups = c("DAs", "Buffer"),
                     options = layersControlOptions(collapsed = FALSE))
}

Comparing 1 km and 5 km Buffers

compare_buffers <- function(lat, lon, region = list(CMA = "59933")) {
  site <- st_sfc(st_point(c(lon, lat)), crs = 4326)
  buf1 <- st_transform(st_buffer(st_transform(site, 3347), 1000), 4326)
  buf5 <- st_transform(st_buffer(st_transform(site, 3347), 5000), 4326)

  data <- get_census(dataset = "CA21",
                     regions = region,
                     vectors = c("v_CA21_1"),
                     level = "DA",
                     geo_format = "sf")

  leaflet() %>%
    addTiles() %>%
    addPolygons(data = data, color = "#888888", weight = 1, fillOpacity = 0.3, group = "DAs") %>%
    addPolygons(data = buf1, color = "#E31837", weight = 2, fillOpacity = 0.25, group = "1 km buffer") %>%
    addPolygons(data = buf5, color = "#003DA5", weight = 2, fillOpacity = 0.15, group = "5 km buffer") %>%
    addMarkers(lng = lon, lat = lat, popup = "Target Site") %>%
    addLegend(position = "bottomright",
              colors = c("#E31837", "#003DA5", "#888888"),
              labels = c("1 km buffer", "5 km buffer", "Dissemination areas"),
              opacity = 0.8,
              title = "Map legend: buffers") %>%
    addLayersControl(overlayGroups = c("DAs", "1 km buffer", "5 km buffer"),
                     options = layersControlOptions(collapsed = FALSE))
}

Color and Symbol Notes

  • Red (#E31837) is the 1 km buffer. It represents the immediate trade area.
  • Blue (#003DA5) is the 5 km buffer. It represents the extended trade area.
  • Gray (#888888) outlines dissemination areas. This is the smallest census geography.

Exporting Maps as Images

Create static images for reports or slide decks with mapview::mapshot() and webshot2.

install.packages(c("mapview", "webshot2"))
library(mapview)
library(webshot2)

m <- compare_buffers(lat = 43.6532, lon = -79.3832)
mapshot(m, file = "Toronto_Buffer_Comparison.png", vwidth = 1600, vheight = 1000, zoom = 1)

Appendix: Formatting Exported Maps

In Quarto

![](Toronto_Buffer_Comparison.png){fig-cap="Toronto 1 km and 5 km trade area buffers" width=85% fig-align="center"}

In PowerPoint

  • Export at 16:9, for example 1600 by 900 pixels.
  • Add a thin border if your slide background is dark.
  • Use concise legend labels, for example 1 km Trade Area and 5 km Trade Area.
  • Add a small credit below the image: Statistics Canada (2021) and CensusMapper.

Troubleshooting and Common Issues

API key or authorization

Error: Error: Unauthorized (HTTP 401)

Fix:

Sys.setenv(CANCENSUS_API_KEY = "your_api_key_here")
options(cancensus.api_key = Sys.getenv("CANCENSUS_API_KEY"))

CRS mismatch

Error: Error in st_intersection(): arguments have different crs

Confirm and align CRS before intersections. Buffer in a projected CRS and convert back for mapping.

st_crs(data)
st_crs(buffer)
# Align and buffer in EPSG 3347, then transform back
buffer <- st_transform(st_buffer(st_transform(site, 3347), 1000), 4326)

Empty or partial returns

Message: no data for the selected region

Check that the region code exists in CA21. Use helper functions to explore regions:

list_census_regions("CA21")

Leaflet rendering

If basemap tiles fail, try a stable provider.

leaflet() %>% addProviderTiles(providers$CartoDB.Positron)

webshot2 timeouts

If exports are blank or time out, ensure a modern browser is available. Increase delay if needed.

mapshot(m, file = "map.png", delay = 5)

Implementation Checklist

  1. Obtain and store a CensusMapper API key in your environment.
  2. Choose working geography, most often DA for fine detail or CT for city scale.
  3. Download boundary files if you plan to cache or pre process geographies.
  4. Build buffers with a projected CRS and intersect with census units.
  5. Aggregate results with population based weights for demand oriented metrics.
  6. Produce maps and export static images for deliverables.
  7. Cite Statistics Canada and CensusMapper in outputs.

References

  • Statistics Canada. 2021 Census of Population Profile.
  • CensusMapper API documentation.
  • cancensus R package documentation.
  • mapview and webshot2 documentation.

This guide provides a complete Canadian workflow that mirrors your U.S. process, including buffer logic, aggregation, interactive maps, static exports, and practical troubleshooting. It ships with a banner for visual identity and a PDF footer for attribution when you render a handout.