Oxygen concentrations

## ===============================================================================
## 3-D scatter plot of nutrient and oxygen concentration of deep water
## data from the world Ocean Database
## implemented by Karline Soetaert
## ===============================================================================
library(OceanView)

## Load if not yet loaded
library(RNetCDF)
Directory <- "ftp://ftp.nodc.noaa.gov/pub/data.nodc/woa/WOA05nc/annual/"
url <- Directory
ftplist <- c("i00an1", "n00an1", "O00an1", 
             "p00an1", "s00an1", "t00an1")

path <- getwd()
Listfiles <- list.files(path)
for (i in 1:length(ftplist)) {
  File  <- paste(ftplist[i], ".nc", sep = "") 
  if (!File %in% Listfiles) {
    Url       <- paste(Directory, File, sep = "")
    download.file(Url, File, method =  "auto", quiet = FALSE, 
      mode = "wb", cacheOK = TRUE)   
  }
}

# read data and put in list

readWOA <- function(input) {
  D.nc <- open.nc(paste(path,"/",input, ".nc", sep = ""))
  lon      <- var.get.nc(D.nc, 'lon')
  lat      <- var.get.nc(D.nc, 'lat')     
  time     <- var.get.nc(D.nc, 'time')     
  depth    <- var.get.nc(D.nc, 'depth')     
  value    <- var.get.nc(D.nc, input)
  name     <- att.get.nc(D.nc, input, "long_name")
  units    <- att.get.nc(D.nc, input, "units")
 # create list with data
  return (list(lon = lon, lat = lat, depth = depth, time= time,
      value = value, name = name, units = units))
}
silicate  <- readWOA("i00an1") 
nitrate   <- readWOA("n00an1") 
oxygen    <- readWOA("O00an1") 
phosphate <- readWOA("p00an1")

# Scatterplot
par(mfrow = c(1, 1))
D <- 28
Depth <- phosphate$depth[D]
scatter3D(x = phosphate$value[,, D], y = nitrate$value[,, D], 
          z = silicate$value[,, D], colvar = oxygen$value[,, D],
          xlab = "Phosphate", ylab = "Nitrate", zlab = "Silicate",
          clab = c("Oxygen", "%"), pch = ".", 
          ticktype = "detailed", phi = 10, 
          main = paste("concentrations at ", Depth, " m"))