# Convert sample data to sf object
sample_sf <- sf::st_as_sf(sample_points, 
                          coords = c("lon", "lat"), 
                          crs = 4326)

# Basic point map
basic_map <- create_spatial_map(
  spatial_data = sample_sf,
  fill_variable = "ndvi",
  map_type = "points",
  title = "Sample NDVI Points",
  verbose = TRUE
)

print(basic_map)

# NDVI-specific colors
ndvi_map <- create_spatial_map(
  spatial_data = sample_sf,
  fill_variable = "ndvi",
  color_scheme = "ndvi",
  title = "NDVI with Specialized Colors"
)

# Terrain colors for elevation data
# terrain_map <- create_spatial_map(
#   spatial_data = elevation_data,
#   color_scheme = "terrain",
#   title = "Elevation Map"
# )

# Plasma colors for general data
plasma_map <- create_spatial_map(
  spatial_data = sample_sf,
  fill_variable = "yield",
  color_scheme = "plasma",
  title = "Yield with Plasma Colors"
)

print(plasma_map)

# Create sample raster for demonstration
sample_raster <- terra::rast(nrows = 50, ncols = 50, 
                             xmin = -84, xmax = -82, 
                             ymin = 39, ymax = 41)
terra::values(sample_raster) <- runif(2500, 0.2, 0.9)
names(sample_raster) <- "NDVI"

# Fast raster plot
plot_raster_fast(
  raster_data = sample_raster,
  title = "Sample NDVI Raster",
  color_scheme = "ndvi"
)

# With custom breaks
plot_raster_fast(
  raster_data = sample_raster,
  title = "NDVI with Custom Classes",
  color_scheme = "ndvi",
  breaks = c(0, 0.3, 0.5, 0.7, 1.0)
)

# Create sample multi-band raster
red_band <- terra::rast(nrows = 30, ncols = 30, 
                        xmin = -84, xmax = -82, 
                        ymin = 39, ymax = 41)
terra::values(red_band) <- runif(900, 0.1, 0.3)

green_band <- red_band
terra::values(green_band) <- runif(900, 0.2, 0.4)

blue_band <- red_band
terra::values(blue_band) <- runif(900, 0.05, 0.15)

# Stack bands
rgb_stack <- c(red_band, green_band, blue_band)
names(rgb_stack) <- c("Red", "Green", "Blue")

# RGB composite plot
plot_rgb_raster(
  raster_data = rgb_stack,
  r = 1, g = 2, b = 3,
  stretch = "lin",
  title = "RGB Composite"
)

# False color composite
plot_rgb_raster(
  raster_data = rgb_stack,
  r = 2, g = 1, b = 3,  # Green-Red-Blue
  stretch = "hist",
  title = "False Color Composite"
)

# Simulate Ohio boundary for demonstration
ohio_boundary <- sf::st_polygon(list(matrix(c(
  -84.5, 38.5, -80.5, 38.5, -80.5, 42.0, -84.5, 42.0, -84.5, 38.5
), ncol = 2, byrow = TRUE)))
ohio_sf <- sf::st_sf(geometry = sf::st_sfc(ohio_boundary, crs = 4326))

# Map with region boundary
map_with_boundary <- create_spatial_map(
  spatial_data = sample_sf,
  fill_variable = "ndvi",
  region_boundary = ohio_sf,  # Would normally use "Ohio"
  title = "NDVI in Ohio",
  color_scheme = "ndvi"
)

print(map_with_boundary)

# The package supports many boundary types:
# create_spatial_map(data, region_boundary = "Ohio")          # US State
# create_spatial_map(data, region_boundary = "CONUS")         # Continental US
# create_spatial_map(data, region_boundary = "Ohio:Franklin") # State:County
# create_spatial_map(data, region_boundary = c(-84, 39, -82, 41)) # Bounding box

# Get available color schemes
color_schemes <- c("viridis", "plasma", "ndvi", "terrain", "water", "categorical")

# Apply different schemes to the same data
for (scheme in color_schemes[1:3]) {
  map <- create_spatial_map(
    spatial_data = sample_sf,
    fill_variable = "ndvi",
    color_scheme = scheme,
    title = paste("NDVI with", scheme, "colors"),
    point_size = 4
  )
  
  print(paste("Created map with", scheme, "color scheme"))
}

# Customize point appearance
styled_map <- create_spatial_map(
  spatial_data = sample_sf,
  fill_variable = "yield",
  map_type = "points",
  point_size = 6,
  color_scheme = "plasma",
  title = "Customized Point Map"
)

print(styled_map)

# Map with transparent points
# (Demonstrated conceptually - actual implementation may vary)

# Create high-resolution map for publication
publication_map <- create_spatial_map(
  spatial_data = sample_sf,
  fill_variable = "ndvi",
  color_scheme = "ndvi",
  title = "NDVI Distribution in Study Area",
  output_file = "publication_ndvi_map.png"
)

# Customize for journal specifications
journal_map <- create_spatial_map(
  spatial_data = sample_sf,
  fill_variable = "yield",
  color_scheme = "viridis",
  title = "",  # No title for journal figure
  output_file = "figure_1.png"
)

# The package automatically handles legends and layouts
# Maps include appropriate legends, scale bars, and formatting

# For complex layouts, combine with other packages:
# library(patchwork)  # For multi-panel figures
# library(ggplot2)    # For additional customization

# Test data validation
tryCatch({
  # This should work
  valid_map <- create_spatial_map(sample_sf, fill_variable = "ndvi")
  print("Valid map created successfully")
}, error = function(e) {
  print(paste("Error:", e$message))
})

# Handle missing data
sample_with_na <- sample_sf
sample_with_na$ndvi[1:2] <- NA

na_map <- create_spatial_map(
  spatial_data = sample_with_na,
  fill_variable = "ndvi",
  title = "Data with Missing Values"
)

print("Map with NA values handled automatically")

# Quick diagnostic check
diagnostic_result <- quick_diagnostic()

# Test specific mapping functions
mapping_test <- tryCatch({
  test_map <- create_spatial_map(sample_sf, fill_variable = "ndvi")
  "Mapping functions working"
}, error = function(e) {
  paste("Mapping error:", e$message)
})

print(mapping_test)

# Always check your data first
print("Sample data structure:")
print(head(sample_sf))

# Verify coordinate reference system
print(paste("CRS:", sf::st_crs(sample_sf)$input))

# Check for valid geometries
valid_geoms <- sf::st_is_valid(sample_sf)
print(paste("Valid geometries:", all(valid_geoms)))

# Start simple, then add complexity
simple_map <- quick_map(sample_sf)

# Add customization progressively
enhanced_map <- create_spatial_map(
  spatial_data = sample_sf,
  fill_variable = "ndvi",
  color_scheme = "ndvi",
  title = "Enhanced NDVI Map",
  point_size = 5
)

# Add interactivity if needed
# interactive_version <- create_interactive_map(sample_sf, fill_variable = "ndvi")

# For large datasets, consider:
# 1. Simplifying geometries
# 2. Reducing point density
# 3. Using raster instead of vector for very dense data

# Example: Check data size
print(paste("Number of features:", nrow(sample_sf)))
print(paste("Number of variables:", ncol(sample_sf) - 1))  # Minus geometry

# For large rasters, use terra's efficient plotting
if (exists("sample_raster")) {
  print(paste("Raster dimensions:", paste(dim(sample_raster), collapse = " x ")))
}

# GeoSpatialSuite maps work well with ggplot2
library(ggplot2)

# Extract ggplot object for further customization
base_map <- create_spatial_map(sample_sf, fill_variable = "ndvi")

# Customize with ggplot2 (if the map is a ggplot object)
if (inherits(base_map, "ggplot")) {
  enhanced_ggplot <- base_map +
    theme_minimal() +
    labs(caption = "Data source: Field measurements") +
    theme(
      plot.title = element_text(size = 16, face = "bold"),
      legend.position = "bottom"
    )
}

# For interactive web maps
if (requireNamespace("leaflet", quietly = TRUE)) {
  web_map <- create_interactive_map(
    spatial_data = sample_sf,
    fill_variable = "ndvi",
    popup_vars = c("ndvi", "yield"),
    basemap = "terrain"
  )
  
  # Further customize with leaflet functions
  enhanced_web_map <- web_map %>%
    leaflet::addMiniMap() %>%
    leaflet::addScaleBar()
}

# Create sample NDVI raster
ndvi_raster <- terra::rast(nrows = 40, ncols = 40, 
                           xmin = -84, xmax = -82, 
                           ymin = 39, ymax = 41)
terra::values(ndvi_raster) <- runif(1600, 0.1, 0.9)
names(ndvi_raster) <- "NDVI"

# Specialized NDVI visualization
ndvi_map <- create_ndvi_map(
  ndvi_data = ndvi_raster,
  title = "NDVI Analysis",
  ndvi_classes = "none"  # Use continuous colors
)

print("NDVI map created with specialized colors")

# Create sample water quality data
water_points <- data.frame(
  lon = c(-83.0, -83.1, -83.2, -82.9, -82.8),
  lat = c(40.0, 40.1, 40.2, 39.9, 40.3),
  dissolved_oxygen = c(8.2, 7.5, 6.8, 8.9, 7.1),
  temperature = c(18.5, 19.2, 20.1, 17.8, 19.5)
)

water_sf <- sf::st_as_sf(water_points, 
                         coords = c("lon", "lat"), 
                         crs = 4326)

# Water quality visualization
water_map <- create_water_quality_plot(
  water_data = water_sf,
  variable = "dissolved_oxygen",
  title = "Dissolved Oxygen Levels"
)

# High-resolution PNG export
create_spatial_map(
  spatial_data = sample_sf,
  fill_variable = "ndvi",
  color_scheme = "ndvi",
  title = "NDVI Distribution",
  output_file = "high_res_ndvi.png"
)

# PDF export for vector graphics
create_spatial_map(
  spatial_data = sample_sf,
  fill_variable = "yield",
  color_scheme = "viridis",
  title = "Yield Distribution",
  output_file = "yield_map.pdf"
)

# Export interactive map as HTML
if (requireNamespace("leaflet", quietly = TRUE)) {
  interactive_map <- create_interactive_map(
    spatial_data = sample_sf,
    fill_variable = "ndvi",
    title = "Interactive NDVI Map"
  )
  
  # Save as HTML file
  if (requireNamespace("htmlwidgets", quietly = TRUE)) {
    htmlwidgets::saveWidget(
      interactive_map,
      "interactive_ndvi_map.html",
      selfcontained = TRUE
    )
  }
}

# Create comparison visualization
create_comparison_map(
  data1 = sample_raster,
  data2 = sample_raster * 1.2,  # Simulated change
  comparison_type = "side_by_side",
  titles = c("Year 1", "Year 2"),
  color_scheme = "ndvi"
)

# The package automatically detects appropriate mapping approaches

# Point data -> scatter plot with colors
point_auto <- create_spatial_map(
  spatial_data = sample_sf,
  fill_variable = "ndvi",
  map_type = "auto"  # Auto-detects as points
)

# Raster data -> raster plot
raster_auto <- create_spatial_map(
  spatial_data = sample_raster,
  map_type = "auto"  # Auto-detects as raster
)

print("Auto-detection completed successfully")

# 1. Use terra plotting for speed
plot_raster_fast(sample_raster, title = "Fast Plotting")

# 2. Simplify data when appropriate
# simplified_sf <- sf::st_simplify(complex_sf, dTolerance = 100)

# 3. Use appropriate map types
# For very dense points, consider raster interpolation
# For large rasters, consider aggregation

# 4. Monitor memory usage
print(paste("Sample data memory usage:", 
            format(object.size(sample_sf), units = "KB")))

# Clean up large objects when done
# rm(large_raster)
# gc()  # Garbage collection

# Use temporary files for intermediate results
temp_file <- tempfile(fileext = ".tif")
print(paste("Temporary file:", temp_file))

# For very large analyses, process in chunks
# chunk_size <- 1000  # Adjust based on available memory

# Check CRS compatibility
sample_utm <- sf::st_transform(sample_sf, crs = 3857)  # Web Mercator

# The package handles CRS automatically in most cases
mixed_crs_map <- create_spatial_map(
  spatial_data = sample_utm,  # UTM coordinates
  fill_variable = "ndvi",
  title = "Map with Different CRS"
)

print("CRS handling completed automatically")

# Test with different data formats
formats_test <- list(
  sf_object = sample_sf,
  data_frame = sample_points,
  raster_object = sample_raster
)

for (format_name in names(formats_test)) {
  tryCatch({
    test_map <- quick_map(formats_test[[format_name]])
    print(paste(format_name, "format: OK"))
  }, error = function(e) {
    print(paste(format_name, "format error:", e$message))
  })
}