#
# Below are the R packages needed for this routine, you can get them from CRAN
#
#  This R routine reads the temperature-depth grid from CSIRO ".nc" files
#  It then produces a matrix of years/depth/temperature
#
# Data originally from: CSIRO
#
# Andy May 3/16/2026
#
# Required packages 
library(RNetCDF)
library(terra)
library(tidyterra)
library(ggplot2)
library(dplyr)
library(readr)
library(lubridate)
library(sf)
library(rnaturalearth)
library(RColorBrewer)

# ────────────────────────────────────────────────────────────────
# Function compute latitude weights for Earth's ellipsoid
# ────────────────────────────────────────────────────────────────

lat_weights_ellipsoid <- function(r) {
  # Extract latitude centers from the raster
  lats <- yFromRow(r, 1:nrow(r))
  
  # Grid spacing (HadSST and ERSST are 5°)
  dg <- abs(lats[2] - lats[1])
  
  # WGS-84 ellipsoid parameters
  a  <- 6378.137                # semi-major axis (km)
  e2 <- 0.00669437999014        # eccentricity squared
  
  # Ellipsoidal area integrand
  integrand <- function(phi) {
    2 * pi * a^2 * (1 - e2) * cos(phi) / (1 - e2 * sin(phi)^2)^2
  }
  
  # Integrate area for each latitude band
  nlat <- length(lats)
  Area_lat <- numeric(nlat)
  
  for (j in 1:nlat) {
    phi1 <- (lats[j] - dg/2) * pi/180
    phi2 <- (lats[j] + dg/2) * pi/180
    # Divide by number of longitude cells (360/dg) to get per‑cell area
    Area_lat[j] <- integrate(integrand, phi1, phi2)$value / (360 / dg)
  }
  
  # Normalize by equatorial band (constant factor cancels in weighted.mean)
  eq_idx <- which.min(abs(lats))
  Wt_lat <- Area_lat / Area_lat[eq_idx]
  
  return(Wt_lat)
}
# ────────────────────────────────────────────────────────────────

setwd("d:/Climate_Change/Ocean_Warming/CSIRO")
getwd()
#
# ────────────────────────────────────────────────────────────────
# Plot themes
# ────────────────────────────────────────────────────────────────

 theme_white <- theme(
  plot.background   = element_rect(fill = "white", colour = NA),
  panel.background  = element_rect(fill = "white", colour = NA),
  panel.grid.major  = element_line(colour = "grey92"),
  panel.grid.minor  = element_line(colour = "grey96"),
  axis.text         = element_text(colour = "black"),
  axis.title        = element_text(colour = "black"),
  plot.title        = element_text(colour = "black", hjust = 0.5)
 )

 theme_map <- theme(
  text = element_text(color = "blue"),

  # Make title blue
  plot.title = element_text(
    hjust = 0.5,
    size = 14,
    color = "blue"
  ),

  # Remove axis labels ("x" and "y")
  axis.title = element_blank(),

  # Remove axis tick labels
  axis.text = element_blank(),
  axis.ticks = element_blank(),

  # Background
  plot.background = element_rect(fill = "lightgray", color = NA),
  panel.background = element_rect(fill = "lightgray", color = NA),

  # Remove default grid
  panel.grid = element_blank(),

  # Legend text blue
  legend.text = element_text(color = "blue"),
  legend.title = element_text(color = "blue"),
  legend.position = "bottom",
  legend.key.width = unit(2, "cm")
 )

#
# fname_bot <- "temperature_cars2009_bot.nc"
fname_a <- "temperature_cars2009a.nc"
nc=open.nc(fname_a)
print.nc(nc)
# nc2=open.nc(fname_bot)
# print.nc(nc2)
#
#
# SST is the sea-surface temperature
# Create SpatRaster
# SST <- rast(fname_a)

# mean is the mean temperature array extracted from the NCDF file
# lon, lat, depth
mean_t = var.get.nc(nc,"mean",start=c(1,1,1),count=c(721,331,79))
var.inq.nc(nc, "mean")
depths = var.get.nc(nc,"depth")
lat = var.get.nc(nc,"lat")
lon = var.get.nc(nc,"lon")
str(depths)
str(lat)
str(long)
str(mean_t)
#
mean_t <- mean_t*0.00057986052828346 + 13

  nlon <- length(lon)
  nlat <- length(lat)

# CSIRO longitudes are 0–360; convert to -180–180 and reorder
  lon_shift <- ifelse(lon > 180, lon - 360, lon)
  ord_lon   <- order(lon_shift)
  lon2      <- lon_shift[ord_lon]

  #Latitudes are already south to north

# Reorder longitude dimension
mean_t2 <- mean_t[ord_lon, , ]   # apply same ordering as lon2

# Build a dummy raster to feed your weighting function
r_dummy <- rast(nrows = length(lat), ncols = length(lon2),
                xmin = min(lon2), xmax = max(lon2),
                ymin = min(lat), ymax = max(lat))

Wt_lat <- lat_weights_ellipsoid(r_dummy)

global_profile <- numeric(length(depths))

for (k in seq_along(depths)) {
  T_k <- mean_t2[ , , k]   # lon × lat slice
  
  # Weight by latitude only (lon cells equal width)
  w <- matrix(Wt_lat, nrow = length(lon2), ncol = length(lat), byrow = TRUE)
  
  global_profile[k] <- weighted.mean(T_k, w, na.rm = TRUE)
}

df_profile <- data.frame(
  depth = depths,
  temp  = global_profile
)
#
# ────────────────────────────────────────────────────────────────
# Plot depth-temperature profile
# ────────────────────────────────────────────────────────────────

p_profile <- ggplot(df_profile, aes(temp, depth)) +
                geom_path(color = "blue", size = 1.2) +
                scale_y_reverse(limits = c(6000, 0),
                       breaks = seq(0, 6000, 1000)) +
                scale_x_continuous(breaks = seq(0, 18, by = 1),
                       limits = c(0,18)) +
                labs(x = "Temperature (°C)", y = "Depth (m)",
                title = "Global Area-Weighted Temperature Profile") +
                theme_white +
                theme(panel.border = element_rect(colour = "black", 
                      fill = NA, linewidth = 0.8))
print(p_profile)
ggsave("CSIRO_area_weighted_depth_temp_profile.png", p_profile,
       width = 9, height = 5, dpi = 300, bg = "white")

#
# ────────────────────────────────────────────────────────────────
# Map temperature at 4500 meters
# ────────────────────────────────────────────────────────────────

# Robinson projection
rob <- "+proj=robin +lon_0=0 +datum=WGS84 +units=m +no_defs"

# Load and transform world map
world <- ne_countries(scale = "medium", returnclass = "sf")
world_rob <- st_transform(world, crs = rob)

# Build graticule in geographic coords, then transform
grat <- st_graticule(
  lat = seq(-60, 60, 15),
  lon = seq(-180, 180, 30)
) |> st_transform(rob)

# Custom longitude labels (bottom only)
lon_labs <- data.frame(
  lon = seq(-180, 180, 30),
  lat = -80,                     # place labels below map
  lab = paste0(seq(-180, 180, 30), "°")
) |>
  st_as_sf(coords = c("lon", "lat"), crs = 4326) |>
  st_transform(rob)

lat_labs <- data.frame(
  lon = -170,                     # place labels left of the map
  lat = seq(-60, 60, 15),
  lab = paste0(seq(-60, 60, 15), "°")
) |>
  st_as_sf(coords = c("lon", "lat"), crs = 4326) |>
  st_transform(rob)


# Extract the depth index closest to 4500 m
k4500 <- which.min(abs(depths - 4500))
T4500 <- mean_t2[ , , k4500]

# Build a data frame for mapping
df4500 <- expand.grid(
  lon = lon2,
  lat = lat
)

df4500$temp <- as.vector(T4500)
df4500 <- df4500 %>% filter(!is.na(temp))

# Convert to sf points
sf4500 <- st_as_sf(df4500, coords = c("lon", "lat"), crs = 4326)
sf4500_rob <- st_transform(sf4500, crs = rob)

# Plot map
p_map <- ggplot() +
            geom_sf(data = world_rob, fill = "grey80", color = "grey40") +
            geom_sf(data = sf4500_rob, aes(color = temp), size = 0.6) +
             scale_color_gradientn(
                colours = c("blue", "green", "yellow"),
                values  = scales::rescale(c(-1, 2, 5)),
                limits  = c(-1, 5),
                name    = "°C"
             ) +
            geom_sf(data = grat, color = "gray", linewidth = 0.2, alpha = 0.5) +

            # Longitude labels
              geom_sf_text(
                data = lon_labs,
                aes(label = lab),
                color = "blue", size = 3
              ) +

            # Latitude labels
              geom_sf_text(
                data = lat_labs,
                aes(label = lab),
                color = "blue",
                size = 3
              ) +
            labs(title = "Temperature at ~4500 m (CARS 2009)") +
            theme_map
print(p_map)
  ggsave(paste0("Temp at 4500m (CSIRO)", ".png"), p_map, width = 10, height = 6, 
                        dpi = 300, bg = "white")

#
# ────────────────────────────────────────────────────────────────
# Mapping error v depth
# ────────────────────────────────────────────────────────────────

map_err <- var.get.nc(nc, "map_error")
map_err <- map_err * 0.00057986052828346   # scale factor
map_err2 <- map_err[ord_lon, , ]
err_profile <- numeric(length(depths))

# map-error is all null.