Westerschelde bathymetry

## ============================================================================
## A transparent bathymetry on google maps
## the bathymetry data is too large to share
## implemented by Karline Soetaert
## ============================================================================

require(RgoogleMaps)
require(OceanView)
require(RNetCDF)
Nc  <- open.nc("WSdiepte2012_LL.nc")
lat <- var.get.nc(Nc, "lat")
lon <- var.get.nc(Nc, "lon")
depth <- var.get.nc(Nc, "diepte2012_m")

destfile = tempfile(fileext=".png")
MetaInfo <- GetMap(center = c( 51.38, 4.18), zoom= 13, maptype = "satellite", destfile = destfile)

require(png)
img <- as.raster(readPNG(destfile))
par (mar = c(4, 4, 4, 6))
x <- seq(MetaInfo$BBOX$ll[2], MetaInfo$BBOX$ur[2], length.out = MetaInfo$size[2])
y <- seq(MetaInfo$BBOX$ll[1], MetaInfo$BBOX$ur[1], length.out = MetaInfo$size[1])
image2D(x, y, z = img[nrow(img):1,], 
       ylab = expression(""^o~N), xlab = expression(""^o~E), 
       main = "Westerschelde bathymetry",
       xaxs = "i", yaxs = "i")

ix <- which(lon > MetaInfo$BBOX$ll[2]-0.1 & lon <=  MetaInfo$BBOX$ur[2]+0.1)
iy <- which(lat > MetaInfo$BBOX$ll[1]-0.1 & lat <=  MetaInfo$BBOX$ur[1]+0.1)

DD <- depth[ix,iy]

image2D(lon[ix], lat[iy], z = DD, rasterImage = TRUE, 
  add = TRUE, NAcol = "transparent", clim = c(-33, 0),
   axes = FALSE, xlab = "", ylab = "", ylim = range(lat), shade = 0.3,
   colkey = list(length = 0.3, dist = 0.1, width = 0.5), alpha = 0.3,
   contour = list(level = -2, drawlabels = FALSE, col = "darkgrey"), clab = "m")