diff --git a/LIST_GROUPS.pro b/LIST_GROUPS.pro index 1398b5cef6f2fe860f92a682703341f8234f6af6..f54746e0401eaa1b6aeef2f7a117405127eb7826 100644 --- a/LIST_GROUPS.pro +++ b/LIST_GROUPS.pro @@ -38,7 +38,7 @@ ; 'UK Met-Office'] ; groups ngrp = n_elements(grpList) - if verbose then print, '... defined ',ngrp,' data groups' + if self.verbose then print, '... defined ',ngrp,' data groups' if (n_elements(grpList) ne n_elements(grpLong)) then print, '*** Warning: group-list and group-list-long dont have the same size.' diff --git a/cws_read__create_group_list.pro b/cws_read__create_group_list.pro index 6802819d7c171535fd7319500a0c9f12f62249e0..3d3c5c6679ea575ffe5bb11f7d0fb38fcdf6d63b 100644 --- a/cws_read__create_group_list.pro +++ b/cws_read__create_group_list.pro @@ -1,6 +1,9 @@ -function cws_read::create_group_list, product=product, single=single, group=group, guteInd=guteInd, no_FUB=no_FUB - COMMON FEEDBACK, quiet, verbose, debug +;; + +function cws_read::create_group_list, product=product, single=single, group=group, guteInd=guteInd, no_FUB=no_FUB + + ; set product if n_elements(product) ne 0 then begin product_org = self.product @@ -12,12 +15,17 @@ function cws_read::create_group_list, product=product, single=single, group=grou ; read gutInd for the product @cws_stats_params - + + + ; create group list, of those groups we wish plots if not keyword_set(single) then begin ; for multi plot, plot for all good groups group = 'all' - group_list = (*(self.grpList))[guteInd] + + + group_list = (self.coef['prd',product,'groups']).where(1) + endif else begin ; create group list if n_elements(group) eq 0 then begin @@ -39,12 +47,16 @@ function cws_read::create_group_list, product=product, single=single, group=grou guteInd = guteInd [where(group_list ne 'FUB')] endif - if verbose then $ + if self.verbose then $ for gg = 0, n_elements(group_list)-1 do print, ' '+group_list[gg] ; restore original product if n_elements(product) ne 0 then self->set_product, product_org - + + + + +stop return, group_list end diff --git a/cws_read__define.pro b/cws_read__define.pro index 9d107add850039a5306424954e1c63e646f9473f..b36cb35d154bc1565eae4d39c54aba8904fb95da 100755 --- a/cws_read__define.pro +++ b/cws_read__define.pro @@ -1056,6 +1056,9 @@ function cws_read::init $ self.fileMetaList = ptr_new(fileMetaList) self.mask_path = cws_mask_path self.code_source = cws_source_path + self.verbose = "true" + self.quiet = "false" + self.debug = "false" return,1 @@ -1223,7 +1226,10 @@ pro cws_read__define nodata_value_raw : 0 , $ nodata_value_physical : -1. , $ scale : 0. , $ - currentWorkingDirOrig : '' } + currentWorkingDirOrig : '', $ + verbose : Boolean(1) $ + , quiet : Boolean(1) $ + , debug : Boolean(1) } end diff --git a/cws_stats__all_one_image.pro b/cws_stats__all_one_image.pro index 682880c68edb78cf77c4a0fb801dc73fc2e55ba3..fbc4047bfa821adfee1d48eb12903e90b21e72a0 100644 --- a/cws_stats__all_one_image.pro +++ b/cws_stats__all_one_image.pro @@ -1,3 +1,4 @@ +; pro cws_stats::all_one_image $ , product = product $ , cut = cut diff --git a/cws_stats__constants.pro b/cws_stats__constants.pro index 965d779a44370e2f48360b53ff73bf807d02f2a4..bcb3fc977b900abc6b2c55ce3031771617b68bf8 100644 --- a/cws_stats__constants.pro +++ b/cws_stats__constants.pro @@ -10,14 +10,15 @@ grp = orderedhash() prd = orderedhash() reg = orderedhash() ; group list - grpList =['CMS','EUM','OCA','MPF','FUB','DLR','MFR','RMB','AWG','UKM','GSF','LAR','COX','OCA2','KNM','TPS','SUI','CLV'] - has_ctp = [1,1,1,1,0,0,1,0,1,1,0,0,0,0,0,1,1,1] + grpList =['CMS','EUM','OCA','MPF','FUB','DLR','MFR','RMB','AWG','UKM','GSF','LAR','COX','OCA2','KNM','TPS','SUI','CLV','LARN'] + has_ctp = [1,1,1,1,0,0,1,0,1,1,0,0,0,0,0,1,1,1,1] + has_cmb = [1,1,1,1,0,0,1,0,1,1,0,0,0,0,0,1,1,1,1] colList=['Green', 'Red', 'Blue', 'Cyan', 'Orange',$ 'Magenta', 'Dark Green', 'Cornflowerblue', 'Brown', 'Yellow' ,$ 'Dark Orchid','Tan', 'LightSeaGreen', 'Violet', 'Khaki',$ - 'LightCoral', 'Hot Pink', 'Yellow', 'Olive Drab', 'Dark Gray', 'Black','Green'] + 'DARK_OLIVE_GREEN', 'Hot Pink', 'Yellow', 'Olive Drab', 'Dark Gray', 'Black','Green','Tan'] grpLong = ['Climate Monitoring - SAF', $ @@ -37,7 +38,8 @@ colList=['Green', 'Red', 'Blue', 'Cyan', 'Orange',$ , 'KNMI' $ , 'Tropos Leipzig' $ , 'Meteo Suisse' $ - , 'CLAVR-x'] ; groups + , 'CLAVR-x' $ + ,'NASA, Langley Research Center'] ; groups for i=0,n_elements(grpList) - 1 do grp[grpList[i],'longname'] = grpLong[i] @@ -93,9 +95,11 @@ for ii = 0,n_elements(prdList) -1 do begin endfor for ii = 0,n_elements(has_ctp) -1 do begin - prd['ctp','groups',grp_idx] = has_ctp[ii] + prd['ctp','groups',grplist[ii]] = has_ctp[ii] + print,ii,grpList[ii] endfor prd['ctp','pmulti'] = [0,3,4] +prd['cmb','pmulti'] = [0,3,4] all_hash = hash() diff --git a/cws_stats__define.pro b/cws_stats__define.pro index 35149aa42c0816e286da123ecf4722db97eb64e4..44b845105bc23e49fab11ee60e423ea5c7d9b19d 100755 --- a/cws_stats__define.pro +++ b/cws_stats__define.pro @@ -971,117 +971,6 @@ end ;- -pro do_it $ - ,cut=cut $ - ,product=product $ - ,com=com $ - ,kind=kind $ - ,hour = hour $ - ,overpass = overpass $ - , group = group $ - , month = month $ - , day = day $ - ,_extra=_extra - - error_n=0 - - if n_elements(hour) eq 0 and n_elements(cut) eq 1 then begin - case cut of - 2: hour = 1230 - 3: hour = 1345 - 4: hour = 1345 - 5: hour = 1345 - 6: hour = 1400 - else: hour=1200 - endcase - endif - - default, cut, 0 - default, product, 'cod' - default, kind, 1 - - default, day, 13 - default, month, 6 - default, year, 2008 - default, hour, 1200 - - ;if total(product eq ['cmb','ovw','rgb'] ) eq 1 then myCom=0 else myCom = com - o=obj_new('cws_stats', cut=cut , com=com, _extra=_extra) - if (n_elements(group) eq 0) then group='AWG' - - IF product eq 'rgb' and ( not between(kind,3,4)) THEN RETURN - IF product eq 'cmb' and ( not between(kind,3,4)) THEN RETURN - IF product eq 'cph' and ( not between(kind,3,4)) THEN RETURN - ; CATCH,error - ; IF error ne 0 THEN BEGIN - ; error_n ++ - ; print,'Fehleranzahl: ',error_n - ; RETURN - ; ENDIF - - hr = strmid(string(hour,format='(i4)'),0,2) - mi = strmid(string(hour,format='(i4)'),2,2) - if total(mi eq ['00','15','30','45']) ne 1 then mi = '00' - - o->set_date,year,month,day,hr,mi - - thisX=!x - thisY=!y - thisP=!p - print,'kind: ',kind - CASE kind OF - - 1: o->make_corThumb, product=product, _extra=_extra, hour=hour ; in cws_stats__define - 2: o->make_corSingle, product=product, hour=hour - 3: o->make_image, product=product, hour=hour - 4: o->make_image, /single, product=product, hour=hour - 5: o->hist1d, product=product, hour=hour - 6: o->diff_maps, group='AWG', product=product, hour=hour - 7: o->make_tables, product=product, hour=hour - 8: o->calipso_cpr, product=product, overpass=overpass - 9: begin - if group eq 'ALL' then $ - group =['CMS','EUM','OCA','MPF','FUB','DLR','MFR','AWG','UKM','GSF','LAR'] - for tt=0,n_elements(group) -1 do $ - o->cpr_profile,group=group[tt],cut=cut,_extra=_extra - end - 10: o->cod_ctp , group = group - 11: o->cod_ctp , group = group - 12: begin - if group eq 'ALL' then begin - group =['CMS','EUM','OCA','MPF','FUB','DLR','MFR','AWG','UKM','GSF','LAR'] - for tt=0,n_elements(group) -1 do $ - o->caliop_profile, group=group[tt], cut=cut, _extra=_extra - endif else o->caliop_profile, cut=cut, _extra=_extra - end - 13: o->cth_all, cut=cut - 14: o->lwp_all, cut=cut - 15: o->cod_all, cut=cut - 16: o->ref_all, cut=cut - 17: o->brian_profile, cut=cut - 18: o->brian_lwp, cut=cut - 19: begin - if group eq 'ALL' then $ - group =['CMS','EUM','OCA','MPF','FUB','DLR','MFR','AWG','UKM','GSF','LAR'] - for tt=0,n_elements(group) -1 do $ - o->cpr_profile_stats,group=group[tt],cut=cut,_extra=_extra - end - 20: o->cph_all,cut =cut - 21: o->all_one_image,product=product,cut = cut - 22: o->cloudsat_profile_lonlat,cut=cut - 23: o->taylor, group='AWG', product='CTH', atrainSensor='CPR' ; atrainSensor='AMSR'/'CALIOP'/'CPR'/'MODIS' ('ECMWF') - else: begin - print, '*** Error in do_it (cws_stats__define.pro), unknown kind ', kind,' of plot (1-23)' - end - - ENDCASE - - !x=thisX - !y=thisY - !p=thisP - catch,/cancel - -end ;--------------------------------------------------------- diff --git a/cws_stats__hist1d.pro b/cws_stats__hist1d.pro index cc60cd56c7996552414705fa08a48fdc6bc13cdc..569042e457448fa3cfa861d2edc91340abd69067 100644 --- a/cws_stats__hist1d.pro +++ b/cws_stats__hist1d.pro @@ -26,7 +26,7 @@ pro cws_stats::hist1d, hour=hour, product=product, mask=mask, surf=surf, phase=p save_data=save_data, read_data=read_data, xrange=xrange, paper=paper, no_FUB=no_FUB compile_opt idl2 - COMMON FEEDBACK, quiet, verbose, debug + self.kind = 'his' ; choose histogram as kind @@ -34,7 +34,7 @@ pro cws_stats::hist1d, hour=hour, product=product, mask=mask, surf=surf, phase=p if n_elements(surf) eq 0 then surf = 'all' if n_elements(phase) eq 0 then phase = 'all' - if not quiet then print, '... start hist1d, surf='+surf+', phase='+phase+' (cws_start__define.pro)' + if not self.quiet then print, '... start hist1d, surf='+surf+', phase='+phase+' (cws_start__define.pro)' ; replace product, if user wants another product IF keyword_set(product) then self->set_product,product @@ -85,12 +85,12 @@ pro cws_stats::hist1d, hour=hour, product=product, mask=mask, surf=surf, phase=p ;--------------------------- first=0 if (mask eq 'COM' or mask eq 'both') then begin - - FOR gg = 0, n_elements(guteInd) -1 DO BEGIN - + + ;FOR gg = 0, n_elements(guteInd) -1 DO BEGIN + foreach group, groups do begin ; set group ;self->set_group, (*(self.grpList))[guteInd[gg mod n_elements(guteInd)]] - self->set_group, groups[gg] + self->set_group, group ; read data img = self->get_disks(surf=surf, phase=phase, n_disks=n_disks) if self.product eq 'cmb' then img = img-1. @@ -103,7 +103,7 @@ pro cws_stats::hist1d, hour=hour, product=product, mask=mask, surf=surf, phase=p common_mask = (img ge 0) first = 1 endif else common_mask = common_mask * (img ge 0) - ENDFOR + ENDFOREACH endif endif @@ -130,13 +130,13 @@ pro cws_stats::hist1d, hour=hour, product=product, mask=mask, surf=surf, phase=p if ~keyword_set(read_data) then begin - if verbose and im eq 0 then print, '... calculate histogram for individual mask ', im - if verbose and im eq 1 then print, '... calculate histogram for common mask ', im + if self.verbose and im eq 0 then print, '... calculate histogram for individual mask ', im + if self.verbose and im eq 1 then print, '... calculate histogram for common mask ', im print, 'n_grp', n_grp ; FOR gg = 0, 1 -1 DO BEGIN FOR gg = 0, n_grp -1 DO BEGIN - + ; set group ; self->set_group,(*(self.grpList))[guteInd[gg]] self->set_group, groups[gg] @@ -195,7 +195,7 @@ pro cws_stats::hist1d, hour=hour, product=product, mask=mask, surf=surf, phase=p hist[im,gg,*] = histo [0:n_hist-1] - if verbose then begin + if self.verbose then begin print, ' '+(*(self.grpList))[guteInd[gg]], reform(hist[im,gg,*]) endif @@ -313,6 +313,7 @@ pro cws_stats::hist1d, hour=hour, product=product, mask=mask, surf=surf, phase=p p.ytitle = xtitle p.position = [0.23, 0.12, 0.79, 0.85] + stop endelse p_list = list() @@ -326,7 +327,7 @@ pro cws_stats::hist1d, hour=hour, product=product, mask=mask, surf=surf, phase=p thick=1., color=fsc_color(colList[guteInd[gg mod n_elements(guteInd)]+2]) endif else begin ; p0 = plot(hist[im,gg,*], x_axis_center,/overplot, color = fsc_color(colList[guteInd[gg mod n_elements(guteInd)]+2]),name=groups[gg] ) - p0 = plot(hist[im,gg,*], x_axis_center,/overplot, color = (self.constants())['grp',groups[gg],'color'] ,name=groups[gg] ) + p0 = plot(hist[im,gg,*], x_axis_center,/overplot, color = (self.coef)['grp',groups[gg],'color'] ,name=groups[gg] ) p_list.add,p0 endelse @@ -336,10 +337,11 @@ pro cws_stats::hist1d, hour=hour, product=product, mask=mask, surf=surf, phase=p if self.product ne 'ctp' then begin ; devide legend in two parts if more than 6 groups + IF n_elements(guteInd) GE 6 THEN BEGIN n_elem_gute = n_elements(guteInd) ; first half of the legend - + stop cws_legend, groups[0:((n_elem_gute/2) -1)] $ ,linestyle=0 $ ,cws_color=guteInd[0:((n_elem_gute/2)-1)]+2 $ diff --git a/cws_stats__make_image.pro b/cws_stats__make_image.pro index 00ef962f996de5d543deb48e27ddea4115cec72a..ffeac960d33a780a772a268a6efc871ac5acdea5 100644 --- a/cws_stats__make_image.pro +++ b/cws_stats__make_image.pro @@ -458,249 +458,8 @@ pro cws_stats::make_image, group=group, hour=hour, product=product, single=singl ;----------------------------------------------------- 'cmb': begin - - x_def = !x - y_def = !y - - ; adjust the size of the diagram / and open device - if keyword_set(single) then begin - ; single plot - !p.multi=0 - !x.margin = 2.0 - !y.margin = [0.0,0.2] - !y.omargin = [0.2,0.2] - charsize = 0.8 - endif else begin - !p.multi=pmultImg - self->startDevice, xscale=0.2+(2./3.*float(pmultimg[1])), yscale=2./3.* float(pmultimg[1]), addn=addName - !x.omargin=1. - !y.omargin=[2.,6.] - !x.margin=0.9 - !y.margin=[0.2,2.2] - ;self->startDevice, xscale=2.1, yscale=1.8 - endelse - - ; read land sea mask - ls = self->land_sea() - - ; get longitude - lon = self->lon() - - ; initialize counter for disagreements - cld = ls * 0 - obs = cld - n_data = 0 - - ; loop over groups, take only the groups with data (guteInd) - ;FOR gg = 0, n_elements(guteInd) -1 DO BEGIN - FOR gg = 0, n_elements(group_list)-1 DO BEGIN - - ; skip the COX group, if phase is not ice - if group_list[gg] eq 'COX' and phase ne 'ice' then continue - - ; set the group - ;group_str = (*(self.grpList))[guteInd[gg]] - group_str = group_list[gg] - if not quiet then print, ' create cmb plot for ', group_str - - ;; read the data (store it in array imp) - ;if group_str eq 'OCA' then begin - ; ; special MPF has same cmb as OCA - ; self->set_group,'MPF' - ; img = self->get_data() - ; self->set_group,'OCA' - ;endif else begin - ;; read data (normal case) -; self->set_group,(*(self.grpList))[guteInd[gg]] - self->set_group, group_str - img = self->get_data(phase=phase, surf=surf ) - ;endelse - - ; skip plot if single plot is wanted and there is no data - if keyword_set(single) and n_elements(img) eq 1 then return - - ; for single plot open device for each group - if keyword_set(single) then self->startDevice, xscale=1.2, yscale=1.2, /single, n1 = group_str, addn=addName - - ; first normal data, then no data - if n_elements(img) ne 1 then begin - - if group ne 'AVG' then begin - ; modify cloud data, cloud -> 3, land -> 1, sea -> 0 - dum = (0 * img) $ - + 3*(img gt 1.5) $ - + 1*(between(img,0.5,1.5) And (ls eq 255) ) $ - + 2*(between(img,0.5,1.5) And (ls eq 127) ) - ; specify color table - color_table = 14 - mini=1 - maxi=3 - ; calculate fraction of cloudy pixels for this scene - frc = round(100.* total(img gt 1.5)/float(total(img gt 0.0)) ) - print, ' ', group_str + ' (', frc, '%)', min(img), max(img) - color_bar = 0 - - ; plot image - view2d, dum $ - , color_bar = color_bar $ - , no_data_idx=where(img le 0) $ - , coast_vec=size(coast,/type) eq 10 ? coast : 0 $ - , title = group_str + ' ('+makes(frc,frc,1,di=2) + '%)' $ - , aspect=1 $ - , coast_col=0 $ - , mini=mini, maxi=maxi $ - , ticklength =0. $ - , /cool $ - , /no_axes $ - , charsize=charsize $ - , space_idx = where(lon lt -200) $ - , space_col= !p.background $ - , col_table = color_table, position=[0.05,0.0,1.00,1.00] ;, pos2= [1200, 300, 7500, 7500] - delvarx, crop - endif else begin - ; do not modify the multi algorithm average - dum = img - ; specify color table - color_table = 13 - mini=0 - maxi=1 - ; calculate fraction of cloudy pixels for this scene - frc = round(100.* total(img gt 0.5)/float(total(img gt 0.0)) ) - print, ' ', group_str + ' (', frc, '%)', min(img), max(img) - bartitle = 'cloud cover' - - ; !!! cheating !!! - img[where(img eq max(img))] = 1.0 - - ;; plot image - view2d, img $ - ,/cool, color_bar = 1 $ - ,no_data_idx = where(img eq -1.) $ - ,mini=mini, maxi=maxi $ - ;,mini=0, maxi=2.3 $ - ,title = 'algorithm average' + ' ('+makes(frc,frc,1,di=2) + '%)' $ - ,coast_vec=size(coast,/type) eq 10 ? coast : 0 $ - ,coast_thick = 1.1 $ - ,coast_color = fsc_color('black') $ - ,/no_axes $ - ,barformat=Imgformat $ - ,aspect=1 $ - ,barcharsize=1. $ - ,bartitle = bartitle $ - ,Col_Table = 17 $ - ,brewer = 1 $ - ,space_idx = Where(self->lon() Lt -200.) $ - ,space_color = !P.background $ - ,barsize=0.18 $ - ,bardist = -0.14 $ - ,charsize = charsize, position=[0.05,0.18,0.97,0.97] - crop = '725x610+30+0' - endelse - - - ; do statistics for disagreement plot - ; add 1 to each pixel with cloud (for gg=0 initialize) - cld = gg eq 0 ? img gt 1.5 : temporary(cld) + (img gt 1.5) - ; add 1 to each pixel with obersavtion (for gg=0 initialize) - obs = gg eq 0 ? img ne 0 : temporary(obs) + (img ne 0) - ; add 1 to number of good data sets - n_data++ - - endif else begin - - dum = fltarr( 742, 742) - dum = + 1*(ls eq 255) + 2*(ls eq 127) - - ; if no data is available, plot empty disc - view2d,dum $ - ,no_data_idx=where(img le 0) $ - ,coast_vec=size(coast,/type) eq 10 ? coast : 0 $ - ,title = group_str $ - ,aspect=1 $ - ,coast_col=0 $ - ,mini=1 $ - ,maxi=3 $ - ,ticklength =0. $ - ,/cool $ - ,/no_axes $ - ,charsize=charsize $ - ,space_idx = where(lon lt -200) $ - ,space_col= !p.background $ - ,col_table = 14, position=[0.05,0.0,1.00,1.00] - ;,/nodata - - ; write 'NO DATA' over the plot - xyouts,((size(dum,/dim))[0])/2.,((size(dum,/dim))[0])/2. $ - ,'NO DATA',orientation = 45,align = 0.5 ,charsize=2. $ - , color =fsc_color('gold') - - endelse - - - ;; write marker of plot, that OCA copies mask from MPF - ;if group_str eq 'OCA' then begin - ; xyouts,((size(dum,/dim))[0])/2.,((size(dum,/dim))[0])/2.,'USES MPF !C cloud mask'$ - ; , orientation = 45,align = 0.5 ,charsize=1.4 $ - ; , color =fsc_color('red') - ;endif - - ; single plot for each group ... - if keyword_set(single) then begin - ; write date and region above the title - xyouts, 0.5, 0.95, 'Cloud mask ' + ' ' + self->get_date(/nice) $ - +' Cut ' + makes(self.cut,self.cut,1,di=2), $ - color=fsc_color('black'), /normal, $ - alignment=0.5, charsize=1.0 - - ; close file for each group - self->endDevice - if ~keyword_set(ps) then self->makeJPG, crop=crop - endif - - ENDFOR - - ; multiplot for all groups - if not keyword_set(single) then begin $ - ; write write date and region - xyouts, 0.5, 0.95, 'Cloud mask ' + ' ' + self->get_date(/nice) $ - +' Cut ' + makes(self.cut,self.cut,1,di=2) + addtitle, $ - color=fsc_color('black'), /normal, $ - alignment=0.5 ,charsize=1.5 - - ; close file - self->endDevice - if ~keyword_set(ps) then self->makeJPG - endif - - ; ----------- - ; create plot for disagreement - if group eq 'all' or group eq 'ALL' then begin - if not quiet then print,'... create agreement plot' - !p.multi=0 - !x.margin = 0.5 - !y.margin = [0.05,0.1] - !y.omargin = [0.2,0.2] - charsize = 0.97 - self->startDevice, xscale=1.3, yscale=1.3, n1 = 'AGR', addn=addName ;, addname='_MATCH_PERC' - - ; maximal number of disagreements - n_max = float(floor(0.5*obs)) - - ; plot image - view2d, float((obs-cld) < cld) / n_max, /cool, /colo, no_data_idx=where(cld eq 0),$ - title='Cloud mask disagreement', $ - coast_vec=size(coast,/type) eq 10 ? coast : 0, coast_col=0, aspect=1, $ - /no_axes, min=0., max=1., bar_nlev=5,$ - bardist=-0.15, /discrete, barformat='(F3.1)',charsize=charsize, bartitle='disagreement', $ - space_idx = where(self->lon() lt -200), space_col = !p.background, $ - col_table = 15, position=[0.1, 0.0, 1.0, 1.0] ; col_table = 16, - - - self->endDevice - if ~keyword_set(ps) then self->makeJPG, addn='_MATCH_DISAGR', crop='780x660+130+160' - endif - ; ----------- - + self.make_image_cmb,single = single, addtitle = addtitle, group = group + end ;----------------------------------------------------- diff --git a/cws_stats__make_image_cmb.pro b/cws_stats__make_image_cmb.pro new file mode 100644 index 0000000000000000000000000000000000000000..949215096230a3d62c0b143cff01486f98035b61 --- /dev/null +++ b/cws_stats__make_image_cmb.pro @@ -0,0 +1,260 @@ + +pro cws_stats::make_image_cmb $ + , single = single $ + , addtitle=addtitle $ + , group = group + + + +COMMON FEEDBACK, quiet, verbose, debug + + + + x_def = !x + y_def = !y + group_list = self->create_group_list(product=product, group=group, guteInd=guteInd, single=single, no_FUB=no_FUB) + stop + + ; adjust the size of the diagram / and open device + if keyword_set(single) then begin + ; single plot + !p.multi=0 + !x.margin = 2.0 + !y.margin = [0.0,0.2] + !y.omargin = [0.2,0.2] + charsize = 0.8 + endif else begin + !p.multi=self.coef['prd','cmb','pmulti'] + ;self->startDevice, xscale=0.2+(2./3.*float(pmultimg[1])), yscale=2./3.* float(pmultimg[1]), addn=addName + window,1,xsize=500,ysize=1000 + !x.omargin=1. + !y.omargin=[2.,6.] + !x.margin=0.9 + !y.margin=[0.2,2.2] + ;self->startDevice, xscale=2.1, yscale=1.8 + endelse + + ; read land sea mask + ls = self->land_sea() + + ; get longitude + lon = self->lon() + + ; initialize counter for disagreements + cld = ls * 0 + obs = cld + n_data = 0 + stop + ; loop over groups, take only the groups with data (guteInd) + ;FOR gg = 0, n_elements(guteInd) -1 DO BEGIN + FOR gg = 0, n_elements(group_list)-1 DO BEGIN + + ; skip the COX group, if phase is not ice + if group_list[gg] eq 'COX' and phase ne 'ice' then continue + + ; set the group + ;group_str = (*(self.grpList))[guteInd[gg]] + group_str = group_list[gg] + if not quiet then print, ' create cmb plot for ', group_str + + ;; read the data (store it in array imp) + ;if group_str eq 'OCA' then begin + ; ; special MPF has same cmb as OCA + ; self->set_group,'MPF' + ; img = self->get_data() + ; self->set_group,'OCA' + ;endif else begin + ;; read data (normal case) +; self->set_group,(*(self.grpList))[guteInd[gg]] + self->set_group, group_str + img = self->get_data(phase=phase, surf=surf ) + ;endelse + + ; skip plot if single plot is wanted and there is no data + if keyword_set(single) and n_elements(img) eq 1 then return + + ; for single plot open device for each group + if keyword_set(single) then self->startDevice, xscale=1.2, yscale=1.2, /single, n1 = group_str, addn=addName + + ; first normal data, then no data + if n_elements(img) ne 1 then begin + + if group ne 'AVG' then begin + ; modify cloud data, cloud -> 3, land -> 1, sea -> 0 + dum = (0 * img) $ + + 3*(img gt 1.5) $ + + 1*(between(img,0.5,1.5) And (ls eq 255) ) $ + + 2*(between(img,0.5,1.5) And (ls eq 127) ) + ; specify color table + color_table = 14 + mini=1 + maxi=3 + ; calculate fraction of cloudy pixels for this scene + frc = round(100.* total(img gt 1.5)/float(total(img gt 0.0)) ) + print, ' ', group_str + ' (', frc, '%)', min(img), max(img) + color_bar = 0 + + ; plot image + view2d, dum $ + , color_bar = color_bar $ + , no_data_idx=where(img le 0) $ + , coast_vec=size(coast,/type) eq 10 ? coast : 0 $ + , title = group_str + ' ('+makes(frc,frc,1,di=2) + '%)' $ + , aspect=1 $ + , coast_col=0 $ + , mini=mini, maxi=maxi $ + , ticklength =0. $ + , /cool $ + , /no_axes $ + , charsize=charsize $ + , space_idx = where(lon lt -200) $ + , space_col= !p.background $ + , col_table = color_table, position=[0.05,0.0,1.00,1.00] ;, pos2= [1200, 300, 7500, 7500] + delvarx, crop + endif else begin + ; do not modify the multi algorithm average + dum = img + ; specify color table + color_table = 13 + mini=0 + maxi=1 + ; calculate fraction of cloudy pixels for this scene + frc = round(100.* total(img gt 0.5)/float(total(img gt 0.0)) ) + print, ' ', group_str + ' (', frc, '%)', min(img), max(img) + bartitle = 'cloud cover' + + ; !!! cheating !!! + img[where(img eq max(img))] = 1.0 + + ;; plot image + view2d, img $ + ,/cool, color_bar = 1 $ + ,no_data_idx = where(img eq -1.) $ + ,mini=mini, maxi=maxi $ + ;,mini=0, maxi=2.3 $ + ,title = 'algorithm average' + ' ('+makes(frc,frc,1,di=2) + '%)' $ + ,coast_vec=size(coast,/type) eq 10 ? coast : 0 $ + ,coast_thick = 1.1 $ + ,coast_color = fsc_color('black') $ + ,/no_axes $ + ,barformat=Imgformat $ + ,aspect=1 $ + ,barcharsize=1. $ + ,bartitle = bartitle $ + ,Col_Table = 17 $ + ,brewer = 1 $ + ,space_idx = Where(self->lon() Lt -200.) $ + ,space_color = !P.background $ + ,barsize=0.18 $ + ,bardist = -0.14 $ + ,charsize = charsize, position=[0.05,0.18,0.97,0.97] + crop = '725x610+30+0' + endelse + + + ; do statistics for disagreement plot + ; add 1 to each pixel with cloud (for gg=0 initialize) + cld = gg eq 0 ? img gt 1.5 : temporary(cld) + (img gt 1.5) + ; add 1 to each pixel with obersavtion (for gg=0 initialize) + obs = gg eq 0 ? img ne 0 : temporary(obs) + (img ne 0) + ; add 1 to number of good data sets + n_data++ + + endif else begin + + dum = fltarr( 742, 742) + dum = + 1*(ls eq 255) + 2*(ls eq 127) + + ; if no data is available, plot empty disc + view2d,dum $ + ,no_data_idx=where(img le 0) $ + ,coast_vec=size(coast,/type) eq 10 ? coast : 0 $ + ,title = group_str $ + ,aspect=1 $ + ,coast_col=0 $ + ,mini=1 $ + ,maxi=3 $ + ,ticklength =0. $ + ,/cool $ + ,/no_axes $ + ,charsize=charsize $ + ,space_idx = where(lon lt -200) $ + ,space_col= !p.background $ + ,col_table = 14, position=[0.05,0.0,1.00,1.00] + ;,/nodata + + ; write 'NO DATA' over the plot + xyouts,((size(dum,/dim))[0])/2.,((size(dum,/dim))[0])/2. $ + ,'NO DATA',orientation = 45,align = 0.5 ,charsize=2. $ + , color =fsc_color('gold') + + endelse + + + ;; write marker of plot, that OCA copies mask from MPF + ;if group_str eq 'OCA' then begin + ; xyouts,((size(dum,/dim))[0])/2.,((size(dum,/dim))[0])/2.,'USES MPF !C cloud mask'$ + ; , orientation = 45,align = 0.5 ,charsize=1.4 $ + ; , color =fsc_color('red') + ;endif + + ; single plot for each group ... + if keyword_set(single) then begin + ; write date and region above the title + xyouts, 0.5, 0.95, 'Cloud mask ' + ' ' + self->get_date(/nice) $ + +' Cut ' + makes(self.cut,self.cut,1,di=2), $ + color=fsc_color('black'), /normal, $ + alignment=0.5, charsize=1.0 + + ; close file for each group + self->endDevice + if ~keyword_set(ps) then self->makeJPG, crop=crop + endif + + ENDFOR + + ; multiplot for all groups + if not keyword_set(single) then begin $ + ; write write date and region + xyouts, 0.5, 0.95, 'Cloud mask ' + ' ' + self->get_date(/nice) $ + +' Cut ' + makes(self.cut,self.cut,1,di=2) + addtitle, $ + color=fsc_color('black'), /normal, $ + alignment=0.5 ,charsize=1.5 + + ; close file + self->endDevice + if ~keyword_set(ps) then self->makeJPG + endif + + ; ----------- + ; create plot for disagreement + if group eq 'all' or group eq 'ALL' then begin + if not quiet then print,'... create agreement plot' + !p.multi=0 + !x.margin = 0.5 + !y.margin = [0.05,0.1] + !y.omargin = [0.2,0.2] + charsize = 0.97 + ;self->startDevice, xscale=1.3, yscale=1.3, n1 = 'AGR', addn=addName ;, addname='_MATCH_PERC' + + ; maximal number of disagreements + n_max = float(floor(0.5*obs)) + + ; plot image + view2d, float((obs-cld) < cld) / n_max, /cool, /colo, no_data_idx=where(cld eq 0),$ + title='Cloud mask disagreement', $ + coast_vec=size(coast,/type) eq 10 ? coast : 0, coast_col=0, aspect=1, $ + /no_axes, min=0., max=1., bar_nlev=5,$ + bardist=-0.15, /discrete, barformat='(F3.1)',charsize=charsize, bartitle='disagreement', $ + space_idx = where(self->lon() lt -200), space_col = !p.background, $ + col_table = 15, position=[0.1, 0.0, 1.0, 1.0] ; col_table = 16, + + + ; self->endDevice + ; if ~keyword_set(ps) then self->makeJPG, addn='_MATCH_DISAGR', crop='780x660+130+160' + endif + ; ----------- + + + +end diff --git a/cws_stats_do_it.pro b/cws_stats_do_it.pro new file mode 100644 index 0000000000000000000000000000000000000000..3eefa953886e41550123f8c9b08c4e7582f352de --- /dev/null +++ b/cws_stats_do_it.pro @@ -0,0 +1,111 @@ +pro cws_stats_do_it $ + ,cut=cut $ + ,product=product $ + ,com=com $ + ,kind=kind $ + ,hour = hour $ + ,overpass = overpass $ + , group = group $ + , month = month $ + , day = day $ + ,_extra=_extra + + error_n=0 + + if n_elements(hour) eq 0 and n_elements(cut) eq 1 then begin + case cut of + 2: hour = 1230 + 3: hour = 1345 + 4: hour = 1345 + 5: hour = 1345 + 6: hour = 1400 + else: hour=1200 + endcase + endif + + default, cut, 0 + default, product, 'cod' + default, kind, 1 + + default, day, 13 + default, month, 6 + default, year, 2008 + default, hour, 1200 + + ;if total(product eq ['cmb','ovw','rgb'] ) eq 1 then myCom=0 else myCom = com + o=obj_new('cws_stats', cut=cut , com=com, _extra=_extra) + if (n_elements(group) eq 0) then group='AWG' + + IF product eq 'rgb' and ( not between(kind,3,4)) THEN RETURN + IF product eq 'cmb' and ( not between(kind,3,4)) THEN RETURN + IF product eq 'cph' and ( not between(kind,3,4)) THEN RETURN + ; CATCH,error + ; IF error ne 0 THEN BEGIN + ; error_n ++ + ; print,'Fehleranzahl: ',error_n + ; RETURN + ; ENDIF + + hr = strmid(string(hour,format='(i4)'),0,2) + mi = strmid(string(hour,format='(i4)'),2,2) + if total(mi eq ['00','15','30','45']) ne 1 then mi = '00' + + o->set_date,year,month,day,hr,mi + + thisX=!x + thisY=!y + thisP=!p + print,'kind: ',kind + CASE kind OF + + 1: o->make_corThumb, product=product, _extra=_extra, hour=hour ; in cws_stats__define + 2: o->make_corSingle, product=product, hour=hour + 3: o->make_image, product=product, hour=hour + 4: o->make_image, /single, product=product, hour=hour + 5: o->hist1d, product=product, hour=hour + 6: o->diff_maps, group='AWG', product=product, hour=hour + 7: o->make_tables, product=product, hour=hour + 8: o->calipso_cpr, product=product, overpass=overpass + 9: begin + if group eq 'ALL' then $ + group =['CMS','EUM','OCA','MPF','FUB','DLR','MFR','AWG','UKM','GSF','LAR'] + for tt=0,n_elements(group) -1 do $ + o->cpr_profile,group=group[tt],cut=cut,_extra=_extra + end + 10: o->cod_ctp , group = group + 11: o->cod_ctp , group = group + 12: begin + if group eq 'ALL' then begin + group =['CMS','EUM','OCA','MPF','FUB','DLR','MFR','AWG','UKM','GSF','LAR'] + for tt=0,n_elements(group) -1 do $ + o->caliop_profile, group=group[tt], cut=cut, _extra=_extra + endif else o->caliop_profile, cut=cut, _extra=_extra + end + 13: o->cth_all, cut=cut + 14: o->lwp_all, cut=cut + 15: o->cod_all, cut=cut + 16: o->ref_all, cut=cut + 17: o->brian_profile, cut=cut + 18: o->brian_lwp, cut=cut + 19: begin + if group eq 'ALL' then $ + group =['CMS','EUM','OCA','MPF','FUB','DLR','MFR','AWG','UKM','GSF','LAR'] + for tt=0,n_elements(group) -1 do $ + o->cpr_profile_stats,group=group[tt],cut=cut,_extra=_extra + end + 20: o->cph_all,cut =cut + 21: o->all_one_image,product=product,cut = cut + 22: o->cloudsat_profile_lonlat,cut=cut + 23: o->taylor, group='AWG', product='CTH', atrainSensor='CPR' ; atrainSensor='AMSR'/'CALIOP'/'CPR'/'MODIS' ('ECMWF') + else: begin + print, '*** Error in do_it (cws_stats__define.pro), unknown kind ', kind,' of plot (1-23)' + end + + ENDCASE + + !x=thisX + !y=thisY + !p=thisP + catch,/cancel + +end diff --git a/cws_stats_params.pro b/cws_stats_params.pro index 5338f258ef5637d18964e5fed46456e1283eac42..105964f30eb39c9a99a139274bcee5778260e84d 100755 --- a/cws_stats_params.pro +++ b/cws_stats_params.pro @@ -394,6 +394,9 @@ 'ctp2': guteInd =[ 2 ] endcase + + + nightInd = [2,3,5,8,9] ;pmultImg = [0,4,3] ; org @@ -497,7 +500,7 @@ if not keyword_set(add_COX) then begin if strlowcase(self.phase) ne 'ice' and strlowcase(self.product) ne 'rgb' and strlowcase(self.product) ne 'ovw' then begin if not keyword_set(quiet) and not keyword_set(quiet2) then $ - print, '... remove COX group for non icy plots, phase=', self.phase + i_COX = 12 ; with FUB ; i_COX = 11 ; without FUB if n_elements(guteInd) gt 0 then $ diff --git a/icwg2_wrap_images.pro b/icwg2_wrap_images.pro index 2ff760308a6655191626a3e84eef8b5e8ecdd6f1..988139e6745870d1c0ff345ad710c519867d5cdc 100644 --- a/icwg2_wrap_images.pro +++ b/icwg2_wrap_images.pro @@ -1,9 +1,12 @@ pro icwg2_wrap_images o = cws_stats() + group_list = self->create_group_list(product='cmb', group=group, guteInd=guteInd, single=single, no_FUB=no_FUB) - o.make_image,product='ctp' + + + o.make_image,product='cmb',/single stop - o.make_corthumb,product='ctp' - + o.make_corthumb,product='cmb' +stop end