Oxygen saturation

## ====================================================================================
## World hypoxia, from the Global Ocean DataBase
## implemented by Karline Soetaert
## ====================================================================================
require(RNetCDF)
require(OceanView)
# The directory and name of the file with oxygen saturation data:
Directory <- "ftp://ftp.nodc.noaa.gov/pub/data.nodc/woa/WOA05nc/annual/"
File <- "O00an1.nc"
url <- paste(Directory, File, sep="")

# download only if not yet done
if (!File %in% list.files())
download.file(url, File, method = "auto", quiet = FALSE,
mode = "wb", cacheOK = TRUE)

# get data
Nc <- open.nc(File)
print.nc(Nc)

lat <- var.get.nc(Nc, "lat")
lon <- var.get.nc(Nc, "lon")
depth <- var.get.nc(Nc, "depth")
oxsat <- var.get.nc(Nc, "O00an1")

# convert to Atlantic-central
lon <- c(lon[181:360] - 360, lon[1:180])
oxsat <- oxsat[c(181:360, 1:180), , ]
par(bg = grey(0.5))
par(mfrow = c(1, 1))
par(mar = c(2, 2, 2, 2), oma = c(2, 0, 0, 0))

ii <- 1:360 #which (lon > -120 & lon < 150)
jj <- which (lat > -70 & lat < 90)

osat <- oxsat[ii,jj,]

x <- lon[ii]
y <- lat[jj]

xs <- c(x[1], x[length(x)]) # E boundary
ys <- y[seq(1,length(y), len = 8)] # S boundary

slice3D(colvar = osat, x = x,
y = y, z = -depth, #lighting = TRUE,
NAcol = "grey", xs = xs, ys = ys, zs = NULL,
theta = 10, phi = 60, colkey = list(length = 0.5),
expand = 0.5, ticktype = "detailed",
clab = "%", main = "Oxygen saturation",
xlab = "longitude", ylab = "latitude",
zlab = "depth", plot = FALSE)

image3D(col = "grey", x = x, y = y, z = -5000, #alpha = 0.05,
NAcol = "black", add = TRUE, plot = FALSE)

sea <- osat[,,1]
sea[! is.na(sea)] <- 1

image3D(colvar = sea, x = x, y = y, z = 50, #alpha = 0.05,
col = "transparent", NAcol = "black", add = TRUE, colkey = FALSE)