combine.R

# combine.R
# combine scaled and cropped geotiffs together, in a single plot, applying a consistent CRS.

# Load required libraries
library(terra)
library(sf)
library(ggplot2)
library(maps)
library(dplyr)
library(tidyr)
library(terrainr)

# Create images directory if it doesn't exist
if (!dir.exists("images")) {
  dir.create("images")
}

# Get all _scaled_cropped.tif files from the output folder
tif_files <- list.files("output", pattern = "_scaled_cropped\\.tif$", full.names = TRUE)

# Check if any files were found
if (length(tif_files) == 0) {
  stop("No _scaled_cropped.tif files found in the output folder")
}

cat("Found", length(tif_files), "raster files\n")

# Load all raster files
raster_list <- lapply(tif_files, rast)

# Extract file names for labels
file_labels <- basename(tif_files)
file_labels <- gsub("_scaled_cropped\\.tif$", "", file_labels)

# Set CRS to UTM Zone 11N for all rasters (appropriate for California)
utm_crs <- "EPSG:32611"

cat("Setting CRS to UTM Zone 11N for all rasters...\n")
for (i in 1:length(raster_list)) {
  crs(raster_list[[i]]) <- utm_crs
}

# Convert to WGS84 for mapping
wgs84_crs <- "EPSG:4326"

cat("Projecting all rasters to WGS84...\n")
projected_rasters <- list()
for (i in 1:length(raster_list)) {
  cat("Projecting", file_labels[i], "...\n")
  projected_rasters[[i]] <- project(raster_list[[i]], wgs84_crs, method = "near")
}

# Convert rasters to data frames for ggplot
cat("Converting rasters to data frames...\n")
raster_dfs <- list()
is_rgb_list <- logical(length(projected_rasters))

for (i in 1:length(projected_rasters)) {
  cat("Converting", file_labels[i], "...\n")
  
  # Convert raster to data frame
  df <- as.data.frame(projected_rasters[[i]], xy = TRUE)
  
  # Check the number of bands and set appropriate column names
  num_bands <- ncol(df) - 2  # subtract x and y columns
  cat("  Number of bands:", num_bands, "\n")
  
  if (num_bands == 3) {
    # RGB data
    colnames(df) <- c("x", "y", "r", "g", "b")
    is_rgb_list[i] <- TRUE
    
    # Remove NA values for RGB data
    df <- df[!is.na(df$r) & !is.na(df$g) & !is.na(df$b), ]
  } else if (num_bands == 1) {
    # Single band data
    colnames(df) <- c("x", "y", "value")
    is_rgb_list[i] <- FALSE
    
    # Remove NA values
    df <- df[!is.na(df$value), ]
  } else {
    # Multiple bands but not RGB - use first band
    colnames(df)[1:2] <- c("x", "y")
    colnames(df)[3] <- "value"
    is_rgb_list[i] <- FALSE
    
    # Remove NA values
    df <- df[!is.na(df$value), ]
    cat("  Warning: Multiple bands detected, using first band for visualization\n")
  }
  
  # Add dataset name
  df$dataset <- file_labels[i]
  
  raster_dfs[[i]] <- df
}

# No coastline background needed

# Calculate the combined extent of all bounding boxes for zooming
all_extents <- lapply(projected_rasters, ext)
combined_extent <- ext(
  min(sapply(all_extents, function(x) x[1])),  # xmin (longitude)
  max(sapply(all_extents, function(x) x[2])),  # xmax (longitude)
  min(sapply(all_extents, function(x) x[3])),  # ymin (latitude)
  max(sapply(all_extents, function(x) x[4]))   # ymax (latitude)
)

# Add some buffer around the bounding boxes for better visualization
buffer <- 0.02  # degrees
xlim <- c(combined_extent[1] - buffer, combined_extent[2] + buffer)
ylim <- c(combined_extent[3] - buffer, combined_extent[4] + buffer)

# Create the plot
cat("Creating plot with separate raster layers...\n")
p <- ggplot()

# Add each raster as a separate layer
for (i in 1:length(raster_dfs)) {
  cat("Adding raster layer for", file_labels[i], "...\n")
  
  if (is_rgb_list[i]) {
    # RGB data - use geom_spatial_rgb
    p <- p + geom_spatial_rgb(data = raster_dfs[[i]], 
                              mapping = aes(x = x, y = y, r = r, g = g, b = b))
  } else {
    # Single band data - use geom_tile
    p <- p + geom_tile(data = raster_dfs[[i]], 
                       aes(x = x, y = y, fill = value))
  }
}

# Continue building the plot
p <- p +
  # Add color scale only if we have non-RGB data
  {if (any(!is_rgb_list)) {
    scale_fill_gradientn(colors = terrain.colors(100),
                         name = "Raster\nValue",
                         na.value = "transparent")
  }} +
  
  # Set coordinate system and zoom to Santa Barbara region
  coord_sf(crs = st_crs(4326), 
           xlim = xlim, 
           ylim = ylim) +
  
  # Add titles and labels
  labs(title = "Combined Dibblee Maps -- Santa Barbara Region",
       x = "Longitude",
       y = "Latitude") +
  
  # Theme
  theme_minimal() +
  theme(panel.grid = element_line(color = "white", linewidth = 0.2),
        panel.background = element_rect(fill = "white"),
        legend.position = "bottom",
        plot.title = element_text(face = "bold", size = 14),
        plot.subtitle = element_text(size = 12))

# Save the plot
ggsave("images/combined.png", plot = p, width = 12, height = 10, dpi = 300)

# Print summary information
cat("Bounding box summary:\n")
cat("====================\n")
for (i in 1:length(projected_rasters)) {
  e <- ext(projected_rasters[[i]])
  cat("Dataset:", file_labels[i], "\n")
  cat("  Longitude:", round(e[1], 4), "to", round(e[2], 4), "\n")
  cat("  Latitude:", round(e[3], 4), "to", round(e[4], 4), "\n")
  cat("\n")
}

cat("Multiple raster layers plot saved to: images/combined.png\n")
cat("Successfully processed", length(raster_list), "raster files as separate layers\n")