Introduction to SSIMmap

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).

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:

The function returns a list with three elements:

args(ssim_bandwidth)
## function (shape, map1, map2, max_bandwidth, transform = c("normal_score", 
##     "percentile", "none", "minmax"), option = "midpoint") 
## NULL

How to execute

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

S2 <- ssim_bandwidth(
  shape         = Toronto_SSIM,
  map1          = "CIMD",
  map2          = "Commute",
  max_bandwidth = 100,
  transform     = "percentile",
  option        = "midpoint"
)
S2$bandwidth
## [1] 54
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:

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:

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)
## function (shape, map1, map2, global = TRUE, bandwidth = NULL, 
##     transform = c("normal_score", "percentile", "none", "minmax"), 
##     k1 = NULL, k2 = NULL, do_test = FALSE, R = 1000, fdr = TRUE, 
##     alpha = 0.05, seed = NULL) 
## NULL

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
)
## 
## 
## |Statistic |      SSIM|       SIM|       SIV|       SIP|
## |:---------|---------:|---------:|---------:|---------:|
## |Mean      | 0.8557107| 0.9968253| 0.9965763| 0.8614715|
## |Min       | 0.7107253| 0.9779483| 0.9739425| 0.7121337|
## |Max       | 0.9397111| 0.9999999| 0.9999998| 0.9405333|
## |SD        | 0.0507395| 0.0046724| 0.0039565| 0.0521172|
## 
## Global permutation p-values (two-sided):
## 
## 
## |Component | Global_Mean| P_value|
## |:---------|-----------:|-------:|
## |SSIM      |      0.8557|  0.0010|
## |SIM       |      0.9968|  0.7413|
## |SIV       |      0.9966|  0.0699|
## |SIP       |      0.8615|  0.0010|
# Global means and permutation p-values
S1_global$p_global
##   Component Global_Mean     P_value
## 1      SSIM   0.8557107 0.000999001
## 2       SIM   0.9968253 0.741258741
## 3       SIV   0.9965763 0.069930070
## 4       SIP   0.8614715 0.000999001
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
)
## 
## 
## |Statistic |       SSIM|       SIM|       SIV|        SIP|
## |:---------|----------:|---------:|---------:|----------:|
## |Mean      | -0.2732415| 0.8240093| 0.9498834| -0.3225624|
## |Min       | -0.6975540| 0.4251059| 0.7150255| -0.7138978|
## |Max       |  0.0544271| 0.9999914| 0.9999978|  0.1675890|
## |SD        |  0.1997007| 0.1790784| 0.0730944|  0.2186079|
## 
## Global permutation p-values (two-sided):
## 
## 
## |Component | Global_Mean| P_value|
## |:---------|-----------:|-------:|
## |SSIM      |     -0.2732|   0.001|
## |SIM       |      0.8240|   0.001|
## |SIV       |      0.9499|   0.001|
## |SIP       |     -0.3226|   0.001|
S2_global$p_global
##   Component Global_Mean     P_value
## 1      SSIM  -0.2732415 0.000999001
## 2       SIM   0.8240093 0.000999001
## 3       SIV   0.9498834 0.000999001
## 4       SIP  -0.3225624 0.000999001

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")])
## Simple feature collection with 6 features and 7 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 7228684 ymin: 946321.9 xmax: 7243310 ymax: 954707.4
## Projected CRS: NAD83 / Statistics Canada Lambert
##        SSIM       SIM       SIV       SIP     p_value     q_value  sig
## 1 0.8591248 0.9933423 0.9967759 0.8676804 0.000999001 0.000999001 TRUE
## 2 0.8573882 0.9908127 0.9959857 0.8688260 0.000999001 0.000999001 TRUE
## 3 0.8638841 0.9867431 0.9946716 0.8801804 0.000999001 0.000999001 TRUE
## 4 0.8582752 0.9934988 0.9965025 0.8669235 0.000999001 0.000999001 TRUE
## 5 0.8406321 0.9969807 0.9993703 0.8437092 0.000999001 0.000999001 TRUE
## 6 0.8732341 0.9840033 0.9942157 0.8925931 0.000999001 0.000999001 TRUE
##                         geometry
## 1 MULTIPOLYGON (((7238905 952...
## 2 MULTIPOLYGON (((7233381 952...
## 3 MULTIPOLYGON (((7232286 949...
## 4 MULTIPOLYGON (((7235821 950...
## 5 MULTIPOLYGON (((7240141 950...
## 6 MULTIPOLYGON (((7229747 949...
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.

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.

args(ssim_raster)
## function (map1, map2, global = TRUE, w = 1, transform = c("normal_score", 
##     "percentile", "none", "minmax"), k1 = NULL, k2 = NULL, do_test = FALSE, 
##     local_test = FALSE, R = 1000, fdr = TRUE, alpha = 0.05, seed = NULL) 
## NULL

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)
## $summary
##   metric      mean         min       max         sd
## 1   SSIM 0.5380046 -0.98968521 0.9999790 0.58855978
## 2    SIM 0.5507247 -0.99427575 1.0000000 0.60271292
## 3    SIV 0.9822690  0.40999905 1.0000000 0.03216969
## 4    SIP 0.9924066  0.03063862 0.9999998 0.02029212
## 
## $settings
## $settings$global
## [1] TRUE
## 
## $settings$w
## [1] 1
## 
## $settings$transform
## [1] "normal_score"
## 
## $settings$do_test
## [1] FALSE
## 
## $settings$local_test
## [1] FALSE
ssim_raster(fwi_0816_bc, fwi_1101_bc)   # summer vs fall   (lower similarity expected)
## $summary
##   metric      mean        min       max         sd
## 1   SSIM 0.2909122 -0.9914894 0.9998444 0.61660743
## 2    SIM 0.3063775 -0.9964454 1.0000000 0.65287253
## 3    SIV 0.9561909  0.2553303 1.0000000 0.06145141
## 4    SIP 0.9845342 -0.1710329 0.9999987 0.04173728
## 
## $settings
## $settings$global
## [1] TRUE
## 
## $settings$w
## [1] 1
## 
## $settings$transform
## [1] "normal_score"
## 
## $settings$do_test
## [1] FALSE
## 
## $settings$local_test
## [1] FALSE

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