R/plot_zg_scatter.R

In dendRoAnalyst: A Tool for Processing and Analyzing Dendrometer Data

# Internal scatter plot: each row in one selected table is one observation.
.dm_zg_scatter <- function(x, source, xvar, yvar, fit, Year, DOY, point_size) {
  dat <- x[[source]]
  if (!is.data.frame(dat)) stop(source, " is unavailable. Re-run phase.zg with drought_indices = TRUE for drought data.", call. = FALSE)
  valid <- function(v) is.character(v) && length(v) == 1L && v %in% names(dat) && is.numeric(dat[[v]])
  if (!valid(xvar) || !valid(yvar)) stop("xvar and yvar must name numeric columns in ", source, call. = FALSE)
  dat <- as.data.frame(dat)
  if (!is.null(Year) || !is.null(DOY)) {
    time_col <- intersect(c("Start", "DATE", "Date", "date", "TIME", "Time"), names(dat))
    if (!length(time_col)) stop("No date column available for filtering.", call. = FALSE)
    time <- as.POSIXct(dat[[time_col[1]]], tz = "UTC")
    keep <- !is.na(time)
    if (!is.null(Year)) keep <- keep & as.integer(format(time, "%Y")) %in% Year
    if (!is.null(DOY)) {
      if (length(DOY) != 2L || any(!is.finite(DOY)) || any(DOY < 1 | DOY > 366)) stop("DOY must contain two values between 1 and 366.", call. = FALSE)
      day <- as.integer(format(time, "%j"))
      keep <- keep & if (DOY[1] <= DOY[2]) day >= DOY[1] & day <= DOY[2] else day >= DOY[1] | day <= DOY[2]
    }
    dat <- dat[keep, , drop = FALSE]
  }
  omitted <- sum(!is.finite(dat[[xvar]]) | !is.finite(dat[[yvar]]))
  dat <- dat[is.finite(dat[[xvar]]) & is.finite(dat[[yvar]]), , drop = FALSE]
  if (!nrow(dat)) stop("No finite paired observations remain.", call. = FALSE)
  dat$.scatter_x <- dat[[xvar]]
  dat$.scatter_y <- dat[[yvar]]
  p <- ggplot2::ggplot(dat, ggplot2::aes(x = .data$.scatter_x, y = .data$.scatter_y)) +
    ggplot2::geom_point(size = point_size, alpha = 0.65, colour = "#795238") +
    ggplot2::theme_minimal() + ggplot2::labs(x = xvar, y = yvar, title = paste(source, "relationship"),
      subtitle = paste(nrow(dat), "paired observations;", omitted, "non-finite pairs excluded"))
  model <- NULL
  statistics <- data.frame(n = nrow(dat), excluded = omitted)
  if (fit == "linear") {
    if (nrow(dat) < 3L || length(unique(dat$.scatter_x)) < 2L || length(unique(dat$.scatter_y)) < 2L)
      stop("Linear fit requires at least three pairs and variation in both variables.", call. = FALSE)
    model <- stats::lm(.scatter_y ~ .scatter_x, data = dat)
    sm <- summary(model)
    statistics$intercept <- unname(stats::coef(model)[1])
    statistics$slope <- unname(stats::coef(model)[2])
    statistics$r_squared <- sm$r.squared
    statistics$p_value <- sm$coefficients[2, 4]
    p <- p + ggplot2::geom_smooth(method = "lm", formula = y ~ x, se = TRUE, colour = "#795238") +
      ggplot2::labs(subtitle = sprintf("n = %d; slope = %.4g; R-squared = %.3f; p = %.3g", nrow(dat), statistics$slope, statistics$r_squared, statistics$p_value),
        caption = "Pooled linear fit with 95% confidence band. Inference assumes independent observations; repeated events/days may violate this.")
  }
  attr(p, "scatter_stats") <- statistics
  attr(p, "scatter_model") <- model
  p
}