--- title: "Introduction to SSIMmap" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Introduction to SSIMmap} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- The goal of the *SSIMmap* package is to extend [the classical SSIM method](https://ieeexplore.ieee.org/stamp/stamp.jsp?arnumber=1284395) 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](https://www.sciencedirect.com/science/article/pii/S0198971501000096?via%3Dihub) 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*](https://CRAN.R-project.org/package=tmap), [*ggplot2*](https://CRAN.R-project.org/package=ggplot2), [*patchwork*](https://CRAN.R-project.org/package=patchwork), and [*tidyterra*](https://CRAN.R-project.org/package=tidyterra)). ```{r setup, echo = FALSE} 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. ```{r, message = FALSE} 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. ```{r} args(ssim_bandwidth) ``` ### How to execute ```{r, fig.width = 6, fig.height = 4} S1 <- ssim_bandwidth( shape = Toronto_SSIM, map1 = "CIMD", map2 = "PP", max_bandwidth = 100, transform = "percentile", option = "midpoint" ) S1$bandwidth S1$plot ``` ```{r, fig.width = 6, fig.height = 4} 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`**: ```{r, fig.width = 10, fig.height = 5} # 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. ```{r} args(ssim_polygon) ``` ### Global SSIM with permutation test ```{r} 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 ``` ```{r} 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`. ```{r} 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")]) ``` ```{r} 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` ```{r, fig.show = "hold", out.width = "45%"} 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`). ```{r} 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): ```{r} 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. ```{r, eval = FALSE} # 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: ```{r, eval = FALSE} 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: ```{r, eval = FALSE} # 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: ```{r, eval = FALSE} 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. * 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. * 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. * 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. * 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.