inst/figures/render_figures.R

In SSIMmap: The Structural Similarity Index Measure for Maps

# Portable manuscript artwork from freshly reproduced numerical outputs.
# source(system.file("figures", "render_figures.R", package = "SSIMmap"))
# render_manuscript_figures("reproduction-output")
#
# The default output is 600-dpi CMYK TIFF with LZW compression. These R layouts preserve
# scientific values, classes, panel order and palette; they are not pixel copies
# of the separately prepared CMYK journal print masters. Figure 3 uses a
# simplified provincial boundary for display only; the raster mask is unchanged.
# Figures 3 and 4 are reprojected for display only (nearest and bilinear,
# respectively); all scientific computations use the original raster grids.
# Figure 7 requires the full fixed-seed permutation results, never legacy q.

render_manuscript_figures <- function(reproduction_dir,
                                      output_dir = file.path(reproduction_dir, "figures"),
                                      figures = 3:7, dpi = 600, overwrite = FALSE,
                                      colourspace = c("CMYK", "RGB"),
                                      boundary_file = system.file("figures", "BC_display_boundary.rds", package="SSIMmap")) {
  colourspace <- match.arg(colourspace)
  needed <- c("terra", "sf", "ggplot2", "patchwork", "tidyterra")
  if (colourspace == "CMYK") needed <- c(needed,"magick")
  missing <- needed[!vapply(needed, requireNamespace, logical(1), quietly = TRUE)]
  if (length(missing)) stop("Install figure dependencies: ", paste(missing, collapse = ", "))
  if (!capabilities("tiff")) stop("This R installation has no TIFF device.")
  if (any(!figures %in% 3:7)) stop("figures must contain numbers 3 through 7.")
  if (length(dpi) != 1L || !is.finite(dpi) || dpi <= 0) stop("dpi must be positive.")
  reproduction_dir <- normalizePath(reproduction_dir, mustWork = TRUE)
  fd <- file.path(reproduction_dir, "figure_data")
  paths <- list(
    `3` = file.path(fd, paste0("Figure3_", c("20230816", "20230818", "20231111"), ".tif")),
    `4` = file.path(fd, paste0("Figure4_", c("summer_summer", "summer_fall"), "_native.tif")),
    `5` = file.path(fd, "Toronto.gpkg"),
    `6` = file.path(fd, c("Figure6_curves.csv", "Figure6_bandwidth.csv")),
    `7` = file.path(reproduction_dir, "Toronto_tests",
                    c("cimd_msdi_b87_local.rds", "cimd_commute_b91_local.rds")))
  required <- unlist(paths[as.character(figures)], use.names = FALSE)
  absent <- required[!file.exists(required)]
  if (length(absent)) stop("Missing regenerated input(s). Run reproduce_manuscript(..., full=TRUE) for Figure 7:\n",
                           paste(absent, collapse = "\n"))
  dir.create(output_dir, recursive = TRUE, showWarnings = FALSE)
  targets <- file.path(output_dir, sprintf("Figure_%02d.tiff", figures))
  if (!overwrite && any(file.exists(targets))) stop("Output exists; use overwrite=TRUE to replace it.")

  g <- asNamespace("ggplot2")
  # Local aliases keep this source helper independent of attached packages.
  ggplot <- g$ggplot; aes <- g$aes; annotate <- g$annotate
  theme <- g$theme; element_text <- g$element_text; element_rect <- g$element_rect
  colors5 <- c("#440154", "#3B528B", "#21908D", "#5DC963", "#FDE725")
  qcolors <- c("#CA0020", "#F4A582", "#FFFFFF", "#BABABA", "#404040")
  toronto_crs <- "+proj=aea +lat_0=40 +lon_0=-96 +lat_1=50 +lat_2=70 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs"
  audit <- list()
  map_theme <- function() g$theme_void(base_size = 8, base_family = "sans") +
    theme(plot.tag = element_text(face = "plain", size = 10),
          plot.tag.position = c(0.02, 0.98), legend.title = element_text(size = 8, face = "plain"),
          legend.text = element_text(size = 7),
          plot.background = element_rect(fill = "white", colour = NA),
          panel.background = element_rect(fill = "white", colour = NA),
          plot.margin = g$margin(3, 3, 3, 3))
  save_plot <- function(p, number, height) {
    f <- file.path(output_dir, sprintf("Figure_%02d.tiff", number))
    args <- list(filename = f, width = round(6.8*dpi), height = round(height*dpi), units = "px",
                 res = dpi, compression = "lzw", bg = "white")
    if (capabilities("cairo")) args$type <- "cairo"
    do.call(grDevices::tiff, args)
    device <- grDevices::dev.cur()
    tryCatch(print(p),finally=grDevices::dev.off(device))
    if (colourspace == "CMYK") {
      im <- magick::image_read(f)
      im <- magick::image_background(im,"white",flatten=TRUE)
      im <- magick::image_convert(im,colorspace="cmyk",depth=8,matte=FALSE)
      magick::image_write(im,path=f,format="tiff",depth=8,
                          density=paste0(dpi,"x",dpi),compression="lzw")
      if (!identical(magick::image_info(im)$colorspace,"CMYK")) stop("CMYK conversion failed.")
    }
    invisible(f)
  }
  # A scale bar and two-tone north arrow in the map's meter coordinates.
  furniture <- function(p, x, y, length_m, arrow_x, arrow_h, text_size = 2.5, middle = TRUE) {
    aw <- arrow_h * 0.55
    tick <- if (middle) c(0, .5, 1) else c(0, 1)
    p + annotate("segment", x=x, xend=x+length_m, y=y, yend=y, linewidth=.25) +
      annotate("segment", x=x+tick*length_m,
               xend=x+tick*length_m, y=y, yend=y-arrow_h*.16, linewidth=.25) +
      annotate("text", x=x+tick*length_m, y=y+arrow_h*.25,
               label=format(tick*length_m/1000, trim=TRUE),
               size=text_size, vjust=0) +
      annotate("text", x=x+length_m*if(middle)1.3 else 1.65, y=y+arrow_h*.25,
               label="km", size=text_size, vjust=0, hjust=0) +
      annotate("polygon", x=c(arrow_x, arrow_x, arrow_x-aw/2),
               y=c(y+arrow_h, y+arrow_h*.34, y), fill="grey75", colour=NA) +
      annotate("polygon", x=c(arrow_x, arrow_x+aw/2, arrow_x),
               y=c(y+arrow_h, y, y+arrow_h*.34), fill="white", colour=NA) +
      annotate("path", x=c(arrow_x, arrow_x+aw/2, arrow_x, arrow_x-aw/2, arrow_x),
               y=c(y+arrow_h, y, y+arrow_h*.34, y, y+arrow_h), linewidth=.35)
  }
  if (3 %in% figures) {
    if (!nzchar(boundary_file) || !file.exists(boundary_file)) stop("The packaged BC display boundary was not found.")
    boundary <- sf::st_transform(readRDS(boundary_file),"EPSG:3005")
    raw <- lapply(paths$`3`, terra::rast)
    shown <- lapply(raw, terra::project, y="EPSG:3005", method="near")
    reds <- c("#fff5f0", "#fee0d2", "#fcbba1", "#fc9272", "#fb6a4a", "#ef3b2c", "#cb181d", "#a50f15", "#67000d")
    # QGIS interpolates this ramp in RGB. Dense RGB samples also preserve that
    # interpolation through ggplot2's Lab interpolation between adjacent stops.
    ramp_at <- c(0,.13,.26,.39,.52,.65,.78,.90,1)
    rgb_channels <- grDevices::col2rgb(reds)
    dense <- vapply(1:3,function(i)stats::approx(ramp_at,rgb_channels[i,],xout=seq(0,1,length.out=1101))$y,numeric(1101))
    dense_reds <- grDevices::rgb(dense[,1],dense[,2],dense[,3],maxColorValue=255)
    plots <- lapply(seq_along(shown), function(i) {
      p <- ggplot() + g$geom_sf(data=boundary,fill=NA,colour="#232323",linewidth=.15) +
        tidyterra::geom_spatraster(data=shown[[i]], maxcell=1e8) +
        g$scale_fill_gradientn(colours=dense_reds, values=seq(0,1,length.out=1101),
          space="Lab", limits=c(0,110), breaks=c(0,110), oob=scales::squish,
          name="FWI", na.value="transparent") +
        g$coord_sf(crs="EPSG:3005", xlim=c(145422.206660677,2097700.61868637),
          ylim=c(-595489.229186362,2755693.19199417), datum=NA, expand=FALSE, clip="off") +
        g$labs(tag=paste0(letters[i], ")")) + map_theme()
      if (i == 3) p <- furniture(p, 6e5, -1.5e5, 2e5, 3.5e5, 2.6e5, 2.25, middle=FALSE)
      p
    })
    p <- patchwork::wrap_plots(plots, nrow=1, guides="collect") & theme(legend.position="right")
    save_plot(p, 3, 3.8)
    audit$Figure3 <- data.frame(panel=c("a)","b)","c)"), date=c("2023-08-16","2023-08-18","2023-11-11"),
      cells=vapply(raw,terra::ncell,numeric(1)), valid=vapply(raw,function(r)sum(!is.na(terra::values(r))),numeric(1)),
      legend_min=0,legend_max=110, display_resampling="nearest")
  }
  if (4 %in% figures) {
    raw <- lapply(paths$`4`, function(f) terra::rast(f)[["SSIM"]])
    shown <- lapply(raw, terra::project, y="EPSG:3005", method="bilinear")
    xs <- range(vapply(shown,terra::xmin,numeric(1)),vapply(shown,terra::xmax,numeric(1)))
    ys <- range(vapply(shown,terra::ymin,numeric(1)),vapply(shown,terra::ymax,numeric(1)))
    plots <- lapply(seq_along(shown), function(i) {
      p <- ggplot() + tidyterra::geom_spatraster(data=shown[[i]],maxcell=1e8) +
        g$scale_fill_viridis_c(limits=c(-1,1),breaks=c(-1,-.5,0,.5,1),na.value="transparent",name="SSIM") +
        g$coord_sf(crs="EPSG:3005",xlim=xs,ylim=ys,datum=NA,expand=FALSE,clip="off") +
        g$labs(tag=paste0(letters[i],")")) + map_theme()
      if (i == 2) p <- furniture(p,xs[1]+.23*diff(xs),ys[1]+.06*diff(ys),3e5,
                                 xs[1]+.16*diff(xs),.08*diff(ys),2.5)
      p
    })
    p <- patchwork::wrap_plots(plots,nrow=1,guides="collect") & theme(legend.position="right")
    save_plot(p,4,6.8*7/12)
    audit$Figure4 <- data.frame(panel=c("a)","b)"),
      minimum=vapply(raw,function(r)terra::global(r,"min",na.rm=TRUE)[1,1],numeric(1)),
      maximum=vapply(raw,function(r)terra::global(r,"max",na.rm=TRUE)[1,1],numeric(1)),
      legend_min=-1,legend_max=1,display_resampling="bilinear")
  }
  toronto_map <- function(x, classes, colors, tag, xlim=c(1320000,1372000), ylim=c(550000,593000)) {
    x$.figure_class <- as.integer(classes)
    p <- ggplot(x) + g$geom_sf(aes(fill=factor(.figure_class)),colour="#00000099",linewidth=.12,
                             show.legend=FALSE) +
      g$scale_fill_manual(values=setNames(colors,seq_along(colors)),drop=FALSE) +
      g$coord_sf(crs=toronto_crs,xlim=xlim,ylim=ylim,datum=NA,expand=FALSE,clip="off") +
      g$labs(tag=tag) + map_theme()
    p
  }
  if (5 %in% figures) {
    x <- sf::st_transform(sf::st_read(paths$`5`,layer="input_values",quiet=TRUE),toronto_crs)
    variables <- c("CIMD","PP","Commute")
    breaks <- lapply(variables,function(v)as.numeric(stats::quantile(x[[v]],seq(0,1,.2),type=7)))
    labels <- list(c("-0.66 - 0.08","0.08 - 0.31","0.31 - 0.6","0.6 - 1.09","1.09 - 2.75"),
                   c("-0.05 - -0.02","-0.02 - 0","0 - 0.01","0.01 - 0.02","0.02 - 0.06"),
                   c("1.18 - 2.23","2.23 - 3.62","3.62 - 5.16","5.16 - 6.7","6.7 - 12.61"))
    classes <- lapply(seq_along(variables),function(i)cut(x[[variables[i]]],breaks[[i]],include.lowest=TRUE,labels=FALSE))
    maps <- lapply(seq_along(variables),function(i)toronto_map(x,classes[[i]],colors5,paste0(letters[i],")")))
    lp <- ggplot() + g$xlim(0,1) + g$ylim(0,1) + map_theme()
    titles <- c("a)  CIMD Index","b)  MSDI Index","c)  Derived commuting\n     measure (%)")
    for (i in 1:3) {
      xx <- c(.01,.53,.53)[i]; yy <- c(.89,.89,.43)[i]
      lp <- lp + annotate("text",x=xx,y=yy+if(i==3).025 else 0,label=titles[i],size=2.6,hjust=0)
      for(j in 1:5) {
        yyj <- yy-j*.066
        lp <- lp + annotate("rect",xmin=xx,xmax=xx+.05,ymin=yyj-.020,ymax=yyj+.020,
                            fill=rev(colors5)[j],colour="grey50",linewidth=.1) +
          annotate("text",x=xx+.07,y=yyj,label=rev(labels[[i]])[j],size=2.6,hjust=0)
      }
    }
    # Scale is drawn in panel c in projected meters, not in the legend panel.
    maps[[3]] <- furniture(maps[[3]],1330000,550800,10000,1324000,4500,2.4)
    save_plot(patchwork::wrap_plots(c(maps,list(lp)),ncol=2),5,6.8*210/297)
    audit$Figure5 <- do.call(rbind,lapply(seq_along(variables),function(i)
      data.frame(variable=variables[i],class=1:5,lower=breaks[[i]][1:5],upper=breaks[[i]][2:6],
                 count=tabulate(classes[[i]],nbins=5),color=colors5,legend=labels[[i]])))
  }
  if (6 %in% figures) {
    curves <- utils::read.csv(paths$`6`[1],check.names=FALSE,stringsAsFactors=FALSE)
    bw <- utils::read.csv(paths$`6`[2],check.names=FALSE,stringsAsFactors=FALSE)
    stopifnot(all(bw$search_max==200),identical(as.integer(bw$chosen),c(87L,91L)))
    plots <- lapply(seq_len(nrow(bw)),function(i) {
      z <- curves[curves$panel==bw$panel[i],]; a <- z[z$variable=="A",]; b <- z[z$variable=="B",]
      ggplot() + g$geom_line(data=a,aes(Bandwidth,Bias,color="Variable A"),linewidth=.55) +
        g$geom_line(data=a,aes(Bandwidth,Variance,color="Variable A"),linetype="dashed",linewidth=.55) +
        g$geom_line(data=b,aes(Bandwidth,Bias,color="Variable B"),linewidth=.55) +
        g$geom_line(data=b,aes(Bandwidth,Variance,color="Variable B"),linetype="dashed",linewidth=.55) +
        annotate("rect",xmin=min(bw$k_A[i],bw$k_B[i]),xmax=max(bw$k_A[i],bw$k_B[i]),
                 ymin=0,ymax=1,fill="grey70",alpha=.3) +
        g$geom_vline(xintercept=bw$square_root[i],linewidth=.35,color="grey40") +
        g$geom_vline(xintercept=c(bw$k_A[i],bw$k_B[i]),color="red",linewidth=.2,alpha=.6) +
        g$geom_vline(xintercept=bw$chosen[i],color="purple",linewidth=.55) +
        annotate("text",x=bw$square_root[i]+1.5,y=1.02,label=paste0("sqrt(n)=",bw$square_root[i]),
                 size=2.4,color="grey40",hjust=0) +
        g$scale_color_manual(name="Variables",values=c("Variable A"="forestgreen","Variable B"="royalblue")) +
        g$labs(x="Bandwidth (k-nearest neighbors)",y="Standardized bias/variance",tag=bw$panel[i],
          subtitle=paste0("Range [",bw$k_A[i],", ",bw$k_B[i],"]; Chosen = ",bw$chosen[i])) +
        g$theme_minimal(base_size=8) + theme(plot.tag=element_text(face="plain",size=10))
    })
    p <- patchwork::wrap_plots(plots,nrow=1,guides="collect") & theme(legend.position="bottom")
    p <- p + patchwork::plot_annotation(caption="* Solid line: Bias / Dashed line: Variance",
      theme=theme(plot.caption=element_text(size=7,hjust=.5)))
    save_plot(p,6,3.4); audit$Figure6 <- bw
  }
  if (7 %in% figures) {
    x <- lapply(paths$`7`,function(f)sf::st_transform(readRDS(f),toronto_crs))
    if (any(!vapply(x,function(z)inherits(z,"sf") && nrow(z)==279L &&
                        all(c("ADAUID","SSIM","q_value") %in% names(z)),logical(1)))) stop("Invalid Figure 7 result data.")
    breaks <- list(c(-Inf,.810055722167976,.843590429700674,.883828895813985,.901766020412953,Inf),
                   c(-Inf,-.490662829749321,-.328723313386720,-.172218333963502,-.082529548642821,Inf))
    labels <- list(c("0.71 - 0.81","0.81 - 0.84","0.84 - 0.88","0.88 - 0.9","0.9 - 0.94"),
                   c("-0.7 - -0.49","-0.49 - -0.33","-0.33 - -0.17","-0.17 - -0.08","-0.08 - 0.05"))
    qbreaks <- c(-Inf,.01,.05,.1,.5,Inf)
    qlabels <- c("0 - 0.01","0.01 - 0.05","0.05 - 0.1","0.1 - 0.5","0.5 - 1")
    plots <- list(); audit_rows <- list()
    for(i in 1:2) for(j in 1:2) {
      v <- if(j==1) x[[i]]$SSIM else x[[i]]$q_value
      cls <- cut(v,if(j==1)breaks[[i]] else qbreaks,include.lowest=TRUE,labels=FALSE)
      cols <- if(j==1)colors5 else qcolors; labs <- if(j==1)labels[[i]] else qlabels
      tag <- paste0(letters[(i-1)*2+j],")")
      p <- toronto_map(x[[i]],cls,cols,tag,xlim=c(1316500,1395000),ylim=c(549000,594000))
      title <- paste0("Scenario ",i,if(j==1)" SSIM" else " q-value")
      p <- p + annotate("text",x=1372000,y=566000,label=title,size=2.6,hjust=0)
      for(k in 1:5) p <- p + annotate("rect",xmin=1372000,xmax=1375000,
        ymin=566000-k*2100-650,ymax=566000-k*2100+650,fill=rev(cols)[k],colour="grey50",linewidth=.1) +
        annotate("text",x=1375800,y=566000-k*2100,label=rev(labs)[k],size=2.6,hjust=0)
      if(i==2 && j==1) p <- furniture(p,1326000,550000,10000,1320000,4500,2.4)
      plots[[length(plots)+1L]] <- p
      audit_rows[[length(audit_rows)+1L]] <- data.frame(panel=tag,ADAUID=x[[i]]$ADAUID,value=v,class=cls,color=cols[cls])
    }
    save_plot(patchwork::wrap_plots(plots,ncol=2),7,6.8*210/297)
    audit$Figure7 <- do.call(rbind,audit_rows)
    audit$Figure7_significance <- data.frame(scenario=1:2,n=279,
      q_lt_01=vapply(x,function(z)sum(z$q_value<.01),integer(1)),
      source=c("cimd_msdi_b87_local.rds","cimd_commute_b91_local.rds"))
  }
  for (name in names(audit)) utils::write.csv(audit[[name]],file.path(output_dir,paste0(name,"_render_audit.csv")),row.names=FALSE)
  writeLines(c(paste0("Portable R-generated manuscript figures: ",colourspace," TIFF, LZW compression."),
    paste("Requested resolution:",dpi,"dpi; width: 6.8 inches."),
    "Scientific values/classes are regenerated inputs; typography and layout are portable equivalents.",
    "Default CMYK conversion uses ImageMagick without an assigned ICC print profile; RGB is optional.",
    "Separate journal print masters preserve their original ICC profile and are not claimed to be pixel-identical.",
    "Figure 3 uses the packaged simplified BC boundary for display only; raster masks are unchanged.",
    "Figure 3 display uses nearest-neighbor EPSG:3005 reprojection; Figure 4 uses bilinear display reprojection.",
    "Figure 5 uses exact type-7 quintiles; printed legend labels are rounded.",
    "Figure 7 uses saved full-run fixed-seed q-values and the original unrounded SSIM class thresholds."),
    file.path(output_dir,"figure_render_notes.txt"))
  invisible(targets)
}