Wadden Sea RWS sediment atlas

## =============================================================================
## Sediment properties from the WaddenSea
## implemented by Karline Soetaert
## =============================================================================
url <- "http://opendap.deltares.nl/thredds/catalog/opendap/rijkswaterstaat/sedimentatlas_waddenzee/"

# names can be created based on what is in a catalog
#korrel <- paste(url, "korrel.nc", sep = "")
#korrel <- "http://opendap.deltares.nl/thredds/dodsC/opendap/rijkswaterstaat/sedimentatlas_waddenzee/korrel.nc"
#download.file(korrel, "korrel.nc", method = "auto", 
#  quiet = FALSE, mode = "wb", cacheOK = TRUE) 
require(OceanView)
require(shape) #function A4
require(RNetCDF)
Nc <- open.nc("korrel.nc")
print.nc(Nc)
lon    <- var.get.nc(Nc, "lon")
lat    <- var.get.nc(Nc, "lat")
phi    <- var.get.nc(Nc, "phi")
cumphi <- var.get.nc(Nc, "cumphi")
dia    <- var.get.nc(Nc, "diameter")   # micrometer

# Four derived parameters
Mdphi <- apply (cumphi, MARGIN = 2, 
  FUN = function(x) approx(x, y = dia, xout = 50)$y)
Meanphi <- apply (phi, MARGIN = 2, 
  FUN = function(x) mean(x * dia /100))
gravel <- 100 - cumphi[dia == 1000, ]
silt   <- cumphi[dia == 63, ]

# 4 figures
require(mapdata)
mapNl <- function() 
  map("worldHires","Netherlands", xlim = range(lon), ylim = range(lat), 
  add = TRUE)

A4()
par(mfrow = c(4, 1), oma = c(4, 0, 2, 1), mar = c(1, 4, 3, 2))
scatter2D(lon, lat, colvar = Mdphi, pch = ".", cex = 4, log = "c", asp = TRUE,
  main = "md phi", clab = "uM", xlab = "", ylab = "", axes = FALSE, frame.plot = TRUE)
mapNl() 
scatter2D(lon, lat, colvar = Meanphi, pch = ".", cex = 4, log = "c", asp = TRUE,
  main = "mean phi", clab = "uM", xlab = "", ylab = "", axes = FALSE, frame.plot = TRUE)
mapNl() 
scatter2D(lon, lat, colvar = gravel+1e-2, pch = ".", cex = 4, log = "c", asp = TRUE, 
  main = "Gravel", clab = "%", clim = c(0.01,100), xlab = "", ylab = "", axes = FALSE, frame.plot = TRUE)
mapNl() 
scatter2D(lon, lat, colvar = silt, pch = ".", cex = 4, log = "c", asp = TRUE,
  main = "Silt-clay", clab = "%", xlab = "", ylab = "", axes = FALSE, frame.plot = TRUE)
mapNl() 

mtext(outer = TRUE, side = 3, text = "RWS sediment atlas", cex = 1.2)
AddSource(url)

# and a google maps figure - the mapping here is approximate...

require(RgoogleMaps)
destfile <- "NNL.png"
MetaInfo <- GetMap(center = c(53.3, 6.0), zoom = 8, 
             maptype = "satellite", destfile = destfile)
require(png)

windows()
par(mfrow = c(1, 1), oma  = c(2, 0, 0, 3))
scatter2D(c(MetaInfo$BBOX$ll[2], MetaInfo$BBOX$ur[2]), 
     c(MetaInfo$BBOX$ll[1], MetaInfo$BBOX$ur[1]), type = "n", 
     xlab = expression(""^o~E), ylab = expression(""^o~N), 
     xaxs = "i", yaxs = "i", cex.lab = 0.8,
     main = "RWS silt-clay", cex.main = 1.2)
img <- as.raster(readPNG(destfile))
rasterImage(img, xleft = MetaInfo$BBOX$ll[2], ybottom = MetaInfo$BBOX$ll[1], 
            xright = MetaInfo$BBOX$ur[2], ytop = MetaInfo$BBOX$ur[1])
box()
scatter2D(lon, lat, colvar = silt, pch = ".", cex = 3, log = "c", add = TRUE, 
  clab = "%", colkey = list(length = 0.5, width = 0.5, dist = 0.1))