Introduction to SSIMmap

The goal of the SSIMmap package extend the classical SSIM method for irregular lattice-based maps and raster images. The original SSIM method was developed to compare two image mimicking the human visual system. The SSIMmap package applies this method to two types of maps (polygon and raster). The geographical SSIM method incorporates well-developed geographically weighted summary statistics with an adaptive bandwidth kernel function for irregular lattice-based maps. This package includes four key functions: ssim_bandwidth, ssim_constant, ssim_polygon,ssim_raster. Users who want to compare two maps and quantify their similarities can utilize this package and visualize the results using other R packages (e.g.,tmap and ggplot).

Data: Toronto, groups2nm.tif, and single2nm.tif

The package has three example data: 1) Toronto, 2) groups2nm (terra raster object), and 3) single2nm (terra raster object).

First, Toronto(an object of sf class) is an example data, which stands for Toronto, ON. This contains neighborhood deprivation indices(the Canadian Index of Multiple Deprivation and Pampalon index) and the census variable(the percentage of households who commute within the census subdivision of residence) for 2016.

library(SSIMmap)
library(sf)
library(terra)
data("Toronto")
plot(Toronto, border=NA)

The second and third data are terra raster objects. These raster objects represent the geographical location from 38 to 41◦N and 0.5–5◦E, centering on the islands of Ibiza, Mallorca and Menorca. The raster objects, labeled as ‘groups2nm’ and ‘single2nm’, contain information on sperm whales. They represent the predicted probability of presence for groups and singletons at a spatial resolution of 2 nautical miles(NM) on a regular grid. For further information about ‘groups2nm’ and ‘single2nm’ please refer to the paper “Novel application of a quantitative spatial comparison tool to species distribution data”.

single<-system.file("/ex/single2nm.tif", package="SSIMmap")
group<-system.file("/ex/groups2nm.tif", package="SSIMmap")
whale_single<-terra::rast(single)
whale_groups<-terra::rast(group)
plot(whale_single)

plot(whale_groups)

Functions: ssim_bandwidth

This function calculates the bandwidth size for the computation of the Structural Similarity Index (SSIM) on polygon maps. A user can decide whether or not to standardize the variables of maps. If the maps have different ranges of variables or include negative values, it is necessary to use standardize =TRUE. This ensures that the conditions for SSIM(symmetry and boundedness (a range from -1 to 1)) are not violated. If the maps share a similar range, there is no need to use standardize= TRUE.The function provides two options for selecting the bandwidth size, which will be displayed in the console window:

  1. The square root of N.
  2. The trade-off between the bias and variance of the two maps.

The first method for the bandwidth selection calculates the square root of the sample size n, which is based on the relative rate of convergence that could minimize the mean integrated squared error (Abramson,1982).

The second method aims to balance local variances and bias in spatial data (geographically weighted mean), with an intent to minimize the difference between them. The primary goal of this method is to find a bandwidth size that balances the bias-variance trade-off and therefore effectively captures the scale of the underlying spatial structure. More information about the bias-variance trade-off can be found in the paper “On the measurement of bias in geographically weighted regression models”. By default for the second method, the function chooses the mid-point of the optimal trade-off between the bias and variance for map1 and map2. A user can choose the upper (the larger one) or the lower (the smaller one) value between the two of the optimal trade-off between the bias and variance.

The function takes as input an object of sf class, which includes columns necessary for SSIM calculation. It outputs two possible selections for the bandwidth size. Additionally, the function generates a plot illustrating the relationship between bias, variance, and bandwidth size with a vertical line of the square root of N result.

args(ssim_bandwidth)
## function (shape, map1, map2, max_bandwidth = max_bandwidth, standardize = TRUE, 
##     option = "midpoint") 
## NULL

How to execute

# ssim_bandwidth(Toronto,"CIMD_SDD","PP_SDD",max_bandwidth= 200)
# ssim_bandwidth(Toronto,"CIMD_SDD","P_commute",max_bandwidth= 200)

Functions: ssim_constant

This function calculates constants (k1 and k2) for the computation of the SSIM on polygon maps based on the ranges of values. It considers a range of an object of sf class including columns for the SSIM calculation and returns the rescaled constants on the console window.

args(ssim_constant)
## function (shape, map1, map2, standardize = TRUE) 
## NULL

How to execute

ssim_constant(Toronto,"CIMD_SDD","PP_SDD")
## Rescaled K1: 0.00031 Rescaled K2: 0.00094

Functions: ssim_polygon

This function calculates the SSIM index for a given polygon. It takes an object of the sf class as input, which includes the necessary columns for SSIM calculation. The function then returns either the global SSIM values (global = TRUE) or the local SSIM values for each given polygon (global = FALSE). By default, the bandwidth size is determined by the square root of N. A user has the flexibility to select their own bandwidth size or use the suggestions from the ssim_bandwidth function.

args(ssim_polygon)
## function (shape, map1, map2, standardize = TRUE, bandwidth = NULL, 
##     k1 = NULL, k2 = NULL, global = TRUE) 
## NULL

How to execute

ssim_polygon(Toronto,"CIMD_SDD","PP_SDD") 
## 
## 
## |Statistic |      SSIM|       SIM|       SIV|       SIP|
## |:---------|---------:|---------:|---------:|---------:|
## |Mean      | 0.7660945| 0.9985874| 0.9929325| 0.7724228|
## |Min       | 0.3922861| 0.9923626| 0.8943716| 0.3995261|
## |Max       | 0.9036843| 1.0000000| 1.0000000| 0.9081441|
## |SD        | 0.0913415| 0.0015052| 0.0132319| 0.0897773|
ssim_polygon(Toronto,"CIMD_SDD","P_commute") 
## 
## 
## |Statistic |       SSIM|       SIM|       SIV|        SIP|
## |:---------|----------:|---------:|---------:|----------:|
## |Mean      | -0.0280161| 0.9914652| 0.9796289| -0.0279894|
## |Min       | -0.4718162| 0.9016800| 0.9104065| -0.4738823|
## |Max       |  0.5389244| 1.0000000| 1.0000000|  0.5901537|
## |SD        |  0.2187775| 0.0139315| 0.0211432|  0.2263333|
df<-ssim_polygon(Toronto,"CIMD_SDD","PP_SDD",global = FALSE)
df_2<-ssim_polygon(Toronto,"CIMD_SDD","P_commute",global = FALSE) 

plot(df,border=NA)

plot(df_2,border=NA)

Functions: ssim_raster

This function calculates the SSIM index for raster images. It takes as input a image file importing from the terra and returns either the global SSIM values (global=TRUE) or the SSIM values for each given cell as the local SSIM (global=FALSE). Default of the window size is 3*3 and a user can use own window size.

args(ssim_raster)
## function (img1, img2, global = TRUE, w = 3, k1 = NULL, k2 = NULL) 
## NULL

How to execute

ssim_raster(whale_single,whale_groups)
## SSIM: 0.01788 SIM: 0.3739 SIV: 0.49631 SIP: 0.02351
result_raster<-ssim_raster(whale_single,whale_groups,global = FALSE)
plot(result_raster)

Example of visualization: the local SSIM for irregular lattice maps using ggplot2

library(ggplot2)
library(RColorBrewer)

#Visualize the SSIM result (the CIMD vs. Pampalon) on the map
ggplot() +
  geom_sf(data = df, aes(fill = SSIM), color = NA) +
  scale_fill_gradientn(colors = brewer.pal(8, "PRGn"), limits = c(-1, 1)) +
  theme_void() +
  theme(
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    axis.text.x = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks = element_blank()
  )

#Visualize the SSIM result (the CIMD vs. commuting pattern) on the map
ggplot() +
  geom_sf(data = df_2, aes(fill = SSIM), color = NA) +
  scale_fill_gradientn(colors =brewer.pal(8, "PRGn"), limits = c(-1, 1)) +
  theme_void() +
  theme(
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    axis.text.x = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks = element_blank()
  )

Reference