Skip to content
Snippets Groups Projects
Commit 5440c0bc authored by Andi Walther's avatar Andi Walther :speech_balloon:
Browse files

update tools

parent b18eb50f
No related branches found
No related tags found
No related merge requests found
pro aw_view2d, data_in ,_extra = _extra $
, no_data_idx = no_data_idx $
, space_idx = space_idx $
, jpg = jpg $
, bartitle = bartitle $
, xrange = xrange $
......@@ -8,19 +9,25 @@ pro aw_view2d, data_in ,_extra = _extra $
, coast_vec = coast_vec $
, buffer = buffer $
, layout = layout $
, no_new = no_new $
, no_save = no_save $
, nr_images = nr_images ; msg coaast point
, color_bar = color_bar $
, cmb = cmb $
, noax = noax $
, w =w ; msg coast point
default,nr_images,1
default,bartitle,''
default, no_new, 'false'
default, no_save,'false'
default,position, [0.25,0.1,0.95,0.9]
default,position, [0.2,0.1,0.95,0.9]
data = data_in
if n_elements(no_data_idx) gt 0 then begin
data[no_data_idx] = !VALUES.F_NAN
endif
endif
if n_elements(color_bar) eq 1 then begin
endif else begin
position = [0.05,0.05,0.95,0.9]
endelse
n_dim = size(data_in,/dim)
xra = findgen(n_dim[0])
......@@ -40,15 +47,11 @@ pro aw_view2d, data_in ,_extra = _extra $
a_Y = yrange[0]
b_Y = (yrange[1]-yrange[0])/double(n_dim[1])
endif
aspect_ratio = n_dim[0]/float(n_dim[1])
ct = COLORTABLE(74, /REVERSE)
; - get color table. I will later add a call to a bw_to_color function
ct = [ [ 000, 000, 125, 000 ], $ ;peter cool von blau nach rot
[ 000, 000, 255, 025 ], $
......@@ -69,27 +72,32 @@ pro aw_view2d, data_in ,_extra = _extra $
int_x = reform(ct[3, *])
n_col = 255l
col_v = 255. * findgen(n_col) / (float(n_col) - 1.)
int_r = (round(interpol(int_r, int_x, col_v)))
int_g = (round(interpol(int_g, int_x, col_v)))
int_b = (round(interpol(int_b, int_x, col_v)))
ctt = [[int_r],[int_g],[int_b]]
col_v = 255. * findgen(n_col) / (float(n_col) - 1.)
int_r = (round(interpol(int_r, int_x, col_v)))
int_g = (round(interpol(int_g, int_x, col_v)))
int_b = (round(interpol(int_b, int_x, col_v)))
ctt = [[int_r],[int_g],[int_b]]
w = window(dimensions=[600,400],buffer=buffer)
if keyword_set(cmb) then begin
ctt = ['white','green','blue','ivory']
data = bytscl(data,top=3,min=0,max=3)
endif
;ct[0,*] = cgcolor('gray',/triple)
if ~no_new then w = window(dimensions=[600,400],buffer=buffer)
undefine,position
print,'fff',layout,no_new
p = image(data,rgb_table=ctt,_extra = _extra, position = position $
, font_size = 13.,background_color='gray',aspect_ratio=aspect_ratio,/current,layout = layout )
p.axis_style = 3
xax = AXIS('X', LOCATION='bottom', TICKDIR=0, MINOR=0,coord_TRANS=[a_x,b_x],target = p)
yax = AXIS('Y', LOCATION='left', TICKDIR=0, MINOR=1,coord_TRANS=[a_y,b_y],target = p)
print,'ggg'
p = image(data,rgb_table=ctt,_extra = _extra, position = position $
, font_size = 13.,background_color='gray',aspect_ratio=aspect_ratio,/current,layout = layout,min_value=0 )
if ~keyword_set(noax) then begin
p.axis_style = 3
xax = AXIS('X', LOCATION='bottom', TICKDIR=0, MINOR=0,coord_TRANS=[a_x,b_x],target = p)
yax = AXIS('Y', LOCATION='left', TICKDIR=0, MINOR=1,coord_TRANS=[a_y,b_y],target = p)
endif
; make coast
IF keyword_set(coast_vec) THEN BEGIN
if keyword_set(coast_vec) then begin
coast_col = 'black'
FOR i = 0l, n_elements(coast_vec) -1l do begin
print,i
......@@ -99,32 +107,22 @@ pro aw_view2d, data_in ,_extra = _extra $
endif
; compute colorbar position
;
bar_size_x = (position[2] -position[0])/10.
bar_size_y = (position[3] -position[1])
bar_pos_x = (position[0] - bar_size_x - 0.1) > 0.01
bar_pos_y = position[1] > 0.01
pos_bar = [bar_pos_x,bar_pos_y,bar_pos_x+bar_size_x, bar_pos_y+bar_size_y ]
print,pos_bar
if n_elements(color_bar) eq 1 then begin
; compute colorbar position
bar_size_x = (position[2] -position[0])/10.
bar_size_y = (position[3] -position[1])
bar_pos_x = (position[0] - bar_size_x - 0.1) > 0.1
bar_pos_y = position[1] > 0.01
pos_bar = [bar_pos_x,bar_pos_y,bar_pos_x+bar_size_x, bar_pos_y+bar_size_y ]
;
c = colorbar ( target = p , ORIENTATION=1, $
POSITION=pos_bar, $
TITLE= bartitle)
;p.xrange = [0,30]
if n_elements(jpg) eq 1 then p.save,jpg,/JPG
;p.title='KKK'
;p.colorbar
print,'===> ', no_save
if no_save eq 'false' then begin
w.save,'test.png',resolution = 200
print,'save'
endif
c = colorbar ( target = p , ORIENTATION=1, $
POSITION=pos_bar, $
TITLE= bartitle)
endif
w.save,'test.png',resolution = 200
end
......@@ -149,10 +147,22 @@ o=msg_data_cl(xrange=[1600,2400],yrange=[2900,3500])
coast=o->coast_line(/vec)
lon = o.lon()
;view2d_max,lon,coast_v = coast,no_data_idx = where(lon lt -200),/cool,/colo
;aw_view2d,lon,coast_v = coast;,no_data_idx = where(lon lt -200)
aw_view2d,lon,layout=[3,3,1],/no_save
aw_view2d,lon,layout=[3,3,2],no_data_idx = where(lon lt -0),/no_new,/no_save
aw_view2d,lon,layout=[3,3,3],no_data_idx = where(lon lt -10),/no_new,/no_save
aw_view2d,lon,layout=[3,3,4],no_data_idx = where(lon lt 10),/no_new
;aw_view2d,lon,coast_v = coast,no_data_idx = where(lon lt -200)
; make multiplot 3x3
ww =window(dimens=[500,500])
for ii =1,9 do begin
aw_view2d,lon,no_data_idx = where(lon lt (ii*2)), w = w,/buffer
rr = w.CopyWindow(BORDER=0)
ww.setCurrent
pp = image(rr,layout = [3,3,ii], /current )
endfor
stop
end
......@@ -111,7 +111,7 @@ pro cws_stats::constants
prd['ctp','groups',grplist[ii]] = has_ctp[ii]
prd['cth','groups',grplist[ii]] = has_ctp[ii]
prd['ctt','groups',grplist[ii]] = has_ctt[ii]
print,ii,grpList[ii]
; print,ii,grpList[ii]
endfor
for ii = 0,n_elements(has_cmb) -1 do begin
......@@ -119,7 +119,7 @@ pro cws_stats::constants
print,ii,grpList[ii], has_cmb[ii]
endfor
print,(prd['cmb','groups',grplist]).where(1)
; print,(prd['cmb','groups',grplist]).where(1)
n_ctp = total(has_ctp)
n_cmb = total(has_cmb)
......
......@@ -167,16 +167,19 @@ pro cws_stats::hist1d, hour=hour, product=product, mask=mask, surf=surf, phase=p
case self.product of
'ctp': begin
histo = ( histogram(img, binsize=binsHist[cut], min=minhist[cut], max=maxhist[cut], LOCATIONS=locations)) / float(n_elements(img))
histo = ( histogram(img, binsize=binsHist[cut], min=minhist[cut] $
, max=maxhist[cut], LOCATIONS=locations)) / float(n_elements(img))
end
'cod': begin
histo = (log_histogram(img, binsize=binsHist[cut], min=minhist[cut], max=maxhist[cut], LOCATIONS=locations)) / float(n_elements(img))
histo = (log_histogram(img, binsize=binsHist[cut], min=minhist[cut] $
, max=maxhist[cut], LOCATIONS=locations)) / float(n_elements(img))
end
else: begin
histo = (histogram(img, binsize=binsHist[cut], min=minhist[cut], max=maxhist[cut], LOCATIONS=locations)) / float(n_elements(img))
histo = (histogram(img, binsize=binsHist[cut], min=minhist[cut] $
, max=maxhist[cut], LOCATIONS=locations)) / float(n_elements(img))
end
endcase
......
......@@ -11,9 +11,10 @@ pro cws_stats::make_corThumb, product=product, hour=hour, phase=phase, surf=surf
if n_elements(phase) eq 0 then phase = 'all'
if n_elements(surf) eq 0 then surf = 'all'
cd, !PROJECTS.cws_path+'tools/', current =current
@./cws_stats_params
cd,current
print,minCorr[self.cut]
print,maxCorr[self.cut]
cut = self.cut
self.kind = 'cor'
......@@ -23,7 +24,7 @@ pro cws_stats::make_corThumb, product=product, hour=hour, phase=phase, surf=surf
minCorr[cut] = 10.^minCorr[cut]
binCorr[cut] = 10.^binCorr[cut]
endif
print, "BULLSHIT COD: ", minCorr[cut], maxCorr[cut], binCorr[cut]
print, "B COD: ", minCorr[cut], maxCorr[cut], binCorr[cut]
valueRange = (maxCorr[cut]-minCorr[cut])
......@@ -91,7 +92,7 @@ pro cws_stats::make_corThumb, product=product, hour=hour, phase=phase, surf=surf
!p.multi=[0,n_elements(group_list),n_elements(group_list)]
Window,3,xsize=1400,ysize=1400
;Window,3,xsize=1400,ysize=1400
p_preserve=!p
......@@ -100,8 +101,9 @@ pro cws_stats::make_corThumb, product=product, hour=hour, phase=phase, surf=surf
flag_axes = 0
dumFlag = 0
group_list = self.create_group_list(product = product )
n_groups = n_elements(group_list)
w00 = window(DIMENSION=[900,900],/buffer)
; w = window(xsize=800,ysize= 900)
foreach group1,group_list, gg do begin
; read the data
......@@ -117,8 +119,8 @@ pro cws_stats::make_corThumb, product=product, hour=hour, phase=phase, surf=surf
foreach group2,group_list, tt do begin
print, ' reference group: ', group1, gg, n_elements(group_list) -1
print, ' comparison group: ', group2, tt, n_elements(group_list) -1
print, ' reference / comparison group: ', group1, group2, gg, tt
; read data of the 2nd group
img2 = self->get_data(group=group2, phase=phase, surf=surf)
if self.product eq 'cod' then begin
......@@ -146,24 +148,30 @@ pro cws_stats::make_corThumb, product=product, hour=hour, phase=phase, surf=surf
dum = make_array((dum=(maxCorr[self.cut]-minCorr[self.cut])/binCorr[self.cut]),dum,value=0.)
dum = dum*0
dum[0]=1
view2d,100.*dum,/cool,no_data_idx=where(dum eq 0) $
, range=[minCorr[cut],maxCorr[cut],minCorr[cut],maxCorr[cut]] $
, /no_axes ,aspect=1,charthick=2.3,charsize=thumb_charsize[cut] $
, min=0,max=maxCorrProc[cut],no_data_col = fsc_color('khaki')
xyouts, minCorr[cut] + 0.5*valueRange, minCorr[cut] +0.85 *valueRange $
, ' '+ group2 +' ' $
+ group1 +'!X' $
,charsize=1.2*thumb_charsize[cut] ,color=fsc_color('black') $
,align = 0.5
dum =minCorr[self.cut] + ((maxCorr[self.cut]-minCorr[self.cut])/binCorr[self.cut])/2
xyouts,dum ,dum,'NO DATA',charsize=1.0,orientation=45,align=0.5 $
, color =fsc_color('gold')
imgFormat='(f5.1)'
aw_view2d,100.*dum,no_data_idx=where(dum eq 0), w = w_dum,/noax
t11 = text( 0.26,0.8 , ' '+ group2 +' '+group1+' !X',font_size=25)
t12 = text( 0.26,0.5 , ' NO DATA ',font_size=35)
rr = w_dum.CopyWindow(BORDER=0)
w00.setCurrent
pp = image(rr,layout = [n_elements(group_list),n_elements(group_list),gg+1], /current )
w_dum.close
continue
ENDIF
; calculate statistics for existing data
; if self.product eq 'ctp' then begin
; temp_max = maxCorr[cut]
; temp_min = minCorr[cut]
; maxCorr[cut] = temp_min
; minCorr[cut] = temp_max
; endif
dum = hist2d(img2[ind],img1[ind],min=minCorr[self.cut],max=maxCorr[self.cut],bin=binCorr[self.cut])
corrFct = correlate(img2[ind],img1[ind])
biasFct = bias(img2[ind],img1[ind])
......@@ -180,40 +188,21 @@ pro cws_stats::make_corThumb, product=product, hour=hour, phase=phase, surf=surf
dum[0]=1
; plot the scatter plot
; work later on obect based view2d
imgFormat='(f5.1)'
view2d, 100.*dum,/cool,no_data_idx=where(dum eq 0) $
, range=[minCorr[cut],maxCorr[cut],minCorr[cut],maxCorr[cut]] $
, /no_axes ,aspect=1,charthick=2.3,charsize=thumb_charsize[cut] $
, min=0,max=maxCorrProc[cut],no_data_col = fsc_color('khaki')
;
xyouts, minCorr[cut] + 0.5*valueRange, minCorr[cut] +0.85 *valueRange $
, ' '+ group2+' ' $
+ group1 +'!X' $
,charsize=1.2*thumb_charsize[cut] ,color=fsc_color('black') $
,align = 0.5
xyouts, minCorr[cut] + 0.06*valueRange, minCorr[cut] +0.64 *valueRange $
, 'bias: ' + string(biasFct, format = '(f6.1)') ,charsize=thumb_charsize[cut] $
, color = abs(biasFct) gt 0.2 * valueRange ? cgcolor('red') : cgcolor('black')
xyouts, minCorr[cut] + 0.06*valueRange, minCorr[cut] +0.64 *valueRange $
, 'bias: ' ,charsize=thumb_charsize[cut] ,color = cgcolor('black')
xyouts, minCorr[cut] + 0.06*valueRange, minCorr[cut] +0.46 *valueRange $
, 'rmse: ' + string(rmseFct, format = '(f6.1)') ,charsize=thumb_charsize[cut] ,color = cgcolor('black')
xyouts, minCorr[cut] + 0.06*valueRange, minCorr[cut] +0.28 *valueRange $
, 'corr: ' + string(corrFct, format = '(f6.2)') ,charsize=thumb_charsize[cut] $
, color = corrFct gt 0.6 ? fsc_color('black') : fsc_color('red')
xyouts, minCorr[cut] + 0.06*valueRange, minCorr[cut] +0.28 *valueRange $
, 'corr: ' ,charsize=thumb_charsize[cut] ,color='black'
xyouts, minCorr[cut] + 0.06*valueRange, minCorr[cut] +0.10 *valueRange $
, '#: ' + string(smplFct,format='(i7)') ,charsize=thumb_charsize[cut] ,color = cgcolor('black')
aw_view2d,100.*dum,no_data_idx=where(dum eq 0), w = w_dum,/noax,/buffer
t11 = text( 0.26,0.9 , ' '+ group2 +' '+group1+' !X',font_size=38,/buffer)
t12 = text( 0.26,0.7,'bias: ' + string(biasFct, format = imgFormat), font_size=38)
t13 = text( 0.26,0.5,'rmse: ' + string(rmseFct, format = imgFormat), font_size=38)
t14 = text( 0.26,0.3,'corr: ' + string(corrFct, format = '(f6.2)'), font_size=38)
t15 = text( 0.26,0.1,'#: ' + string(smplFct, format = '(i7)'), font_size=38)
rr = w_dum.CopyWindow(BORDER=0)
w00.setCurrent
pp = image(rr,layout = [n_elements(group_list),n_elements(group_list),n_groups * gg + tt+ 1], /current )
w_dum.close
ENDIF
......@@ -221,35 +210,22 @@ pro cws_stats::make_corThumb, product=product, hour=hour, phase=phase, surf=surf
if tt eq gg then begin
dum = dum*0
dum[0]=1
view2d,100.*dum,/cool,no_data_idx=where(dum eq 0) $
, range=[minCorr[cut],maxCorr[cut],minCorr[cut],maxCorr[cut]] $
, /no_axes ,aspect=1,charthick=2.3,charsize=thumb_charsize[cut] $
, min=0,max=maxCorrProc[cut],no_data_col = fsc_color('orange')
plots, [minCorr[cut],minCorr[cut],maxCorr[cut],maxCorr[cut],minCorr[cut]] $
, [minCorr[cut],maxCorr[cut],maxCorr[cut],minCorr[cut],minCorr[cut]] ,color=fsc_color('blue')$
, psym=-3,thick=2.5
;
xyouts, minCorr[cut] + 0.5*valueRange, minCorr[cut] +0.85 *valueRange $
, ' '+ group2 +'!X' $
,charsize=1.2* thumb_charsize[cut] ,color=fsc_color('black') ,align=0.5
xyouts, minCorr[cut] + 0.06*valueRange, minCorr[cut] +0.64 *valueRange $
, 'mean: ' + string(meanFct, format = imgFormat) $
,charsize= thumb_charsize[cut] ,color=cgcolor('black')
xyouts, minCorr[cut] + 0.06*valueRange, minCorr[cut] +0.46 *valueRange $
, 'median: ' + string(mediFct, format = imgFormat) $
,charsize= thumb_charsize[cut],color=cgcolor('black')
xyouts, minCorr[cut] + 0.06*valueRange, minCorr[cut] +0.28 *valueRange $
, 'stdev: ' + string(stdevFct, format = imgFormat) $
,charsize= thumb_charsize[cut],color=cgcolor('black')
xyouts, minCorr[cut] + 0.06*valueRange, minCorr[cut] +0.10 *valueRange $
, '#: ' + string(smplFct,format='(i7)') ,charsize=thumb_charsize[cut] ,color=cgcolor('black')
aw_view2d,100.*dum,no_data_idx=where(dum eq 0), w = w_dum,/noax,/buffer
t11 = text( 0.22,0.9 , ' '+ group2 +'!X',font_size=38)
imgFormat='(f5.1)'
t12 = text( 0.22,0.7,'mean: ' + string(meanFct, format = imgFormat), font_size=38)
t13 = text( 0.22,0.5,'median: ' + string(mediFct, format = imgFormat), font_size=38)
t14 = text( 0.22,0.3,'stdev: ' + string(stdevFct, format = imgFormat), font_size=38)
t15 = text( 0.22,0.1,'#: ' + string(smplFct, format = '(i7)'), font_size=38)
rr = w_dum.CopyWindow(BORDER=0)
w00.setCurrent
pp = image(rr,layout = [n_elements(group_list),n_elements(group_list),n_groups * gg + tt+ 1], /current )
w_dum.close
endif
......@@ -262,34 +238,47 @@ pro cws_stats::make_corThumb, product=product, hour=hour, phase=phase, surf=surf
minCorr[cut] = temp_max
endif
if tt ne 0 and gg ne 0 and dumFlag eq 0 then flag_axes = 0 else flag_axes = 1
view2d,alog(100.*dum/float(total(dum))),/cool,no_data_idx=where(dum eq 0) $
; ,title= string(corrFct,format='(f5.2)') $
; + ' '+string(rmseFct,format='(f6.1)') $
, range=[minCorr[cut],maxCorr[cut],minCorr[cut],maxCorr[cut]] $
, no_axes = flag_axes $
, aspect=1,charsize = 1.5*thumb_charsize[cut] $
, min=alog(100./float(total(dum))),max=alog(maxCorrProc[cut])
plots,[0,0],[1000,1000],color=fsc_color('white')
if dumFlag eq 0 and flag_axes eq 0 then dumFlag = 1
if self.product eq 'ctp' then begin
xyouts, minCorr[cut] - 0.01*valueRange $
, maxCorr[cut] + 0.25*valueRange , $
' '+group2+' ' $
+ group1 +'!X' $
,charsize=1.2*thumb_charsize[cut],color=fsc_color('white') ,/data
endif else begin
xyouts, minCorr[cut] + 0.01*valueRange,minCorr[cut] +0.85 *valueRange, $
' '+group2+' ' $
+ group1 +'!X' $
,charsize=1.2*thumb_charsize[cut],color=fsc_color('white')
; IF self.comStyle eq 'IND' then xyouts, minCorr[cut] + 0.55*valueRange,minCorr[cut] +0.05 *valueRange $
; ,'# '+string(total(dum),format='(i7)'),charsize=0.5
aw_view2d,alog(100.*dum/float(total(dum))),no_data_idx=where(dum eq 0), w = w_dum,noax = (tt ne 2) or (gg ne 2) , /buffer $
, xrange=[minCorr[cut],maxCorr[cut]], yrange = [minCorr[cut],maxCorr[cut]] $
, min_value = alog(100./float(total(dum))) , max_value = alog(maxCorrProc[cut])
t11 = text ( 0.26, 0.7 , ' '+group2+' '+ group1 +' ',color= 'white', font_size = 38)
rr = w_dum.CopyWindow(BORDER=0)
w00.setCurrent
pp = image(rr,layout = [n_elements(group_list),n_elements(group_list),n_groups * gg + tt+ 1], /current )
w_dum.close
endelse
; view2d,alog(100.*dum/float(total(dum))),/cool,no_data_idx=where(dum eq 0) $
; ; ,title= string(corrFct,format='(f5.2)') $
; ; + ' '+string(rmseFct,format='(f6.1)') $
; , range=[minCorr[cut],maxCorr[cut],minCorr[cut],maxCorr[cut]] $
; , no_axes = flag_axes $
; , aspect=1,charsize = 1.5*thumb_charsize[cut] $
; , min=alog(100./float(total(dum))),max=alog(maxCorrProc[cut])
;
; plots,[0,0],[1000,1000],color=fsc_color('white')
;
; if dumFlag eq 0 and flag_axes eq 0 then dumFlag = 1
;
; if self.product eq 'ctp' then begin
; xyouts, minCorr[cut] - 0.01*valueRange $
; , maxCorr[cut] + 0.25*valueRange , $
; ' '+group2+' ' $
; + group1 +'!X' $
; ,charsize=1.2*thumb_charsize[cut],color=fsc_color('white') ,/data
; endif else begin
; xyouts, minCorr[cut] + 0.01*valueRange,minCorr[cut] +0.85 *valueRange, $
; ' '+group2+' ' $
; + group1 +'!X' $
; ,charsize=1.2*thumb_charsize[cut],color=fsc_color('white')
; ; IF self.comStyle eq 'IND' then xyouts, minCorr[cut] + 0.55*valueRange,minCorr[cut] +0.05 *valueRange $
; ; ,'# '+string(total(dum),format='(i7)'),charsize=0.5
;
; endelse
if self.product eq 'ctp' then begin
maxCorr[cut] = temp_max
......@@ -301,16 +290,19 @@ pro cws_stats::make_corThumb, product=product, hour=hour, phase=phase, surf=surf
idx_prd = where(product eq *self.prdList)
; write title to graph
xyouts,0.5,0.98,' '+self.coef['prd',product,'longname']+', '+self->get_date(/nice)+ $
addtitle+'!X',align=0.5,/Normal $
, charsize =1.5*thumb_charsize[cut]
w00.setCurrent
t33 = text(0.5,0.98,' '+self.coef['prd',product,'longname']+', '+self->get_date(/nice)+ $
addtitle+'!X' , align = 0.5 )
outname = self->build_result_filename(_extra=_extra, single=single, n_disks=n_disks)+addname
t12 = text(0.02,0.01,'created on '+systime()+' AW',font_size=7)
w00.save,outname+'.png'
print,'saved in',outname+'.png'
;xyouts,0.99,0.97,self->get_date(/nice),/Normal,charsize=1.5*thumb_charsize[cut],align=1
if self.cut ne 0 then xyouts,0.05,0.97,'Cut:'+strtrim(cut,1),/Normal,charsize=1.5*thumb_charsize[cut]
;if self.cut ne 0 then xyouts,0.05,0.97,'Cut:'+strtrim(cut,1),/Normal,charsize=1.5*thumb_charsize[cut]
; fobj=obj_new('msg_data_cl')
; fobj->set_date,self.year, self.month, self.day, self.hour, self.minute
......@@ -318,15 +310,7 @@ pro cws_stats::make_corThumb, product=product, hour=hour, phase=phase, surf=surf
; obj_destroy,fObj
; stop
!p=p_preserve
!x =x_def
!y = y_def
!p.multi=0
write_png,self->build_result_filename(_extra=_extra,single=single,n_disks=n_disks)+addname+'.png', tvrd(/true)
; cphFlag ++
; if phseFlag and (cphFlag lt 3) then begin
......
......@@ -5,244 +5,159 @@ pro cws_stats::make_image_cmb $
, group = group $
, phase = phase
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)
n_groups = n_elements(group_list)
; 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 = 1.4
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=200 * !p.multi[1] ,ysize=200 * !p.multi[2]
!x.omargin=1.
!y.omargin=[2.,6.]
!x.margin=0.9
!y.margin=[0.2,2.2]
charsize=2.2
;self->startDevice, xscale=2.1, yscale=1.8
endelse
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 = 1.4
endif else begin
!P.multi=self.coef['prd','cmb','pmulti']
w1 = window ( dimensions = [200 * ceil((n_groups)/3.),206*3] , location =[0.1,0.1],/buffer )
endelse
; read land sea mask
ls = self->land_sea()
; read land sea mask
ls = self->land_sea()
; get longitude
lon = self->lon()
; get longitude
lon = self->lon()
coast=self->coast_line(/vec)
; initialize counter for disagreements
cld = ls * 0
obs = cld
n_data = 0
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
for gg = 0, n_groups -1 do begin
; set the group
group_str = group_list[gg]
if group_str ne 'AVG' then begin
self->set_group, group_str
img = self->get_data(phase=phase, surf=surf )
endif
; skip plot if single plot is wanted and there is no data
if keyword_set(single) and n_elements(img) eq 1 then return
; set the group
;group_str = (*(self.grpList))[guteInd[gg]]
group_str = group_list[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
; 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
if n_elements(img) ne 1 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) )
mini=0
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
aw_view2d,dum, w = w ,/buffer $
, title = group_str + ' ('+makes(frc,frc,1,di=2) + '%)' $
, min_value = mini, font_size = 18,/cmb,/noax
rr = w.CopyWindow(BORDER=0)
w1.setCurrent
pp = image(rr,layout = [ceil((n_groups)/3.),3,gg+1], /current )
w.close
; 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)
; 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
; 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
aw_view2d,dum,w=w,/buffer $
, title = group_str+' ( NO DATA)' $
, min_value = mini, font_size = 18,/cmb,/noax
t1 = text(300, 300, 'NO DATA' $
, orientation = 45. , alignment = 0.5 $
, font_size = 40. ,color='red',/data)
rr = w.CopyWindow(BORDER=0)
w1.setCurrent
pp = image(rr,layout = [ceil((n_groups)/3.),3,gg+1], /current )
endelse
; 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
; write to png
endif
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
outname = self->build_result_filename(_extra=_extra, single=single, n_disks=n_disks)
png_transparent_bg,name=outname+'.png',tvrd(/true)
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
window ,2,xsize=600,ysize=600
;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,
ENDFOR
; multiplot for all groups
if not keyword_set(single) then begin
t11 = text ( 0.5,0.96,'Cloud Mask '+ self->get_date(/nice) $
+' Subregion ' + makes(self.cut,self.cut,1,di=2) + addtitle $
, color='black',/normal,alignment=0.5, font_size=16 )
t12 = text(0.02,0.01,'created on '+systime(),font_size=7)
outname = self->build_result_filename(_extra=_extra, single=single, n_disks=n_disks)
w1.save, outname+'.png',/transparency
print,'file saved: ',outname+'.png'
w1.close
endif
; -----------
; create plot for disagreement
if group eq 'all' or group eq 'ALL' then begin
if not self.quiet then print,'... create agreement plot'
n_max = float(floor(0.5*obs))
outname = self->build_result_filename(_extra=_extra, single=single, n_disks=n_disks,addn='_MATCH_PERC')
png_transparent_bg,name=outname+'.png',tvrd(/true)
; self->endDevice
; if ~keyword_set(ps) then self->makeJPG, addn='_MATCH_DISAGR', crop='780x660+130+160'
endif
; -----------
aw_view2d,float((obs-cld) < cld) / n_max, no_data_idx=where(cld eq 0),bartitle = 'disagreement' $
, w = w_dis, /color_bar,/buffer
t11 = text ( 0.5,0.96,'Cloud Mask '+ self->get_date(/nice) $
+ ' Subregion ' + makes(self.cut,self.cut,1,di=2) + addtitle $
, color='black',/normal,alignment=0.5, font_size=12 )
t12 = text(0.02,0.01,'created on '+systime(),font_size=7)
outname = self->build_result_filename(_extra=_extra, single=single, n_disks=n_disks,addn='_MATCH_PERC')
w_dis.save,outname+'.png'
print,'file saved: ',outname+'.png'
w_dis.close
endif
......
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)
for cc=0,6 do begin
if cc eq 1 then continue
for hh =0,23 do begin
o.make_image,product='cmb',/single
stop
o.make_corthumb,product='cmb'
for mi =00,45,15 do begin
hour = string(hh,mi,format='(2(i2.2))')
CWS_STATS_DO_IT,product='cmb',kind=3,cut=cc,hour = hour
endfor
endfor
endfor
stop
end
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment