Fixed tide representation and a low resolution mode was added

elbe-2d-visualization.R

@@ -19,9 +19,10 @@ library(stringr)
 setwd('.')
 where = 'FLEXOUT' ## dir where the results of getm-elbe-2d were merged\
 where.to.save = './figures' ## dir where the results of getm-elbe-2d were merged
-file.name = 'out' ## name of file to plot
+file.name = 'out_ergom' ## name of file to plot
 prefix = c('nn_','nut_','phy_','zoo_','det_') ## model prefix used in the fabm.yaml file TODO: read from file
 ftype = 'png' ## output file type. Only 'pdf' or 'png'
+hi.res = T ## high.res uses Barnes interpolation (slow)
 
 hor.slices = 50 ## number of slices along the river in grid interpolation 
 ver.slices = 50 ## number of slices along the depth profile in grid interpolation
@@ -39,6 +40,7 @@ dates.to.plot = c('2020-01-15','2020-02-15','2020-03-15','2020-04-15',
 ########
 out = nc_open(paste0(where,'/',file.name,'.nc'))
 depth = ncvar_get(out, "bathymetry")
+height = ncvar_get(out, "D")
 hn = ncvar_get(out, "hn")
 long = ncvar_get(out, "lonc")
 lati = ncvar_get(out, "latc")
@@ -90,7 +92,7 @@ for(variable.name in variables.to.plot){ ## variable.name is the short name of v
   if(variable.name=='nn_oxy') palette(cmocean('oxy')(n.levels))
   if(variable.name=='salt') palette(cmocean('haline')(n.levels))
   if(variable.name=='temp') palette(rev(hcl.colors(n.levels,'Heat',rev=F)))
-  if(variable.name%in%c('nn_dia','nn_fla','nn_cya')) palette(colorRampPalette(c('tan','yellowgreen','forestgreen'))(n.levels))
+  if(variable.name%in%c('nn_dia','nn_fla','nn_cya')) palette(colorRampPalette(c('lightyellow','yellowgreen','forestgreen'))(n.levels))
   if(variable.name%in%c('light_par','nn_PAR')) palette(cmocean('solar')(n.levels))
   
   ##########
@@ -131,15 +133,25 @@ for(variable.name in variables.to.plot){ ## variable.name is the short name of v
   ##########
   for(i in time.list){
     ti=t[,,i]
-    
-    for(j in xx)for(z in 1:levels){
-      h = hn[j,,z]
+    D=height[,i]
+    for(j in xx){
+      h = hn[j,,i]
       dh[j,] = -depth[j]+cumsum(rev(h))
     }
     
-    tii = interpBarnes(rep(xx,levels),dh,ti,xgl=hor.slices, ygl=ver.slices,iterations = 1)
+    if(hi.res) tii = interpBarnes(rep(xx,levels),dh,ti,xgl=hor.slices, ygl=ver.slices,iterations = 1)
+    if(!hi.res){
+      tii = binMean2D(rep(xx,levels),dh,ti,xbreaks=pretty(xx,500), pretty(-40:5,150)) ## average over large regions
+      hor.slices = length(tii$xmids) ## resolutions are overwritten      
+      ver.slices = length(tii$ymids) ## resolutions are overwritten
+      tii$xg = tii$xmids ## renaming for using same code
+      tii$yg = tii$ymids
+      tii$zg = tii$result
+    } 
     censored = 1 + tii$zg * 0  
     for(k in 1:hor.slices)for(j in 1:ver.slices) if(tii$yg[j] < -depth[tii$xg[k]]) censored[k,j]=NA
+    for(k in 1:hor.slices)for(j in 1:ver.slices) if((tii$yg[j] + depth[tii$xg[k]]) > D[tii$xg[k]]) censored[k,j]=NA
+    
     tii$zg=tii$zg*censored
     local.lim = range(tii$zg,na.rm=T)
     image(tii$xg,tii$yg,tii$zg,ylim=c(-42,5),xaxt='n',yaxt='n',
@@ -147,6 +159,7 @@ for(variable.name in variables.to.plot){ ## variable.name is the short name of v
           col = palette()[(1+(local.lim[1]-abs.lim[1])/d.lim):(n.levels+n.levels*(local.lim[2]-abs.lim[2])/d.lim)])
     if(has.contours)contour(tii$xg,tii$yg,tii$zg,levels=scale.levels,add=T,col='white',drawlabels = F)
     lines(-depth,col=alpha('black',1.5),lwd=1,lty=1)
+    lines(D-depth,col=alpha('black',1.5),lwd=1,lty=1)
     abline(h = 0,col=alpha('black',0.5),lwd=1,lty=2)
     text(max(tii$xg),min(tii$yg),time[i],adj=c(1.1,0),font=2,cex=1.5)
     axis(side=2, at = seq(5,-40,by=-5),line=0,tcl=.25,labels=NA)