plot.R

#!/usr/bin/env Rscript

# plot.R - Plot a single GeoTIFF on Santa Barbara coastline background

# Load required libraries
library(terra)
library(sf)
library(ggplot2)
library(maps)
library(terrainr)
# Get command line arguments
args <- commandArgs(trailingOnly = TRUE)

# Check if GeoTIFF path is provided
if (length(args) == 0) {
  stop("Usage: Rscript plot.R <path_to_geotiff>")
}

geotiff_path <- args[1]

# Check if file exists
if (!file.exists(geotiff_path)) {
  stop("GeoTIFF file not found: ", geotiff_path)
}

# Extract base filename for output
base_name <- tools::file_path_sans_ext(basename(geotiff_path))

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

# Load the GeoTIFF
cat("Loading GeoTIFF:", geotiff_path, "\n")
raster_data <- rast(geotiff_path)

# Check if raster needs CRS setting
if (is.na(crs(raster_data)) || crs(raster_data) == "") {
  cat("Setting CRS to UTM Zone 11N...\n")
  crs(raster_data) <- "EPSG:32611"
}

# Project to WGS84 for mapping with coastline data
cat("Projecting raster to WGS84...\n")
wgs84_crs <- "EPSG:4326"
projected_raster <- project(raster_data, wgs84_crs, method = "near")

# Convert raster to data frame for ggplot
cat("Converting raster to data frame...\n")
raster_df <- as.data.frame(projected_raster, xy = TRUE)

# Check the number of bands and set appropriate column names
num_bands <- ncol(raster_df) - 2  # subtract x and y columns
cat("Number of bands:", num_bands, "\n")

if (num_bands == 3) {
  # RGB data
  colnames(raster_df) <- c("x", "y", "r", "g", "b")
  is_rgb <- TRUE
} else if (num_bands == 1) {
  # Single band data
  colnames(raster_df) <- c("x", "y", "value")
  is_rgb <- FALSE
} else {
  # Multiple bands but not RGB
  colnames(raster_df)[1:2] <- c("x", "y")
  colnames(raster_df)[3] <- "value"
  is_rgb <- FALSE
  cat("Warning: Multiple bands detected, using first band for visualization\n")
}

# Remove NA values
if (is_rgb) {
  raster_df <- raster_df[!is.na(raster_df$r) & !is.na(raster_df$g) & !is.na(raster_df$b), ]
} else {
  raster_df <- raster_df[!is.na(raster_df$value), ]
}

# Get Santa Barbara area coastline data
cat("Getting Santa Barbara coastline data...\n")

# Get California state outline and county data for detailed coastline
california <- map_data("state", region = "california")
california_county <- map_data("county", region = "california")

# Filter to Santa Barbara area counties for more detailed coastline
sb_counties <- c("santa barbara", "ventura", "kern", "san luis obispo")
santa_barbara_region <- california_county[california_county$subregion %in% sb_counties, ]

# Calculate plot extent based on raster extent with buffer
raster_extent <- ext(projected_raster)
buffer <- 0.02  # degrees
xlim <- c(raster_extent[1] - buffer, raster_extent[2] + buffer)
ylim <- c(raster_extent[3] - buffer, raster_extent[4] + buffer)

# Create the plot
cat("Creating plot...\n")
p <- ggplot() +
  # Add California outline for context
  geom_polygon(data = california, 
               aes(x = long, y = lat, group = group),
               fill = "lightgray", 
               color = "darkgray", 
               linewidth = 0.3,
               alpha = 0.5) +
  
  # Add detailed coastline from county boundaries
  geom_polygon(data = santa_barbara_region,
               aes(x = long, y = lat, group = group),
               fill = "lightblue",
               color = "darkblue",
               linewidth = 0.8,
               alpha = 0.3) +
  
  # Add raster data (conditional based on data type)
  {if (is_rgb) {
    geom_spatial_rgb(data = raster_df, 
                     mapping = aes(x=x, y=y, r=r, g=g, b=b),
                     alpha = 0.8)
  } else {
    geom_tile(data = raster_df, 
              aes(x = x, y = y, fill = value),
              alpha = 0.8)
  }} +
  
  # Set color scale for single band raster values
  {if (!is_rgb) {
    scale_fill_gradientn(colors = terrain.colors(100),
                         name = "Value",
                         na.value = "transparent")
  }} +
  
  # Set coordinate system and zoom to raster extent
  coord_sf(crs = st_crs(4326), 
           xlim = xlim, 
           ylim = ylim) +
  
  # Add titles and labels
  labs(title = paste("GeoTIFF Plot:", base_name),
       subtitle = "Santa Barbara Coastline Region",
       x = "Longitude",
       y = "Latitude") +
  
  # Theme
  theme_minimal() +
  theme(panel.grid = element_line(color = "white", linewidth = 0.2),
        panel.background = element_rect(fill = "aliceblue"),
        legend.position = "right",
        plot.title = element_text(face = "bold", size = 14),
        plot.subtitle = element_text(size = 12))

# Save the plot
output_file <- file.path("images", paste0(base_name, ".png"))
cat("Saving plot to:", output_file, "\n")
ggsave(output_file, plot = p, width = 12, height = 8, dpi = 300)

# Print summary
cat("\nPlot creation complete!\n")
cat("Input file:", geotiff_path, "\n")
cat("Output file:", output_file, "\n")
cat("Raster dimensions:", dim(raster_data), "\n")
cat("Raster extent (WGS84):", as.vector(raster_extent), "\n")