-
Andi Walther authoredAndi Walther authored
cws_stats__lwp_all.pro 4.52 KiB
pro cws_stats::lwp_all $
, start $
, nr $
,overpass=overpass $
,cut = cut $
,noPS = noPS
orbt_nr = time_2_orbit( $
julday(self.month,self.day,self.year,self.hour,self.minute))
default,overpass,orbt_nr
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)
pxl = geo_2_msg(lon,lat)
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
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
o = obj_new('cws2',overpass=overpass)
o->atrain::set_overpass,overpass,/allow_nonexi
ind = indgen(nr_idx)+start
lwp_ams = o->get_data(pro='LWP',sens='AMSR',/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)
lwp_ams = lwp_ams[ind]
loadct,39
outFile = !PROJECTS.cws_path+'/results/atrain/cloudsat/C'+string(cut,fo='(i2.2)')+'_lwp_' $
+'all_' +string(overpass,format='(i5.5)')
file_mkdir,file_dirname(outFile)
@LIST_GROUPS.pro
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
idx=where(lwp_ams le 0,c)
if c gt 0 then lwp_ams[idx] = !VALUES.F_NAN
plot,lat[ind],lwp_ams*1000.,/xstyle,xtitle='Latitude',charsize=0.8,/nodata $
, yTitle = 'lwp [g/m!U2!N]',yRange=[0,500] $
, title='Liquid Water Path'
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(lwp_ams))
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]
if size(lwp,/TNAME) eq 'STRING' or (max(lwp) le 0) then begin
lwp = (0.833*0.65 * ref*(cod<50)) * (round(cph) eq 1) $
+ (((cod<50)^(1/0.84))/0.065 ) * (round(cph) ne 1)
lwp >=0
endif
if ii eq 0 then ice = where(round(cph) ne 1,c)
if c gt 0 then begin
lwp[ice] = !values.f_nan
endif
erg[ii,*] = lwp
;oplot,lat[ind],cth_affi,psym=1,color=fsc_color(colors[ii]),symsize=0.5
endfor
; win,1
;plot,cth_cpr
meani = make_array(n_elements(lwp_ams),value=!values.f_nan)
for uu=1,n_elements(lwp_ams)-2 do begin
idx = where(erg[*,uu] gt 0,cc)
if cc eq 0 then continue
mini = min(erg[idx,uu])
maxi = max(erg[idx,uu])
sddevi = cc gt 1 ? stddev(erg[idx,uu]) : 0
meani[uu] = median(erg[idx,uu],/even)
percentiles = percentiles(erg[idx,uu],value=[0.33,0.66])
plots,[lat[ind[uu-1]],lat[ind[uu]],lat[ind[uu]],lat[ind[uu-1]],lat[ind[uu-1]]] $
,[mini,mini,maxi,maxi,mini],color=t_c_i(240,240,240),thick=5.,noclip=0.
plots,[lat[ind[uu-1]],lat[ind[uu]],lat[ind[uu]],lat[ind[uu-1]],lat[ind[uu-1]]] $
,[percentiles[0],percentiles[0],percentiles[1],percentiles[1],percentiles[0]] $
,color=t_c_i(201,201,201),thick=5.,noclip=0.
endfor
oplot,lat[ind],meani,color=fsc_color('black')
lwp_ams[ice] = !values.f_nan
oplot,lat[ind],lwp_ams*1000.,linestyle=0,color=fsc_color('red'),thick=2.5
legend,['AMSR-E','SEVIRI median']$
,color=[fsc_color('red'),fsc_color('black')] $
,linestyle=0,box=0,pos=[0.1,0.8],/normal,charsize=0.7
xyouts,0.18,0.63,'MIN/MAX',charthick=5.,color =t_c_i(241,241,241),/normal,charsize=0.8 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
end