library(terra)
library(sf)
library(dplyr)

# 1. SET PATHS AND PARAMETERS -------------------------------------------------
input_raster_path <- "C:/Project/Adjusting Carbon maps/Layers/Change_2023/Change_ALL_2023/Change_All_2023.tif"
shapefile_path <- "C:/Project/Adjusting Carbon maps/Layers/Stock_2023/Sweden_Län/Sweden_Lan.shp"
output_dir <- "C:/Project/Adjusting Carbon maps/Layers/Change_2023/Change_ALL_2023/Change_All_2023_Lan"

# Create output directory if it doesn't exist
if (!dir.exists(output_dir)) {
  dir.create(output_dir, recursive = TRUE)
}

# 2. LOAD DATA ----------------------------------------------------------------
# Load raster
raster_data <- rast(input_raster_path)

# Load shapefile
vector_data <- st_read(shapefile_path, quiet = TRUE)

# Check if OBJEKT_ID column exists
if (!"OBJEKT_ID" %in% names(vector_data)) {
  stop("The shapefile does not contain an 'OBJEKT_ID' column. Please check your data.")
}

# Check coordinate systems
cat("Raster CRS:", crs(raster_data), "\n")
cat("Vector CRS:", st_crs(vector_data)$wkt, "\n")

# Reproject vector to match raster if needed
if (!identical(crs(raster_data), st_crs(vector_data)$wkt)) {
  vector_data <- st_transform(vector_data, crs(raster_data))
  cat("Reprojected vector to match raster CRS\n")
}

# 3. CLIP RASTER FOR EACH FEATURE ---------------------------------------------
# Get raster base name for output naming
raster_basename <- tools::file_path_sans_ext(basename(input_raster_path))

# Loop through each feature
for (i in 1:nrow(vector_data)) {
  
  # Get individual feature
  single_feature <- vector_data[i, ]
  
  # Get OBJEKT_ID value
  objekt_id <- as.character(single_feature$OBJEKT_ID)
  
  # Check if OBJEKT_ID is valid
  if (is.na(objekt_id) || objekt_id == "") {
    warning(paste("Feature", i, "has missing or empty OBJEKT_ID. Using index instead."))
    objekt_id <- paste0("feature_", i)
  }
  
  # Clean the OBJEKT_ID for filename (remove invalid characters)
  clean_id <- gsub("[^[:alnum:]_\\-]", "_", objekt_id)
  
  cat("Processing OBJEKT_ID:", objekt_id, "\n")
  
  # Convert sf to SpatVector for terra
  feature_vect <- vect(single_feature)
  
  # Crop and mask the raster
  cropped_raster <- crop(raster_data, feature_vect)
  masked_raster <- mask(cropped_raster, feature_vect)
  
  # Optional: Check if the clipped raster has data
  if (global(masked_raster, "notNA")[[1]] == 0) {
    warning(paste("OBJEKT_ID", objekt_id, "does not intersect with the raster. Skipping."))
    next
  }
  
  # 4. CREATE OUTPUT FILENAME WITH OBJEKT_ID ----------------------------------
  # Create output filename with format: rastername_OBJEKT_ID.tif
  output_filename <- file.path(output_dir, 
                               paste0(raster_basename, "_", clean_id, ".tif"))
  
  # Alternative naming options (uncomment if preferred):
  # 1. Only OBJEKT_ID: output_filename <- file.path(output_dir, paste0(clean_id, ".tif"))
  # 2. With prefix: output_filename <- file.path(output_dir, paste0("clipped_", clean_id, ".tif"))
  
  # 5. SAVE OUTPUT -------------------------------------------------------------
  # Write raster to file
  writeRaster(masked_raster, 
              filename = output_filename,
              overwrite = TRUE,
              # Optional compression for smaller files
              gdal = c("COMPRESS=DEFLATE", "PREDICTOR=2", "ZLEVEL=9"))
  
  cat("Saved:", basename(output_filename), "\n")
}

cat("\nProcessing complete!\n")
cat("Output saved to:", output_dir, "\n")
cat("Total features processed:", nrow(vector_data), "\n")