Introduction to SSIMmap

In SSIMmap: The Structural Similarity Index Measure for Maps

The goal of the SSIMmap package is to extend the classical SSIM method to irregular lattice-based maps and raster images. The original SSIM method was developed to compare two images by mimicking the human visual system. The SSIMmap package applies this method to two types of maps (polygon and raster). For irregular lattice-based maps, the geographical SSIM method incorporates geographically weighted summary statistics with an adaptive k-nearest-neighbour (k-NN) Gaussian kernel.

This package includes three key functions: ssim_bandwidth, ssim_polygon, and ssim_raster. In addition to point estimates of SSIM and its components (SIM, SIV, SIP), the polygon and raster functions also support permutation tests for global and local significance, with optional Benjamini–Hochberg FDR correction for local p-values.

Users who want to compare two maps and quantify their similarities can use this package and visualize the results with other R packages (e.g., tmap, ggplot2, patchwork, and tidyterra).

suppressPackageStartupMessages(library(SSIMmap))
knitr::opts_chunk$set(warning = FALSE, message = FALSE) 

Example data

The package ships with the following example data sets:

1. Toronto_SSIM — an sf polygon object for Toronto, ON, containing two neighbourhood deprivation indices (the Canadian Index of Multiple Deprivation, CIMD, and the Pampalon index, PP) and a census variable (Commute: percentage of households commuting within the census subdivision of residence) for 2016.

2. Daily Canadian Fire Weather Index (FWI) maps for British Columbia (BC) — three packed terra raster objects at a 2 km spatial resolution, derived from the Canadian Forest Fire Weather Index System. FWI is a numeric rating that is driven by meteorological inputs such as temperature, wind speed, humidity and precipitation, and typically ranges from near 0 (minimal fire danger) to values exceeding 100 (extreme fire-weather conditions).

The three dates were chosen to represent contrasting fire-weather conditions:

• fwi_0816_bc — 16 August 2023 (peak summer fire season)
• fwi_0818_bc — 18 August 2023 (peak summer fire season)
• fwi_1101_bc — 1 November 2023 (late fall)

Two days within the summer season are expected to be highly similar in space, whereas a summer–fall comparison should show much lower spatial similarity, because fall conditions in western Canada are typically cooler and wetter, with reduced fire danger across most of the province.

Because terra raster objects rely on external pointers, the FWI rasters are stored as PackedSpatRaster objects (created with terra::wrap()) and need to be converted back to SpatRaster with terra::unwrap() before use.

library(SSIMmap)
library(sf)
library(terra)
library(ggplot2)
library(patchwork)
library(tidyterra)

# Polygon data
data("Toronto_SSIM")
plot(Toronto_SSIM["CIMD"], border = NA, main = "CIMD")

# Raster data (FWI for British Columbia)
# Stored as PackedSpatRaster — unwrap to SpatRaster for use with terra
data("fwi_0816_bc"); fwi_0816_bc <- terra::unwrap(fwi_0816_bc)
data("fwi_0818_bc"); fwi_0818_bc <- terra::unwrap(fwi_0818_bc)
data("fwi_1101_bc"); fwi_1101_bc <- terra::unwrap(fwi_1101_bc)

plot(fwi_0816_bc, main = "FWI – 16 August 2023")
plot(fwi_0818_bc, main = "FWI – 18 August 2023")
plot(fwi_1101_bc, main = "FWI – 1 November 2023")

Function: ssim_bandwidth

ssim_bandwidth() selects a bandwidth (number of nearest neighbours, k) for the adaptive k-NN Gaussian kernel used in ssim_polygon(), by computing bias–variance trade-off curves for each of the two input variables. Note that this function does not compute SSIM itself; it is a helper for choosing a suitable bandwidth.

Key arguments:

• shape: an sf polygon object.
• map1, map2: column names in shape for the two variables.
• max_bandwidth: upper bound of k (must be $\geq 12$ and not exceed the number of polygons).
• transform: one of "normal_score" (Blom normal scores), "percentile" (empirical percentiles), "none", or "minmax" (min–max normalization to [0, 1]). Use a transformation when the two variables have very different scales or units.
• option: how to choose a single bandwidth from the suggested range — "midpoint" (default), "lower", or "upper".

The function returns a list with three elements:

• plot: a ggplot object showing standardized bias (solid line) and variance (dashed line) against bandwidth, with the suggested range highlighted.
• bandwidth: the chosen bandwidth k_star.
• tradeoff: the underlying trade-off data and sqrt(n) as a reference.

args(ssim_bandwidth)

How to execute

S1 <- ssim_bandwidth(
  shape         = Toronto_SSIM,
  map1          = "CIMD",
  map2          = "PP",
  max_bandwidth = 100,
  transform     = "percentile",
  option        = "midpoint"
)
S1$bandwidth
S1$plot

S2 <- ssim_bandwidth(
  shape         = Toronto_SSIM,
  map1          = "CIMD",
  map2          = "Commute",
  max_bandwidth = 100,
  transform     = "percentile",
  option        = "midpoint"
)
S2$bandwidth
S2$plot

Combining two bandwidth plots with patchwork

When comparing several pairs of maps, you can place the trade-off plots side-by-side using patchwork:

# Tag each panel and remove individual captions
p1 <- S1$plot +
  labs(tag = "a", caption = NULL) +
  theme(plot.tag = element_text(face = "bold", size = 14))

p2 <- S2$plot +
  labs(tag = "b", caption = NULL) +
  theme(plot.tag = element_text(face = "bold", size = 14))

# Combine side-by-side with a shared legend at the bottom
combined_bandwidth_plot <- (p1 + p2) +
  plot_layout(guides = "collect") &
  theme(legend.position = "bottom")

# Add a single global caption
combined_bandwidth_plot +
  plot_annotation(
    caption = "* Solid line: Bias / Dashed line: Variance",
    theme   = theme(plot.caption = element_text(size = 10, hjust = 0.5))
  )

Function: ssim_polygon

ssim_polygon() computes SSIM and its components (SIM, SIV, SIP) for two numeric variables stored in a polygon sf object. Local statistics are calculated using an adaptive k-nearest-neighbour Gaussian kernel based on polygon centroids. The number of neighbours is controlled by bandwidth; if bandwidth = NULL, the default is ceiling(sqrt(n)), where n is the number of polygons.

Depending on the arguments, the function can return:

• a global summary table for SSIM, SIM, SIV, and SIP when global = TRUE, or
• an sf object with local SSIM, SIM, SIV, and SIP values for each polygon when global = FALSE.

The argument transform controls preprocessing of the two input variables. Available options are "normal_score", "percentile", "minmax", and "none".

When do_test = TRUE, the function performs permutation-based inference:

• Global test (global = TRUE): returns two-sided permutation p-values for global SSIM, SIM, SIV, and SIP.
• Local test (global = FALSE): returns per-polygon SSIM p-values in p_value. If fdr = TRUE, Benjamini-Hochberg FDR-adjusted values are returned in q_value; otherwise, q_value is equal to p_value. The logical column sig flags polygons with q_value < alpha.

Under the null hypothesis, both variables are randomly permuted at each iteration, breaking both spatial structure and association between the two,variables.

args(ssim_polygon)

Global SSIM with permutation test

S1_global <- ssim_polygon(
  shape     = Toronto_SSIM,
  map1      = "CIMD",
  map2      = "PP",
  global    = TRUE,
  bandwidth = 87,
  transform = "percentile",
  k1 = NULL, k2 = NULL,
  do_test = TRUE, R = 1000, fdr = TRUE, alpha = 0.05,
  seed = 1
)

# Global means and permutation p-values
S1_global$p_global

S2_global <- ssim_polygon(
  shape     = Toronto_SSIM,
  map1      = "CIMD",
  map2      = "Commute",
  global    = TRUE,
  bandwidth = 91,
  transform = "percentile",
  do_test = TRUE, R = 1000, fdr = TRUE, alpha = 0.05,
  seed = 1
)
S2_global$p_global

Local SSIM with permutation test and FDR correction

When global = FALSE, the function returns an sf object that can be passed directly to mapping packages such as ggplot2 or tmap.

S1_local <- ssim_polygon(
  shape     = Toronto_SSIM,
  map1      = "CIMD",
  map2      = "PP",
  global    = FALSE,
  bandwidth = 87,
  transform = "percentile",
  do_test = TRUE, R = 1000, fdr = TRUE, alpha = 0.05,
  seed = 1
)
head(S1_local[, c("SSIM", "SIM", "SIV", "SIP", "p_value", "q_value", "sig")])

S2_local <- ssim_polygon(
  shape     = Toronto_SSIM,
  map1      = "CIMD",
  map2      = "Commute",
  global    = FALSE,
  bandwidth = 91,
  transform = "percentile",
  do_test = TRUE, R = 1000, fdr = TRUE, alpha = 0.05,
  seed = 1
)

Visualizing local SSIM with ggplot2

library(RColorBrewer)

# CIMD vs. Pampalon
ggplot() +
  geom_sf(data = S1_local, aes(fill = SSIM), color = NA) +
  scale_fill_gradientn(colors = brewer.pal(8, "PRGn"), limits = c(-1, 1)) +
  theme_void()

# CIMD vs. Commute
ggplot() +
  geom_sf(data = S2_local, aes(fill = SSIM), color = NA) +
  scale_fill_gradientn(colors = brewer.pal(8, "PRGn"), limits = c(-1, 1)) +
  theme_void()

You can also overlay FDR-significant polygons (sig = TRUE) on top of the local SSIM map by adding a separate layer with geom_sf(data = subset(S1_local, sig), ...).

Function: ssim_raster

ssim_raster() computes the Structural Similarity Index Measure (SSIM) between two raster images. It also returns the three SSIM components: SIM, SIV, and SIP.

The function uses a moving window of size $(2w + 1) \times (2w + 1)$.For example, the default w = 1 uses a 3 x 3 window. Inputs must be single-layer terra::SpatRaster objects with the same geometry.

• If global = TRUE, the function returns a list containing a global summary table for SSIM, SIM, SIV, and SIP.
• If global = FALSE, the function returns a multi-layer terra::SpatRaster with per-cell SSIM, SIM, SIV, and SIP.

The argument transform controls the preprocessing applied to both raster maps before SSIM calculation. Available options are "normal_score", "percentile", "minmax", and "none".

When do_test = TRUE, permutation-based statistical inference is performed.

• If global = TRUE, the function returns global p-values and, when fdr = TRUE, FDR-adjusted q-values.
• If global = FALSE and local_test = TRUE, the function returns local p-value layers (SSIM_p, SIM_p, SIV_p, SIP_p), FDR-adjusted q-value layers (SSIM_q, SIM_q, SIV_q, SIP_q) when fdr = TRUE, and significance layers (SSIM_sig, SIM_sig, SIV_sig, SIP_sig).

args(ssim_raster)

Global SSIM

A quick global comparison between the two summer FWI maps (16 vs. 18 August 2023) and between summer and fall (16 August vs. 1 November 2023):

ssim_raster(fwi_0816_bc, fwi_0818_bc)   # summer vs summer (high similarity expected)
ssim_raster(fwi_0816_bc, fwi_1101_bc)   # summer vs fall   (lower similarity expected)

Local SSIM with a permutation test

In the example below we use w = 1 (a $3 \times 3$ moving window), do_test = TRUE and local_test = TRUE to obtain both the local SSIM surfaces and per-cell p-values. The chunks are set to eval = FALSE to keep the vignette build time short — when running interactively you can simply remove that option.

# Comparison 1: summer vs summer
p1_l <- ssim_raster(
  fwi_0816_bc, fwi_0818_bc,
  global     = FALSE,
  w          = 1,
  do_test    = TRUE,
  local_test = TRUE,
    fdr= TRUE,
  R          = 999
)

# Comparison 2: summer vs fall
p2_l <- ssim_raster(
  fwi_0816_bc, fwi_1101_bc,
  global     = FALSE,
  w          = 1,
  do_test    = TRUE,
  local_test = TRUE,
  fdr= TRUE,
  R          = 999
)

Visualizing the SSIM components with tidyterra

We re-project the SSIM result to an Albers Equal-Area projection (EPSG:3005, the standard projection for Statistics Canada / British Columbia) for area-faithful display, then use tidyterra::geom_spatraster() to facet over the four SSIM layers:

crs_albers <- "EPSG:3005"

# Re-project local SSIM result (continuous values → bilinear)
p1_l_alb <- terra::project(p1_l, crs_albers, method = "bilinear")

# Select SSIM + components and assign panel labels
ssim_maps_alb <- p1_l_alb[[c("SSIM", "SIM", "SIV", "SIP")]]
facet_labs <- c(SSIM = "a", SIM = "b", SIV = "c", SIP = "d")

gg <- ggplot() +
  geom_spatraster(data = ssim_maps_alb) +
  facet_wrap(~lyr, ncol = 2, labeller = labeller(lyr = facet_labs)) +
  scale_fill_viridis_c(
    limits   = c(-1, 1),
    na.value = "transparent",
    name     = "Values"
  ) +
  coord_sf(crs = crs_albers, expand = FALSE) +
  theme_minimal() +
  theme(
    strip.text       = element_text(size = 14, face = "bold", hjust = 0),
    strip.background = element_blank(),
    legend.title     = element_text(size = 13, face = "bold"),
    legend.text      = element_text(size = 11),
    panel.background = element_blank(),
    panel.grid       = element_blank()
  )
gg

Comparing two SSIM rasters side-by-side with patchwork

A common workflow is to compare a baseline image against several follow-up images and place the resulting SSIM maps in a single figure. Below we build one panel for each comparison (summer vs summer; summer vs fall) and combine them with a shared legend:

# Re-project the second comparison and keep only SSIM
p2_l_alb         <- terra::project(p2_l, crs_albers, method = "bilinear")
ssim_maps_alb_p2 <- p2_l_alb[[c("SSIM")]]

ssim_maps_alb_p2 <- p2_l_alb[[c("SSIM", "SIM", "SIV", "SIP")]]
facet_labs <- c(SSIM = "a", SIM = "b", SIV = "c", SIP = "d")

gg2 <- ggplot() +
  geom_spatraster(data = ssim_maps_alb_p2) +
  facet_wrap(~lyr, ncol = 2, labeller = labeller(lyr = facet_labs)) +
  scale_fill_viridis_c(
    limits   = c(-1, 1),
    na.value = "transparent",
    name     = "Values"
  ) +
  coord_sf(crs = crs_albers, expand = FALSE) +
  theme_minimal() +
  theme(
    strip.text       = element_text(size = 14, face = "bold", hjust = 0),
    strip.background = element_blank(),
    legend.title     = element_text(size = 13, face = "bold"),
    legend.text      = element_text(size = 11),
    panel.background = element_blank(),
    panel.grid       = element_blank())
gg2

We expect the summer–summer comparison (panel a) to show high SSIM values across most of British Columbia, whereas the summer–fall comparison (panel b) should show substantially lower similarity, reflecting the seasonal shift from dry summer fire weather to cooler, wetter fall conditions.

Per-cell p-value maps

When do_test = TRUE and local_test = TRUE, the returned SpatRaster also contains SSIM_p, SIM_p, SIV_p, and SIV_p layers, which can be plotted directly:

p_maps <- p1_l[[c("SSIM_p", "SIM_p", "SIV_p", "SIV_p")]]
plot(p_maps,
     main = c("p-value: SSIM", "p-value: SIM",
              "p-value: SIV",  "p-value: SIP"))

References

• Brunsdon, C., Fotheringham, A.S. and Charlton, M. (2002). Geographically weighted summary statistics — a framework for localised exploratory data analysis. Computers, Environment and Urban Systems, 26(6), pp. 501–524. https://doi.org/10.1016/S0198-9715(01)00009-6
• Dimitrakopoulos, A.P., Bemmerzouk, A.M. and Mitsopoulos, I.D. (2011). Evaluation of the Canadian fire weather index system in an eastern Mediterranean environment. Meteorological Applications, 18(1), pp. 83–93. https://doi.org/10.1002/met.214
• Lu, B., Harris, P., Charlton, M. and Brunsdon, C. (2014). The GWmodel R package: further topics for exploring spatial heterogeneity using geographically weighted models. Geo-spatial Information Science, 17(2), pp. 85–101. https://doi.org/10.1080/10095020.2014.917453
• Wang, Z., Bovik, A.C., Sheikh, H.R. and Simoncelli, E.P. (2004). Image quality assessment: from error visibility to structural similarity. IEEE Transactions on Image Processing, 13(4), pp. 600–612. https://doi.org/10.1109/TIP.2003.819861
• Benjamini, Y. and Hochberg, Y. (1995). Controlling the false discovery rate: a practical and powerful approach to multiple testing. Journal of the Royal Statistical Society: Series B, 57(1), pp. 289–300. https://doi.org/10.1111/j.2517-6161.1995.tb02031.x

SSIMmap documentation built on May 8, 2026, 1:07 a.m.