source: trunk/SRC/ToBeReviewed/TRIANGULATION/drawsectionbottom.pro @ 76

;------------------------------------------------------------
;------------------------------------------------------------
;------------------------------------------------------------
;+
; NAME:drawsectionbottom
;
; PURPOSE:fill and draw the bottom continents for a real section.
;
; CATEGORY:
;
; CALLING SEQUENCE:
; 
; INPUTS:
;
; KEYWORD PARAMETERS:
;
;        COAST_COLOR: the color of the coastline.
;                     defaut value is 0 => black
;
;        COAST_THICK: the thickness of the coastline.
;                     defaut value is 1
;
;        CONT_COLOR: the color of the continent. defaut value is
;        (!d.n_colors - 1) < 255 => white
;
; OUTPUTS:
;
; COMMON BLOCKS:common.pro
;
; SIDE EFFECTS:
;
; RESTRICTIONS:simple way to fill continents for a section (using the
; fact that continents are wider at the bottom than at the top). 
;
; EXAMPLE:
;
; MODIFICATION HISTORY:Sebastien Masson (smasson@lodyc.jussieu.fr)
;                      June 14, 2002
;-
;------------------------------------------------------------
;------------------------------------------------------------
;------------------------------------------------------------
PRO drawsectionbottom, maskin, xxaxisin, depthsin $
                       , COAST_COLOR = coast_color, COAST_THICK = coast_thick $
                       , CONT_COLOR = cont_color, CONT_NOFILL = cont_nofill $
                       , OVERPLOT = overplot, _extra = ex
;---------------------------------------------------------
@cm_general
  IF NOT keyword_set(key_forgetold) THEN BEGIN
@updatekwd
  ENDIF
;--------------------------------------------------------- 
  if keyword_set(overplot) then return
;
; mask is from bottom to top
;----------
; boundaries conditions:
  nx = (size(maskin))[1]
  nz = (size(maskin))[2]+1
;
  IF size(xxaxisin, /n_dimensions) EQ 1 THEN $
    xxaxisin = temporary(xxaxisin)#replicate(1, nz)
  IF size(depthsin, /n_dimensions) EQ 1 THEN $
    depthsin = replicate(1, nx)#temporary(depthsin)
; for the mask : we add ocean at the top then it is always possible to
; find one ocean point on each water column.
  mask = [[maskin], [replicate(1, nx)]]
; for x axis we also add one level
;  xxaxis = [[xxaxisin], [xxaxisin[*, 0]]]
; x axis must cover nx+1 points because we will draw the edge of the
; mask.
; if it was mot possible in decoupeterre.pro to extend the xxaxis we
; do it now by hand
  xxaxis = xxaxisin[*, 0]
;  if (size(xxaxis))[1] EQ nx then begin
  if n_elements(xxaxis) EQ nx then begin
    deltax = abs(!x.range[1]-!x.range[0])/10.
    xxaxis = [xxaxis[0]-deltax, xxaxis ]
;    x0 = xxaxis[0]-deltax
;    xxaxis = [replicate(x0, 1, nz), xxaxis ]
  ENDIF
; for the depth,
  usepartial = total(depthsin, 2)
  usepartial = total(usepartial NE usepartial[0]) GE 1
  depths = depthsin[0:nx-1,*]
; we add one level according to the ocean level we had to the mask. 
  deltaz = abs(!y.range[1]-!y.range[0])/10.
  zmax = max(depthsin)+deltaz
;   depths = [[depths], [replicate(zmax, nx+1)]]
  depths = [[depths], [replicate(zmax, nx)]]
; if min(depths) gt -1 then we must add one line at the bottom (this
; appens when the bottom limit is defined between T[k] and W[k]
; points)
  IF min(depthsin) GT -1 THEN BEGIN
    zmin = min(!y.range)-deltaz
    depths = [[replicate(zmin, nx)], [depths]]
    mask = [[replicate(0, nx)], [mask]]
    nz = nz+1
  ENDIF
;
;----------
  xleft = xxaxis[0:nx-1]
  xright = xxaxis[1:nx]
; looking for the position of the bottom of the ocean.
  pos = nz-1-total(mask, 2)
  depths = depths[lindgen(nx)+nx*pos]
;
  xx = transpose([[xleft], [xright]])
  xx = !x.range[0] > xx[*] < !x.range[1]
  zz = transpose([[depths], [depths]])
  out = where(zz GE max(depthsin) OR zz LT min(!y.range), countout)
  zz = min(!y.range) > zz[*] < max(!y.range)
;----------
  xx = float(xx)
  zz = float(zz)
;----------
; filling of the continents
  IF NOT keyword_set(cont_color) THEN cont_color = (!d.n_colors - 1) < 255
;   if NOT keyword_set(cont_nofill) then $
;     polyfill, [min(xx,  max = maxx), xx[*], maxx] $
;               , [min(zz,  max = mazz), zz[*], mazz] $
;               , color = cont_color
   if NOT keyword_set(cont_nofill) then $
     polyfill, [min(xx, max = maxx), xx[*], maxx] $
               , [!y.range[0], zz[*], !y.range[0]], color = cont_color
;----------
; drawing of the coast (bottom) line...
; we could have plot directly xx and yy but if countout ne 0, doing
; this will draw an non-existing bottom line along !y.range values
; which is not so good...
; we thus do this ugly for/if loops to make sure that we don't draw
; these lines but we keep all vertical lines...
   IF countout NE 0 THEN BEGIN
     FOR i = 0, countout -1 DO BEGIN
       CASE 1 OF
         out[i] EQ 0:BEGIN 
; if we start with a out point
           xxx = !values.f_nan
           zzz = !values.f_nan           
         END 
         i EQ 0:BEGIN 
; i eq 0 but out[i] ne 0
           xxx = [xx[0:out[i]], !values.f_nan]
           zzz = [zz[0:out[i]], !values.f_nan]
         END
         ELSE:BEGIN
; two consecutive out values at the same depth: we just keep
; !values.f_nan values until the next change of depth
           IF (out[i-1] EQ out[i]-1) AND (zz[out[i-1]] EQ zz[out[i]]) THEN BEGIN 
             xxx = [xxx, !values.f_nan]
             zzz = [zzz, !values.f_nan]             
           ENDIF ELSE BEGIN 
; we keep everything inbetween the out values (including themselves
; for the vertical lines) but we had !values.f_nan to remove the horizontal lines
             xxx = [xxx, xx[out[i-1]:out[i]], !values.f_nan]
             zzz = [zzz, zz[out[i-1]:out[i]], !values.f_nan]             
           ENDELSE
         END
       ENDCASE
       IF i EQ countout-1 AND out[i] NE n_elements(xx)-1 THEN BEGIN
         xxx = [xxx, xx[out[i]:*]]
         zzz = [zzz, zz[out[i]:*]]             
       ENDIF
     ENDFOR
     plots, xxx, zzz, color = coast_color, thick = coast_thick, _extra = ex
   ENDIF ELSE plots, xx, zz, color = coast_color, thick = coast_thick, _extra = ex
;
  return
end