Update vignette with improved hexagon utility functions and workflow documentation

Author: MaxMLang

vignettes/hexsmoothR-complete-guide.Rmd

@@ -94,7 +94,7 @@ Spatial topology determines smoothing algorithm behavior. Each hexagon receives
 
 ## Step 1: Library Loading and Data Preparation
 
-```{r load-libraries}
+```{r load-libraries, eval = TRUE}
 library(hexsmoothR)
 library(sf)
 library(terra)
@@ -151,11 +151,9 @@ if (file.exists(sentinel_file)) {
   true_coastal_vegetation[true_coastal_vegetation > 0.8] <- 0.8
 }
 
-
 # Use original data for hexagonal smoothing demonstration
 observed_coastal_vegetation <- true_coastal_vegetation
 
-
 # Data summary
 if (file.exists(sentinel_file)) {
   cat("Sentinel-2 L2A NDVI data loaded\n")
@@ -168,7 +166,6 @@ cat("CRS:", crs(observed_coastal_vegetation, proj = TRUE), "\n")
 missing_pct <- round(global(is.na(observed_coastal_vegetation), "mean")$mean * 100, 1)
 cat("Missing values:", missing_pct, "%\n")
 
-
 # Create color palette for NDVI visualization
 plot_colors <- colorRampPalette(c("darkred", "red", "orange", "yellow", "green", "darkgreen"))(100)
 cat("Color palette: 100 colors (red to green)\n")
@@ -192,7 +189,7 @@ if (file.exists(sentinel_file)) {
   plot_title_2 <- "Observed Coastal Vegetation Data (NDVI-like)\nSimulated data"
 }
 
-plot(observed_coastal_vegetation$sentinel2_ndvi_2,
+plot(observed_coastal_vegetation$sentinel2_ndvi_3,
      main = plot_title_2,
      col = plot_colors,
      axes = FALSE)
@@ -218,7 +215,7 @@ Use projected CRS (UTM) for accurate measurements.
 
 Calculate cell size to match target hexagon count and raster area.
 
-```{r create-grid}
+```{r create-grid, eval = TRUE}
 # Create study area polygon from raster extent
 study_area_wgs <- st_sf(geometry = st_sfc(
   st_polygon(list(matrix(
@@ -272,7 +269,7 @@ plot(st_geometry(hex_grid_wgs84), main = paste("Hexagonal Grid -", nrow(hex_grid
 
 ### Helper Function for Optimal Cell Size
 
-```{r optimal-cell-size, eval =FALSE}
+```{r optimal-cell-size, eval = TRUE}
 # Package suggests optimal cell size for target hexagon count
 target_cells <- 5000
 optimal_cell_size <- find_hex_cell_size_for_target_cells(
@@ -302,7 +299,7 @@ hex_grid <- existing_grid
 
 Extract raster values into hexagon cells using existing grid.
 
-```{r extract-raster}
+```{r extract-raster, eval = TRUE}
 # Extract raster values into hexagon cells
 # Function handles CRS transformations and pixel aggregation
 
@@ -442,28 +439,26 @@ extracted_auto <- extract_raster_data(
 
 ### Standard 2-Order Topology (Backward Compatible)
 
-```{r configure-topology, eval = FALSE}
+```{r configure-topology, eval = TRUE}
 # Compute spatial topology for smoothing (default: 2 orders)
-# Note: This example assumes you have a hex_grid object from previous steps
-# topology <- compute_topology(hex_grid)
+topology <- compute_topology(hex_grid)
 
-# Example output structure (not executable in vignette):
-cat("Topology computed: 100 cells\n")
-cat("Average distance: 15000 m\n")
-cat("Sigma (bandwidth): 7500 m\n")
+cat("Topology computed:", length(topology$neighbors[[1]]), "cells\n")
+cat("Average distance:", round(topology$avg_distance, 0), "m\n")
+cat("Sigma (bandwidth):", round(topology$sigma, 0), "m\n")
 
 cat("Gaussian weights:\n")
-cat("Center: 0.667\n")
-cat("First-order: 0.222\n") 
-cat("Second-order: 0.111\n")
+cat("Center:", round(topology$weights$center_weight, 3), "\n")
+cat("First-order:", round(topology$weights$neighbor_weights[1], 3), "\n") 
+cat("Second-order:", round(topology$weights$neighbor_weights[2], 3), "\n")
 
 # Weights applied to individual neighbors, not averages
 cat("Weights applied to each neighbor individually\n")
 
 # Neighbor examples
-cat("Cell 1: 6 1st, 12 2nd neighbors\n")
-cat("Cell 2: 6 1st, 12 2nd neighbors\n")
-cat("Cell 3: 6 1st, 12 2nd neighbors\n")
+cat("Cell 1:", length(topology$neighbors[[1]][[1]]), "1st,", length(topology$neighbors[[2]][[1]]), "2nd neighbors\n")
+cat("Cell 2:", length(topology$neighbors[[1]][[2]]), "1st,", length(topology$neighbors[[2]][[2]]), "2nd neighbors\n")
+cat("Cell 3:", length(topology$neighbors[[1]][[3]]), "1st,", length(topology$neighbors[[2]][[3]]), "2nd neighbors\n")
 
 # Note: Visualization code would go here but requires actual grid objects
 cat("Neighborhood visualization would show:\n")
@@ -474,7 +469,7 @@ cat("- 2nd order neighbors (yellow)\n")
 
 ### Advanced: N-Order Topology for Enhanced Smoothing
 
-```{r configure-n-order-topology, eval = FALSE}
+```{r configure-n-order-topology, eval = TRUE}
 # Compute topology with 3 neighbor orders for more extensive smoothing
 topology_3 <- compute_topology(hex_grid, neighbor_orders = 3)
 
@@ -505,7 +500,7 @@ for (i in 1:min(3, length(topology_3$neighbors[[1]]))) {
 topology_custom <- compute_topology(
   hex_grid, 
   neighbor_orders = 4,
-  neighbor_weights = list(0.5, 0.3, 0.15, 0.05)  # Custom decay pattern
+  neighbor_weights_param = list(0.5, 0.3, 0.15, 0.05)  # Custom decay pattern
 )
 
 cat("Custom 4-order topology with manual weights\n")
@@ -528,28 +523,31 @@ Higher orders provide more extensive smoothing but increase computation time and
 
 ### Standard 2-Order Smoothing (Backward Compatible)
 
-```{r apply-smoothing, eval = FALSE}
+```{r apply-smoothing, eval = TRUE}
 # Apply spatial smoothing using C++-optimized function (2 orders)
-# Note: This example assumes you have topology and extracted_values from previous steps
-# smoothing_results <- smooth_variables(
-#   variable_values = list(coastal_vegetation = extracted_values),
-#   neighbors = topology$neighbors,
-#   weights = topology$weights,
-#   var_names = c("coastal_vegetation")
-# )
+smoothing_results <- smooth_variables(
+  variable_values = list(coastal_vegetation = extracted_values),
+  neighbors = topology$neighbors,
+  weights = topology$weights,
+  var_names = c("coastal_vegetation")
+)
 
 cat("Spatial smoothing complete\n")
-cat("Results: raw, neighbors_1st, neighbors_2nd, weighted_combined\n")
+cat("Results:", paste(names(smoothing_results$coastal_vegetation), collapse = ", "), "\n")
 
 # Example smoothing effectiveness analysis:
 cat("Smoothing effectiveness:\n")
-cat("Original - Mean: 0.456 SD: 0.234\n")
+cat("Original - Mean:", round(mean(extracted_values, na.rm = TRUE), 3), "SD:", round(sd(extracted_values, na.rm = TRUE), 3), "\n")
 
 cat("Smoothed coastal vegetation values:\n")  
-cat("  Mean: 0.478\n")
-cat("  SD: 0.189\n")
+cat("  Mean:", round(mean(smoothing_results$coastal_vegetation$weighted_combined, na.rm = TRUE), 3), "\n")
+cat("  SD:", round(sd(smoothing_results$coastal_vegetation$weighted_combined, na.rm = TRUE), 3), "\n")
+
+noise_reduction <- (sd(extracted_values, na.rm = TRUE) - 
+                   sd(smoothing_results$coastal_vegetation$weighted_combined, na.rm = TRUE)) / 
+                   sd(extracted_values, na.rm = TRUE) * 100
 
-cat("Coastal vegetation noise reduction: 19.2%\n")
+cat("Coastal vegetation noise reduction:", round(noise_reduction, 1), "%\n")
 cat("(Lower standard deviation = more consistent coastal vegetation patterns)\n")
 
 # The smoothing results contain:
@@ -563,7 +561,7 @@ cat("\n★ weighted_combined is typically the best result for noise reduction\n"
 
 ### Advanced: N-Order Smoothing for Enhanced Results
 
-```{r apply-n-order-smoothing, eval = FALSE}
+```{r apply-n-order-smoothing, eval = TRUE}
 # Apply 3-order smoothing using the new N-order system
 smoothing_results_3 <- smooth_variables(
   variable_values = list(coastal_vegetation = extracted_values),
@@ -595,9 +593,9 @@ cat("- neighbors_2nd:", round(mean(smoothing_results_3$coastal_vegetation$neighb
 cat("- neighbors_3rd:", round(mean(smoothing_results_3$coastal_vegetation$neighbors_3rd, na.rm = TRUE), 3), "\n")
 
 # Calculate additional noise reduction from 3-order smoothing
-noise_reduction_3 <- (sd(raw_values, na.rm = TRUE) - 
+noise_reduction_3 <- (sd(extracted_values, na.rm = TRUE) - 
                      sd(smoothing_results_3$coastal_vegetation$weighted_combined, na.rm = TRUE)) / 
-                     sd(raw_values, na.rm = TRUE) * 100
+                     sd(extracted_values, na.rm = TRUE) * 100
 
 cat("3-Order noise reduction:", round(noise_reduction_3, 1), "%\n")
 cat("Additional reduction over 2-order:", round(noise_reduction_3 - noise_reduction, 1), "%\n")
@@ -612,20 +610,19 @@ cat("Additional reduction over 2-order:", round(noise_reduction_3 - noise_reduct
 
 ## Step 6: Combine Results and Analyze
 
-```{r combine-results, eval = FALSE}
+```{r combine-results, eval = TRUE}
 # Combine the hexagonal grid with smoothing results
-# Note: This example assumes you have hex_grid and smoothing_results from previous steps
-# hex_grid_with_results <- hex_grid
-# hex_grid_with_results$coastal_vegetation_raw <- smoothing_results$coastal_vegetation$raw
-# hex_grid_with_results$coastal_vegetation_neighbors_1st <- smoothing_results$coastal_vegetation$neighbors_1st
-# hex_grid_with_results$coastal_vegetation_neighbors_2nd <- smoothing_results$coastal_vegetation$neighbors_2nd
-# hex_grid_with_results$coastal_vegetation_weighted_combined <- smoothing_results$coastal_vegetation$weighted_combined
+hex_grid_with_results <- hex_grid
+hex_grid_with_results$coastal_vegetation_raw <- smoothing_results$coastal_vegetation$raw
+hex_grid_with_results$coastal_vegetation_neighbors_1st <- smoothing_results$coastal_vegetation$neighbors_1st
+hex_grid_with_results$coastal_vegetation_neighbors_2nd <- smoothing_results$coastal_vegetation$neighbors_2nd
+hex_grid_with_results$coastal_vegetation_weighted_combined <- smoothing_results$coastal_vegetation$weighted_combined
 
 # Transform to WGS84 for plotting
-# hex_grid_with_results_wgs84 <- st_transform(hex_grid_with_results, crs = 4326)
+hex_grid_with_results_wgs84 <- st_transform(hex_grid_with_results, crs = 4326)
 
 cat("Results combined with grid successfully\n")
-cat("Available columns: grid_id, grid_index, coastal_vegetation_raw, coastal_vegetation_neighbors_1st, coastal_vegetation_neighbors_2nd, coastal_vegetation_weighted_combined\n")
+cat("Available columns:", paste(names(hex_grid_with_results), collapse = ", "), "\n")
 
 # Analyze the smoothing effect
 cat("\n=== SMOOTHING ANALYSIS ===\n")
@@ -636,37 +633,53 @@ cat("- neighbors_2nd: Mean of second-order neighbors only\n")
 cat("- weighted_combined: Weighted average (center + all neighbors, weights applied to each neighbor)\n\n")
 
 # Example analysis output:
-cat("Original values - Count: 100 Mean: 0.4562 SD: 0.2341\n")
-cat("Smoothed values - Count: 100 Mean: 0.4789 SD: 0.1893\n")
-cat("Variance reduction: 19.1%\n")
+cat("Original values - Count:", length(extracted_values), "Mean:", round(mean(extracted_values, na.rm = TRUE), 4), "SD:", round(sd(extracted_values, na.rm = TRUE), 4), "\n")
+cat("Smoothed values - Count:", length(smoothing_results$coastal_vegetation$weighted_combined), "Mean:", round(mean(smoothing_results$coastal_vegetation$weighted_combined, na.rm = TRUE), 4), "SD:", round(sd(smoothing_results$coastal_vegetation$weighted_combined, na.rm = TRUE), 4), "\n")
+cat("Variance reduction:", round(noise_reduction, 1), "%\n")
 ```
 
 ## Step 7: Visualize the Improvement
 
-```{r visualize-results, eval = FALSE}
+```{r visualize-results, eval = TRUE}
 # Create spectacular before/after visualization showing the power of spatial smoothing
-# Note: This example assumes you have hex_grid_wgs84 and smoothing_results from previous steps
-
 cat("=== CREATING BEFORE/AFTER VISUALIZATION ===\n")
 
-# Example visualization workflow (not executable in vignette):
-cat("Visualization would show:\n")
-cat("1. Original noisy coastal vegetation raster (NDVI-like)\n")
-cat("2. True underlying coastal vegetation pattern (for reference)\n")
-cat("3. Spatially smoothed coastal vegetation result\n")
-cat("4. Color legend with viridis palette\n")
+# Create visualization plots
+par(mfrow = c(2, 2))
 
-cat("Progression of smoothing would display:\n")
-cat("- Raw data\n")
-cat("- Smoothed (1st order)\n")
-cat("- Smoothed (All orders)\n")
+# Plot 1: Original extracted values
+plot(st_geometry(hex_grid_wgs84), 
+     col = plot_colors[cut(extracted_values, 100)],
+     main = "Original Extracted Values",
+     border = "white", lwd = 0.01)
 
-cat("All visualization plots would be created successfully!\n")
-```
+# Plot 2: Smoothed values (1st order)
+plot(st_geometry(hex_grid_wgs84), 
+     col = plot_colors[cut(smoothing_results$coastal_vegetation$neighbors_1st, 100)],
+     main = "Smoothed (1st Order Neighbors)",
+     border = "white", lwd = 0.01)
 
+# Plot 3: Smoothed values (2nd order)
+plot(st_geometry(hex_grid_wgs84), 
+     col = plot_colors[cut(smoothing_results$coastal_vegetation$neighbors_2nd, 100)],
+     main = "Smoothed (2nd Order Neighbors)",
+     border = "white", lwd = 0.01)
 
+# Plot 4: Final weighted combined result
+plot(st_geometry(hex_grid_wgs84), 
+     col = plot_colors[cut(smoothing_results$coastal_vegetation$weighted_combined, 100)],
+     main = "Final Weighted Combined Result",
+     border = "white", lwd = 0.01)
 
-## Step 9: Working with Multiple Variables
+cat("All visualization plots created successfully!\n")
+cat("Progression shows:\n")
+cat("- Top left: Raw extracted data\n")
+cat("- Top right: 1st order neighbor smoothing\n")
+cat("- Bottom left: 2nd order neighbor smoothing\n")
+cat("- Bottom right: Final weighted combined result\n")
+```
+
+## Step 8: Working with Multiple Variables
 
 ```{r multiple-variables, eval = FALSE}
 # You can work with multiple raster files and variables
@@ -701,9 +714,9 @@ smoothing_multi <- smooth_variables(
 cat("Variables processed:", paste(names(smoothing_multi), collapse = ", "), "\n")
 ```
 
-## Step 10: Hexagon Measurement Utilities
+## Step 9: Hexagon Measurement Utilities
 
-```{r hexagon-utilities}
+```{r hexagon-utilities, eval = TRUE}
 # The package provides helper functions for hexagon measurements
 # IMPORTANT: cell_size_flat is in METERS (UTM coordinates)
 cell_size_flat <- 20000  # 20km flat-to-flat distance in meters
@@ -759,16 +772,19 @@ This vignette demonstrates the complete hexsmoothR workflow:
 ### Understanding Units and CRS
 
 **Coordinate Reference Systems:**
+
 - **UTM coordinates**: Use meters for cell sizes (e.g., 20000 for 20km hexagons)
 - **WGS84 coordinates**: Use degrees for cell sizes (e.g., 0.1 for ~11km hexagons)
 - **Functions automatically convert** between CRS when needed
 
 **Cell Size Units:**
+
 - **`cell_size` parameter**: Always in the units of your study area's CRS
 - **UTM zones**: cell_size in meters (recommended for real-world analysis)
 - **Geographic coordinates**: cell_size in degrees (use with caution for large areas)
 
 **Smoothing Weights:**
+
 - **Weights apply to each individual neighbor**, not to neighbor means
 - **Center weight**: Influence of the cell itself
 - **First-order weight**: Influence of each immediate neighbor