;=========================================================================== ;comb_MODISdust_v6.3.4_GRIB.pro ;author: Jian Zeng (ERT) ;date: 2012-01-25 ;purpose: generate hourly MODIS dust AOD and concentration in GRIB format ;============================================================================ pro disp_aod, data, lats, lons,outfile,title=title,bartitle=bartitle,$ limits=limits,latdel=latdel,londel=londel ;sz=size(lats,/dim) ;x_size=sz(0) ;y_size=sz(1) map_limit=[min(lats),min(lons),max(lats),max(lons)] ; projection as MODIS rgb IMAGE ;center=fltarr(2) ;center[0]=lats(fix(x_size/2),fix(y_size/2)) ;center[1]=lons(fix(x_size/2),fix(y_size/2)) if not keyword_set(limits) then limits=map_limit MISVAL=-99999.0 xmag=8 ymag=8 sz=size(lats,/dim) data_=congrid(data,sz(0)*xmag,sz(1)*ymag, /center) lat_=congrid(lats,sz(0)*xmag,sz(1)*ymag, cubic=-0.5,/center,/interp) lon_=congrid(lons,sz(0)*xmag,sz(1)*ymag, cubic=-0.5,/center,/interp) bot=30 top=250 range=[0,1.2] ncolors=top-bot+1 img=bytscl(data_,min=range(0),max=range(1),top=ncolors-1)+bot ; liang if min(data(where(data gt -99998.0))) lt -201.0 then begin if min(data(where(data gt -999999.0))) lt -201.0 then begin ind=where(data_ lt 0.,cp) if cp gt 0 then img(ind)=254 ind00=where(lat_ gt -90 and lat_ lt 90. and lon_ gt -180 and lon_ lt 180.) endif else begin ;non-dust ; liang ind=where(data_ lt -100. and data_ gt -99998.0,cp) ind=where(data_ lt -100. and data_ gt -999999.0,cp) if cp gt 0 then img(ind)=254 ;dust flag with no AOD ind=where(data_ lt -99. and data_ gt -101.0,cp) if cp gt 0 then img(ind)=254 ;253 ind00=where(data_ gt -99998.0) endelse ;no coverage(missing) ; liang ind=where(data_ lt -99998.0,cp) ind=where(data_ lt -999999.0,cp) if cp gt 0 then img(ind)=0 set_plot,'Z' device,set_resolution=[900,750],z_buffer=0,SET_FONT='Times Bold',/TT_FONT !p.charthick=3 !p.charsize=1 loadct,39,/silent tvlct,r,g,b,/get r(255)=0 g(255)=0 b(255)=0 r(0)=255 g(0)=255 b(0)=255 r(254)=180 g(254)=180 b(254)=180 r(253)=250 g(253)=100 b(253)=250 tvlct,r,g,b ind00=where(lat_ gt -90 and lat_ lt 90 and lon_ gt -180 and lon_ lt 180) imagemap_unscl,img(ind00),lat_(ind00),lon_(ind00),$ /noborder,position=[0.05,0.18,0.95,0.9],$ limit=limits map_continents,/countries,/noerase,color=255 map_continents,/usa,/HIRES,/cont,/noerase,color=255 map_grid,color=255,/label,/box,latdel=latdel,londel=londel,chars=1.5 xyouts,0.05,0.95,title,charsize=2.8,color=255,charthick=4.0,/normal xyouts,0.83,0.011,'>',/norm,charsize=2,color=255 colorbar,position=[0.2,0.06,0.8,0.09],range=range,$ format='(f3.1)',div=6,bottom=bot,minor=1,$ chars=2,title=bartitle,ncolors=ncolors,color=255 make_png,outfile+'.png' set_plot,'x' print,'creating----',outfile+'.png' end ;-------------------------------------------------------------------------- pro disp_conc,data,lats,lons,outfile,title=title,bartitle=bartitle,$ limits=limits,latdel=latdel,londel=londel levels=[2,10,20,30,50] cols=[3,5,6,7,8] sunglt_col=4 cld_col=1 clr_col=14 sz=size(data,/dim) xsz=sz(0) ysz=sz(1) ind=where(data gt 0.,nd) if nd gt 0 then data(ind)=data(ind)*1.e9 img=bytarr(xsz,ysz)+0B ;default color=white FOR ilev=0,n_elements(levels)-1 DO BEGIN IF ilev EQ n_elements(levels)-1 THEN ind=where(data gt levels(ilev),ndx) ELSE ind=where(data gt levels(ilev) and data le levels(ilev+1),ndx) IF ndx GT 0 THEN img(ind)=cols(ilev) ENDFOR ;sunglint idx=where(abs(data+2) lt 1.e-6,ndx) if ndx gt 0 then img(idx)=sunglt_col ;cloud idx=where(abs(data+1) lt 1.e-6,ndx) if ndx gt 0 then img(idx)=cld_col ;clear idx=where(data lt -9998.0,ndx) if ndx gt 0 then img(idx)=clr_col xmag=4 ymag=4 img_=congrid(img,sz(0)*xmag,sz(1)*ymag, /center) lat_=congrid(lats,sz(0)*xmag,sz(1)*ymag, cubic=-0.5,/center,/interp) lon_=congrid(lons,sz(0)*xmag,sz(1)*ymag, cubic=-0.5,/center,/interp) set_plot,'Z' device,set_resolution=[900,750],z_buffer=0,SET_FONT='Times Bold',/TT_FONT !p.charthick=3 !p.charsize=1 TEK_color TVLCT, 0, 0, 0, 255 TVLCT, 255, 255, 255, 0 TVLCT, 220, 220, 220, 1 ind00=where(lat_ gt -90 and lat_ lt 90 and lon_ gt -180 and lon_ lt 180) imagemap_unscl,img_(ind00),lat_(ind00),lon_(ind00),$ /noborder,position=[0.05,0.18,0.95,0.9],$ limit=limits map_continents,/countries,/noerase,color=255 map_continents,/usa,/HIRES,/cont,/noerase,color=255 map_grid,color=255,/label,/box,latdel=latdel,londel=londel,chars=1.5 xyouts,0.05,0.95,title,charsize=2.8,color=255,charthick=4.0,/normal col_all=[clr_col,sunglt_col,cld_col,cols] bar_name=[' other',' sunglint',' cloud',string(levels,format='(i3)'),'>'] chars=1.5 & bar_w=0.1 hori_color_bar, nbar=n_elements(col_all),bar_w=bar_w,bar_h=0.02,$ ystart=0.07,bar_name=bar_name,chars=chars,$ bcolors=col_all, bar_title=bartitle,title_color=255,$ bartit_start=0.4,index='range' make_png,outfile+'.png' set_plot,'x' print,'creating----',outfile+'.png' END ;---------------------------------------------------------- @read_binary_all.pro @calc_latlon.pro @convert_julianday.pro @write_binary_all.pro @imagemap_unscl.pro @colorbar.pro @make_png.pro @hori_color_bar.pro PRO comb_MODISdust_V634_GRIB,Year,Day,SAT,ver=ver,plot_dust=plot_dust IF N_ELEMENTS(ver) EQ 0 THEN ver='v6.3.4' IF N_ELEMENTS(plot_dust) EQ 0 THEN plot_dust=1 print, 'in comb_MODISdust_V634_GRIB' ; dir0='/net/orbit184l/disk2/pub/zeng/MODIS_dust/DB/' dir0='/data/test/Dust/output/' convert_julianday,Year,Day,mon,day1 month=string(mon,format='(i2.2)') nday=string(day1,format='(i2.2)') nlat=251 & nlon=601 clat=25. & clon=-125. dlat=0.1 & dlon=0.1 calc_latlon,nlon,nlat,clat,clon,dlat,dlon,lon_center,lat_center,xl=xl,xr=xr,yb=yb,yt=yt yr=string(year,format='(i4.4)') jday=string(day,format='(i3.3)') convert_julianday,yr,jday,mon,day month=string(mon,format='(i2.2)') nday=string(day,format='(i2.2)') date=yr+month+nday dir=dir0+yr+month+nday+'/Output/' ; ftp_dir='/net/www/aftp/pub/smcd/spb/shobha/DUST/GRIB/' print, 'dirin com='+dir files=file_search(dir+SAT+'.?'+yr+jday+'*.dat',count=nfile) print,'files6666='+files IF nfile EQ 0 THEN print, 'no input files found for combining!' ELSE BEGIN pos=strpos(files,ver+'.mv1.dat') IF min(pos) EQ -1 THEN print,'wrong files!!!' ELSE BEGIN hrs=strmid(files,pos(0)-5,2) FOR ihr=16,22 DO BEGIN ind=where(hrs eq ihr,nhr) IF nhr gt 0 THEN BEGIN stamp=date+'.hr'+string(ihr,format='(i2.2)') IF nhr EQ 1 THEN BEGIN infile=files(ind(0)) read_binary_all,infile,nx,ny,lon,lat,aod,conc combaod=aod combconc=conc ENDIF ELSE BEGIN FOR ifl=0,nhr-1 DO BEGIN infile=files(ind(ifl)) print,infile read_binary_all,infile,nx,ny,lon,lat,aod,conc ;for aod ;non-dust pixels within the granule set to -200 ;domain not covered by the granule set to -99999.0 ;dust aod > 0 ;dust pixel but no aod retreival set to -100 ;for conc ;cloud within the granule set to -1 ;sunglint within the granule set to -2 ;domain not covered by the granule set to -3 ;non dust set to -9999.0 ;dust conc > 0 IF ifl EQ 0 THEN BEGIN combaod=aod combconc=conc ENDIF ELSE BEGIN ;overlapped scans for aod idx0=where(combaod gt -99998.0 and aod gt -99998.0,ndx0) if ndx0 gt 0 then begin idx=where(combaod gt 0.0 and aod gt 0.0,ndx) if ndx gt 0 then begin temp=combaod(idx)+aod(idx) combaod(idx)=temp/2 endif idx=where(combaod le 0.0 and combaod gt -99998.0 and aod gt 0.0,ndx) if ndx gt 0 then begin combaod(idx)=aod(idx) endif endif ;new scans for aod idx=where(combaod lt -99998.0 and aod gt -99998.0,ndx) if ndx gt 0 then combaod(idx)=aod(idx) ;overlapped scans for conc idx0=where(abs(combconc+3) gt 0.5 and abs(conc+3) gt 0.5,ndx0) if ndx0 gt 0 then begin idx=where(combconc gt 0.0 and conc gt 0.0,ndx) if ndx gt 0 then begin temp=combconc(idx)+conc(idx) combconc(idx)=temp/2 endif idx=where(combconc le 0.0 and abs(combconc+3) gt 0.5 and conc gt 0.0,ndx) if ndx gt 0 then begin combconc(idx)=conc(idx) endif idx=where((combconc lt -9998.0) and ((abs(conc+2) lt 1.e-6) or (abs(conc+1) lt 1.e-6)),ndx) if ndx gt 0 then begin combconc(idx)=conc(idx) endif endif ;new scans for conc idx=where(abs(combconc+3) lt 1.e-6 and abs(conc+3) gt 0.5,ndx) if ndx gt 0 then combconc(idx)=conc(idx) ENDELSE ;IF ifl EQ 0 ENDFOR ;FOR ifl=0,nhr ENDELSE ;IF nhr EQ 1 outfile=dir+SAT+'.comb.'+stamp+'.'+ver+'.dat' print, 'outfile='+outfile write_binary_all,outfile,nlon,nlat,lon_center,lat_center,combaod,combconc outfl_grib=dir0+yr+month+nday+'/GRIB/'+SAT+'dust.aod_conc.'+ver+'.'+date+'.hr'+string(ihr,format='(i2.2)')+'.grib' openw,1,'infile0' printf,1,outfile printf,1,outfl_grib printf,1,strmid(date,2,6)+string(ihr,format='(i2.2)')+'00' close,1 print,'before grib_main_dus' ; spawn,'/net/orbit184l/home/jzeng/DUST/Package/GRIB/grib_main_dust' spawn,'/home/oper/Dust/OSPO_Package_Dust_GRIB/IDL_code/grib_main_dust' ; spawn,'rm -f infile0' print,'createdsss '+outfl_grib ; file_copy, outfl_grib,ftp_dir,/overwrite ; if plot_dust eq 1 then begin outimg=dir0+yr+month+nday+'/AOD/'+SAT+'dustaod.'+ver+'.'+stamp bartit='Dust Aerosol Optical Thickness' disp_aod, combaod,lat_center, lon_center,outimg,title=SAT+'_Dust_'+ver+' '+yr+'/'+month+'/'+nday+' '+string(ihr,format='(i2.2)')+'Z-'+string(ihr+1,format='(i2.2)')+'Z',bartit=bartit,latdel=5,londel=10 outimg=dir0+yr+month+nday+'/AOD/'+SAT+'dustconc.'+ver+'.'+stamp bartit='Dust Column Concentration(!4l!3g/m!u3!n)' disp_conc,combconc,lat_center,lon_center,outimg,title=SAT+'_Dust_'+ver+' '+yr+'/'+month+'/'+nday+' '+string(ihr,format='(i2.2)')+'Z-'+string(ihr+1,format='(i2.2)')+'Z',bartit=bartit,latdel=5,londel=10 ; endif ENDIF ELSE print,'no MODIS granules at hour of',ihr ENDFOR ENDELSE ENDELSE END