pro find_gpm_overpass,forbit ;************************************* ; Input ;************************************* ;path_prefix='/uufs/chpc.utah.edu/common/home/' path_prefix='/Volumes/' ; Choose which ship track to use region='investigator_2020' ;region='investigator_2019' ;region='investigator_2018' ; This is the region covering the entire ship track of the experiment if region eq 'investigator_2019' then begin ulat=-42.0 & llat=-67.0 & llon=137.0 & rlon=157.0;153.0 ;big region for 2019 shipdir=path_prefix+'u0079358/gpm/' shipfilename='investigator_ship_track.20190118.210915.20190304.211315.1min.cdf' ship_search='2019' endif else if region eq 'investigator_2018' then begin ulat=-42.0 & llat=-66.0 & llon=130.0 & rlon=152.0 ;big region for 2018 shipdir=path_prefix+'mace-group5/field_programs/socrates/ship_track/' shipfilename='socrates_ship_track.20180110.232300.20180221.052900.cdf' ship_search='2018' endif else if region eq 'investigator_2020' then begin ulat=-42.0 & llat=-66.0 & llon=130.0 & rlon=152.0 ;big region for 2018 shipdir=path_prefix+'mace-group5/field_programs/investigator/in2020_v08/' shipfilename='in2020_v08_ship_track.20201204.124200.20210115.060900.cdf' ship_search='2020' endif ; If you put in a gpm orbit, this is the only one that will run ; 2018 investigator overpasses ;forbit='021988' ;20180111.00 ;no cband ;forbit='022028' ;20180113.15 ;forbit='022136' ;20180120.14 ;forbit='022167' ;20180122.14 ;forbit='022198' ;20180124.13 ;no cband ;forbit='022229' ;20180126.13 ;forbit='022260' ;20180128.13 ;no cband ;forbit='022353' ;20180203.13 ;no cband ;forbit='022368' ;20180204.12 ; ;forbit='022370' ;20180204.13 ; ;forbit='022383' ;20180205.11 ;forbit='022398' ;20180206.10 ; ;forbit='022448' ;20180209.14 ;**** ;forbit='022460' ;20180210.10 ; ;forbit='022541' ;20180215.13 ;forbit='022567' ;20180217.17 ; ;forbit='022618' ;20180220.12 ; 2019 investigator overpasses ;forbit='027858' ;20190123.09 ;forbit='027869' ;20190124.03 ;forbit='027873' ;20190124.08 ;forbit='027885' ;20190125.04 ;forbit='027888' ;20190125.07 ;forbit='027903' ;20190126.06 ;forbit='027932' ;20190128.04 ;only one with ka ;forbit='027947' ;20190129.03 ;forbit='027949' ;20190129.05 ;forbit='027964' ;20190130.04 ;forbit='028085' ;20190207.00 ;forbit='028087' ;20190207.02 ;forbit='028102' ;20190208.01 ;error don't have ka file ;forbit='028146' ;20190210.23 ;forbit='028162' ;20190211.23 ;error don't have ka file ;forbit='028193' ;20190213.23 ;error don't have ka file ;forbit='028225' ;20190215.23 ;forbit='028392' ;20190226.18 ; 2020 investigator overpasses ;forbit='038478' ;20201206.022621 203.34 ;no gpm in small box ;forbit='038503' ;20201207.181326 5.31 ;forbit='038580' ;20201212.165924 168.08 ;no gpm in small box ;forbit='038586' ;20201213.010054 264.13 ;no gpm in small box ;forbit='038657' ;20201217.154550 353.78 ;no gpm in small box ;forbit='038724' ;20201221.215229 43.17 forbit='038734' ;20201222.143534 151.47 ;2021 01 11 14 38 41 39046 141.98 -55.90 287.95 ;2021 01 12 15 22 32 39062 144.27 -51.49 143.24 ;2021 01 14 07 12 04 39087 145.85 -46.34 257.11 ;2021 01 14 15 14 27 39093 146.67 -44.91 243.34 ; Output file directories output_dir='' ;imagedir=path_prefix+'u0079358/gpm/' imagedir='./' ;imagedir='/uufs/chpc.utah.edu/common/home/u0079358/gpm/' ;imagedir='/uufs/chpc.utah.edu/common/home/u0079358/public_html/gpm/overpasses_power/' ;************************************* ; Constants ;************************************* ; KU Radar ku_name='KuPR' ku_freq='13.6 GHz' ku_swath='245 km' ku_dir=path_prefix+'mace-group5/gpm/GPM_L2/GPM_2AKu.06/' ku_sub_pt=24 ;sub satellite point ku_deltah=250 ;m for NS ;ku_sfc_idx=175 ; bottom index original one ku_sfc_idx=167 ; KA Radar ka_name='KaPR' ka_freq='35.5 GHz' ka_swath='120 km' ka_dir=path_prefix+'mace-group5/gpm/GPM_L2/GPM_2AKa.06/' ka_sub_pt=12 ;sub satellite point ;ka_deltah=500 ;m for HS ka_deltah=250 ;m for MS ;ka_sfc_idx=175;87 ;surface plane ka_sfc_idx=167 ; GMI ***choose the channel** 10-183ghz gmi_dir=path_prefix+'mace-group5/gpm/GPM_L1B/GPM_1BGMI.05/' ; Raw power directory pow_dir=path_prefix+'/mace-group5/gpm/GPM_L1B/GPM_PRL1KU.05/' ; Common dbz min and max dmin_dbz=0 dmax_dbz=35 ;***************************** ; Read in the ship track file ;***************************** fid=ncdf_open(shipdir+shipfilename) xid=ncdf_varid(fid,'julian_day') & ncdf_varget,fid,xid,julian_day_ship xid=ncdf_varid(fid,'lat') & ncdf_varget,fid,xid,lat_ship xid=ncdf_varid(fid,'lon') & ncdf_varget,fid,xid,lon_ship ncdf_close,fid ;****************************************************************************************** ; Calculate co-lon coordinates ;****************************************************************************************** lcolon=llon if lcolon lt 0 then lcolon=lcolon+360.0 rcolon=rlon if rcolon lt 0 then rcolon=rcolon+360.0 print,ulat,llat print,llon,lcolon,rlon,rcolon print,llon,rlon,rlon,llon,llon,'*',ulat,ulat,llat,llat,ulat ;********************** ; Read in Ku to start working with ;********************** ; Read Ku read_gpm_2aku_v2,ku_dir,forbit,lat_ku,lon_ku,colon_ku,julian_day_ku,$ dbz_ku,dbz_measured_ku,lat_ku_sub,lon_ku_sub,colon_ku_sub,nbin_ku ;************************************************ ; Check to see if this orbit track is in our big region during the experiment ;************************************************* ; Check for big region overpass idx_profiles=where(lat_ku_sub ge llat and lat_ku_sub le ulat and $ colon_ku_sub ge lcolon and colon_ku_sub le rcolon and $ julian_day_ku ge julian_day_ship[0] and $ julian_day_ku le julian_day_ship[-1],count_profiles) print,count_profiles,' profiles in big region during experiment' if count_profiles eq 0 then print,' overpass not in big region' if count_profiles gt 1 then begin print,'overpass in big region',count_profiles,'count_profiles' ;************************************************ ; Find the ships location at this time ; Calculate an average position and time for the ship location during overpass ; julian_day must be monotonic increasing ;************************************************ result_ship=where(julian_day_ship ge julian_day_ku[idx_profiles[0]] and $ julian_day_ship le julian_day_ku[idx_profiles[-1]],count_ship) ; Ship was out at this time if count_ship gt 0 then begin print,count_ship,'ship times in gpm overpass time interval' ; Average time and location of ship during gpm overpass julian_day_ship_ovp=total(julian_day_ship[result_ship])/float(count_ship) lat_ship_ovp=total(lat_ship[result_ship])/float(count_ship) lon_ship_ovp=total(lon_ship[result_ship])/float(count_ship) colon_ship_ovp=lon_ship_ovp if colon_ship_ovp lt 0 then colon_ship_ovp=colon_ship_ovp+360.0 ; Check the times caldat,julian_day_ku[0],smm,sdd,syy,shh,smi,sss print,syy,smm,sdd,shh,smi,sss,'gpm start' caldat,julian_day_ship_ovp,omm,odd,oyy,ohh,omi,oss print,oyy,omm,odd,ohh,omi,oss,'ship' caldat,julian_day_ku[-1],emm,edd,eyy,ehh,emi,ess print,eyy,emm,edd,ehh,emi,ess,'gpm end' ; Find the closest distance between ship track and gpm sub point delta_dist=make_array(n_elements(lon_ku_sub),/float,value=1e7) for ii=0,n_elements(lon_ku_sub)-1 do begin if lon_ku_sub[ii] gt -9999 and lat_ku_sub[ii] gt -9999 then begin delta_dist[ii]=map_2points(lon_ship_ovp,lat_ship_ovp,$ lon_ku_sub[ii],lat_ku_sub[ii],/meters) endif endfor ; Pick out the time and location of gpm closest approach min_dist=min(delta_dist,didx) julian_day_ovp=julian_day_ku[didx] lat_ovp=lat_ku_sub[didx] lon_ovp=lon_ku_sub[didx] colon_ovp=colon_ku_sub[didx] print,'min_dist between satillite and ship track ',min_dist/1000.0,'km' print,'min_dist idx',didx ; Make a date string for the output grid file caldat,julian_day_ovp,omm,odd,oyy,ohh,omi,oss odate_str=string(oyy,format='(I4)')+string(omm,format='(I02)')+string(odd,format='(I02)') otime_str=string(ohh,format='(I02)')+':'+string(omi,format='(I02)')+':'+string(oss,format='(I02)') otime_str2=string(ohh,format='(I02)')+string(omi,format='(I02)')+string(oss,format='(I02)') ; The satellite needs to pass close to the ship to be in cband range ; Calculate ovp box to be delta_degrees aroud the ship delta_degree=0.5 ulat_ovp=lat_ship_ovp+delta_degree & llat_ovp=lat_ship_ovp-delta_degree llon_ovp=lon_ship_ovp-2.0*delta_degree & rlon_ovp=lon_ship_ovp+2.0*delta_degree lcolon_ovp=colon_ship_ovp-2.0*delta_degree & rcolon_ovp=colon_ship_ovp+2.0*delta_degree ; Check again to see if closest approach is close enough to ship idx_ovp=where(lat_ku_sub ge llat_ovp and lat_ku_sub le ulat_ovp and $ colon_ku_sub ge lcolon_ovp and colon_ku_sub le rcolon_ovp,count_ovp) if count_ovp eq 0 then begin print,'no satellite in small ship ovp box' endif if count_ovp gt 1 then begin print,'satellite data in small ship ovp box',count_ovp,' count_ovp' ; Calculate plotting ovp box to be delta_degrees around the ship delta_degree=2.0 ulat_ovp=lat_ship_ovp+delta_degree & llat_ovp=lat_ship_ovp-delta_degree llon_ovp=lon_ship_ovp-2.0*delta_degree & rlon_ovp=lon_ship_ovp+2.0*delta_degree lcolon_ovp=colon_ship_ovp-2.0*delta_degree & rcolon_ovp=colon_ship_ovp+2.0*delta_degree idx_ovp=where(lat_ku_sub ge llat_ovp and lat_ku_sub le ulat_ovp and $ colon_ku_sub ge lcolon_ovp and colon_ku_sub le rcolon_ovp,count_ovp) ; calculate range rings around ship location (25km) lon_cir=fltarr(4,360) lat_cir=fltarr(4,360) for r=1,4 do begin arc_dist=(float(r)*25000.0)/(6371D*1000D) ;m for s=0,359 do begin result=ll_arc_distance([lon_ship_ovp,lat_ship_ovp],arc_dist,s,/degrees) lon_cir[r-1,s]=result[0] lat_cir[r-1,s]=result[1] endfor endfor ; Calculate regular grid array of the 2degree plotting overpass box delta=0.01 ynum=fix((ulat_ovp-llat_ovp)/delta) grid_lat=findgen(ynum)*delta+llat_ovp xnum=fix((rcolon_ovp-lcolon_ovp)/delta) grid_colon=findgen(xnum)*delta+lcolon_ovp grid_lon=grid_colon result=where(grid_lon gt 180,count) if count gt 0 then grid_lon[result]=grid_lon[result]-360.0 ; gpm height array height_ku=(findgen(nbin_ku)*ku_deltah)/1000.0 height_ku_reverse=reverse(height_ku) ;************************ ; Set up the plot ;************************ ; Set up plot size pxdim=1500 & pydim=900 ; Position the plots xl=0.05 & xr=0.95 yb=0.08 & yt=0.90 sx=0.05 sy=0.12 numplots_x=6 numplots_y=2 position_plots,xl,xr,yb,yt,sx,sy,numplots_x,numplots_y,pos ; Colorbar position cbpos=pos ;cbpos[*,0]=pos[*,2]+0.08 ;vertical ;cbpos[*,2]=cbpos[*,0]+0.01 cbpos[*,1]=pos[*,3]+0.019 ;horizontal cbpos[*,3]=cbpos[*,1]+0.01 ; Top is the last color to scale 256 colors, 0-255 top_color=252 ; Colortable 0-252 253=white mytable=colortable(39,ncolors=254) ;254=hot pink ;gray=255 mytable=[mytable,transpose([238,18,137]),transpose([200,200,200])] mycbtable=mytable[0:top_color,*] ; Strings dx=0.01 & dy=0.01 fs1=12 x_above_cb=3.5 ;Clear plot space pnum=0 p0=plot([0,1],[0,1],position=pos[pnum,*],/buffer,dimensions=[pxdim,pydim],axis_style=4,/nodata) d=p0.convertcoord([pos[pnum,0],pos[pnum,2]],[pos[pnum,1],pos[pnum,3]],/normal,/to_device) isx=(d[0,1,0]-d[0,0,0]) isy=(d[1,1,0]-d[1,0,0]) ;***************************** ; Plot map - Ku - z factor corrected ;***************************** pnum=6 ;dmin_dbz=10 ;dmax_dbz=35 ku_sfc_idx1=ku_sfc_idx plot_map_gpm_dbz,pnum,dbz_ku,ku_sfc_idx1,lat_ku,lon_ku,colon_ku,grid_colon,grid_lat,$ dmin_dbz,dmax_dbz,top_color,llat_ovp,lcolon_ovp,ulat_ovp,rcolon_ovp,$ pxdim,pydim,mytable,pos,lon_cir,lat_cir,lon_ku_sub,lat_ku_sub,idx_ovp,height_ku_reverse,$ mycbtable,cbpos,fs1,p0,p1,'ku_zFactorCorrected.'+odate_str+'.'+otime_str2+'.'+forbit+'.cdf' ;t1=text(pos[pnum,0]+0*dx,pos[pnum,3]+x_above_cb*dy,'ZFactorCorrected (dBZ)',font_size=fs1) t3=text(pos[pnum,0]+0*dx,pos[pnum,3]+x_above_cb*dy,ku_name+' Reflectivity Factor with attenuation correction (dBZ)',font_size=fs1) t3=text(pos[pnum,0]+0*dx,pos[pnum,3]-2*dy,'surface height='+string(height_ku_reverse[ku_sfc_idx1],format='(F5.2)')+'km',font_size=fs1) ;***************************** ; Plot curtain - Ku - z factor corrected ;***************************** pnum=7 plot_curtain_gpm_dbz,pnum,dbz_ku,ku_sub_pt,idx_ovp,top_color,dmin_dbz,dmax_dbz,$ height_ku,isx,isy,pos,mytable,lon_ku_sub,fs1,p0,p1 t1=text(pos[pnum,0]+0*dx,pos[pnum,3]+1.5*dy,'Sub Satellite point',font_size=fs1) ;************************** ; Plot map - Ku - z factor measured ;sfc_idx-15, ku_sfc_idx-14 ;************************** pnum=8 dmin_dbz_measured=dmin_dbz;-10 dmax_dbz_measured=dmax_dbz;50 ku_sfc_idx1=ku_sfc_idx plot_map_gpm_dbz,pnum,dbz_measured_ku,ku_sfc_idx1,lat_ku,lon_ku,colon_ku,grid_colon,grid_lat,$ dmin_dbz_measured,dmax_dbz_measured,top_color,llat_ovp,lcolon_ovp,ulat_ovp,rcolon_ovp,$ pxdim,pydim,mytable,pos,lon_cir,lat_cir,lon_ku_sub,lat_ku_sub,idx_ovp,height_ku_reverse,$ mycbtable,cbpos,fs1,p0,p1,'ku_zFactorMeasured.'+odate_str+'.'+otime_str2+'.'+forbit+'.cdf' t3=text(pos[pnum,0]+0*dx,pos[pnum,3]+x_above_cb*dy,ku_name+' Reflectivity Factor as measured (dBZ)',font_size=fs1) t3=text(pos[pnum,0]+0*dx,pos[pnum,3]-2*dy,'surface height='+string(height_ku_reverse[ku_sfc_idx1],format='(F5.2)')+'km',font_size=fs1) ;***************************** ; Plot curtain - Ku - z factor measured ;***************************** pnum=9 plot_curtain_gpm_dbz,pnum,dbz_measured_ku,ku_sub_pt,idx_ovp,top_color,$ dmin_dbz_measured,dmax_dbz_measured,$ height_ku,isx,isy,pos,mytable,lon_ku_sub,fs1,p0,p1 t1=text(pos[pnum,0]+0*dx,pos[pnum,3]+1.5*dy,'Sub Satellite point',font_size=fs1) ;***************************** ; Plot map - KA - z factor corrected ;***************************** pnum=0 ;dmin_dbz=10 ;dmax_dbz=35 read_gpm_2aka,ka_dir,forbit,lat_ka,lon_ka,colon_ka,julian_day_ka,$ dbz_ka,bin_echo_bottom_ka,dbz_measured_ka,bin_real_surface_ka,$ lat_ka_sub,lon_ka_sub,colon_ka_sub,nbin_ka height_ka=(findgen(nbin_ka)*ka_deltah)/1000.0 height_ka_reverse=reverse(height_ka) idx_ovp_ka=where(lat_ka_sub ge llat_ovp and lat_ka_sub le ulat_ovp and $ colon_ka_sub ge lcolon_ovp and colon_ka_sub le rcolon_ovp,count_ovp_ka) ka_sfc_idx1=ka_sfc_idx plot_map_gpm_dbz,pnum,dbz_ka,ka_sfc_idx1,lat_ka,lon_ka,colon_ka,grid_colon,grid_lat,$ dmin_dbz,dmax_dbz,top_color,llat_ovp,lcolon_ovp,ulat_ovp,rcolon_ovp,$ pxdim,pydim,mytable,pos,lon_cir,lat_cir,lon_ka_sub,lat_ka_sub,idx_ovp_ka,height_ka_reverse,$ mycbtable,cbpos,fs1,p0,p1,'ka_zFactorCorrected.'+odate_str+'.'+otime_str2+'.'+forbit+'.cdf' ;t1=text(pos[pnum,0]+0*dx,pos[pnum,3]+x_above_cb*dy,'ZFactorCorrected (dBZ)',font_size=fs1) t3=text(pos[pnum,0]+0*dx,pos[pnum,3]+x_above_cb*dy,ka_name+' Reflectivity Factor with attenuation correction (dBZ)',font_size=fs1) t3=text(pos[pnum,0]+0*dx,pos[pnum,3]-2*dy,'surface height='+string(height_ka_reverse[ka_sfc_idx1],format='(F5.2)')+'km',font_size=fs1) ;***************************** ; Plot curtain - Ka - z factor corrected ;***************************** pnum=1 plot_curtain_gpm_dbz,pnum,dbz_ka,ka_sub_pt,idx_ovp_ka,top_color,dmin_dbz,dmax_dbz,$ height_ka,isx,isy,pos,mytable,lon_ka_sub,fs1,p0,p1 t1=text(pos[pnum,0]+0*dx,pos[pnum,3]+1.5*dy,'Sub Satellite point',font_size=fs1) ;************************** ; Plot map - Ka - z factor measured ;sfc_idx-15 ;************************** pnum=2 dmin_dbz_measured=dmin_dbz;-10 dmax_dbz_measured=dmax_dbz;50 plot_map_gpm_dbz,pnum,dbz_measured_ka,ka_sfc_idx1,lat_ka,lon_ka,colon_ka,grid_colon,grid_lat,$ dmin_dbz_measured,dmax_dbz_measured,top_color,llat_ovp,lcolon_ovp,ulat_ovp,rcolon_ovp,$ pxdim,pydim,mytable,pos,lon_cir,lat_cir,lon_ka_sub,lat_ka_sub,idx_ovp_ka,height_ka_reverse,$ mycbtable,cbpos,fs1,p0,p1,'ka_zFactorMeasured.'+odate_str+'.'+otime_str2+'.'+forbit+'.cdf' t3=text(pos[pnum,0]+0*dx,pos[pnum,3]+x_above_cb*dy,ka_name+' Reflectivity Factor as measured (dBZ)',font_size=fs1) t3=text(pos[pnum,0]+0*dx,pos[pnum,3]-2*dy,'surface height='+string(height_ka_reverse[ka_sfc_idx1],format='(F5.2)')+'km',font_size=fs1) ;***************************** ; Plot curtain - Ka - z factor measured ;***************************** pnum=3 plot_curtain_gpm_dbz,pnum,dbz_measured_ka,ka_sub_pt,idx_ovp_ka,top_color,$ dmin_dbz_measured,dmax_dbz_measured,$ height_ka,isx,isy,pos,mytable,lon_ka_sub,fs1,p0,p1 t1=text(pos[pnum,0]+0*dx,pos[pnum,3]+1.5*dy,'Sub Satellite point',font_size=fs1) ; ;************** ; ; ku power map ; ;************** ; pnum=2 ; ; read_gpm_prl1ku,pow_dir,syear,sdoy3,ulat_ovp,llat_ovp,llon_ovp,rlon_ovp,$ ; lcolon_ovp,rcolon_ovp,$ ; colon_pow,lat_pow,lon_pow,pow,forbit,bin_size_pow,nbin_pow,bin_echo_peak,$ ; lat_sub_pow,colon_sub_pow,lon_sub_pow ; ; ;idx_ovp_pow=where(lat_sub_pow ge llat_ovp and lat_sub_pow le ulat_ovp and $ ; ; colon_sub_pow ge lcolon_ovp and colon_sub_pow le rcolon_ovp,count_ovp_pow) ; idx_ovp_pow=where(lat_sub_pow ge min(lat_sub[idx_ovp]) and $ ; lat_sub_pow le max(lat_sub[idx_ovp]) and $ ; lon_sub_pow ge min(lon_sub[idx_ovp]) and $ ; lon_sub_pow le max(lon_sub[idx_ovp]),count_ovp_pow) ; ; ; pull out surface ; bin_echo_peak_sub=reform(bin_echo_peak[sub_pt,*]) ; sfc_echo_idx=min(bin_echo_peak_sub[idx_ovp_pow]) ; ; dmin_pow=-114 ; dmax_pow=-105;-95 ; plot_map_gpm_dbz,pnum,pow,sfc_echo_idx-10,lat_pow,lon_pow,colon_pow,$ ; grid_colon,grid_lat,$ ; dmin_pow,dmax_pow,top_color,llat_ovp,lcolon_ovp,ulat_ovp,rcolon_ovp,$ ; pxdim,pydim,mytable,pos,lon_cir,lat_cir,lon_sub_pow,lat_sub_pow,idx_ovp_pow,$ ; mycbtable,cbpos,fs1,p0,p1 ; ; t3=text(pos[pnum,0],pos[pnum,3]+x_above_cb*dy,'L1 Echo Power (dBm)',font_size=fs1) ; pow_sub1=reform(pow[*,sub_pt,*]) ; p8=image(pow_sub1,rgb_table=mytable,max_value=dmax_pow,min_value=dmin_pow,$ ; image_dimensions=[400,400]) ; pow_sub2=pow_sub1[0:sfc_echo_idx,*] ; p9=image(pow_sub2,rgb_table=mytable,max_value=dmax_pow,min_value=dmin_pow,$ ; image_dimensions=[400,400]) ; ;************** ; ; ku power curtain ; ;************** ; pnum=3 ; ; ; Subset and make height array for pow ; pow_sub=pow[0:sfc_echo_idx,*,*] ; s=size(pow_sub,/dimensions) ; nbin_sfc=s[0] ; height_pow=(findgen(nbin_sfc)*deltah)/1000.0 ; ; plot_curtain_gpm_dbz,pnum,pow_sub,sub_pt,idx_ovp_pow,top_color,$ ; dmin_pow,dmax_pow,$ ; height_pow,isx,isy,pos,mytable,lon_sub_pow,fs1,p0,p1 ; t1=text(pos[pnum,0]+0*dx,pos[pnum,3]+1.5*dy,'Sub Satellite point',font_size=fs1) ; ;************************************************************************* ;****************** ; Get the cband radar ;************************************************************************* pnum=10 read_cband_pyart,julian_day_ship_ovp,$ ulat_ovp,llat_ovp,llon_ovp,rlon_ovp,dbz_cband,$ lat_cband,lon_cband,alt_cband,julian_day_cband,colon_cband,path_prefix if dbz_cband ne !null then begin sidx=4 ;cband surface dmin_cband=dmin_dbz;-10 dmax_cband=dmax_dbz;25 data_var=reform(dbz_cband[*,*,sidx]) lat_var=reform(lat_cband[*,*,sidx]) lon_var=reform(lon_cband[*,*,sidx]) data_image=bytscl(data_var,top=top_color,min=dmin_cband,max=dmax_cband) result=where(data_var lt -9000 or data_var gt 1e36,count) if count gt 0 then data_image[result]=255 ;gray ;p2=image(data_var,lon_var,lat_var,$ p2=image(data_image,lon_var,lat_var,$ limit=[llat_ovp,llon_ovp,ulat_ovp,rlon_ovp],$ /current,dimensions=[pxdim,pydim],grid_units='degrees',$ ;center_longitude=(rlon_ovp-llon_ovp)/2.0,$ rgb_table=mytable,$;min_value=dmin_cband,max_value=dmax_cband,$ ;label_format='MapGrid_Labels',$ /box_axes,label_position=0,label_show=1,irregular=0,$ map_projection='Mercator',position=pos[pnum,*],font_size=fs1) ; ;p0.save,'test_test_3.png',height=pydim ; ;stop p3=mapcontinents(color='black',/hires) p2.mapgrid.label_color='black' p2.mapgrid.linestyle='dot' for r=0,3 do begin p7=plot(lon_cir[r,*],lat_cir[r,*],/current,/overplot) endfor p8=plot(lon_ku_sub[idx_ovp],lat_ku_sub[idx_ovp],/overplot,color='gray',linestyle=2,thick=2) c2=colorbar(orientation=0,rgb_table=mycbtable,range=[dmin_cband,dmax_cband],$ position=reform(cbpos[pnum,*]),font_size=10,tickdir=1,border_on=1) caldat,julian_day_cband[0],mm,dd,yy,hh,mi,ss date_time=string(yy,format='(I04)')+string(mm,format='(I02)')+string(dd,format='(I02)')+$ ' '+string(hh,format='(I02)')+':'+string(mi,format='(I02)')+':'+string(ss,format='(I02)') t1=text(pos[pnum,0]+0*dx,pos[pnum,3]+x_above_cb*dy,'CBand '+date_time,font_size=fs1) t3=text(pos[pnum,0]+0*dx,pos[pnum,3]-2*dy,'surface height='+string(alt_cband[0,0,sidx]/1000.0,format='(F5.2)')+'km',font_size=fs1) if 1 eq 0 then begin ; Write this out to a netcdf file output_var=data_var result=where(data_var lt -9000 or data_var gt 1e36,count) if count gt 0 then output_var[result]=-9999 s=size(output_var,/dimensions) outfilename='cband_dbz.'+odate_str+'.'+otime_str2+'.'+forbit+'.cdf' fid=ncdf_create(outfilename,/clobber) nlons_did=ncdf_dimdef(fid,'numlons',s[0]) nlats_did=ncdf_dimdef(fid,'numlats',s[1]) vid=ncdf_vardef(fid,'lon_ship_ovp') vid=ncdf_vardef(fid,'lat_ship_ovp') vid=ncdf_vardef(fid,'grid_lat',[nlats_did,nlons_did]) vid=ncdf_vardef(fid,'grid_lon',[nlats_did,nlons_did]) vid=ncdf_vardef(fid,'grid_var',[nlats_did,nlons_did],/float) vid=ncdf_vardef(fid,'surface_height_km') ncdf_control,fid,/endef ncdf_varput,fid,'lon_ship_ovp',lon_ship_ovp ncdf_varput,fid,'lat_ship_ovp',lat_ship_ovp ncdf_varput,fid,'grid_lat',lat_var ncdf_varput,fid,'grid_lon',lon_var ncdf_varput,fid,'grid_var',output_var ncdf_varput,fid,'surface_height_km',alt_cband[0,0,sidx]/1000.0 ncdf_close,fid endif ;************************** ; Plot cband curtain ;************************** if 1 eq 1 then begin pnum=11 ; Pull data out along gpm subsatellite point height_cband=alt_cband[0,0,*]/1000.0 lat1=reform(lat_cband[0,*,0]) lon1=reform(lon_cband[*,0,0]) data_var=make_array(count_ovp,n_elements(height_cband),value=-8888,/float) for k=0,count_ovp-1 do begin closest_lat=min(abs(lat_ku_sub[idx_ovp[k]]-lat1),idx_lat) closest_lon=min(abs(lon_ku_sub[idx_ovp[k]]-lon1),idx_lon) ;print,closest_lat,closest_lon if closest_lon lt 0.1 and closest_lat le 0.1 then begin data_var[k,*]=dbz_cband[idx_lon,idx_lat,*] endif endfor ; Plot the curtain hidx=where(height_cband lt 10) data_image=bytscl(data_var,top=top_color,min=dmin_cband,max=dmax_cband) result=where(data_var eq -9999 or data_var gt 1e36,count) if count gt 0 then data_image[result]=255 ;gray result=where(data_var eq -8888,count) if count gt 0 then data_image[result]=253 ;white p11=image(data_image[*,hidx],/current,image_dimensions=[isx,isy],$ position=pos[pnum,*],rgb_table=mytable) c1=contour(data_var[*,hidx],lon_ku_sub[idx_ovp],height_cband[hidx],$ ytitle='height (km)',xtitle='Longitude',irregular=0,xtickdir=1,ytickdir=1,$ xstyle=1,ystyle=1,/nodata,/current,position=pos[pnum,*],font_size=fs1,$ axis_style=2,$ yrange=[height_cband[hidx[0]],height_cband[hidx[-1]]]);,xrange=[lat_var[0,0],lat_var[-1,0]]) ;xaxis=axis('X',location=-2,title='Latitude',target=c1,$ ; coord_transform=[intercept,slope_m],tickfont_size=fs1) p8=plot(lon_ku_sub[idx_ovp],make_array(count_ovp,value=height_cband[sidx]),$ /overplot,color='gray',thick=2,linestyle=2) ;t1=text(pos[pnum,0]+0*dx,pos[pnum,3]+x_above_cb*dy,string(alt_cband[0,0,sidx]/1000.0,format='(F4.2)')+'km'+$ ; ' CBand '+date_time,font_size=fs1) endif endif else begin ;cband_dbz array exists pnum=10 t1=text(pos[pnum,0]+dx,pos[pnum,3]-2*dy,'no cband file',font_size=fs1) endelse ;************************************************************************* ; GMI ;Swath S1 has channels 1-9: 9channels ;10V 10H 19V 19H 23V 37V 37H 89V 89H. ;Swath S2 has channels 10-13: 4channels ;166V 166H 183+/-3V 183+/-8V. ;************************************************************************* ; Read in all the data read_gpm_1bgmi_all_chan,gmi_dir,forbit,$ lat_s1,lon_s1,colon_s1,julian_day_s1,tb_s1,$ lat_s2,lon_s2,colon_s2,julian_day_s2,tb_s2 channels_s1=['10V','10H','19V','19H','23V','37V','37H','89V','89H'] channels_s2=['166V','166H','183+/-3V','183+/-8V'] ; Channel 8 pnum=4 dmin_gmi=180 & dmax_gmi=265 chan8=reform(tb_s1[*,*,8]) idx_ovp_gmi=where(lat_s1 ge llat_ovp and lat_s1 le ulat_ovp and $ lon_s1 ge lcolon_ovp and lon_s1 le rcolon_ovp,count_ovp_gmi) data_var=chan8[idx_ovp_gmi] lon_var=lon_s1[idx_ovp_gmi] lat_var=lat_s1[idx_ovp_gmi] p5=image(data_var,lon_var,lat_var,$ limit=[llat_ovp,lcolon_ovp,ulat_ovp,rcolon_ovp],$ /current,dimensions=[pxdim,pydim],grid_units='degrees',$ /box_axes,label_position=0,$ ;label_format='MapGrid_Labels',$ rgb_table=mytable,min_value=dmin_gmi,max_value=dmax_gmi,$ map_projection='Mercator',position=pos[pnum,*],font_size=fs1) p6=mapcontinents(color='black',/hires) p5.mapgrid.linestyle='dot' for r=0,3 do begin p7=plot(lon_cir[r,*],lat_cir[r,*],/current,/overplot) endfor p8=plot(lon_ku_sub[idx_ovp],lat_ku_sub[idx_ovp],/overplot,color='white',linestyle=2,thick=2) c0=colorbar(rgb_table=mycbtable,range=[dmin_gmi,dmax_gmi],orientation=0,$ position= reform(cbpos[pnum,*]),font_size=10,tickdir=1,border_on=1) t1=text(pos[pnum,0]+0*dx,pos[pnum,3]+x_above_cb*dy,$ 'GMI Channel: '+channels_s1[8]+' Kelvin',font_size=fs1) ; Channel 166H pnum=5 dmin_gmi=230 & dmax_gmi=270 chan2=reform(tb_s2[*,*,1]) idx_ovp_gmi=where(lat_s2 ge llat_ovp and lat_s2 le ulat_ovp and $ lon_s2 ge lcolon_ovp and lon_s2 le rcolon_ovp,count_ovp_gmi) data_var=chan2[idx_ovp_gmi] lon_var=lon_s2[idx_ovp_gmi] lat_var=lat_s2[idx_ovp_gmi] ;dmin_gmi=min(data_var) & dmax_gmi=max(data_var) p5=image(data_var,lon_var,lat_var,$ limit=[llat_ovp,lcolon_ovp,ulat_ovp,rcolon_ovp],$ /current,dimensions=[pxdim,pydim],grid_units='degrees',$ /box_axes,label_position=0,$ ;label_format='MapGrid_Labels',$ rgb_table=mytable,min_value=dmin_gmi,max_value=dmax_gmi,$ map_projection='Mercator',position=pos[pnum,*],font_size=fs1) p6=mapcontinents(color='black',/hires) p5.mapgrid.linestyle='dot' for r=0,3 do begin p7=plot(lon_cir[r,*],lat_cir[r,*],/current,/overplot) endfor p8=plot(lon_ku_sub[idx_ovp],lat_ku_sub[idx_ovp],/overplot,color='white',linestyle=2,thick=2) c0=colorbar(rgb_table=mycbtable,range=[dmin_gmi,dmax_gmi],orientation=0,$ position= reform(cbpos[pnum,*]),font_size=10,tickdir=1,border_on=1) t1=text(pos[pnum,0]+0*dx,pos[pnum,3]+x_above_cb*dy,$ 'GMI Channel: '+channels_s2[1]+' Kelvin',font_size=fs1) ; Subset the data to write to a cdf file37 GHz, 89 GHz, 166 GHz idx_ovp_gmi=where(lat_s1 ge llat_ovp and lat_s1 le ulat_ovp and $ lon_s1 ge lcolon_ovp and lon_s1 le rcolon_ovp,count_ovp_gmi) lon_var=lon_s1[idx_ovp_gmi] lat_var=lat_s1[idx_ovp_gmi] chan37=reform(tb_s1[*,*,6]) chan37=chan37[idx_ovp_gmi] chan89=reform(tb_s1[*,*,8]) chan89=chan89[idx_ovp_gmi] chan166=reform(tb_s2[*,*,1]) chan166=chan166[idx_ovp_gmi] s=size(chan89,/dimensions) outfilename='gmi_89H.'+odate_str+'.'+otime_str2+'.'+forbit+'.cdf' fid=ncdf_create(outfilename,/clobber) num_did=ncdf_dimdef(fid,'num',s[0]) ovp_did=ncdf_dimdef(fid,'ovp',count_ovp) num_ring_did=ncdf_dimdef(fid,'num_ring',4) deg_did=ncdf_dimdef(fid,'degree',360) vid=ncdf_vardef(fid,'ulat_ovp') vid=ncdf_vardef(fid,'llat_ovp') vid=ncdf_vardef(fid,'rcolon_ovp') vid=ncdf_vardef(fid,'lcolon_ovp') vid=ncdf_vardef(fid,'lon_sub_satellite',ovp_did) vid=ncdf_vardef(fid,'lat_sub_satellite',ovp_did) vid=ncdf_vardef(fid,'lon_range_rings',[num_ring_did,deg_did]) vid=ncdf_vardef(fid,'lat_range_rings',[num_ring_did,deg_did]) vid=ncdf_vardef(fid,'lat',num_did) vid=ncdf_vardef(fid,'lon',num_did) vid=ncdf_vardef(fid,'chan37',num_did,/float) vid=ncdf_vardef(fid,'chan89',num_did,/float) vid=ncdf_vardef(fid,'chan166',num_did,/float) ;vid=ncdf_vardef(fid,'channel',/string) ncdf_control,fid,/endef ncdf_varput,fid,'ulat_ovp',ulat_ovp ncdf_varput,fid,'llat_ovp',llat_ovp ncdf_varput,fid,'rcolon_ovp',rcolon_ovp ncdf_varput,fid,'lcolon_ovp',lcolon_ovp ncdf_varput,fid,'lon_sub_satellite',lon_ku_sub[idx_ovp] ncdf_varput,fid,'lat_sub_satellite',lat_ku_sub[idx_ovp] ncdf_varput,fid,'lon_range_rings',lon_cir ncdf_varput,fid,'lat_range_rings',lat_cir ncdf_varput,fid,'lat',lat_var ncdf_varput,fid,'lon',lon_var ncdf_varput,fid,'chan37',chan37 ncdf_varput,fid,'chan89',chan89 ncdf_varput,fid,'chan166',chan166 ncdf_close,fid ; ;********************* ; ; Map of big region and ship track ; ;********************* ; ; pnum=9 ; p7=map('mercator',limit=[llat,llon-2.0,ulat,rlon+2.0],$ ; position=pos[pnum,*],/current,/box_axes,label_position=0,font_size=fs1,$ ; linestyle='dot') ; p1=mapcontinents(color='black',/hires) ; s0=symbol(lon_ovp,lat_ovp,'star',/data,/sym_filled,sym_color='red',sym_size=1) ; p2=plot(lon_ship,lat_ship,/overplot,color='gray') ;************** ; Image strings ;************** stime=min(julian_day_ku) etime=max(julian_day_ku) caldat,stime,smm,sdd,syy,shh,smi,sss caldat,etime,emm,edd,eyy,ehh,emi,ess sdate_str=string(syy,format='(I4)')+string(smm,format='(I02)')+string(sdd,format='(I02)') stime_str=string(shh,format='(I02)')+string(smi,format='(I02)')+string(sss,format='(I02)') edate_str=string(eyy,format='(I4)')+string(emm,format='(I02)')+string(edd,format='(I02)') etime_str=string(ehh,format='(I02)')+string(emi,format='(I02)')+string(ess,format='(I02)') stime_cstr=string(shh,format='(I02)')+':'+string(smi,format='(I02)')+':'+string(sss,format='(I02)') etime_cstr=string(ehh,format='(I02)')+':'+string(emi,format='(I02)')+':'+string(ess,format='(I02)') ;datastream=file_basename(filesd[i]) ;datastream=strmid(datastream,0,9) caldat,julian_day_ovp,omm,odd,oyy,ohh,omi,oss odate_str=string(oyy,format='(I4)')+string(omm,format='(I02)')+string(odd,format='(I02)') otime_str=string(ohh,format='(I02)')+':'+string(omi,format='(I02)')+':'+string(oss,format='(I02)') otime_str2=string(ohh,format='(I02)')+string(omi,format='(I02)')+string(oss,format='(I02)') ;t1=text(0.1,0.98-0*dy,odate_str+' '+otime_str+' '+datastream+' '+$ ; radarname+' '+radarfreq+' '+radarswath+' swath',font_size=fs1) t1=text(0.1,0.98-0*dy,odate_str+' '+otime_str+' orbit:'+forbit,font_size=fs1) ;t2=text(0.1,0.93,'Dual-Frequency Precipitation Radar (DPR)') ;t2=text(0.1,0.98-2.*dy,radarname+' '+radarfreq+' '+radarswath,font_size=fs1) ; Imagename imagename='ovp.ku_ka_gmi.dbz.'+odate_str+'.'+otime_str2+'.'+forbit imagename=imagename+'.png' print,imagedir+imagename p0.save,imagedir+imagename,height=pydim ;*************************** ; Write out ovp to a file ;*************************** ;openw,1,imagedir+'ship_gpm_overpasses_'+ship_search+'.txt',/append ;printf,1,oyy,omm,odd,ohh,omi,oss,forbit,lon_ship_ovp,lat_ship_ovp,min_dist/1000.0,$ ; format='(I4,5(1X,I02),1X,I6,2X,F8.2,2X,F8.2,2X,F8.2)' ;close,1 endif ;gpm profiles in small ship overpass box endif ;ship was out at this time endif ; gpm profiles in big region ; endfor ;end of loop through gpm radar files files ;endif ;files exist end