cws_stats__caliop_profile.pro 12 KB
Newer Older
Andi Walther's avatar
Andi Walther committed
1
;
Andi Walther's avatar
Andi Walther committed
2
;  Oct 29 2018: AW  changed to direct graphic
Andi Walther's avatar
Andi Walther committed
3
;
Andi Walther's avatar
Andi Walther committed
4
5
6
7
pro cws_stats::caliop_profile, start , nr $
    ,overpass=overpass $ 
    ,cut = cut $
    ,noPS = noPS $
Andi Walther's avatar
Andi Walther committed
8
9
    , group = group $
    , cpr = cpr
Andi Walther's avatar
Andi Walther committed
10
    
Andi Walther's avatar
Andi Walther committed
11
12
  
  if not self.quiet then print, '... start caliop_profile'	
Andi Walther's avatar
Andi Walther committed
13

Andi Walther's avatar
Andi Walther committed
14
  default, noPs, 1
Andi Walther's avatar
Andi Walther committed
15
16
  ;noPs=0 
  default, group, 'AWG'
Andi Walther's avatar
Andi Walther committed
17
  if not self.quiet then print, '... group = ', group
Andi Walther's avatar
Andi Walther committed
18
  default, cut, 5
Andi Walther's avatar
Andi Walther committed
19
  if not self.quiet then print, '... cut = ', cut
Andi Walther's avatar
Andi Walther committed
20
21
22
23
24
25
    
  ; convert date to orbit number 
  orbt_nr = time_2_orbit( $
            julday(self.month,self.day,self.year,self.hour,self.minute))

  default, overpass, orbt_nr
Andi Walther's avatar
Andi Walther committed
26
  if not self.quiet then print, '... ATRAIN orbit number = ', orbt_nr
Andi Walther's avatar
Andi Walther committed
27
28
29
30
  
  ; create cws2 object
  oCwsAtrain = obj_new('cws2', overpass=overpass)

Andi Walther's avatar
Andi Walther committed
31
  if self.verbose then print, '... read lon, lat, ydt'
Andi Walther's avatar
Andi Walther committed
32
33
34
35
  lon = oCwsAtrain->get_data(prod='LON',/data)
  lat = oCwsAtrain->get_data(prod='LAT',/data)
  ydt = oCwsAtrain->get_data(pr='YDT',/data)

Andi Walther's avatar
Andi Walther committed
36
  if self.verbose then print, '... convert (lon,lat) to pixel numbers'
Andi Walther's avatar
Andi Walther committed
37
38
39
  pxl = geo_2_msg(lon,lat)
        
  lidar  = oCwsAtrain->get_data(sen='CALIOP', pro='TAB_0532', /DATA)
Andi Walther's avatar
Andi Walther committed
40
  
Andi Walther's avatar
Andi Walther committed
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
        
  ; create object msg_data_cl  (.../tools/subtools/msg/msg_data_cl__define.pro) 
  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
  ; grpList = ['CMS','EUM','OCA','MPF','FUB','DLR','MFR','AWG','UKM','GSF','LAR'] ;  groups

  ind = indgen(nr_idx)+start

Andi Walther's avatar
Andi Walther committed
63
  radar = o->get_data(pro='RFL',sens='CPR',/data)
Andi Walther's avatar
Andi Walther committed
64
65
  lidar = o->get_data(sen='CALIOP'     ,pro='TAB_0532',/DATA)
  height = o->get_data(sen='CALIOP'     ,pro='HGT',/DATA)
Andi Walther's avatar
Andi Walther committed
66
  sfcbin = o->get_data(sen='CALIOP'     ,pro='SHB',/DATA)
Andi Walther's avatar
Andi Walther committed
67
  sfcbin_cpr = o->get_data(sen='CPR'     ,pro='SHB',/DATA)
Andi Walther's avatar
Andi Walther committed
68
69
  time = o->get_data(pro='YDT',/DATA)
  cth_cal =  o->get_data(sen='CALIOP', prod='CTH',/data)
Andi Walther's avatar
Andi Walther committed
70
71
  cth_cpr =  o->get_data(sen='CPR', prod='CTH',/data)
  
Andi Walther's avatar
Andi Walther committed
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
  time = time[ind]

  lon  = o->get_data( pro='lon',/DATA)
  lat  = o->get_data( pro='lat',/DATA)

  cod_modis = o->get_data(sen='MODIS', prod='COD',/data)
  cod_modis = cod_modis[ind]

  ref_modis = o->get_data(sen='MODIS', prod='CPS_REF',/data)
  ref_modis = ref_modis[ind]

  lwp_amsr = o->get_data(sensor='AMSR', prod = 'LWP',/data)
  lwp_amsr = lwp_amsr[ind]



  cth_cal  = cth_cal[ind]
  lidar  = lidar [*,ind]
Andi Walther's avatar
Andi Walther committed
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
  
  
 
  
 
  
  
  
  
  
  
  if keyword_set(cpr) then begin
    radar  = radar [*,ind]
    sfcbin_cpr = sfcbin_cpr[ind]
    h = REVERSE(FINDGEN(99)*0.25) + 0.5 ; this just for quicklook height bins starting 500 m above surfcae
    l = LINDGEN(N_ELEMENTS(ind))

   
    zind =(99-FINDGEN(99))  #  (FLTARR(N_ELEMENTS(sfcbin_cpr))+1)
    zind =(FLTARR(99)+1) #  (sfcbin_cpr-2)  - 2 - zind
    ; horizontal index
    lind =(FLTARR(99)+1) # FINDGEN(N_ELEMENTS(sfcbin_cpr)) 

    r = TRANSPOSE(radar[zind,lind])
    
    
    values = findgen(70)-28.5
    values[0] = -90.
    values[1] = -29.
    values[2] = -28.
  levels = n_elements(values)

  ctable = COLORTABLE([[085,026,139],[135,206,255],[0,255,0]], NCOLORS = levels, /TRANSPOSE)
  ctable[*,0] = fsc_color('gray',/triple)
  ctable[*,1] = fsc_color('gray',/triple)
  ctable[*,2] = fsc_color('gray',/triple)
  h = reverse(h)
  
  
  
  endif  else begin ; CALIOP
Andi Walther's avatar
Andi Walther committed
131
  height = height[*,ind]
Andi Walther's avatar
Andi Walther committed
132
133
134
135
  
  
  
  
Andi Walther's avatar
Andi Walther committed
136
137
138
139
140
141

    
  h = REVERSE(FINDGEN(99)*0.25) + 0.5 ; this just for quicklook height bins starting 500 m above surfcae
  l = LINDGEN(N_ELEMENTS(ind))
  
  h = height[*,1]
Andi Walther's avatar
Andi Walther committed
142
143
for u=0,n_elements(sfcbin)-1 do  radar[sfcbin[u]:*,u]= -999.
  ;sfcbin = make_array(n_elements(ind),/long)  
Andi Walther's avatar
Andi Walther committed
144
145
146
147
148
149
150

    
   
  zind =(n_elements(h)-FINDGEN(n_elements(h)))  #  (FLTARR(N_ELEMENTS(ind))+1)
  ;zind =(FLTARR(99)+1) #  (sfcbin-2)  - 2 - zind
  ;horizontal index
  lind =(FLTARR(n_elements(h))+1) # FINDGEN(N_ELEMENTS(ind)) 
Andi Walther's avatar
Andi Walther committed
151
152
153
154
155
;  do everything for CALIOP
  r_caliop = TRANSPOSE(lidar[zind,lind])
  
  values = [-0.04,0.004,0.015,0.02,0.03,0.04,0.05,0.1,0.3]
  levels = n_elements(values)
Andi Walther's avatar
Andi Walther committed
156

Andi Walther's avatar
Andi Walther committed
157
158
159
160
161
162
  ctable = COLORTABLE([[085,026,139],[135,206,255],[0,255,0]], NCOLORS = levels, /TRANSPOSE)
  ctable[*,0] = fsc_color('gray',/triple)
   r= r_caliop/1000.
  
  endelse
  
Andi Walther's avatar
Andi Walther committed
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
   
  loadct,39
  outFile = !PROJECTS.cws_path+'/results/atrain/cloudsat/C'+string(cut,fo='(i2.2)')+'_cpr_profile_' $
        +group+'_' +string(overpass,format='(i5.5)')
  file_mkdir,file_dirname(outFile)

  data = o->_get_data(sensor='MSG',affi=group,product='CTP')
  cmb  = o->_get_data(sensor='MSG',affi=group,product='CMB')
  cmb = round(cmb[ind])
  cph  = o->_get_data(sensor='MSG',affi=group,product='CPH')
  cph = round(cph[ind])
  cod = o->_get_data(sensor='MSG',affi=group,product='COD')
  cod = cod[ind]
  ref = o->_get_data(sensor='MSG',affi=group,product='REF')
  ref = ref[ind]
  ctt = o->_get_data(sensor='MSG',affi=group,product='CTT')
  ctt = ctt[ind]
  lwp = o->_get_data(sensor='MSG',affi=group,product='LWP')
  lwp = lwp[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
    add_legend_lwp ='calc. from REF/COD'
  endif else add_legend_lwp = ''
  
  
  
  idx = where(cmb eq 1,cCmb)
  if cCmb gt 0 then cph[idx] = 4  
  
  if n_elements(data) le 1 then begin
    print,'no data'
    return
  endif
  
  idx = where(data le 50,c)
  if c gt 0 then data[idx] = !VALUES.F_NAN
  
  cth_affi = o->atrain::p2z(data[ind,*],ind)
  thisP = !p
  !p.charsize=0.8
  
Andi Walther's avatar
Andi Walther committed
206
207
  ;if ~keyword_set(noPS) then startps,outname=outfile,xs=2,ys=2
  
Andi Walther's avatar
Andi Walther committed
208
  w = window(dimension=[700,800])
Andi Walther's avatar
Andi Walther committed
209
210
  

Andi Walther's avatar
Andi Walther committed
211
  
Andi Walther's avatar
Andi Walther committed
212
213
    
    
Andi Walther's avatar
Andi Walther committed
214
     c= contour(r,l,reverse(h),/current,pos =[0.1,0.68,0.9,0.95], yrange=[0,20] $
Andi Walther's avatar
Andi Walther committed
215
216
217
218
219
220
221
222
223
224
225
226
      , c_value = values $
          ,c_color=ctable,/fill,/xstyle ) 
     
 
     c['axis0'].SHOWTEXT = 0
     c.title = group
     
     t1 = text ( 0.05,0.96 , 'SEVIRI scene: '+self->get_date(/string))
     t2 = text (  0.95,0.96,' CUT: '+string(cut,fo='(i2.2)') $
        +'     OVP: '+string(overpass,format='(i5.5)'),align=1. )
     
    
Andi Walther's avatar
up    
Andi Walther committed
227
228
229
    p01 = plot( cth_affi/1000.,'w+',/overplot,pos = [0.1,0.68,0.9,0.95],/xstyle ,sym_size=0.3)
    
    
Andi Walther's avatar
Andi Walther committed
230
231
232
233
234
235
236
237
238
239
240
241
    cloud_idx = where(cmb gt 1.5,n_idx,compl=compl)
    
    
   ; if n_idx gt 0 then  p00=plot(cloud_idx,make_array(n_idx,value=1.5),'b+',/overplot)
    
    
    p01 = plot(indgen(n_elements(l)),make_array(n_elements(l),value=4),'y+',/current  $
      ,pos =  [0.1,0.64,0.9,0.68],/xstyle,/nodata , yrange=[0,5])
    for ii = 0, 3 do (p01.axes)[ii].HIDE = 1
    
    no_cloud_idx = where(cmb eq 1,cNoCloud) 
    if  cNoCloud gt 0 then begin
242
      p010=plot(compl,make_array(cNoCloud,value=1),'y.',/overplot )  
Andi Walther's avatar
Andi Walther committed
243
244
    endif  
  
Andi Walther's avatar
Andi Walther committed
245
    phase_idx_wat = where(cph eq 1,cWat )
Andi Walther's avatar
Andi Walther committed
246
    if cWat gt 1 then begin
247
      p011=plot(phase_idx_wat,make_array(cWat,value=2),'g.',/overplot  )   
Andi Walther's avatar
Andi Walther committed
248
249
    endif 
    
Andi Walther's avatar
Andi Walther committed
250
    phase_idx_ice = where(cph eq 2,cIce )
Andi Walther's avatar
Andi Walther committed
251
    if  cIce gt 1 then begin
252
       p012=plot(phase_idx_ice,make_array(cIce,value=3),'r.',/overplot  )  
Andi Walther's avatar
Andi Walther committed
253
254
    endif
    
Andi Walther's avatar
Andi Walther committed
255
    phase_idx_mix = where(cph eq 3,cMix )
Andi Walther's avatar
Andi Walther committed
256
    if  cMix gt 1 then begin
257
      p013=plot(phase_idx_mix,make_array(cMix,value=4),'c.',/overplot  ) 
Andi Walther's avatar
Andi Walther committed
258
259
260
    endif
    
    
Andi Walther's avatar
Andi Walther committed
261
262
    thisP = !P
    !p.charsize /= 2.4
Andi Walther's avatar
Andi Walther committed
263
    
264
265
266
267
268
   ; t30 = text ( nr_idx+10, -0. , 'Cl' , color='blue',/data) 
    t31 = text ( 0.93,0.63, 'NoCl' , color='gold',font_size=7) 
    t32 = text ( 0.93,0.645, 'Wat' , color='green',font_size=7) 
    t33 = text ( 0.93,0.66, 'Ice' , color='red',font_size=7) 
   ; t34 = text ( nr_idx+10, -3.2 , 'Mix' , color='purple',/data) 
Andi Walther's avatar
Andi Walther committed
269
270
    
   
Andi Walther's avatar
Andi Walther committed
271
272
273
    !p = thisP
      
    cod_wat = cod * (cph eq 1) - 999.*(cph ne 1)
Andi Walther's avatar
Andi Walther committed
274
275
    
    
Andi Walther's avatar
up    
Andi Walther committed
276
277
    p1 = plot( indgen(n_elements(ind)),cod,/current,symbol ='.',/nodata, yrange=[1,160],/xstyle $
      ,/ylog)
Andi Walther's avatar
Andi Walther committed
278
279
280
281
    p1.position = [0.1,0.46,0.9,0.62]
   
    
    p1['axis1'].Title ='cod []'
Andi Walther's avatar
up    
Andi Walther committed
282
    p1['axis1'].Axis_Range =[0,160]
Andi Walther's avatar
Andi Walther committed
283
284
    p1['axis0'].SHOWTEXT = 0
    
Andi Walther's avatar
up    
Andi Walther committed
285
286
    ; plot COD
    
287
288
289
    p10 = plot ( indgen(n_elements(ind)),cod ,'r+' ,/overplot,sym_size=0.5,name='Ice')
    p11 = plot ( indgen(n_elements(ind)),cod_wat,'g+',/overplot,sym_size=0.5,name = 'Water')
    p12 = plot ( indgen(n_elements(ind)),cod_modis, 'b+',/overplot,sym_size=0.5,name='Modis or AMSR-E')
Andi Walther's avatar
Andi Walther committed
290
291
    
    t12 =text ( 0.5,0.627,'Cloud optical depth' , align = 0.5 )
Andi Walther's avatar
up    
Andi Walther committed
292
   
293
294
295
296
297
298
299
    t13 = text ( 0.1,0.624,'+ Ice ','r',font_size=7)
    t14 = text ( 0.18,0.624,'+ Water ','g',font_size=7)
    t15 = text ( 0.26,0.624,'+ Modis ','b',font_size=7)
    
    ;leg = legend(target=[p10,p11,p12],position= [0.8,0.3],/auto_text_color,linestyle = 6 )
    
   
Andi Walther's avatar
Andi Walther committed
300

Andi Walther's avatar
Andi Walther committed
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
   ; plot,indgen(n_elements(ind)),cod,pos=[0.1,0.46,0.9,0.62],/noerase $
;              , yRange=[0,50],psym=1 $
;              ;, xTickformat='NOLABEL'$   ;!!!!!!!!!!!!!!!!!
;              , yTitle = 'cod []' $
;              , symsize=0.3,/nodata,/xstyle, xrange=[0,n_elements(ind)-1]
;    xyouts,0.5,0.627,'Cloud optical depth',charsize=.6,/normal,align=0.5         
;    oplot,indgen(n_elements(ind)),cod,psym=1 $
;                 ,symsize=0.6,col=fsc_color('red')   
;
;    oplot,indgen(n_elements(ind)),cod_wat,psym=1 $
;                 ,symsize=0.6,col=fsc_color('green')  
;                 
;    oplot,indgen(n_elements(ind)),cod_modis,psym=1 $
;                 ,symsize=0.2,col=fsc_color('blue')              
;              
Andi Walther's avatar
Andi Walther committed
316
317
318
319
320
321
322
323
324
325
    ; !!!!!!!!!!!!!!!!!
    ; legend,['Water', ' Ice', ' MODIS'] $
    ;         , psym = [1,1,1],color=[fsc_color('green')$
    ;         , fsc_color('red'),0], charsize = 0.5    $
    ;         , box = 0, spacing = 0.3 ,/Fill , /clear
 
    ; lwp = (0.833*0.65 * ref*cod)
    if size(lwp_amsr,/tname) eq 'STRING' then lwp_amsr=make_array(size(lwp_amsr,/dim),val=0)

    
Andi Walther's avatar
Andi Walther committed
326
327
328
329
330
331
332
333
334
335
336
337
    wat = where(cph eq 1,n_wat)
    
    
    
    p2 = plot( lat[ind],(lwp_amsr>0)*1000. $
              , pos=[0.1,0.06,0.9,0.22] , /current,/xstyle,/nodata $
                , yrange =[0,300])
    p2['axis0'].TITLE =   'Latitude'
    p2['axis1'].TITLE =   'lwp [g/m!U2!N]' 
    
    
    t20= text (   0.5,0.227,'Liquid water path',align=0.5)    
338
339
    t21 = text ( 0.1,0.224,'+ MSG Water ','g',font_size=7)
    t22 = text ( 0.24,0.224,'+ AMSR-E  ','b',font_size=7)
Andi Walther's avatar
Andi Walther committed
340
341
342
343
344
345
346
347
348
349
350
351
352
353
    
    
   ; ax0 = Axis('x',axis_range= [0,n_elements(ind)-1],style = 0)
    
    if n_wat gt 0 then begin
      p20 = plot ( (indgen(n_elements(ind)))[wat],(lwp[wat]>0) $
          , 'g+' $
          , sym_size=0.5,/current, pos=[0.1,0.06,0.9,0.22] $
          , axis_style = 0 )
          
      p21 = plot((indgen(n_elements(ind)))[wat],(lwp_amsr[wat]>0)*1000. $
          , 'b+', sym_size=0.5, /current , pos=[0.1,0.06,0.9,0.22]  $
          , axis_style = 0 )
    
Andi Walther's avatar
Andi Walther committed
354
    endif
Andi Walther's avatar
Andi Walther committed
355
356
357
358
    
    if add_legend_lwp ne '' then t21 = text ( 0.1,0.01,'+ calculated from COD and REF ',font_size=0.3)
    if add_legend_lwp eq '' then t22 = text ( 0.1,0.01,'+ LWP provided by group ',font_size=0.3)
    
Andi Walther's avatar
up    
Andi Walther committed
359

Andi Walther's avatar
Andi Walther committed
360
361
362
363
364
 
    ref_idx = where(ref le 0,cRef)
    if cRef gt 0 then ref[ref_idx] = -999.
    ref_wat = ref * (cph eq 1) - 999.*(cph ne 1)
    ref_ice = ref * (cph ne 1) - 999.*(cph eq 1) 
Andi Walther's avatar
Andi Walther committed
365
366
367
368
369
370
    
    
    p3 = plot(indgen(n_elements(ind)),ref $
              , pos=[0.1,0.26,0.9,0.42] $
                , /current,/xstyle, yrange=[0,50],/nodata  $
                , xrange = [0,n_elements(ind)-1] )
Andi Walther's avatar
Andi Walther committed
371
       
Andi Walther's avatar
Andi Walther committed
372
373
374
375
376
377
378
379
380
381
382
383
384
    p3['axis1'].TITLE =   "ref [$\mu$m]"
    p3['axis0'].SHOWTEXT = 0             
    t30 =  text (  0.5,0.427,'Cloud effective radius', align = 0.5 ) 
    
    ref_modis = float(ref_modis)
    
    p30 = plot ( indgen(n_elements(ind)),ref_modis,'b+' $
                 ,sym_size=0.5,/overplot )
    p31 = plot ( indgen(n_elements(ind)),ref_ice,'r+' $
                 ,sym_size=0.5,/overplot )
    p32 = plot ( indgen(n_elements(ind)),ref_wat, 'g+' $
                 ,sym_size=0.5,/overplot )                          
              
Andi Walther's avatar
up    
Andi Walther committed
385
      
386
387
388
     
   t4 = text (0.05,0.01,'created on '+systime(),font_size=7)   
      
Andi Walther's avatar
Andi Walther committed
389
390
  w.save,outFile+'.png'
  print,outFile+'.png'
Andi Walther's avatar
Andi Walther committed
391
392
393
394
395
396

end