## ----setup, echo = FALSE------------------------------------------------------
suppressPackageStartupMessages(library(SSIMmap))
knitr::opts_chunk$set(warning = FALSE, message = FALSE) 

## ----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")

## -----------------------------------------------------------------------------
args(ssim_bandwidth)

## ----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

## ----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

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

## -----------------------------------------------------------------------------
args(ssim_polygon)

## -----------------------------------------------------------------------------
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

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

## ----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()

## -----------------------------------------------------------------------------
args(ssim_raster)

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

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

## ----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

## ----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

## ----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"))