R/ssim_polygon.R

In SSIMmap: The Structural Similarity Index Measure for Maps

#' The SSIM index for polygon maps.
#'
#' This function computes the SSIM, a measure of similarity between two polygon maps.
#'
#' @param shape A \code{sf} polygon containing the attributes that can create polygon-based maps.
#' @param map1 The name of the first map to compare as a column in the shape.
#' @param map2 The name of the second map to compare as a column in the shape.
#' @param global If global is TRUE, returning the global average of SSIM, SIM, SIV, and SIP. If the option is FALSE, a \code{sf} polygon containing the SSIM, SIM, SIV, and SIP for each polygon is returned.
#' Default is TRUE.
#' @param k1 The constant used in the SSIM calculation. Default is NULL, in which case it is computed from the maximum value of the variables.
#' @param k2 The constant used in the SSIM calculation. Default is NULL, in which case it is computed from the maximum value of the variables.
#' @param bandwidth The size of bandwidth for the Gaussian kernel weighting used in the SSIM calculation. Default is the square root of N.
#' @param standardize If TRUE, standardize the variables before computing the SSIM. Default is TRUE.
#'
#'
#' @details This function computes the SSIM index for two polygon maps.
#' @return If global is TRUE, returning the global average SSIM, SIM, SIV, and SIP on the console window.
#' If global is FALSE, a \code{sf} object containing the SSIM, SIM, SIV, and SIP for each polygon.
#'
#' @importFrom dplyr select contains
#' @importFrom knitr kable
#' @importFrom stats sd var
#' @examples
#' # Load example sf object Toronto Area with three attributes for maps:
#' # Pampalon Index, CIMD Index,
#' # and percentage of household commuting within the same Census Sub-Division of residence.
#' data("Toronto")
#'
#' # Mapping two attributes
#' plot(Toronto$CIMD_SDD)
#' plot(Toronto$PP_SDD)
#' # Finding global ssim
#'\donttest{
#' ssim_polygon(Toronto,"CIMD_SDD","PP_SDD")
#' df<-ssim_polygon(Toronto,"CIMD_SDD","PP_SDD",global=FALSE)
#' }
#' @export ssim_polygon


# Functions ---------------------------------------------------------------

ssim_polygon<-function(shape,map1,map2,standardize=TRUE,bandwidth=NULL,k1=NULL,k2=NULL,global=TRUE){
  if(map1==map2){
    stop("variables are identical")
  }

  num_rows <- nrow(shape)
  sqrt_num_rows <- round(sqrt(num_rows),0)


  if(standardize){
    shape_df<-as.data.frame(shape)
    shape_df$z_score_map1<-(shape_df[,map1]-mean(shape_df[,map1]))/sd(shape_df[,map1])
    shape_df$z_score_map2<-(shape_df[,map2]-mean(shape_df[,map2]))/sd(shape_df[,map2])
    min<-min(min(shape_df$z_score_map1),min(shape_df$z_score_map2))
    shape_df$z_score_map1<-shape_df$z_score_map1-min
    shape_df$z_score_map2<-shape_df$z_score_map2-min
    z_scores<-as.data.frame(cbind(shape_df$z_score_map1,shape_df$z_score_map2))
    colnames(z_scores)<- paste0("zscore",colnames(z_scores))
    names(z_scores)<- sub("V1",map1,names(z_scores))
    names(z_scores)<- sub("V2",map2,names(z_scores))
    shape_merged<-cbind(shape,z_scores)
    map1<-as.character(colnames(z_scores[1]))
    map2<-as.character(colnames(z_scores[2]))
    if(is.null(bandwidth)){
      gwss_result<-suppressWarnings(gwss_new(shape_merged,vars = c(map1,map2),bw=sqrt_num_rows))

    }
    else{
      gwss_result<-suppressWarnings(gwss_new(shape_merged,vars = c(map1,map2),bw=bandwidth))
    }
    mean<-dplyr::select(gwss_result,contains("LM"))
    sd<-dplyr::select(gwss_result,contains("LSD"))
    cov<-dplyr::select(gwss_result,contains("Cov"))
    merged_df<-cbind(mean,sd,cov,z_scores)
    merged_df<-cbind(merged_df,shape)



  }

  else{

    if(is.null(bandwidth)){
      gwss_result<-gwss_new(shape,vars = c(map1,map2),bw=sqrt_num_rows)

    }else{
      gwss_result<-gwss_new(shape,vars = c(map1,map2),bw=bandwidth)
    }
    mean<-dplyr::select(gwss_result,contains("LM"))
    sd<-dplyr::select(gwss_result,contains("LSD"))
    cov<-dplyr::select(gwss_result,contains("Cov"))
    merged_df<-cbind(mean,sd,cov)
    colnames(merged_df)<-gsub("zscore", "", colnames(merged_df))
    merged_df<-cbind(merged_df,shape)




  }


  shape_df<-as.data.frame(merged_df)
  map1<-as.data.frame(dplyr::select(shape_df,contains(map1)))
  map2<-as.data.frame(dplyr::select(shape_df,contains(map2)))

  globalMin <- min(map1,map2)
  l<-max(map1,map2)
  max<-max(map1,map2)

  map1_LM<-dplyr::select(map1,contains('_LM'))
  map2_LM<-dplyr::select(map2,contains('_LM'))
  map1_LSD<-dplyr::select(map1,contains('_LSD'))
  map2_LSD<-dplyr::select(map2,contains('_LSD'))
  Cov_map1_map2<-dplyr::select(shape_df,contains('Cov_'))
  l <- l-globalMin

  if(is.null(k1)){
    k1<-(0.01*max)/255
  }
  if(is.null(k2)){
    k2<-(0.03*max)/255
  }

  if(k1>1){
    stop("k1 needs to be smaller than 1")
  }

  if(k2>1){
    stop("k2 needs to be smaller than 1")
  }

  k <- c(k1, k2)
  C1 <-(k[1]*l)^2
  C2 <-(k[2]*l)^2
  C3 <-C2/2

  mu1<-map1_LM
  mu2<-map2_LM
  sig1<-map1_LSD
  sig2<-map2_LSD
  sig12<-Cov_map1_map2
  SIM<- ((2*mu1*mu2)+C1) / (mu1^2 + mu2^2 + C1)
  if(any(SIM>1) || any(SIM<0)){
    stop("violating the boundness rule,you need to rescale the attributes for maps")
  }

  SIV <- ((2*sig1*sig2)+C2) / (sig1^2 + sig2^2 + C2)
  if(any(SIV>1) || any(SIV<0)){
    stop("violating the boundness rule,you need to rescale the attributes for maps")
  }
  SIP <- (sig12 + C3) / (sig1 * sig2 + C3)
  if(any(SIP>1)|| any(SIP< (-1))){
    stop("violating the boundness rule,you need to rescale the attributes for maps")
  }
  SSIM <- SIM*SIV*SIP
  df<-as.data.frame(cbind(SSIM, SIM, SIV, SIP))
  colnames(df)<-c("SSIM","SIM","SIV","SIP")

  if(global){
    SSIM_mean <- mean(df$SSIM)
    SIM_mean <- mean(df$SIM)
    SIV_mean <- mean(df$SIV)
    SIP_mean <- mean(df$SIP)

    SSIM_min <-min(df$SSIM)
    SIM_min <- min(df$SIM)
    SIV_min <- min(df$SIV)
    SIP_min <- min(df$SIP)

    SSIM_max <- max(df$SSIM)
    SIM_max <- max(df$SIM)
    SIV_max <- max(df$SIV)
    SIP_max <- max(df$SIP)

    SSIM_sd <- sd(df$SSIM)
    SIM_sd <- sd(df$SIM)
    SIV_sd <- sd(df$SIV)
    SIP_sd <- sd(df$SIP)

    result <- data.frame(
      Statistic = c("Mean", "Min", "Max", "SD"),
      SSIM = c(SSIM_mean, SSIM_min, SSIM_max, SSIM_sd),
      SIM = c(SIM_mean, SIM_min, SIM_max, SIM_sd),
      SIV = c(SIV_mean, SIV_min, SIV_max, SIV_sd),
      SIP = c(SIP_mean, SIP_min, SIP_max, SIP_sd)
    )

    print(knitr::kable(result, row.names = FALSE))
  }
  else{
    SSIM_results<-cbind(shape,df)
    return(SSIM_results)
  }
}