#
# PCRaster Script for creating flood maps
# 
# Analyses Water Surface map and removes non-contgouos sections
# Micha Werner
# Delft Hydraulics 2004
# 
# Input data (all taken from command line)
# 1: Digital elevation model 
# 2: Interpolated water surface map
# 3: delineation of main channel 
# 4: Map sections
# 5: resulting flood map
# 6: ID of current flood map section
# call
# pcrcalc -f pcr_flood_clump.mod <arguments>
# 1_ dem map file name
# 2_ interpolated flood surface map file
# 3_ River axis map file
# 4_ Section map file
# 5_ section id
# 6_ Return flood map name


binding
  dem       = $1;
  surf      = $2;
  river     = $3;
  sectmap   = $4;
  sectid    = $5;
  floodmap  = $6;  

  riverbln  = tmp_pcraxisb.map;
  clumpmap  = tmp_pcrclump.map;
  rivclump  = tmp_pcrclumr.map;
  classmap  = tmp_pcrclass.map;
  depthmap  = tmp_pcrdepth.map;
  sectnom   = tmp_sectmap.map;
  tmpsurf   = tmp_surfmap.map;
  
# there is no dynamic section!  
initial 
  # create the flood map sections as a nominal map
  sectnom  = nominal(sectmap);
  
  # filter out all areas except current id
  tmpsurf  = scalar(if(sectnom eq sectid then surf));
  

  # create a flood map using the clump option of all contigous areas
  depthmap  = scalar(if(surf gt dem then surf-dem));
  #clumpmap  = clump(nominal(if(surf gt dem then 1)));
  clumpmap  = clump(nominal(if(defined(depthmap) then 1)));
  # create map that contains the overlap between clump and river axis
  riverbln  = defined(river) and defined(depthmap);
  rivclump  = nominal(if(defined(river) then clumpmap));
  # create a rivclump map using the majority via the areaclass operation
 
  rivclump  = areamajority(rivclump,clumpmap);
  floodmap  = scalar(if(defined(rivclump)  then depthmap));