7 x 3D Surface Plots

# Load necessary libraries
pacman::p_load(pacman, plotly, dplyr, viridis, akima)

# List of file names with corresponding NCEI Accession numbers
file_list <- list(
  "SAOcean20140607_20140717.csv" = "0157403",
  "SAOcean20170205_20170501.csv" = "0184841",
  "SAOcean20190107_20190122.csv" = "0208364(a)",
  "SAOcean20190223_20190414.csv" = "0208364(b)",
  "SAOcean20190806_20190827.csv" = "0208364(c)",
  "SAOcean20191125_20191220.csv" = "0208364(d)",
  "SAOcean20200108_20200202.csv" = "0225435"
)

# Path to your CSV files
data_path <- "/cloud/project"

# Full file paths and accession numbers
full_file_paths <- file.path(data_path, names(file_list))
accession_numbers <- unlist(file_list)

# Function to load, clean, and plot individual datasets
generate_plot_for_file <- function(file, accession_number) {
  # Load the data
  data <- read.csv(file)
  
  # Select and clean the necessary columns
  variables <- c("LAT_dec_degree", "LONG_dec_degree", "Sal", 
                 "SST..deg.C.", "Tequ..deg.C.", "Pequ..hPa.", "Patm..hPa.")
  data_clean <- data %>% select(all_of(variables)) %>% na.omit()
  
  # Create unique grids for latitude and longitude
  x_unique <- seq(min(data_clean$LAT_dec_degree), max(data_clean$LAT_dec_degree), length.out = 100)
  y_unique <- seq(min(data_clean$LONG_dec_degree), max(data_clean$LONG_dec_degree), length.out = 100)
  
  # Function to interpolate and create grid data for each variable
  interpolate_to_grid <- function(variable) {
    interp_result <- interp(
      x = data_clean$LAT_dec_degree,
      y = data_clean$LONG_dec_degree,
      z = data_clean[[variable]],
      xo = x_unique,
      yo = y_unique,
      duplicate = "mean"
    )
    list(x = interp_result$x, y = interp_result$y, z = interp_result$z)
  }
  
  # Interpolate each variable
  sal_grid <- interpolate_to_grid("Sal")
  sst_grid <- interpolate_to_grid("SST..deg.C.")
  tequ_grid <- interpolate_to_grid("Tequ..deg.C.")
  pequ_grid <- interpolate_to_grid("Pequ..hPa.")
  patm_grid <- interpolate_to_grid("Patm..hPa.")
  
  # Create the combined plot with multiple surfaces
  combined_plot <- plot_ly()
  
  # Add Salinity surface
  combined_plot <- combined_plot %>% add_surface(
    x = sal_grid$x, y = sal_grid$y, z = sal_grid$z,
    colorscale = "Viridis",
    colorbar = list(title = "Salinity", len = 0.2, y = 0.95),
    opacity = 0.7
  )
  
  # Add SST surface
  combined_plot <- combined_plot %>% add_surface(
    x = sst_grid$x, y = sst_grid$y, z = sst_grid$z,
    colorscale = "Cividis",
    colorbar = list(title = "SST (°C)", len = 0.2, y = 0.75),
    opacity = 0.7
  )
  
  # Add Tequ surface
  combined_plot <- combined_plot %>% add_surface(
    x = tequ_grid$x, y = tequ_grid$y, z = tequ_grid$z,
    colorscale = "Plasma",
    colorbar = list(title = "Tequ (°C)", len = 0.2, y = 0.55),
    opacity = 0.7
  )
  
  # Add Pequ surface
  combined_plot <- combined_plot %>% add_surface(
    x = pequ_grid$x, y = pequ_grid$y, z = pequ_grid$z,
    colorscale = "Turbo",
    colorbar = list(title = "Pequ (hPa)", len = 0.2, y = 0.35),
    opacity = 0.7
  )
  
  # Add Patm surface
  combined_plot <- combined_plot %>% add_surface(
    x = patm_grid$x, y = patm_grid$y, z = patm_grid$z,
    colorscale = "Blues",
    colorbar = list(title = "Patm (hPa)", len = 0.2, y = 0.15),
    opacity = 0.7
  )
  
  # Add layout details
  combined_plot <- combined_plot %>% layout(
    title = paste("NCEI Accession:", accession_number, 
                  "Salinity, SST, Tequ, Pequ, Patm, Across Latitude and Longitude"),
    scene = list(
      xaxis = list(title = "Latitude"),
      yaxis = list(title = "Longitude"),
      zaxis = list(title = "Variable Scale")
    )
  )
  
  # Return the plot
  combined_plot
}

# Generate plots for all files
plots <- mapply(
  generate_plot_for_file,
  full_file_paths,
  accession_numbers,
  SIMPLIFY = FALSE
)

# Display all plots
for (plot in plots) {
  print(plot)
}