crop.R

# crop.R: scale-down and crop a geotiff from the Dibblee map collection for
# combining with other maps. 
library(terra)
library(sf)

# Get command line arguments
args <- commandArgs(trailingOnly = TRUE)

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

geotiff_path <- args[1]
txt_path <- paste0(geotiff_path, ".txt")

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

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

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

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

# Rectify the raster data (as suggested by the warning)
raster_data <- rectify(raster_data)

# Scale down to 1/16 size by aggregating with factor of 4
cat("Scaling down to 1/16 size...\n")
scaled_raster <- aggregate(raster_data, fact = 4, fun = "mean")

# Save the scaled-down version with original compression settings
# scaled_output <- file.path("output", paste0(base_name, "_scaled.tif"))
# cat("Saving scaled GeoTIFF with original compression settings:", scaled_output, "\n")
# writeRaster(scaled_raster, scaled_output, overwrite = TRUE, 
#             datatype = "INT1U", 
#             gdal = c("COMPRESS=JPEG", "JPEG_QUALITY=85", "TILED=YES", "BLOCKXSIZE=256", "BLOCKYSIZE=256"))

# Read the coordinate points from the .txt file
cat("Reading coordinate points from text file:", txt_path, "\n")
points_data <- read.table(txt_path, header = FALSE, 
                         col.names = c("longitude", "latitude", "proj_x", "proj_y"))

# Create a polygon from the coordinate points using convex hull
cat("Creating polygon from coordinate points...\n")

# Convert projected coordinates to an sf object
coords_sf <- st_as_sf(points_data, coords = c("proj_x", "proj_y"))

# Create convex hull polygon to enclose all points
hull_polygon <- st_convex_hull(st_union(coords_sf))

# Get the CRS from the scaled raster
raster_crs <- crs(scaled_raster)
if (is.na(raster_crs) || raster_crs == "") {
  cat("Setting raster CRS to UTM Zone 11N...\n")
  crs(scaled_raster) <- "EPSG:32611"
  raster_crs <- crs(scaled_raster)
}

# Set the same CRS for the polygon
st_crs(hull_polygon) <- raster_crs

# Convert polygon to terra SpatVector for compatibility with terra functions
hull_vect <- vect(hull_polygon)

# Get the bounding box of the polygon for initial cropping
polygon_bbox <- ext(hull_vect)

cat("Polygon bounding box:\n")
cat("X range:", polygon_bbox[1], "to", polygon_bbox[2], "\n")
cat("Y range:", polygon_bbox[3], "to", polygon_bbox[4], "\n")

# First crop the raster to the polygon's bounding box
cat("Cropping scaled raster to polygon bounding box...\n")
bbox_cropped <- crop(scaled_raster, polygon_bbox)

# Now mask the raster using the polygon to make areas outside transparent
cat("Masking raster with polygon (making outside areas transparent)...\n")
cropped_raster <- mask(bbox_cropped, hull_vect)

# Save the cropped raster with original compression settings
cropped_output <- file.path("output", paste0(base_name, "_scaled_cropped.tif"))
cat("Saving polygon-masked GeoTIFF with original compression settings:", cropped_output, "\n")
writeRaster(cropped_raster, cropped_output, overwrite = TRUE, 
            datatype = "INT1U", 
            gdal = c("COMPRESS=JPEG", "JPEG_QUALITY=85", "TILED=YES", "BLOCKXSIZE=256", "BLOCKYSIZE=256"))

# Also save the polygon as a shapefile for reference
polygon_output <- file.path("output", paste0(base_name, "_polygon.shp"))
cat("Saving polygon shapefile:", polygon_output, "\n")
st_write(hull_polygon, polygon_output, quiet = TRUE)

# Print summary
cat("\nProcessing complete!\n")
cat("Original dimensions:", dim(raster_data), "\n")
cat("Scaled dimensions:", dim(scaled_raster), "\n")
cat("Cropped dimensions:", dim(cropped_raster), "\n")
cat("Size reduction factor (scaled):", prod(dim(raster_data)[1:2]) / prod(dim(scaled_raster)[1:2]), "\n")
cat("Size reduction factor (cropped):", prod(dim(raster_data)[1:2]) / prod(dim(cropped_raster)[1:2]), "\n")

cat("\nOutput files saved to:\n")
cat("-", cropped_output, "\n")
cat("-", polygon_output, "\n")