diff --git a/cws_read__define.pro b/cws_read__define.pro
index 2722e3ee10f4c4b824e4f2b76b2bf8e53423d8c2..cf1de1c88dc95ea111cf72fdfa65f90e46fa47f2 100755
--- a/cws_read__define.pro
+++ b/cws_read__define.pro
@@ -194,7 +194,7 @@ end
function cws_read::get_disks, group=group, product=product, ice=ice, surf=surf, $
phase=phase, img_file=img_file, log=log, n_disks=n_disks, dates=dates, mask=mask
- COMMON FEEDBACK, quiet, verbose, debug
+
default, n_disks, 1
@@ -211,6 +211,7 @@ function cws_read::get_disks, group=group, product=product, ice=ice, surf=surf,
; read new image
img = self->get_data(group=group, product=product, ice=ice, surf=surf, phase=phase, img_file=img_file, log=log, mask=mask)
+
if n_elements(img) eq 1 then begin
; create dummy with no data (-1) value everywhere
img = make_array(nx,/long,value=-1) # make_array(ny,/long,value=1)
diff --git a/cws_stats__constants.pro b/cws_stats__constants.pro
index 7e6f6772e316415369c4cb8a1435860fd915416c..df67f8e70e39f3f9bb425758f58c6ad9609412ac 100644
--- a/cws_stats__constants.pro
+++ b/cws_stats__constants.pro
@@ -12,6 +12,7 @@ reg = orderedhash()
; group list
grpList =['CMS','EUM','OCA','MPF','FUB','DLR','MFR','RMB','AWG','UKM','GSF','LAR','COX','OCA2','KNM','TPS','SUI','CLV','LARN','MFRN']
has_ctp = [1 ,1 ,1 ,1 ,0 ,1 ,0 ,0 ,1 ,1 ,0 ,1 ,0 ,0, 0 ,1 ,1 ,1 ,1, 1]
+ has_ctt = [ 1 , 1 ,1 ,1 ,0 ,1 ,0 ,0 ,1 ,1 ,0 ,1 ,0 ,0, 0 ,1 ,0 ,1 ,1, 1]
;has_ctp = [0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,1 ,0 ,0 ,1 ,0 ,0, 0 ,0 ,0 ,1 ,1]
has_cmb = [1,1,1,1,0,0,1,0,1,1,0,0,0,0,0,1,1,1,1,1]
has_cod = [1,0,1,0,0,0,0,1,1,1,1,1,0,0,0,1,0,1,1,1]
@@ -99,36 +100,37 @@ endfor
,'Uncertainty of Cloud Phase (2)','Uncertainty of Cloud Water Path (2)']
-for ii = 0,n_elements(prdList) -1 do begin
- prd[prdList[ii],'longname'] = prdListLong[ii]
-
+ for ii = 0,n_elements(prdList) -1 do begin
+ prd[prdList[ii],'longname'] = prdListLong[ii]
+ endfor
-endfor
-for ii = 0,n_elements(version) -1 do prd['ctp','groups',grplist[ii],'version'] = version[ii]
-for ii = 0,n_elements(has_cod) -1 do prd['cod','groups',grplist[ii]] = has_cod[ii]
+ for ii = 0,n_elements(version) -1 do prd['ctp','groups',grplist[ii],'version'] = version[ii]
+ for ii = 0,n_elements(has_cod) -1 do prd['cod','groups',grplist[ii]] = has_cod[ii]
-for ii = 0,n_elements(has_ctp) -1 do begin
- prd['ctp','groups',grplist[ii]] = has_ctp[ii]
- print,ii,grpList[ii]
-endfor
+ for ii = 0,n_elements(has_ctp) -1 do begin
+ 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]
+ endfor
-for ii = 0,n_elements(has_cmb) -1 do begin
- prd['cmb','groups',grplist[ii]] = has_cmb[ii]
- print,ii,grpList[ii], has_cmb[ii]
-endfor
+ for ii = 0,n_elements(has_cmb) -1 do begin
+ prd['cmb','groups',grplist[ii]] = has_cmb[ii]
+ 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)
-prd['ctp','pmulti'] = [0,ceil((n_ctp+2)/3.),3]
-prd['cmb','pmulti'] = [0,ceil((n_ctp+2)/3.),3]
+ n_ctp = total(has_ctp)
+ n_cmb = total(has_cmb)
+ prd['ctp','pmulti'] = [0,ceil((n_ctp+2)/3.),3]
+ prd['cmb','pmulti'] = [0,ceil((n_ctp+2)/3.),3]
-all_hash = hash()
-all_hash['grp'] = grp
-all_hash['prd'] = prd
-all_hash['reg'] = reg
-self.coef = all_hash
+ all_hash = hash()
+ all_hash['grp'] = grp
+ all_hash['prd'] = prd
+ all_hash['reg'] = reg
+ self.coef = all_hash
end
diff --git a/cws_stats__cph_all.pro b/cws_stats__cph_all.pro
index 84387036043ce01c96e0855dd01344a872474223..7576add4a72acb841d2e095371b0deb2b2e5fe33 100644
--- a/cws_stats__cph_all.pro
+++ b/cws_stats__cph_all.pro
@@ -1,3 +1,5 @@
+; Review AW 2018/11/11
+;
pro cws_stats::cph_all $
, start $
, nr $
@@ -5,87 +7,93 @@ pro cws_stats::cph_all $
,cut = cut $
,noPS = noPS
-orbt_nr = time_2_orbit( $
+ orbt_nr = time_2_orbit( $
julday(self.month,self.day,self.year,self.hour,self.minute))
-default,overpass,orbt_nr
+ default,overpass,orbt_nr
-oCwsAtrain = obj_new('cws2',overpass=overpass)
+ oCwsAtrain = obj_new('cws2',overpass=overpass)
-lon = oCwsAtrain->get_data(prod='LON',/data)
-lat = oCwsAtrain->get_data(prod='LAT',/data)
-ydt = oCwsAtrain->get_data(pr='YDT',/data)
+ lon = oCwsAtrain->get_data(prod='LON',/data)
+ lat = oCwsAtrain->get_data(prod='LAT',/data)
+ ydt = oCwsAtrain->get_data(pr='YDT',/data)
-pxl = geo_2_msg(lon,lat)
+ pxl = geo_2_msg(lon,lat)
-radar = oCwsAtrain->get_data(sen='CPR', pro='RFL', /DATA)
+ radar = oCwsAtrain->get_data(sen='CPR', pro='RFL', /DATA)
-ooo = obj_new('msg_data_cl',xrange=[self.x_0,self.x_1],yrange=[self.y_0,self.y_1])
-mer_vec = ooo->meridians(/vec,interval = ceil((self.x_1-self.x_0) /100.))
-obj_destroy,ooo
+ ooo = obj_new('msg_data_cl',xrange=[self.x_0,self.x_1],yrange=[self.y_0,self.y_1])
+ mer_vec = ooo->meridians(/vec,interval = ceil((self.x_1-self.x_0) /100.))
+ obj_destroy,ooo
-idx = where (between((ydt mod (3600D*24D))/3600D $
+ idx = where (between((ydt mod (3600D*24D))/3600D $
,float(self.hour)+float(self.minute)/60. $
,float(self.hour)+(float(self.minute))/60.+0.3) $
and between(pxl.(0),self.x_0,self.x_1) $
and between(pxl.(1),self.y_0,self.y_1) ,cW)
-start = idx[0]
-nr_idx = max(idx)-min(idx)
-if nr_idx le 0 then stop
+ start = idx[0]
+ nr_idx = max(idx)-min(idx)
+ if nr_idx le 0 then stop
-o = obj_new('cws2',overpass=overpass)
-o->atrain::set_overpass,overpass,/allow_nonexi
+ o = obj_new('cws2',overpass=overpass)
+ o->atrain::set_overpass,overpass,/allow_nonexi
-ind = indgen(nr_idx)+start
+ ind = indgen(nr_idx)+start
-cph_mod = o->get_data(pro='CPH',sens='MODIS',/data)
+ cph_mod = o->get_data(pro='CPH',sens='MODIS',/data)
-time = o->get_data(pro='YDT',/DATA)
-time = time[ind]
-lon = o->get_data( pro='lon',/DATA)
-lat = o->get_data( pro='lat',/DATA)
+ time = o->get_data(pro='YDT',/DATA)
+ time = time[ind]
+ lon = o->get_data( pro='lon',/DATA)
+ lat = o->get_data( pro='lat',/DATA)
-cph_mod = cph_mod[ind]
+ cph_mod = cph_mod[ind]
-loadct,39
-outFile = !PROJECTS.cws_path+'/results/atrain/cloudsat/C'+string(cut,fo='(i2.2)')+'_cph_' $
+ loadct,39
+ outFile = !PROJECTS.cws_path+'/results/atrain/cloudsat/C'+string(cut,fo='(i2.2)')+'_cph_' $
+'all_' +string(overpass,format='(i5.5)')
-file_mkdir,file_dirname(outFile)
+ file_mkdir,file_dirname(outFile)
-grpList = ['CMS','EUM','OCA','MPF','FUB','DLR','MFR','AWG','UKM','GSF','LAR'] ; groups
-colors=['red','yellow','blue','khaki','green','brown','orange','olive' $
+ grpList = ['CMS','EUM','OCA','MPF','FUB','DLR','MFR','AWG','UKM','GSF','LAR'] ; groups
+ colors=['red','yellow','blue','khaki','green','brown','orange','olive' $
,'purple','pink','gray']
-if ~keyword_set(noPS) then startps,outname=outfile,xs=3,ys=1.2
-
+ ;if ~keyword_set(noPS) then startps,outname=outfile,xs=3,ys=1.2
+ p = plot(lat[ind],cph_mod)
+ p.title = 'Cloud Phase'
+ p.ytitle ='cloud phase '
+ p.xtitle = 'Latitude'
+
+
+ stop
-plot,lat[ind],cph_mod,/xstyle,xtitle='Latitude',charsize=0.8,/nodata $
+ plot,lat[ind],cph_mod,/xstyle,xtitle='Latitude',charsize=0.8,/nodata $
, yTitle = 'cloud phase --- [g/m!U2!N]',yRange=[0,6] $
, title='Cloud Phase'
-xyouts,0.05,0.96,'SEVIRI scene: '+self->get_date(/string),/normal,charsize=0.7
+ xyouts,0.05,0.96,'SEVIRI scene: '+self->get_date(/string),/normal,charsize=0.7
xyouts, 0.95,0.96,' CUT: '+string(cut,fo='(i2.2)') $
+' OVP: '+string(overpass,format='(i5.5)'),align=1.,/normal,charsize=0.7
-erg = fltarr(n_elements(grpList),n_elements(cph_mod))
+ erg = fltarr(n_elements(grpList),n_elements(cph_mod))
-for ii=0, n_elements(grpList) -1 do begin
+ for ii=0, n_elements(grpList) -1 do begin
- group = grpList[ii]
- cph = round(o->_get_data(sensor='MSG',affi=group,product='CPH'))
- cod = o->_get_data(sensor='MSG',affi=group,product='COD')
- ref = o->_get_data(sensor='MSG',affi=group,product='REF')
- lwp = o->_get_data(sensor='MSG',affi=group,product='LWP')
- lwp = lwp[ind]
- ref = ref[ind]
- cod = cod[ind]
- cph = cph[ind]
- erg[ii,*] = round(cph)
+ group = grpList[ii]
+ cph = round(o->_get_data(sensor='MSG',affi=group,product='CPH'))
+ cod = o->_get_data(sensor='MSG',affi=group,product='COD')
+ ref = o->_get_data(sensor='MSG',affi=group,product='REF')
+ lwp = o->_get_data(sensor='MSG',affi=group,product='LWP')
+ lwp = lwp[ind]
+ ref = ref[ind]
+ cod = cod[ind]
+ cph = cph[ind]
+ erg[ii,*] = round(cph)
;oplot,lat[ind],cth_affi,psym=1,color=fsc_color(colors[ii]),symsize=0.5
-endfor
+ endfor
; stop
@@ -128,13 +136,13 @@ endfor
xyouts,0.18,0.58,'33% - 66% Percentil',charthick=5.,color =t_c_i(201,201,201),/normal ,charsize=0.8
oplot,lat[ind],meani,color=fsc_color('black')
-
- if ~keyword_set(noPS) then begin
- endps
- spawn,'convert -density 200 '+outfile $
- +'.eps '+outfile+'.jpg'
- spawn, 'rm -f '+outfile+'.eps'
- endif
+ stop
+ ; if ~keyword_set(noPS) then begin
+ ; endps
+; spawn,'convert -density 200 '+outfile $
+; +'.eps '+outfile+'.jpg'
+; spawn, 'rm -f '+outfile+'.eps'
+; endif
end
diff --git a/cws_stats__hist1d.pro b/cws_stats__hist1d.pro
index e8b65455108b8d627462775b3154fee320329f8d..9903eea13cd8af21f50917c02b2d734f43d3cb44 100644
--- a/cws_stats__hist1d.pro
+++ b/cws_stats__hist1d.pro
@@ -92,7 +92,7 @@ pro cws_stats::hist1d, hour=hour, product=product, mask=mask, surf=surf, phase=p
;self->set_group, (*(self.grpList))[guteInd[gg mod n_elements(guteInd)]]
self->set_group, group
; read data
- img = self->get_disks(surf=surf, phase=phase, n_disks=n_disks)
+ img = self->get_data(surf=surf, phase=phase, group= group)
if self.product eq 'cmb' then img = img-1.
; skip groups with no data
@@ -107,6 +107,8 @@ pro cws_stats::hist1d, hour=hour, product=product, mask=mask, surf=surf, phase=p
endif
endif
+
+ stop
; ; calculate x_axis of the histograms
; x_axis_borders = binsHist[cut]*indgen((maxHist[cut])/binsHist[cut])
@@ -142,13 +144,13 @@ pro cws_stats::hist1d, hour=hour, product=product, mask=mask, surf=surf, phase=p
self->set_group, group
; read data
- img = self->get_disks(surf=surf, phase=phase, n_disks=n_disks, dates=dates)
+ img = self->get_data(surf=surf, phase=phase)
if self.product eq 'cmb' then img = img-1.
; skip groups with no data
if max(img) eq -1 then begin
- print, '*** Warning, no data for group: ', (*(self.grpList))[guteInd[gg]]
+ print, '*** Warning, no data for group: ', group
continue
endif
@@ -161,7 +163,7 @@ pro cws_stats::hist1d, hour=hour, product=product, mask=mask, surf=surf, phase=p
; !!! histogram returns an array of size nhist for cod, but
; !!! histogram returns an array of size nhist+1 for ctp
; hist[im,gg,*] = (histogram(img,binsize=binsHist[cut],min=minhist[cut],max=maxhist[cut]))/float(n_elements(img))
-
+
case self.product of
'ctp': begin
@@ -291,7 +293,7 @@ pro cws_stats::hist1d, hour=hour, product=product, mask=mask, surf=surf, phase=p
xtitle = self.product eq 'ref' ? textoidl('\mum') : bartitle
xtitle = (*(self.prdListLong))[WHERE(self.product EQ *(self.prdList))] + ' / ' + xtitle
- if self.product ne 'cth' and self.product ne 'ctp' then begin
+ if self.product ne 'cth' and self.product ne 'ctp' and self.product ne 'ctt' then begin
; plot axis (for all products except cth fraction is on y axis)
@@ -300,7 +302,7 @@ pro cws_stats::hist1d, hour=hour, product=product, mask=mask, surf=surf, phase=p
yrange=[0.000001, ylogHist[cut] ? 1 : yhistmax[im]]
p = plot (x_axis_borders, hist[im,0,*])
- p.yrange=yrange
+ p.yrange = yrange
p.xrange = xrange
p.xtitle = 'fraction of pixels'
p.ytitle = xtitle
@@ -312,6 +314,7 @@ pro cws_stats::hist1d, hour=hour, product=product, mask=mask, surf=surf, phase=p
endif else begin
if self.product eq 'ctp' then begin & yrange= [imgMax[cut], 70] & ylog=1 & endif
if self.product eq 'cth' then begin & yrange= [imgMin[cut], imgMax[cut]] & ylog=0 & endif
+ if self.product eq 'ctt' then begin & yrange= [imgMax[cut], imgMin[cut]] & ylog=0 & endif
if n_elements(histmax) eq 0 then histmax = ylogHist[cut] ? 1 : yhistmax[im]
if self.product eq 'ctp' then if im eq m1 then xrange=[0,0.15] else xrange=[0,0.05]
if n_elements(xrange) eq 0 then xrange=[0.000001, histmax]
@@ -348,13 +351,13 @@ pro cws_stats::hist1d, hour=hour, product=product, mask=mask, surf=surf, phase=p
p_list = list()
- ; oplot the histogram data
+
foreach group,group_list,gg DO BEGIN
print, '... oplot with ', group,' ', $
(self.coef)['grp',group,'color']
- if self.product ne 'cth' and self.product ne 'ctp' then begin
-
+ if self.product ne 'cth' and self.product ne 'ctp' and self.product ne 'ctt' then begin
+
p0 = plot(x_axis_center, hist[im,gg,*] ,/overplot, color = (self.coef)['grp',group,'color'] ,name=group )
p_list.add,p0
@@ -377,30 +380,7 @@ pro cws_stats::hist1d, hour=hour, product=product, mask=mask, surf=surf, phase=p
leg.position=[0.9,0.7]
- ; 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 $
-; ,box=0, /norm, pos = legendPos, thick=1.2
-;
-;
-; cws_legend,groups[(n_elem_gute/2):*] $
-; ,linestyle=0 $
-; ,cws_color=guteInd[(n_elem_gute/2):*]+2 $
-; ,box=0, /norm, pos = legendPos - [0.2,0], thick=1.2
-; ENDIF ELSE BEGIN
-;
-;
-;
-; cws_legend, groups $
-; ,linestyle=0 $
-; ,cws_color=guteInd+2 $
-; ,box=0, /norm, pos = legendPos,thick=1.2
-; ENDELSE
-
+
endif else begin
; all groups in one legend
; cws_legend,(*(self.grpList))[guteInd] $
@@ -412,11 +392,7 @@ pro cws_stats::hist1d, hour=hour, product=product, mask=mask, surf=surf, phase=p
endelse
- ; close device and convert to jpg
- ;self->endDevice
- ;self->makeJpg
- ;file_mkdir , self->build_result_filename(_extra=_extra, single=single, n_disks=n_disks)
p.save, self->build_result_filename(_extra=_extra, single=single, n_disks=n_disks)+addname+'.png',/transparent
print,'saved in :',self->build_result_filename(_extra=_extra, single=single, n_disks=n_disks)+addname+'.png'
diff --git a/cws_stats_params.pro b/cws_stats_params.pro
index 50e3933b18777d4693d7c3920dc93cd915d34415..55746c9c770d4d3458f521afac2f9e0e44a0b9c8 100755
--- a/cws_stats_params.pro
+++ b/cws_stats_params.pro
@@ -264,7 +264,7 @@
if n_elements(bartitle) eq 0 then bartitle='km'
maxHist = [18,18,18,18,18,18,18,18,18,18,18,18,18,18,18,18]
minHist = make_array(13,/long,value=0)
- binsHist= [1.0, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25]
+ binsHist= [1.0, 0.25, 0.25, 0.25, 1., 1., 0.25, 0.25]
yhistmax1 = [0.1, 0.7, 0.7, 0.1, 0.3, 0.2, 0.2, 0.6, 0.4, 0.2, 0.2, 0.3, 0.2, 0.2, 0.2] ; common mask
yhistmax2 = [0.1, 0.6, 0.1, 0.001, 0.3, 0.1, 0.1, 0.6, 0.3, 0.15, 0.1, 0.5, 0.1, 0.1, 0.1] ; individual mask
@@ -408,8 +408,10 @@
maxHist = [1000,950,1000,1000,500,1000,1000,1000,1000,1000,1000,1000,1000]
minHist = make_array(13,/double,value=47.863)
- ; binsHist = [50,50,50,50,50,50,50,50,50,50,50,50,50]
- binsHist = make_array(13,/double,value=0.06) ;!!! log scale of ctp!!! -> 0.06 is approx 1km
+ binsHist = [50,50,50,50,50,50,50,50,50,50,50,50,50]
+ binsHist = make_array(13,value=25.)
+ ;binsHist = [50,50,50,50,50,50,50,50,50,50,50,50,50]
+ ; binsHist = make_array(13,/double,value=0.06) ;!!! log scale of ctp!!! -> 0.06 is approx 1km
yhistmax1 = [0.07,0.2,0.6,0.1,0.7,0.3,0.15,0.5,0.65,0.35,0.15,0.6,0.15,0.15,0.15]
yhistmax2 = [0.07,0.2,0.6,0.1,0.7,0.3,0.15,0.5,0.65,0.35,0.15,0.6,0.15,0.15,0.15]
; corr
diff --git a/icwg2_wrap_images.pro b/icwg2_wrap_images.pro
index 988139e6745870d1c0ff345ad710c519867d5cdc..828318fb8d0194a25758083262b3751d92687dc5 100644
--- a/icwg2_wrap_images.pro
+++ b/icwg2_wrap_images.pro
@@ -1,6 +1,6 @@
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)
+ group_list = self->create_group_list(product='cmb', group=group, guteInd=guteInd, single=single, no_FUB=no_FUB)