cws_stats__algorithm_stddev.pro

;+
; function algorithm_stddev
; 
; calculate and plot the algorithm standard deviation of a cloud product
; example call:
; $ std = o->algorithm_stddev (/calc, /jpg, /save)
;
; arguments and keywords
; hour     integer   hour of the satellite picture time slot, e.g. '12'
; product  string    product string, e.g. 'cod'
; relative keyword   calculate relative standard deviation (default
;                    absolute), e.g. /relative 
; no_plot  keyword   scip the plotting
; ps       keyword   produce a ps-file
; jpg      keyword   convert ps to jpg file
; calc     keyword   calculate the standard deviation (default read
;                    from file)
; save     keyword   save the standard deviation in a png file
; mini     float     minimum of the scale in the diagram
; maxi     float     minimum of the scale in the diagram
; mask:    string    specify the mask, possibilities:
;                    'com' for common mask 
;                    (average is only calculated when all SEVIRI retrievals provide data) or
;                    'ind' for individual mask 
;                    (average is calculated for all pixels with at least one SEVIRI retrieval)
;                    default is 'ind'
;
;  initial description UH 2012/04/18
;-

function cws_stats::algorithm_stddev, hour=hour, product=product, relative=relative, no_plot=no_plot, ps=ps, jpg=jpg, $
                                      calc=calc, save=save, mini=mini, maxi=maxi, log=log, mask=mask, no_FUB=no_FUB

  compile_opt idl2
  COMMON FEEDBACK, quiet, verbose, debug

  IF keyword_set(product) then self->set_product, product
  if n_elements(product) eq 0 then product = self.product
  if verbose then print,'    '
  if verbose then print,'... start cws_stats::algorithm_stddev '
  if n_elements(mini) eq 0 then mini = 0
  if n_elements(mask)    eq 0 then mask = 'ind'

  cd,!PROJECTS.cws_path+'/tools/', current=current
  @cws_stats_params
  grpList = self->create_group_list(product=product, guteInd=guteInd, no_FUB=no_FUB)

  if verbose then print, '... scale (1=full resolution): ', self.scale
  if keyword_set(save) and self.scale ne 1.0 then begin
    print, '*** Error in cws_stats::algorithm_stddev'
    print, '    saving results only permitting with a scale factor of 1.0'
    return, -1
  endif

  cd, current
  self.kind='img'
  cut = self.cut
  used_data=''

  IF keyword_set(relative) then ABS_or_REL='REL' else ABS_or_REL='ABS'

  if keyword_set(calc) then begin 
   
    ; get the average 
    if verbose then print,'... read multi algorithm average'
    ;avg = self->algorithm_average(/no_plot, log=log)
    avg = self->get_data(log=log, group='AVG')
    print, '... min/max of average: ', min(avg), max(avg)
   ; if self.product eq 'cod' then begin 
   ;   no_data_index = where(avg eq -1.)
   ;   avg=10.^avg
   ;   avg[no_data_index] = -1.
   ; endif
   ; print, 'average after scaling: ', min(avg), max(avg)
   ; help, avg

    ; define variance array
    var = avg*0.0
  
    ; loog over groups with data and sum up the variance  
    FOR gg = 0 , n_elements(guteInd) -1 DO BEGIN
  
      ; set next algorithm group
      self->set_group,grpList[gg]

      ; skip second cloud layers
      IF strlen(grpList[gg] gt 3) THEN BEGIN 
        print, '... skip data from group ', grpList[gg]
        CONTINUE 
      ENDIF 

      ; read data 
      img = self->get_data(log=log)
      print, '... min/max of group: ', grpList[gg], min(img), max(img)

      ; if no data available, then continue with next group
      IF (n_elements(img) eq 1) THEN CONTINUE

      if self.product eq 'cmb' then begin
        data_idx = img gt 0   ; data exists, where image 1 = cloud free or 2 = cloudy
        img = img-1
      endif else $ 
        data_idx = img ge 0  ; data exists, where image is greater equal 0

      ;; take care of logarithmic quantities  
      ;if self.product eq 'cod' then begin
      ;  no_data_index = where(img eq -1.)
      ;  img=10.^img
      ;  img[no_data_index] = -1.
      ;endif

      ; sum up all values, where data is avaibable 
      if used_data eq '' then begin 
        var = (img-avg)*(img-avg)*data_idx
        obs =  data_idx
      endif else begin 
        var = var + (img-avg)*(img-avg)*data_idx
        obs = obs + data_idx
      endelse

      ; remember used groups 
      if used_data eq '' and keyword_set(log) then used_data = 'log stddev: ' else used_data = 'stddev: '
      used_data = used_data + grpList[gg] + ' '
  
    ENDFOR 
  
    ; calculate the standard deviation (take care of no data)
    stddev = img*0.0
    delvarx, img
    stddev = sqrt(var/(obs-1)) 
    ; remove pixels with only one retrieval
    stddev[where(obs le 1)] = -1.

    ; if stddev is calculated in logarithmic space,
    ; transfer back to normal data
    if keyword_set(log) then begin
      print, '... convert from log10 scale back to normal scale'
      print, '    min/max stddev(log('+self.product+')) = ', min(stddev[where(obs gt 1)]), max(stddev[where(obs gt 1)])
      stddev[where(obs gt 1)] = 10^stddev[where(obs gt 1)]
      print, '    min/max 10^stddev(log('+self.product+')) = ', min(stddev[where(obs gt 1)]), max(stddev[where(obs gt 1)])
    endif 

    ; for common mask remove pixels where not every retrieval has data
    if mask eq 'com' then stddev[where(obs ne max(obs))] = -1.
    IF keyword_set(relative) then begin
      stddev[where(obs gt 1)] = stddev[where(obs gt 1)] / avg[where(obs gt 1)] * 100.
    endif 

  endif else begin 
    ;path = keyword_set(irregular_path) ? (irregular_path) : (self->build_data_path()+self.product_path)
    ;file1 = path + self.product+'-'+'STD-'+ABS_or_REL+'-'+(self->get_date(/string))+'-MSG.sav'
    ;print, '... read standard deviation from IDL save file ', file1
    ;RESTORE, FILENAME=file1
    self->set_group,'STD_'+ABS_or_REL
    stddev = self->get_data()
    ; take care of logarithmic storage of cloud optical depth
    if self.product eq 'cod' then begin 
      no_data_index = where(stddev eq -1.)
      stddev=10.^stddev
      stddev[no_data_index] = -1.
    endif       
    ;help, stddev
    ;stop
  endelse 

  ;; initialize unit and remarks
  self->get_metadata, unit=unit, remarks=bartitle

  ; save the standard deviation somehow
  IF keyword_set(save) then begin 
    ;path = keyword_set(irregular_path) ? (irregular_path) : (self->build_data_path()+self.product_path)
    ;file1 = path + self.product+'-'+'STD-'+ABS_or_REL+'-'+(self->get_date(/string))+'-MSG.sav'
    ;print, '... save standard deviation as an IDL save file ', file1
    ;save, stddev, obs, FILENAME=file1

    ;; BULLSHIT PROGRAMMING!!!
    ;stddev[where(obs gt 1)] = stddev[where(obs gt 1)] / avg[where(obs gt 1)] * 100.
    ;file1 = path + self.product+'-'+'STD-'+'REL'+'-'+(self->get_date(/string))+'-MSG.sav'
    ;print, '... save standard deviation as an IDL save file ', file1
    ;save, stddev, obs, FILENAME=file1

    stddev_scaled=stddev
    self->get_metadata, name=name, version=version, img_type=img_type, img_quality=img_quality, $
                    unit=unit, accuracy=accuracy, slope=slope, offset=offset, $
                    nodata_value_raw=nodata_value_raw, nodata_value_phy=nodata_value_phy, $
                    remarks=remarks

    ; change slope and offset for saving standard deviations
    offset = -1
    slope = (max(stddev_scaled) - offset) / 255.
    print, '*** offset: ', offset, ', slope: ', slope

    print,'min/max (alg_stddev, line 180) = ', min(stddev_scaled), max(stddev_scaled)

    self->scale_data, data=stddev_scaled, slope=slope, offset=offset, $
                    nodata_value_raw=nodata_value_raw, nodata_value_phy=nodata_value_phy

    print,'min/max (alg_stddev, line 184) = ', min(stddev_scaled), max(stddev_scaled)

    file = self.product+'-'+'STD_ABS'+'-'+(self->get_date(/string))+'-MSG'
    if keyword_set(relative) then file = self.product+'-'+'STD_REL'+'-'+(self->get_date(/string))+'-MSG'

    version = systime()

    ;help, stddev_scaled
    ;print,'min/max (save) = ', min(stddev_scaled), max(stddev_scaled)
    ;view2d, stddev, /cool, /colo;, no_data_idx=where(stddev le 0)
    ;stop
    self->save_product, stddev_scaled, name, img_type, img_quality, unit, accuracy, slope, offset, nodata_value_raw, nodata_value_phy,$
  	version, used_data, remarks, irregular_path = path, irregular_filename = file $
;	version, used_data, remarks, irregular_path = path, irregular_filename = file $
	,scale = 1, x_0 = 0, y_0 = 0, x_1 = 3711, y_1 = 3711

    delvarx, stddev_scaled

  ENDIF 


  ; plot the result 
  if not keyword_set(no_plot) then begin 

    if keyword_set(ps) or keyword_set(jpg) then $
      self->startDevice, xscale=1.6, yscale=1.2, n1='STDDEV-'+ABS_or_REL $
    else $ 
      window, 1
  
    coast=self->coast_line(/vec)
  
    self->get_metadata, name=name, version=version, img_type=img_type, img_quality=img_quality, $
                    unit=unit, accuracy=accuracy, slope=slope, offset=offset, $
                    nodata_value_raw=nodata_value_raw, nodata_value_phy=nodata_value_phy, $
                    remarks=remarks

    ; modify maxi 
    if n_elements(maxi) eq 0 then begin 
      IF keyword_set(relative) then maxi=100   $
      else maxi = (imgMax[cut] - imgMin[cut]) / sqrt(2)
      if keyword_set(log) then maxi=5.0
    endif


    ; set format of the plot
    case self.product of
      'cmb': begin 
        bartitle='cloud mask'
        ;mini=0
        if keyword_set(relative) then barformat='(I3)' else barformat='(F3.1)'
        unit='1'
      end
      'cph': begin 
        bartitle='cloud phase'
        ;mini=1
        ;barformat='(I3)'
      end 
      'ctt': begin 
        bartitle='cloud top temperature'
        if not keyword_set(relative) then begin & mini=0 & maxi=40. & endif
        barformat='(I3)'
        unit='K'
      end 
      'cth': begin 
      if not keyword_set(relative) then begin & mini=0 & maxi=6 & endif else begin & mini=0 & maxi=75 & endelse
        barformat='(I3)'
      end 
      'ctp': begin 
        ;mini=1
        barformat='(I3)' 
      end 
      'cod': begin 
        bartitle='cloud optical thickness'
        ;mini=0
        ;maxi=50
        barformat='(F3.1)'
        unit='1'
      end 
      'ref': begin 
        bartitle='effective radius'
        ;maxi=25
        barformat='(I3)'
        unit=textoidl('micro ')+'m'
      end 
      'lwp': begin 
        ;maxi=70
        barformat='(I3)'
      end 
      else: begin 
        print, '*** Error in cws_stats::stddev: plotting details not set for', self.product
        stop
      end
    endcase 

    if keyword_set(log) then barformat='(F4.2)'

    print,'min/max (line 220) = ', min(stddev), max(stddev)

    ; take care of logarithmic plotting
    if self.product eq 'cod' then begin 
      no_data_index = where(stddev eq -1.0)
      stddev = alog10(stddev)
      stddev[no_data_index] = -1.0
    endif 
    print,'min/max (line 226) = ', min(stddev), max(stddev)

    title='standard deviation '+self.product;+', '+ self.get_date(/nice)

    if keyword_set(relative) then title = 'relative '+title

    if keyword_set(log) then bartitle = '10^[ ' else bartitle = ''
    bartitle = bartitle + 'stddev (' 
    if keyword_set(log) then bartitle = bartitle +'log(' 
    bartitle = bartitle + self.product
    if unit ne '1' then bartitle = bartitle + ' / ' + unit
    bartitle = bartitle+')'
    if keyword_set(log) then bartitle = bartitle+') ]'

    ;if self.product eq 'cod' then  bartitle = 'log('+bartitle+')'
    if keyword_set(relative) then bartitle = bartitle+ ' / %'
    view2d, stddev, /cool, /colo, title=title, mini=mini, maxi=maxi, bar_nlev=5, coast_vec=coast, /box, no_data_idx=where(stddev lt 0), $ ;, no_data_idx=where(stddev lt 0)
            bartitle=bartitle, barformat=barformat, bardist=-0.15, aspect=1, position=[0.1, 0.0, 1.0, 0.95], col_table=brewer_tbl, $ ; col_table=16, $
            space_idx = Where(self->lon() Lt -200.), space_color = !P.background
  
    ;view2d, stddev, /cool, /colo, title=title, no_data_idx=where(avg le 0), log=plot_log, coast_vec=coast, aspect=1, $
    ;             /box, bartitle=bartitle, mini=float(imgMin[cut]), maxi=float(imgMax[cut]), bardist=-0.15, barformat=barformat, $
    ;             bar_nlev=5, position=[0.1, 0.0, 1.0, 0.95]

    if keyword_set(ps) or keyword_set(jpg) then begin 
      self->endDevice
      if keyword_set(jpg) then self->makeJPG, crop='970x1000+180+0'  ;, crop='780x650+180+0'
    endif 

  endif 

  return, stddev

end