variogram_random_precip.R

# Load necessary libraries
library(sf)
library(sp)  # For spatial objects
library(gstat)  # For variogram computation and modeling
library(ggplot2)  # For plotting
library(dplyr)  # For data manipulation

# Define file paths
base_dir <- normalizePath("../")  # Adjust base directory path as needed
# observation_shapefile <- file.path(base_dir, "precipitation_compute", "precipitation_points.shp")
# observation_shapefile <- file.path(base_dir, "precipitation_compute/precipitation_points", "sn_ppt_points_20150524_06_09.shp")
# observation_shapefile <- file.path(base_dir, "precipitation_compute/precipitation_points", "sn_ppt_points_20160603_01_04.shp")
# observation_shapefile <- file.path(base_dir, "precipitation_compute/precipitation_points", "gd_ppt_points_20181016_17_20.shp")
observation_shapefile <- file.path(base_dir, "precipitation_compute/precipitation_points", "ln_ppt_points_20181016_17_20.shp")

# Read the shapefile (POINTS shapefile)
observation_points <- st_read(observation_shapefile)

# Check if 'precipitat' column exists
if (!"precipitat" %in% colnames(observation_points)) {
  stop("Error: 'precipitat' column not found in the shapefile.")
}

# Remove the lower 5% of precipitation values
observation_points <- observation_points %>% filter(precipitat >= 0)
# lower_threshold <- quantile(observation_points$precipitat, 0.45, na.rm = TRUE)
lower_threshold <- quantile(observation_points$precipitat, 0.5, na.rm = TRUE) #0.5 for SN 2016
print(lower_threshold)
# lower_threshold <- max(lower_threshold, 5)
# print(lower_threshold)
filtered_points <- observation_points %>% filter(precipitat >= lower_threshold)

# Print the number of filtered points
cat("Number of points after removing lower 5%:", nrow(filtered_points), "\n")

# Number of iterations
num_iterations <- 20
range_values <- numeric(num_iterations)  # Store range values for each iteration
seed_values <- numeric(num_iterations)  # Store seed values

# Run variogram analysis 200 times with different random samples
for (i in 1:num_iterations) {
  
  # Generate a unique seed number
  seed_number <- 40 * i
  set.seed(seed_number)  # Use a different seed for each iteration
  seed_values[i] <- seed_number
  
  # Randomly select 300 points
  selected_points <- filtered_points %>% sample_n(300)
  
  cat("\nIteration", i, "- Seed:", seed_number, "- Number of selected points:", nrow(selected_points), "\n")
  
  # Check CRS and transform if necessary
  if (st_is_longlat(selected_points)) {
    selected_points <- st_transform(selected_points, crs = 32614)  # UTM Zone 14N
  }
  
  # Convert sf to sp object
  selected_points_sp <- as(selected_points, "Spatial")
  selected_points_sp$obs <- selected_points_sp$precipitat
  
  # Automatically determine cutoff distance (2/3 of max pairwise distance)
  coords <- coordinates(selected_points_sp)
  dist_matrix <- as.matrix(dist(coords))
  max_dist <- max(dist_matrix)
  cutoff_auto <- max(max_dist * (2/3))  # Ensures reasonable cutoff
  
  cat("  Cutoff distance:", cutoff_auto, "meters\n")
  
  # Compute empirical variogram
  empirical_variogram <- variogram(obs ~ 1, selected_points_sp, 
                                   cutoff = cutoff_auto, width = cutoff_auto / 25)
  
  # Check if variogram computation is valid
  if (nrow(empirical_variogram) == 0) {
    cat("  Error: Variogram computation failed in iteration", i, "\n")
    next
  }
  
  # Fit a theoretical variogram model (Exponential)
  fitted_variogram <- fit.variogram(empirical_variogram, 
                                    model = vgm(model = "Exp", nugget = 0.01), 
                                    fit.sills = TRUE)
  
  # Store the estimated range (correlation length)
  range_value <- fitted_variogram[fitted_variogram$model == "Exp", "range"]
  range_values[i] <- range_value
  
  cat("  Estimated Variogram Range:", range_value, "meters\n")
}

# Convert results to a data frame
range_df <- data.frame(Iteration = 1:num_iterations, Seed = seed_values, Range = range_values)

# Print summary statistics of the range values
cat("\nSummary Statistics of Variogram Ranges:\n")
summary_stats <- summary(range_df$Range)
print(summary_stats)

# Calculate and print the mean and median
mean_range <- mean(range_df$Range, na.rm = TRUE)
median_range <- median(range_df$Range, na.rm = TRUE)
cat("\nMean Variogram Range:", mean_range, "meters\n")
cat("Median Variogram Range:", median_range, "meters\n")

# Save the range values as a CSV
# write.csv(range_df, "variogram_ranges.csv", row.names = FALSE)

# Print the range values for all iterations
print(range_df)

ggplot(range_df, aes(x = Iteration, y = (Range * 3) / 1000)) +  # Multiply by 3 and convert meters to km
  geom_point(color = "darkgray", size = 2) +  # Plot data points
  geom_line(color = "darkgray", linetype = "dashed") +  # Connect points with a dashed line
  geom_hline(aes(yintercept = (mean_range * 3) / 1000, color = "Mean"), linetype = "dashed", size = 1) +  # Mean line (thicker)
  geom_hline(aes(yintercept = (median_range * 3) / 1000, color = "Median"), linetype = "solid", size = 1) +  # Median line (thicker)
  scale_color_manual(
    values = c("Mean" = "purple", "Median" = "maroon"),  # Define colors for legend
    guide = guide_legend(override.aes = list(fill = NA))  # Remove legend box fill
  ) +
  labs(
    # title = "Variogram Range Across 200 Random Samples",
    x = "Iterations",
    y = "Correlation Length (km)"
  ) +
  theme_minimal(base_size = 14) +  # Increase base font size
  theme(
    plot.title = element_text(size = 16, hjust = 0.5),  # Centered, bold title with larger font
    axis.title = element_text(size = 16),  # Bigger axis titles
    axis.text = element_text(size = 16),  # Bigger axis labels
    legend.text = element_text(size = 14),  # Bigger legend text
    panel.border = element_rect(color = "black", fill = NA, size = 1),  # Single black rectangular boundary
    legend.title = element_blank(),  # Remove the "Statistics" heading
    legend.position = c(0.05, 0.95),  # Move legend to the top-left inside the box
    legend.justification = c(0, 1),  # Align legend inside the top-left corner
    legend.background = element_rect(fill = "white", color = NA, size = 0),  # Transparent background
    panel.grid.major = element_line(color = "gray80"),  # Major grid lines
    panel.grid.minor = element_line(color = "gray90")  # Minor grid lines
  )

# Save the updated plot
# ggsave(
#   filename = file.path(base_dir, "Results_paper/paper_final/precip_sn_2015.pdf"),
#   height = 3.75,  # Rectangular height
#   width = 5,  # Rectangular width
#   units = "in",
#   dpi = 1000,
#   scale = 1
# )