Interactive maps in R with mapgl

Our real Sydney case study

We will use Greater Sydney, as defined by the ABS Greater Capital City Statistical Area (1GSYD). This includes the Central Coast and Blue Mountains; it is broader than the City of Sydney council area. The polygon maps use 373 real Statistical Areas Level 2 (SA2s). An SA2 is a statistical area, not necessarily a suburb or a school catchment.

The school snapshot includes environmental education centres and other government education sites as well as primary and secondary schools. We retain the provider’s types. It covers government education, rather than all schools. The enrolment field is latest reported FTE, not a count measured on the extraction date; the source does not supply an enrolment-year column. Population estimates for 2025 are preliminary. These periods and populations cannot be treated as interchangeable.

The prepared CSV, GeoJSON, GeoPackage, and GeoTIFF files are bundled with the source. The data guide records attribution, original URLs, transformations, exclusions, and checksums. None of the maps uses simulated observations.

You can work through the core sections first, then choose an extension that suits your data. An effective project might need only one well-designed map and a supporting chart. This is an optional teaching resource: follow the assessment instructions for the required data, methods, and submission format.

1. Set up R and Quarto

Use R 4.1 or later, RStudio, and the Quarto CLI.

install.packages(c(
  "mapgl", "sf", "dplyr", "readr", "htmltools", "htmlwidgets",
  "terra", "viridisLite", "knitr", "rmarkdown", "jsonlite"
))

Now load the packages in the report. Unlike installation, this step must run whenever you render it.

library(mapgl)
library(sf)

Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.4.0; sf_use_s2() is TRUE

library(dplyr)

Attaching package: 'dplyr'

The following objects are masked from 'package:stats':

    filter, lag

The following objects are masked from 'package:base':

    intersect, setdiff, setequal, union

library(readr)

# This tutorial uses the mapgl 0.5.2 API.
stopifnot(packageVersion("mapgl") >= "0.5.2")

mapgl creates the maps; sf stores vector geometries and performs spatial operations; dplyr prepares tables; readr imports CSV files. Later, terra handles raster data. Package names and R functions are case-sensitive.

The mapgl documentation describes two interfaces: maplibre() and mapboxgl(). We will use maplibre() throughout. Mapbox requires an access token and its own service terms; it is an optional alternative, rather than a requirement for this lesson.

2. Make your first map

# Start with one actual record from the bundled NSW snapshot.
first_school <- read_csv("data/sydney-public-schools.csv",
                         show_col_types = FALSE) |>
  filter(name == "Glebe Public School")
stopifnot(nrow(first_school) == 1)
first_location <- c(first_school$longitude, first_school$latitude)

first_map <- maplibre(
  style = openfreemap_style("positron"),
  center = first_location,
  zoom = 12,
  projection = "mercator",
  height = 420,
  scrollZoom = FALSE
) |>
  add_markers(
    first_location,
    color = "#0d706f",
    popup = "Glebe Public School · NSW Department of Education"
  ) |>
  add_navigation_control(position = "top-right") |>
  add_fullscreen_control(position = "top-right") |>
  add_scale_control(position = "bottom-left", unit = "metric")

first_map

The final line displays the object. If you finish a chunk with first_map <- ..., R stores the map but does not automatically display it.

A marker is convenient for one or a few locations. For a dataset of points, use a circle layer, which supports styling, filtering, and clustering. See the maplibre() reference and available OpenFreeMap styles.

3. Put geospatial data into the right form

You need either coordinates, geometries, or an identifier that can be joined to boundaries. A place name on its own is not a geometry.

Geocoding is a separate step, with possible ambiguity, service limits, and costs. “Sydney” could mean a city, an administrative area, or a suburb. Document how you resolved it. A map cannot supply geographic precision that your source data does not contain.

schools_raw <- read_csv("data/sydney-public-schools.csv",
  col_types = cols(school_code = col_character()))

schools_raw |>
  select(name, school_type, longitude, latitude, enrolment_fte) |>
  head(4) |>
  knitr::kable(caption = "A few real records from the NSW school snapshot.")

coordinate_ok <- with(
  schools_raw,
  is.finite(longitude) & is.finite(latitude) &
    longitude >= -180 & longitude <= 180 &
    latitude >= -90 & latitude <= 90
)
rejected_schools <- schools_raw[!coordinate_ok, ]
schools_clean <- schools_raw[coordinate_ok, ]
stopifnot(nrow(schools_clean) > 0, !anyDuplicated(schools_clean$school_code))

schools <- st_as_sf(
  schools_clean,
  coords = c("longitude", "latitude"),
  crs = 4326,
  remove = FALSE
)

data.frame(
  rows_imported = nrow(schools_raw),
  rows_excluded = nrow(rejected_schools),
  rows_mapped = nrow(schools)
) |>
  knitr::kable()

st_crs(schools)$input

[1] "EPSG:4326"

unique(st_geometry_type(schools))

[1] POINT
18 Levels: GEOMETRY POINT LINESTRING POLYGON MULTIPOINT ... TRIANGLE

st_bbox(schools)

     xmin      ymin      xmax      ymax 
150.23258 -34.32048 151.58244 -33.13547 

areas <- st_read("data/my-boundaries.geojson", quiet = TRUE)
stopifnot(!is.na(st_crs(areas)))
areas <- st_transform(areas, 4326)

# A GeoPackage can contain several layers.
st_layers("data/my-geography.gpkg")
roads <- st_read("data/my-geography.gpkg", layer = "roads", quiet = TRUE)

# For projected CSV coordinates, use the CRS stated by the data provider.
projected_points <- st_as_sf(my_table, coords = c("easting", "northing"),
                            crs = verified_source_crs)
web_points <- st_transform(projected_points, 4326)

4. Add points, tooltips, and popups

A map contains sources and layers. A source holds the data. A layer tells the browser how to draw those data. Here, providing an sf object directly to source lets mapgl register the source for us.

# Escape text before placing it inside popup HTML.
schools <- schools |>
  mutate(popup_html = paste0(
    "<strong>", htmltools::htmlEscape(name), "</strong><br>",
    htmltools::htmlEscape(school_type), "<br>",
    "Latest reported FTE enrolment: ",
    if_else(is.na(enrolment_fte), "Not reported", as.character(enrolment_fte))
  ))

points_map <- maplibre(
  style = openfreemap_style("positron"),
  bounds = schools,
  fitBoundsOptions = list(padding = 40),
  projection = "mercator",
  height = 440,
  scrollZoom = FALSE
) |>
  add_circle_layer(
    id = "school-points",
    source = schools,
    circle_color = "#0d706f",
    circle_radius = 7,
    circle_opacity = 0.85,
    circle_stroke_color = "#ffffff",
    circle_stroke_width = 1.5,
    tooltip = "name",
    popup = "popup_html"
  ) |>
  add_navigation_control()

points_map

bounds = schools fits the starting view to the data. id is a unique layer identifier, not a column name. tooltip and popup identify columns used for hover and click content. Popup columns can contain HTML, so escape text imported from other sources before adding your own markup.

fitBoundsOptions = list(padding = 40) leaves room around the fitted data so circles at the edges are fully visible. It is a MapLibre option passed through by maplibre().

Read the sources and layers overview and add_circle_layer() reference for more options.

maplibre_view(
  filter(schools, !is.na(enrolment_fte)),
  column = "enrolment_fte",
  palette = viridisLite::viridis,
  style = openfreemap_style("positron"),
  projection = "mercator",
  height = 420,
  scrollZoom = FALSE
) |>
  add_navigation_control()

5. Map categories and quantities

Styling expressions connect an attribute to an appearance. A constant gives every feature the same style; an expression calculates a style from a feature or map property.

school_groups <- c("Primary / infants", "Secondary", "Other government sites")
school_colours <- c("#0072B2", "#D55E00", "#009E73")
schools <- schools |>
  mutate(school_group = case_when(
    school_type %in% c("Primary School", "Infants School") ~ school_groups[1],
    school_type == "Secondary School" ~ school_groups[2],
    TRUE ~ school_groups[3]
  ))

# Missing enrolment is not zero. Keep those locations in the basic map,
# but account for their exclusion from this quantitative encoding.
schools_known <- schools |> filter(!is.na(enrolment_fte))
stopifnot(all(is.finite(schools_known$enrolment_fte)),
          all(schools_known$enrolment_fte >= 0))
schools_known <- schools_known |>
  mutate(radius_px = 6 * sqrt(enrolment_fte / 500))

data.frame(all_locations = nrow(schools),
           missing_enrolment = sum(is.na(schools$enrolment_fte)),
           locations_with_enrolment = nrow(schools_known)) |>
  knitr::kable()

styled_map <- maplibre(
  style = openfreemap_style("positron"), bounds = schools,
  fitBoundsOptions = list(padding = 40), projection = "mercator",
  height = 500, scrollZoom = FALSE
) |>
  add_circle_layer(
    id = "styled-schools", source = schools_known,
    circle_color = match_expr("school_group", values = school_groups,
                              stops = school_colours, default = "#999999"),
    circle_radius = get_column("radius_px"), circle_opacity = 0.8,
    circle_stroke_color = "white", circle_stroke_width = 1,
    tooltip = "name", popup = "popup_html"
  ) |>
  add_categorical_legend(
    legend_title = "School group · click to filter",
    values = school_groups, colors = school_colours,
    patch_shape = "circle", position = "top-left",
    layer_id = "styled-schools", interactive = TRUE,
    filter_column = "school_group"
  ) |>
  add_categorical_legend(
    legend_title = "Latest reported FTE · circle area",
    values = c("500", "1,000", "2,000"), colors = rep("#526478", 3),
    patch_shape = "circle", sizes = 12 * sqrt(c(500, 1000, 2000) / 500),
    position = "bottom-left", add = TRUE
  ) |>
  add_navigation_control()
styled_map

6. Combine real SA2 polygons, lines, and points

Use fill layers for areas, line layers for outlines or paths, and circle layers for points. We will draw actual SA2 polygons, derive line geometries from their boundaries, and overlay the school locations. These lines represent boundaries; they do not imply travel routes.

# Small display copy: simplified by 25 metres in a projected CRS.
sa2 <- st_read("data/greater-sydney-sa2-display.geojson", quiet = TRUE) |>
  st_transform(4326)
stopifnot(!anyDuplicated(sa2$sa2_code),
          !any(st_is_empty(sa2)), all(st_is_valid(sa2)))
sa2_lines <- st_boundary(sa2)

Registering a source once lets several layers reuse it. Add fills first, then lines and points so the points remain visible.

layered_map <- maplibre(
  style = openfreemap_style("positron"), bounds = sa2,
  fitBoundsOptions = list(padding = 40), projection = "mercator",
  height = 480, scrollZoom = FALSE
) |>
  add_source(id = "sa2-data", data = sa2) |>
  add_fill_layer(
    id = "sa2-fill", source = "sa2-data",
    fill_color = "#56B4E9", fill_opacity = 0.18, popup = "sa2_name"
  ) |>
  add_line_layer(
    id = "sa2-outlines", source = sa2_lines,
    line_color = "#526478",
    line_width = interpolate(property = "zoom", values = c(7, 13),
                             stops = c(0.4, 1.5))
  ) |>
  add_circle_layer(
    id = "school-overlay", source = schools,
    circle_color = "#0072B2", circle_radius = 4,
    circle_stroke_color = "white", circle_stroke_width = 0.6,
    popup = "popup_html", tooltip = "name"
  ) |>
  add_layers_control(
    layers = list("SA2 polygons" = "sa2-fill",
                  "SA2 boundary lines" = "sa2-outlines",
                  "Government education locations" = "school-overlay"),
    position = "top-left", collapsible = FALSE
  ) |>
  add_navigation_control()
layered_map

source = "sa2-data" names a source ID. A layer control uses layer IDs. Here the layers are independent so you can inspect fills, outlines, and points separately. To switch fill and outline together, group their IDs in one vector. See add_source() and add_layers_control().

7. Join population statistics to SA2 boundaries

Population is supplied as a table, while geometry is supplied separately. Join them on the nine-character SA2 code, checking uniqueness and unmatched keys in both directions. Keep identifiers as character values; area names can change or be spelled differently.

population <- read_csv("data/sydney-sa2-population.csv",
  col_types = cols(sa2_code = col_character()))
stopifnot(!anyNA(sa2$sa2_code), !anyNA(population$sa2_code),
          !anyDuplicated(sa2$sa2_code), !anyDuplicated(population$sa2_code))

boundaries_without_data <- anti_join(st_drop_geometry(sa2), population,
                                     by = "sa2_code")
data_without_boundaries <- anti_join(population, st_drop_geometry(sa2),
                                     by = "sa2_code")
stopifnot(nrow(boundaries_without_data) == 0,
          nrow(data_without_boundaries) == 0)

sa2_data <- left_join(sa2, select(population, -sa2_name), by = "sa2_code")
stopifnot(nrow(sa2_data) == nrow(sa2), all(sa2_data$area_km2 > 0),
          all(is.finite(sa2_data$population_2025)),
          all(sa2_data$population_2025 >= 0))
sa2_data <- sa2_data |>
  mutate(
    density_2021 = population_2021 / area_km2,
    density_2025 = population_2025 / area_km2,
    popup_html = paste0(
      "<strong>", htmltools::htmlEscape(sa2_name), "</strong><br>",
      "SA2: ", sa2_code, "<br>Population, June 2025: ",
      format(population_2025, big.mark = ",", trim = TRUE),
      "<br>Area: ", sprintf("%.2f km²", area_km2),
      "<br>Density: ", sprintf("%.0f residents/km²", density_2025)
    )
  )
sa2_data |>
  st_drop_geometry() |> arrange(desc(density_2025)) |>
  select(sa2_name, population_2025, area_km2, density_2025) |> head(6) |>
  knitr::kable(digits = c(NA, 0, 4, 0),
               caption = "The six highest-density Sydney SA2s in June 2025.")

The boundary vintage is ASGS 2021, matching the geography used for these population estimates. Using a newer boundary file just because it is newer can break comparability. The bundled snapshot contains population values for all 373 polygons; unmatched keys and missing measurements are still separate problems that you must check in your own data.

8. Build a choropleth of population density

A choropleth shades areas according to a statistic. A sequential colour scale suits a non-negative density. A diverging scale can suit change around zero. Explain the units and the limits of your colour scale.

density_limits <- c(0, 5000, 25000)
density_colours <- c("#f7fbff", "#6baed6", "#08306b")
missing_colour <- "#bdbdbd"
stopifnot(max(sa2_data$density_2025) <= max(density_limits))

density_map <- maplibre(
  style = openfreemap_style("positron"), bounds = sa2_data,
  fitBoundsOptions = list(padding = 40), projection = "mercator",
  height = 480, scrollZoom = FALSE
) |>
  add_source(id = "density-data", data = sa2_data) |>
  add_fill_layer(
    id = "density-fill", source = "density-data",
    fill_color = interpolate("density_2025", values = density_limits,
      stops = density_colours, na_color = missing_colour),
    fill_opacity = 0.85, popup = "popup_html", tooltip = "sa2_name",
    hover_options = list(fill_opacity = 1)
  ) |>
  add_line_layer(
    id = "density-outline", source = "density-data",
    line_color = "#526478", line_width = 0.5
  ) |>
  add_continuous_legend(
    legend_title = "June 2025 · residents per km²",
    values = density_limits, colors = density_colours,
    position = "bottom-left"
  ) |>
  add_navigation_control()
density_map

band_colours <- c("#eff3ff", "#bdd7e7", "#6baed6", "#2171b5", "#08306b")
band_map <- maplibre(
  style = openfreemap_style("positron"), bounds = sa2_data,
  fitBoundsOptions = list(padding = 40), projection = "mercator",
  height = 480, scrollZoom = FALSE
) |>
  add_fill_layer(
    id = "density-bands", source = sa2_data,
    fill_color = step_expr("density_2025", base = band_colours[1],
      values = c(500, 2000, 5000, 8000), stops = band_colours[-1],
      na_color = missing_colour),
    fill_opacity = 0.85, fill_outline_color = "#526478",
    popup = "popup_html"
  ) |>
  add_categorical_legend(
    legend_title = "June 2025 · residents per km²",
    values = c("0 to <500", "500 to <2,000", "2,000 to <5,000",
               "5,000 to <8,000", "8,000 or more"),
    colors = band_colours, position = "top-left"
  ) |>
  add_navigation_control()
band_map

9. Use a spatial join when there is no area code

The school table has coordinates but no SA2 code. A spatial join attaches the code using location. Use the full-resolution analytical boundaries for this operation, not the simplified display copy.

sa2_full <- st_read("data/greater-sydney-sa2.gpkg", layer = "sa2", quiet = TRUE)
# A local projected CRS makes the boundary predicate explicit.
# EPSG:28356 is GDA94 / MGA zone 56; it is not a universal choice.
schools_local <- st_transform(schools, 28356)
sa2_local <- st_transform(sa2_full, 28356)
memberships <- st_within(schools_local, sa2_local)
membership_count <- lengths(memberships)
stopifnot(all(membership_count == 1L))

schools_with_sa2 <- st_join(
  schools_local, select(sa2_local, sa2_code), join = st_within, left = TRUE
)
school_counts <- schools_with_sa2 |>
  st_drop_geometry() |> count(sa2_code, name = "n_locations")
sa2_with_counts <- sa2 |>
  left_join(school_counts, by = "sa2_code") |>
  mutate(n_locations = coalesce(n_locations, 0L))
stopifnot(sum(sa2_with_counts$n_locations) == nrow(schools))
sa2_with_counts |> st_drop_geometry() |> arrange(desc(n_locations)) |>
  select(sa2_name, n_locations) |> head(8) |>
  knitr::kable(caption = "Government education locations assigned by point location.")

Zero means no locations from this government dataset were assigned to the SA2. It does not mean no schools of any kind, no school-age residents, or no access to education. School catchments and student residence can cross SA2 boundaries.

Points on boundaries, outside the study area, or inside overlapping polygons need an explicit rule. st_within() excludes boundary points in this projected example; a different predicate can match more than one polygon. Investigate exceptions before joining instead of dropping duplicate assignments arbitrarily. See st_join().

10. Explore point concentration with a heatmap

A heatmap smooths nearby points across screen pixels. Its appearance depends on zoom, radius, and weighting. We give each government education location the same weight, including those without an enrolment figure.

heat_map <- maplibre(
  style = openfreemap_style("positron"), bounds = schools,
  fitBoundsOptions = list(padding = 40), projection = "mercator",
  height = 440, scrollZoom = FALSE
) |>
  add_heatmap_layer(
    id = "school-heat", source = schools, heatmap_weight = 1,
    heatmap_radius = interpolate(property = "zoom", values = c(7, 14),
                                 stops = c(12, 32)),
    heatmap_intensity = 1,
    heatmap_color = interpolate(property = "heatmap-density",
      values = c(0, 0.25, 0.5, 0.75, 1),
      stops = c("rgba(255,255,255,0)", "#c6dbef", "#6baed6",
                "#2171b5", "#08306b")),
    heatmap_opacity = 0.8
  ) |>
  add_continuous_legend(
    legend_title = "Relative concentration of locations",
    values = c("Lower", "Higher"), colors = c("#c6dbef", "#08306b"),
    position = "bottom-left"
  ) |>
  add_navigation_control()
heat_map

This is a visual concentration map, not a population-adjusted rate or a test of a statistically significant hotspot. For a defined density, calculate a statistic on documented areas or use an appropriate density estimator. See add_heatmap_layer().

11. Cluster crowded school points

Clustering groups nearby points at the current zoom and displays their feature counts. It reduces overplotting; these groups are visual, rather than fixed statistical clusters.

cluster_map <- maplibre(
  style = openfreemap_style("positron"), bounds = schools,
  fitBoundsOptions = list(padding = 40), projection = "mercator",
  height = 470, scrollZoom = FALSE
) |>
  add_circle_layer(
    id = "schools", source = schools, circle_color = "#D55E00",
    circle_radius = 5, popup = "popup_html",
    cluster_options = cluster_options(
      max_zoom = 12, cluster_radius = 45, count_stops = c(0, 20, 100),
      color_stops = c("#56B4E9", "#0072B2", "#172f46"),
      radius_stops = c(18, 25, 34), text_color = "white"
    )
  ) |>
  set_layout_property(layer_id = "schools-cluster-count",
    name = "text-font", value = list("Noto Sans Regular")) |>
  add_categorical_legend(
    legend_title = "Locations in visual cluster",
    values = c("2–19", "20–99", "100 or more", "Individual location"),
    colors = c("#56B4E9", "#0072B2", "#172f46", "#D55E00"),
    patch_shape = "circle", position = "top-left"
  ) |>
  add_navigation_control()
cluster_map

Supply the sf points directly to add_circle_layer() with cluster_options(). The shortcut creates a separate schools-cluster-count symbol layer. We use set_layout_property() to select Noto Sans Regular, available from the OpenFreeMap glyph service. An unavailable font can leave labels blank while circles still render.

12. Add a population year slider without a server

A slider filters features in the browser, so it works on a static website. The filter property must be numeric. We build one polygon record per SA2 and year from the real 2021–2025 population series; the slider shows only one year’s records at a time.

population_long <- read_csv("data/sydney-sa2-population-long.csv",
  col_types = cols(sa2_code = col_character(), year = col_integer()))
stopifnot(!anyDuplicated(paste(population_long$sa2_code, population_long$year)),
          is.numeric(population_long$year),
          nrow(population_long) == nrow(sa2) * 5)
sa2_time <- left_join(sa2, population_long, by = "sa2_code") |>
  mutate(density = population / area_km2,
    year_popup = paste0(htmltools::htmlEscape(sa2_name), "<br>June ", year,
      "<br>Population: ", population,
      "<br>Density: ", sprintf("%.0f residents/km²", density)))
stopifnot(nrow(sa2_time) == nrow(population_long))

time_map <- maplibre(
  style = openfreemap_style("positron"), bounds = sa2,
  fitBoundsOptions = list(padding = 40), projection = "mercator",
  height = 480, scrollZoom = FALSE
) |>
  add_fill_layer(
    id = "population-by-year", source = sa2_time,
    fill_color = interpolate("density", values = density_limits,
                              stops = density_colours),
    fill_opacity = 0.85, fill_outline_color = "#526478", popup = "year_popup"
  ) |>
  add_slider_control(
    layers = "population-by-year", property = "year", values = 2021:2025,
    labels = as.character(2021:2025), mode = "sequential",
    initial_value = 2021, play_button = TRUE,
    title = "Population at 30 June · selected year",
    width = 230, position = "top-left"
  ) |>
  add_continuous_legend(
    legend_title = "All years · residents per km²",
    values = density_limits, colors = density_colours, position = "bottom-left"
  ) |>
  add_navigation_control()
time_map

These are annual population snapshots, so sequential filtering suits the question. Cumulative mode would overlap several estimates of the same area and should not be used to add resident populations across years. The browser reveals existing records; it does not estimate missing years or fit a trend model.

For dates, use as_time_property(date_column, unit = "date") to produce a numeric property: GeoJSON serializes R dates as strings. For large datasets, avoid repeating complex geometries unnecessarily; the slider also supports changing paint expressions on a single geometry source. See add_slider_control().

13. Compare 2021 and 2025 using a shared scale

A swipe comparison helps inspect changes. Both maps must use comparable statistics, units, geography, and colour limits. Here density uses the same official area denominator and the same 0–25,000 scale.

make_year_map <- function(data, year) {
  density_column <- paste0("density_", year)
  stopifnot(density_column %in% names(data))
  maplibre(
    style = openfreemap_style("positron"), bounds = data,
    fitBoundsOptions = list(padding = 40), projection = "mercator",
    scrollZoom = FALSE
  ) |>
    add_fill_layer(
      id = "comparison-density", source = data,
      fill_color = interpolate(density_column, values = density_limits,
                                stops = density_colours, na_color = missing_colour),
      fill_opacity = 0.85, fill_outline_color = "#526478"
    )
}

comparison_map <- compare(
  make_year_map(sa2_data, 2021), make_year_map(sa2_data, 2025),
  mode = "swipe", height = 480
) |>
  add_legend(
    legend_title = "Left: 2021 · Right: 2025 · residents/km²",
    values = density_limits, colors = density_colours,
    type = "continuous", target = "compare", position = "bottom-left"
  )
comparison_map

Use mode = "sync" for synchronized views. Quantify changes in a table or chart as well: colour differences do not establish why population changed. For example, compute population_2025 - population_2021; for percentage change, check for a zero baseline before dividing. See compare().

14. Explore 3D extrusion

Extrusion assigns display heights to polygons. Perspective and occlusion make precise comparisons difficult, so retain a 2D map or table. This example focuses on central Sydney while drawing the real SA2 dataset.

sa2_3d <- sa2_data |> mutate(height_m = density_2025 / 10)
extrusion_map <- maplibre(
  style = openfreemap_style("positron"), center = c(151.20, -33.89),
  zoom = 11.8, projection = "mercator", pitch = 50, bearing = -20,
  height = 480, scrollZoom = FALSE
) |>
  add_fill_extrusion_layer(
    id = "extruded-sa2", source = sa2_3d,
    fill_extrusion_color = interpolate("density_2025", values = density_limits,
                                       stops = density_colours),
    fill_extrusion_height = get_column("height_m"),
    fill_extrusion_opacity = 0.85, popup = "popup_html"
  ) |>
  add_continuous_legend(
    legend_title = "Residents/km² · height = density ÷ 10 m",
    values = density_limits, colors = density_colours, position = "bottom-left"
  ) |>
  add_navigation_control()
extrusion_map

fill_extrusion_height expects metres, so the display conversion is explicit. projection = "mercator" avoids globe-related extrusion artefacts described in the layer overview. See add_fill_extrusion_layer().

15. Display the real ABS population raster

Vector data describe distinct features. Raster data describe values in grid cells. The bundled GeoTIFF is the ABS 2025 population grid cropped and masked to Greater Sydney, retaining its native 1 km cells in EPSG:3577. Its values are modelled residents per km², not observations of individuals at cell centres.

population_grid <- terra::rast("data/sydney-population-grid-2025.tif")
terra::res(population_grid)

[1] 1000 1000

terra::crs(population_grid, proj = TRUE)

[1] "+proj=aea +lat_0=0 +lon_0=132 +lat_1=-18 +lat_2=-36 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs"

terra::global(population_grid, c("min", "max"), na.rm = TRUE)

                        min      max
population_density_2025   0 24015.78

# Classify on the original analytical grid before making a display image.
# 1 = zero, 2 = >0 to <500, 3 = 500 to <2000, ... 6 = >=8000.
grid_classes <- 1 + (population_grid > 0) + (population_grid >= 500) +
  (population_grid >= 2000) + (population_grid >= 5000) +
  (population_grid >= 8000)
grid_colours <- c("#f7f7f7", "#eff3ff", "#bdd7e7", "#6baed6", "#2171b5", "#08306b")
grid_labels <- c("Zero", ">0 to <500", "500 to <2,000", "2,000 to <5,000",
                 "5,000 to <8,000", "8,000 or more")
rgb <- t(grDevices::col2rgb(grid_colours))
terra::coltab(grid_classes) <- data.frame(
  value = 1:6, red = rgb[, 1], green = rgb[, 2], blue = rgb[, 3], alpha = 255
)

raster_map <- maplibre(
  style = openfreemap_style("positron"), bounds = sa2,
  fitBoundsOptions = list(padding = 40), projection = "mercator",
  height = 480, scrollZoom = FALSE
) |>
  add_image_source(id = "population-grid-image", data = grid_classes) |>
  add_raster_layer(
    id = "population-grid", source = "population-grid-image",
    raster_opacity = 0.85, raster_resampling = "nearest"
  ) |>
  add_categorical_legend(
    legend_title = "Modelled residents/km² · June 2025",
    values = grid_labels, colors = grid_colours, position = "top-left"
  ) |>
  add_navigation_control()
raster_map

The colour table makes this a categorical display raster. Mapgl uses nearest-neighbour reprojection for this path; raster_resampling = "nearest" also prevents blended class colours when zooming. Zero is distinct from missing data. The retained analytical GeoTIFF still contains its original numeric density values.

Mapgl turns the raster into an image; this example does not provide cell-value popups. A display image is not a substitute for the original grid in analysis. Cropping by a region and masking edge cells is approximate at the 1 km grid scale, so do not expect the masked grid to reproduce the exact SA2 population total. See add_image_source() and add_raster_layer().

Use terra:: explicitly: terra and mapgl share function names including interpolate() and compare(). If you attach both packages, qualify the intended function with mapgl:: or terra::. For a large raster, crop or reduce a display copy, or consider a tiled source.

16. Put the map in your Quarto report

Interactive htmlwidgets belong in an HTML report. A normal mapgl map returned from an R chunk becomes an interactive widget when Quarto renders the document. No iframe, screenshot, or Shiny server is needed for that workflow. See Quarto’s htmlwidgets guide.

Use this YAML at the top of an independent report:

---
title: "My DATA1X01 Project 2 report"
format:
  html:
    embed-resources: true
    toc: true
    code-copy: true
execute:
  echo: true
  message: false
  error: false
---

Then put your setup, data preparation, and map construction in ordinary executable R chunks. Finish the map chunk with the map object. For example, a report can use this pattern:

```{r}
#| label: fig-schools
#| fig-cap: "Greater Sydney government education locations. Source: NSW Department of Education, 2 October 2026 snapshot."
my_map <- maplibre(
  style = openfreemap_style("positron"), bounds = schools,
  projection = "mercator", height = 480
) |>
  add_circle_layer(
    id = "school-points", source = schools,
    circle_color = "#0072B2", circle_radius = 6, popup = "name"
  )
my_map
```

@fig-schools shows the education locations included in this analysis.

Do not set eval: false, include: false, or output: false on a chunk whose map should appear. Avoid placing maps inside collapsed tabs during initial troubleshooting. A fig-* label and fig-cap provide a caption and cross-reference. Follow your course’s instructions about showing code and submitting supporting files.

quarto preview index.qmd
quarto render index.qmd --to html

# Run once in the Console if needed: install.packages("chromote")
if (knitr::is_html_output()) {
  density_map
} else {
  print_map(density_map, width = 1000, height = 600)
}

# You can also save files outside a report:
save_map(density_map, "density-map.png", width = 1000, height = 600)
htmlwidgets::saveWidget(density_map, "density-map.html", selfcontained = TRUE)

17. Adapt the workflow to Project 2

Start with the analytical question, then choose a geometry and an encoding that help answer it. For example, point locations can support a question about distribution, while area rates can support a comparison between regions. The appearance of the map should follow from that question.

1. Identify the geographic unit. Is each row a location, an event, a route, or an area? Check the data period and coverage.
2. Prepare and validate the data. Read files reproducibly, verify the CRS, inspect invalid or empty geometries, and account for excluded rows.
3. Join only what is needed. Check unique identifiers and unmatched keys. Keep only columns that are appropriate to publish.
4. Choose the encoding. Use category colours for categories, a sequential scale for a rate, or proportional circle areas for suitable counts. Explain the units.
5. Write the finding. Describe a relevant pattern, quantify it, and discuss the source’s limitations. Pair the map with a chart or table when precise comparison matters.
6. Render and check the result. Make sure the HTML opens and the map works in the same format the reader or marker will receive.

18. Troubleshoot common problems

If the basemap service is temporarily unavailable, you can construct a map using style = basemap_style("#f5f5f2") and your local data layers. This provides a plain background for debugging or a deliberately basemap-free design. It lacks place labels and other basemap context, and the mapping library still needs to be available.

Documentation and reproducibility

This lesson uses original examples informed by Kyle Walker’s mapgl documentation. It targets mapgl 0.5.2 and explicitly selects its basemap to avoid relying on a default that may change between versions. Refer to the current documentation when extending the examples.

• mapgl home and getting started
• Sources and layers
• Map design
• Complete function reference
• mapgl release notes
• sf documentation
• terra documentation
• Quarto and htmlwidgets
• OpenFreeMap and basemap attribution
• ABS population methodology
• NSW school data and dictionary

To record your own environment, run sessionInfo() and retain the output with your project. The supplied R/prepare-sydney-data.R rebuilds the teaching files from bundled official-source snapshots. scripts/fetch-sydney-data.R documents how to download those sources again; a later school snapshot may differ. data/README.md documents the fields, periods, licences, and transformations.