;******************************************** ; This program takes the investigator ship track and finds out when ; gpm passed over the ship track. It creates a text table of these ; overpasses. It creates a 4 plot of dbz corrected, dbz measured, ; gmi and L1B power. ; To run it you must have the ship track and gpm KU during the ; ship track time. ;********************************************* pro overpass_table ;************************************* ; Input: date or date & orbit or ship track file ;************************************* ; If there is a date with orbit, only orbit will run ;fdate='20210115' ;& forbit='038657' ;fdate='20180125' & forbit='22218' ;fdate='20180206' ;fdate='20180111' & forbit='021988' ;fdate='20180111' & forbit='021998' ;fdate='20180203' & forbit='022353' ; Choose which ship track to use ;region='investigator_2020' region='investigator_2019' ;region='investigator_2018' ;************************************** ; Directories ;************************************** ; imac or chpc path_prefix='/Volumes/' ;path_prefix='/uufs/chpc.utah.edu/common/home/' ;chpc ; 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 else begin ulat=-42.0 & llat=-67.0 & llon=137.0 & rlon=157.0;153.0 ;big region for 2019 ;region=fdate ;a fix to having no region name? endelse ; Directory output_dir='' ;file output directory imagedir='' ;imagedir='/Users/u0079358/Documents/gpm/' ;imagedir='/uufs/chpc.utah.edu/common/home/u0079358/public_html/gpm/overpasses_power/' ; Create overpass text file and print header ovp_filename='ship_gpm_overpasses_'+region+'.test.txt' openw,1,imagedir+ovp_filename printf,1,'year m d h m s orbit lon lat min distance (km) between ship and sub-satellite point' close,1 ;************************************* ; Constants ;************************************* ; RADAR *** radar='ku' ;only use ku for finding ovp radarname='KuPR' radarfreq='13.6 GHz' radarswath='245 km' gpm_dir='mace-group5/gpm/GPM_L2/GPM_2AKu.06/' ; Sub satellite point sub_pt=24 deltah=250 sfc_idx=175 ; GMI ***choose the channel** gmi_dir=path_prefix+'mace-group5/gpm/GPM_L1B/GPM_1BGMI.05/' chan=8 ;channel 7-89V 8-89H ; Raw power pow_dir=path_prefix+'mace-group5/gpm/GPM_L1B/GPM_PRL1KU.05/' ;***************************** ; 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 ; Start and end time for ship track stime_track=julian_day_ship[0] etime_track=julian_day_ship[-1] ;****************************************************************************************** ; Calculate co-lon coordinates of the big region ;****************************************************************************************** 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 ;****************************************************************************************** ; FIND HOW MANY ORBIT FILES THERE ARE FOR THIS DATE...THEN LOOP THROUGH THE ORBITS! ;****************************************************************************************** ; gpm directory and filename fdir=path_prefix+gpm_dir;+syear+'/'+sdoy3+'/' if forbit eq !NULL and fdate ne !NULL then begin syear=strmid(fdate,0,4) smonth=strmid(fdate,4,2) sday=strmid(fdate,6,2) file_string='2A.GPM.Ku.V8-20180723.'+syear+smonth+sday+'*HDF5' filesd=file_search(fdir+file_string,count=num_files) endif else if forbit ne !NULL and fdate ne !NULL then begin syear=strmid(fdate,0,4) smonth=strmid(fdate,4,2) sday=strmid(fdate,6,2) file_string='2A.GPM.Ku.V8-20180723.'+syear+smonth+sday+'*'+forbit+'*HDF5' filesd=file_search(fdir+file_string,count=num_files) endif else begin caldat,julian_day_ship,smm,sdd,syy,shh,smi,sss fdates=string(syy,format='(I4)')+string(smm,format='(I02)')+string(sdd,format='(I02)') fdates=fdates[uniq(fdates)] for i=0,n_elements(fdates)-1 do begin file_string='2A.GPM.Ku.V8-20180723.'+fdates[i]+'*HDF5' file_list=file_search(fdir+file_string,count=num_files) if filesd eq !NULL then begin filesd=file_list endif else begin filesd=[filesd,file_list] endelse endfor num_files=n_elements(filesd) endelse if num_files eq 0 then print,'NO FILES AT ALL' ;*********************************************** ; Loop through each gpm file ;*********************************************** if num_files gt 0 then begin for i=0,num_files-1 do begin print,filesd[i] parts=strsplit(file_basename(filesd[i]),'.',/extract) forbit=parts[5] ;****************** ; Read data from hdf5 file ;****************** file_id=h5f_open(filesd[i]) if radar eq 'ku' then begin ns_id=h5g_open(file_id,'NS') endif else if radar eq 'ka' then begin ns_id=h5g_open(file_id,'HS') endif dset_id=h5d_open(ns_id,'Latitude') lat=h5d_read(dset_id) h5d_close,dset_id dset_id=h5d_open(ns_id,'Longitude') lon=h5d_read(dset_id) h5d_close,dset_id slv_id=h5g_open(ns_id,'SLV') dset_id=h5d_open(slv_id,'zFactorCorrected') dbz=h5d_read(dset_id) h5d_close,dset_id h5g_close,slv_id time_id=h5g_open(ns_id,'ScanTime') dset_id=h5d_open(time_id,'DayOfYear') doy=h5d_read(dset_id) h5d_close,dset_id dset_id=h5d_open(time_id,'Year') year=h5d_read(dset_id) h5d_close,dset_id dset_id=h5d_open(time_id,'Month') month=h5d_read(dset_id) h5d_close,dset_id dset_id=h5d_open(time_id,'DayOfMonth') day=h5d_read(dset_id) h5d_close,dset_id dset_id=h5d_open(time_id,'Hour') hour=h5d_read(dset_id) h5d_close,dset_id dset_id=h5d_open(time_id,'Minute') minute=h5d_read(dset_id) h5d_close,dset_id dset_id=h5d_open(time_id,'Second') second=h5d_read(dset_id) h5d_close,dset_id dset_id=h5d_open(time_id,'MilliSecond') msecond=h5d_read(dset_id) h5d_close,dset_id h5g_close,time_id pre_id=h5g_open(ns_id,'PRE') dset_id=h5d_open(time_id,'elevation') elevation=h5d_read(dset_id) ;meters dset_id=h5d_open(time_id,'zFactorMeasured') dbz_measured=h5d_read(dset_id) h5d_close,dset_id h5g_close,ns_id h5f_close,file_id s=size(dbz,/dimensions) print,s nbin=s[0] ;176 nray=s[1] ;49 nscan=s[2] ;7937 ; ; loop through the range bins to see where the data is ; for j=0,nbin-1 do begin ; dbz1=reform(dbz[j,*,*]) ; result=where(dbz1 gt -9999,count) ; if count le 0 then begin ; print,j,count ; endif else if count gt 0 then begin ; print,j,count,min(dbz1[result]),max(dbz1[result]) ; endif ; endfor ;***************** ; Caculate colongitude because the pacific crosses the 180 line ;***************** colon=lon result=where(lon lt 0,count) if count gt 0 then colon[result]=colon[result]+360.0 ;***************** ; Pull out subsatellite point ;***************** lat_sub=reform(lat[sub_pt,*]) lon_sub=reform(lon[sub_pt,*]) colon_sub=reform(colon[sub_pt,*]) dbz_sub=reform(dbz[*,sub_pt,*]) ; Only keep valid lat and lon ;r=where(lon_sub ne -9999.90,c) ;if c gt 0 then begin ; lat_sub=lat_sub[r] ; lon_sub=lon_sub[r] ; colon_sub=colon_sub[r] ; dbz_sub=dbz_sub[r] ;endif ;***************** ; Calculate time ;***************** msecond=msecond*1e-3 second=second+msecond julian_day=julday(month,day,year,hour,minute,second) ;************************************************ ; Check to see if this orbit track is in our region during the experiment ;************************************************* ; Check for region overpass idx_profiles=run_pros ;idx_profiles=where(lat_sub ge llat and lat_sub le ulat and $ ; colon_sub ge lcolon and colon_sub le rcolon,count_profiles) idx_profiles=where(lat_sub ge llat and lat_sub le ulat and $ colon_sub ge lcolon and colon_sub le rcolon and $ julian_day ge stime_track and $ julian_day le etime_track,count_profiles) print,count_profiles,'profiles in overpass region during experiment' if count_profiles eq 0 then print,'overpass not in region' if count_profiles gt 1 then begin print,'overpass in region',count_profiles,'count_profiles' ;************************************************ ; Find the ships location at this time ; Calculate an average position and time for the overpass ;************************************************ result_ship=where(julian_day_ship ge julian_day[idx_profiles[0]] and $ julian_day_ship le julian_day[idx_profiles[-1]],count_ship) ; Ship was out at this time if count_ship gt 0 then begin print,count_ship,'ship times in overpass time interval' 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 ;delta_dist=5.0 ;km ;delta_time=60D/(60D*24D) ;1hr ;for ii=0,count_ship-1 do begin ; caldat,julian_day_ship[result_ship[ii]],smm,sdd,syy,shh,smi,sss ; print,syy,smm,sdd,shh,smi,sss ;endfor caldat,julian_day[0],smm,sdd,syy,shh,smi,sss print,syy,smm,sdd,shh,smi,sss caldat,julian_day_ship_ovp,omm,odd,oyy,ohh,omi,oss print,oyy,omm,odd,ohh,omi,oss caldat,julian_day[-1],emm,edd,eyy,ehh,emi,ess print,eyy,emm,edd,ehh,emi,ess ;min_time=min(abs(julian_day_ship-julian_day),tidx) ;min_time=min_time*24D delta_dist=make_array(n_elements(lon_sub),/float,value=1e6) for ii=0,n_elements(lon_sub)-1 do begin if lon_sub[ii] ne -9999.90 then $ delta_dist[ii]=map_2points(lon_ship_ovp,lat_ship_ovp,lon_sub[ii],lat_sub[ii],/meters) endfor min_dist=min(delta_dist,didx) julian_day_ovp=julian_day[didx] lat_ovp=lat_sub[didx] lon_ovp=lon_sub[didx] colon_ovp=colon_sub[didx] ; Calculate ovp box to be delta_degrees aroud 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_sub ge llat_ovp and lat_sub le ulat_ovp and $ colon_sub ge lcolon_ovp and colon_sub le rcolon_ovp,count_ovp) if count_ovp gt 1 then begin print,'satellite data in ship ovp box',count_ovp,'count_ovp' ; plot range rings 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 ;************************ ; Set up the plot ;************************ ; Set up plot size pxdim=1500 & pydim=900 ; Position the plots xl=0.05 & xr=0.90 yb=0.08 & yt=0.90 sx=0.05 sy=0.08 numplots_x=4 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.005 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([150,150,150])] mycbtable=mytable[0:top_color,*] ; Strings dx=0.01 & dy=0.01 fs1=12 ;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]) ; Calculate regular grid array delta=0.05 ynum=(ulat_ovp-llat_ovp)/delta grid_lat=findgen(ynum)*delta+llat_ovp xnum=(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 ; Common dbz min and max dmin_radar=0 dmax_radar=35 ;***************************** ; Plot map - Ku dbz ;***************************** pnum=4 ; pull out surface var=reform(dbz[sfc_idx,*,*]) ;MS 1 ; Get rid of the -9999.90 lats and lons for mask interpolation r=where(lat ne -9999.90 and lon ne -9999.90) good_colon=colon[r] good_lat=lat[r] good_var=var[r] ; Interpolate mask qhull,good_colon,good_lat,triangles,/delaunay,sphere=s grid_mask=griddata(good_colon,good_lat,good_var,xout=grid_colon,yout=grid_lat,$ /grid,/degrees,/sphere,triangles=triangles,$ /kriging,min_points=18,sectors=8,empty_sectors=3,missing=-8888) ; Get rid of the -9999.90 lats and lons and var for cloud interpolation r=where(lat ne -9999.90 and lon ne -9999.90 and var gt -9999,c) if c le 10 then r=where(lat ne -9999.90 and lon ne -9999.90) good_colon=colon[r] good_lat=lat[r] good_var=var[r] ; Triangulate the cloud data qhull,good_colon,good_lat,triangles,/delaunay,sphere=s grid_var=griddata(good_colon,good_lat,good_var,xout=grid_colon,yout=grid_lat,$ /grid,/degrees,/sphere,triangles=triangles,$ /kriging,min_points=8,sectors=8,empty_sectors=3,missing=-9999) ; Bytscal the data data_image=bytscl(grid_var,top=top_color,min=dmin_radar,max=dmax_radar) result=where(grid_var eq -9999,count) if count gt 0 then data_image[result]=255 ;gray ;result=where(grid_var lt 2*dmin_radar and grid_var ne -9999,count) result=where(grid_var lt dmin_radar and grid_var ne -9999,count) if count gt 0 then data_image[result]=255 ;gray result=where(grid_mask eq -8888,count) if count gt 0 then data_image[result]=253 ;white p0=image(data_image,grid_colon,grid_lat,limit=[llat_ovp,lcolon_ovp,ulat_ovp,rcolon_ovp],$ /current,dimensions=[pxdim,pydim],grid_units='degrees',/box_axes,$ center_longitude=(rcolon_ovp-lcolon_ovp)/2.0,rgb_table=mytable,$ map_projection='Mercator',position=pos[pnum,*],font_size=fs1) p1=mapcontinents(color='black') p0.mapgrid.box_color='black' p0.mapgrid.color='black' p0.mapgrid.label_color='black' p0.mapgrid.linestyle='dot' p0.mapgrid.label_position=0 for r=0,3 do begin p7=plot(lon_cir[r,*],lat_cir[r,*],/current,/overplot) endfor c0=colorbar(rgb_table=mycbtable,range=[dmin_radar,dmax_radar],orientation=0,$ position= reform(cbpos[pnum,*]),font_size=10,tickdir=1,border_on=1,title='dBZ') t1=text(pos[pnum,0]+dx,pos[pnum,3]+2*dy,'Surface',font_size=fs1) ; Plot curtain pnum=5 dbz1=reform(dbz[*,sub_pt,idx_ovp]) data_var=rotate(dbz1,3) ;doesn't work 0,2 ;upside down 1 data_image=bytscl(data_var,top=top_color,min=dmin_radar,max=dmax_radar) result=where(data_var lt -9000,count) if count gt 0 then data_image[result]=255 ;gray ; gpm height array height_radar=(findgen(nbin)*deltah)/1000.0 hidx=where(height_radar lt 10) p1=image(data_image[*,hidx],/current,image_dimensions=[isx,isy],position=pos[pnum,*],rgb_table=mytable) c1=contour(data_var[*,hidx],lat_sub[idx_ovp],height_radar[hidx],ytitle='height (km)',xtitle='latitude',$ xstyle=1,ystyle=1,/nodata,/current,position=pos[pnum,*],font_size=fs1,$ irregular=0,xtickdir=1,ytickdir=1,axis_style=2,xrange=[lat_sub[idx_ovp[0]],lat_sub[idx_ovp[-1]]]) ;c2=colorbar(orientation=1,rgb_table=mycbtable,range=[dmin_radar,dmax_radar],title='DbZ',$ ; position=reform(cbpos[pnum,*]),font_size=fs1,tickdir=1,border_on=1) t1=text(pos[pnum,0]+dx,pos[pnum,3]+1*dy,'Sub Satellite point',font_size=fs1) ;************************** ; Plot z factor measured ;************************** pnum=6 ; pull out surface var=reform(dbz_measured[sfc_idx-15,*,*]) ; Get rid of the -9999.90 lats and lons for mask interpolation r=where(lat ne -9999.90 and lon ne -9999.90) good_colon=colon[r] good_lat=lat[r] good_var=var[r] ; Interpolate mask qhull,good_colon,good_lat,triangles,/delaunay,sphere=s grid_mask=griddata(good_colon,good_lat,good_var,xout=grid_colon,yout=grid_lat,$ /grid,/degrees,/sphere,triangles=triangles,$ /kriging,min_points=18,sectors=8,empty_sectors=3,missing=-8888) ; Get rid of the -9999.90 lats and lons and var for cloud interpolation r=where(lat ne -9999.90 and lon ne -9999.90 and var gt -9999,c) if c le 10 then r=where(lat ne -9999.90 and lon ne -9999.90) good_colon=colon[r] good_lat=lat[r] good_var=var[r] ; Triangulate the cloud data qhull,good_colon,good_lat,triangles,/delaunay,sphere=s grid_var=griddata(good_colon,good_lat,good_var,xout=grid_colon,yout=grid_lat,$ /grid,/degrees,/sphere,triangles=triangles,$ /kriging,min_points=8,sectors=8,empty_sectors=3,missing=-9999) ; Bytscal the data data_image=bytscl(grid_var,top=top_color,min=dmin_radar,max=dmax_radar) result=where(grid_var eq -9999,count) if count gt 0 then data_image[result]=255 ;gray ;result=where(grid_var lt 2*dmin_radar and grid_var ne -9999,count) result=where(grid_var lt dmin_radar and grid_var ne -9999,count) if count gt 0 then data_image[result]=255 ;gray result=where(grid_mask eq -8888,count) if count gt 0 then data_image[result]=253 ;white p0=image(data_image,grid_colon,grid_lat,limit=[llat_ovp,lcolon_ovp,ulat_ovp,rcolon_ovp],$ /current,dimensions=[pxdim,pydim],grid_units='degrees',/box_axes,$ center_longitude=(rcolon_ovp-lcolon_ovp)/2.0,rgb_table=mytable,$ map_projection='Mercator',position=pos[pnum,*],font_size=fs1) p1=mapcontinents(color='black') p0.mapgrid.box_color='black' p0.mapgrid.color='black' p0.mapgrid.label_color='black' p0.mapgrid.linestyle='dot' p0.mapgrid.label_position=0 for r=0,3 do begin p7=plot(lon_cir[r,*],lat_cir[r,*],/current,/overplot) endfor c0=colorbar(rgb_table=mycbtable,range=[dmin_radar,dmax_radar],orientation=0,$ position=reform(cbpos[pnum,*]),font_size=10,tickdir=1,border_on=1,title='dBZ') t1=text(pos[pnum,0]+dx,pos[pnum,3]+2*dy,'Surface',font_size=fs1) t3=text(pos[pnum,0],0.98-4.*dy,'Reflectivity Factor as measured (dBZ)',font_size=fs1) ; Plot curtain pnum=7 dbz1=reform(dbz_measured[*,sub_pt,idx_ovp]) data_var=rotate(dbz1,3) ;doesn't work 0,2 ;upside down 1 result=where(data_var gt -9999.0) dmin_var=-10;min(data_var[result]) dmax_var=50;max(data_var[result]) data_image=bytscl(data_var,top=top_color,min=dmin_var,max=dmax_var) result=where(data_var lt -9000,count) if count gt 0 then data_image[result]=255 ;gray ; gpm height array height_radar=(findgen(nbin)*deltah)/1000.0 hidx=where(height_radar lt 10) p1=image(data_image[*,hidx],/current,image_dimensions=[isx,isy],position=pos[pnum,*],rgb_table=mytable) c1=contour(data_var[*,hidx],lat_sub[idx_ovp],height_radar[hidx],ytitle='height (km)',xtitle='latitude',$ xstyle=1,ystyle=1,/nodata,/current,position=pos[pnum,*],font_size=fs1,$ irregular=0,xtickdir=1,ytickdir=1,axis_style=2,xrange=[lat_sub[idx_ovp[0]],lat_sub[idx_ovp[-1]]]) ;c2=colorbar(orientation=1,rgb_table=mycbtable,range=[dmin_var,dmax_var],title='DbZ',$ ; position=reform(cbpos[pnum,*]),font_size=fs1,tickdir=1,border_on=1) t1=text(pos[pnum,0]+dx,pos[pnum,3]+1*dy,'Sub Satellite point',font_size=fs1) ;p0.save,imagedir+'temp_image.png',height=pydim ;stop ;****************** ; Get the cband radar ;****************** do_cband_distance='no' if do_cband_distance eq 'yes' then begin read_cband,julian_day_ship[result_ship[0]],dbz_grid_cband,dbz_mask_cband,x_out,y_out,$ nbins_cband,rstart,rscale,julian_day_cband,elangle_cband ; Plot reflectivity pnum=1 data_var=dbz_grid_cband ;result=where(data_var ne -9999,count) dmin=-20.0;min(data_var[result]) dmax=50.0;max(data_var[result]) data_image=bytscl(data_var,top=top_color,min=dmin,max=dmax) ;result=where(data_var eq -999.,count) result=where(data_var eq -9999,count) if count gt 0 then data_image[result]=253 result=where(dbz_mask_cband eq -8888,count) if count gt 0 then data_image[result]=255 p1=image(data_image,/current,image_dimensions=[isx,isy],position=pos[pnum,*],$ rgb_table=mytable) c1=contour(data_var,x_out/1000.0,y_out/1000.0,ytitle='Y (km)',xtitle='X (km)',$ xstyle=1,ystyle=1,font_size=14,/nodata,/current,position=pos[pnum,*],$ irregular=0,xtickdir=1,ytickdir=1,font_name='Helvetica',axis_style=2) ; plot range rings j=0 for r=0,nbins_cband[j]-1,150 do begin print,'range ring',r ; radius in km range_ring=(rstart[j]+float(r)*rscale[j])/1000.0 result=circle(0,0,range_ring) p2=plot(result,/current,/overplot) endfor c2=colorbar(orientation=1,rgb_table=mycbtable,range=[dmin,dmax],title='Z (dBZ)',$ position=reform(cbpos[pnum,*]),font_size=14,font_name='Helvetica',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)') scan_str=string(elangle_cband[i],format='(F05.2)') t1=text(pos[pnum,0]+dx,pos[pnum,3]+2*dy,'el='+scan_str,font_size=fs1) t1=text(pos[pnum,0]+14*dx,pos[pnum,3]+2*dy,date_time,font_size=fs1) ;p0.save,'temp_map.png',height=pydim endif ;end of do cband do_cband_latlon='no' if do_cband_latlon eq 'yes' and dbz_grid_cband ne !NULL then begin ; Done earlier for range rings ; read_cband_latlon,julian_day_ship[result_ship[0]],dbz_grid_cband,colon_cband,lat_cband,$ ; nbins_cband,rstart,rscale,julian_day_cband,elangle_cband,ulat,llat,llon,rlon,clon,clat ; Plot reflectivity pnum=1 data_var=dbz_grid_cband dmin=dmin_radar;-20.0;min(data_var[result]) dmax=dmax_radar;50.0;max(data_var[result]) data_image=bytscl(data_var,top=top_color,min=dmin,max=dmax) ;result=where(dbz_grid_cband gt -8000 and dbz_grid_cband le -7000,count) result=where(dbz_grid_cband gt -8000 and dbz_grid_cband lt dmin,count) if count gt 0 then data_image[result]=255 ;gray result=where(dbz_grid_cband gt -10000 and dbz_grid_cband le -9000,count) if count gt 0 then data_image[result]=253 ;white p0=image(data_image,colon_cband,lat_cband,limit=[llat,lcolon,ulat,rcolon],$ /current,dimensions=[pxdim,pydim],grid_units='degrees',$ center_longitude=(rcolon-lcolon)/2.0,$ rgb_table=mytable,$;,min_value=dmin,max_value=dmax,$ map_projection='Mercator',$ position=pos[pnum,*],font_size=fs1) p1=mapcontinents(color='black') p0.mapgrid.box_axes=1 p0.mapgrid.box_color='black' p0.mapgrid.color='black' p0.mapgrid.label_color='black' p0.mapgrid.linestyle='dot' p0.mapgrid.label_position=0 ; plot range rings lon_cir=fltarr(360) lat_cir=fltarr(360) for r=1,4 do begin ;radius of earth in meters arc_dist=(float(r)*25000.0)/(6371D*1000D) ;m for s=0,359 do begin result=ll_arc_distance([clon,clat],arc_dist,s,/degrees) lon_cir[s]=result[0] lat_cir[s]=result[1] endfor p7=plot(lon_cir,lat_cir,/current,/overplot) endfor c2=colorbar(orientation=1,rgb_table=mycbtable,range=[dmin,dmax],title='Z (dBZ)',$ position=reform(cbpos[pnum,*]),font_size=14,font_name='Helvetica',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)') scan_str=string(elangle_cband[0],format='(F05.2)') t1=text(pos[pnum,0]+dx,pos[pnum,3]+2*dy,'el='+scan_str,font_size=fs1) t1=text(pos[pnum,0]+14*dx,pos[pnum,3]+2*dy,date_time,font_size=fs1) ;p0.save,'temp_map.png',height=pydim endif else begin ;cband lat-lon pnum=1 t1=text(pos[pnum,0]+dx,pos[pnum,3]-2*dy,'no cband file',font_size=fs1) endelse ;************** ; GMI ;************** do_gmi='yes' if do_gmi eq 'yes' then begin pnum=0 read_gpm_1bgmi,gmi_dir,syear,sdoy3,ulat_ovp,llat_ovp,llon_ovp,rlon_ovp,lcolon_ovp,rcolon_ovp,$ colon_gmi,lat_gmi,tb_gmi,channel_gmi,do_sub,forbit,chan dmin_gmi=180 & dmax_gmi=265 chan8=reform(tb_gmi) idx_ovp_gmi=where(lat_gmi ge llat_ovp and lat_gmi le ulat_ovp and $ colon_gmi ge lcolon_ovp and colon_gmi le rcolon_ovp,count_ovp_gmi) data_var=chan8[idx_ovp_gmi] lon_var=colon_gmi[idx_ovp_gmi] lat_var=lat_gmi[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_sub[idx_ovp],lat_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]+0*dy,$ ; 'GMI Channel: '+channel_gmi+' Kelvin',font_size=fs1) ;; Triangulate the data ;print,'start qhull' ;qhull,colon_gmi,lat_gmi,triangles,/delaunay,sphere=s ;print,'end qhull, start gridding' ;grid_tb=griddata(colon_gmi,lat_gmi,tb_gmi,xout=grid_colon,yout=grid_lat,$ ; /grid,/degrees,/sphere,$ ; /KRIGING,min_points=10,sectors=8,empty_sectors=4,$ ; missing=-9999,$ ; triangles=triangles);,$ ; ;missing=!values.f_nan) ;print,'end gridding' ; ;; Bytscal the data ;;infrared_colortable,ir_table,ir_max,ir_min ;result=where(grid_tb ne -9999) ;dmin=min(grid_tb[result]) ;190 ;dmax=max(grid_tb[result]) ;265 ;data_image=bytscl(grid_tb,top=top_color,min=dmin,max=dmax) ;result=where(grid_tb eq -9999,count) ;if count gt 0 then data_image[result]=253 ;p0=image(data_image,grid_colon,grid_lat,limit=[llat_ovp,lcolon_ovp,ulat_ovp,rcolon_ovp],$ ; /current,dimensions=[pxdim,pydim],grid_units='degrees',$ ; center_longitude=(rcolon_ovp-lcolon_ovp)/2.0,$ ; rgb_table=mytable,$;,min_value=dmin,max_value=dmax,$ ; map_projection='Mercator',$ ; position=pos[pnum,*],font_size=fs1) ; ;p1=mapcontinents(color='black') ;p0.mapgrid.box_axes=1 ;p0.mapgrid.box_color='black' ;p0.mapgrid.color='black' ;p0.mapgrid.label_color='black' ;p0.mapgrid.linestyle='dot' ;p0.mapgrid.label_position=0 ;c0=colorbar($ ; rgb_table=mycbtable,$ ; range=[dmin,dmax],$ ; orientation=0,position= reform(cbpos[pnum,*]),font_size=fs1,$ ; tickdir=1,border_on=1,title='Kelvin') ; ;for r=0,3 do begin ; p7=plot(lon_cir[r,*],lat_cir[r,*],/current,/overplot) ;endfor t3=text(0.1,pos[0,3]+2.0*dy,'GMI Channel: '+channel_gmi,font_size=fs1) endif ;end of do gmi ;p0.save,imagedir+'temp_image.png',height=pydim ;************** ; ku power ;************** do_pku='yes' if do_pku eq 'yes' then begin 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 ; Triangulate the data print,'start qhull' qhull,colon_pow,lat_pow,triangles,/delaunay,sphere=s print,'end qhull, start gridding' ; pull out surface bin_echo_peak_sub=reform(bin_echo_peak[sub_pt,*]) sfc_echo_idx=min(bin_echo_peak_sub[idx_ovp]) print, sfc_echo_idx,'sfc_echo_idx' var=reform(pow[sfc_echo_idx-10,*,*]) grid_var=griddata(colon_pow,lat_pow,var,xout=grid_colon,yout=grid_lat,$ /KRIGING,min_points=10,sectors=8,empty_sectors=4,$ missing=-9999,triangles=triangles,/grid,/degrees,/sphere);,$ ; Bytscal the data result=where(grid_var gt -290,count) dmin=-114;min(var[result]);-120 dmax=-95;max(var[result]);-20 data_image=bytscl(grid_var,top=top_color,min=dmin,max=dmax) result=where(grid_var eq -299.99,count) if count gt 0 then data_image[result]=253 result=where(grid_var eq -9999,count) if count gt 0 then data_image[result]=255 result=where(grid_var eq -300,count) if count gt 0 then data_image[result]=254 p0=image(data_image,grid_colon,grid_lat,limit=[llat_ovp,lcolon_ovp,ulat_ovp,rcolon_ovp],$ /current,dimensions=[pxdim,pydim],grid_units='degrees',$ center_longitude=(rcolon_ovp-lcolon_ovp)/2.0,$ rgb_table=mytable,$;,min_value=dmin,max_value=dmax,$ map_projection='Mercator',$ position=pos[pnum,*],font_size=fs1) p1=mapcontinents(color='black') p0.mapgrid.box_axes=1 p0.mapgrid.box_color='black' p0.mapgrid.color='black' p0.mapgrid.label_color='black' p0.mapgrid.linestyle='dot' p0.mapgrid.label_position=0 c0=colorbar($ rgb_table=mycbtable,$ range=[dmin,dmax],$ orientation=0,position= reform(cbpos[pnum,*]),font_size=fs1,$ tickdir=1,border_on=1,title='dBm') t3=text(pos[pnum,0],pos[pnum,3]+3.*dy,'L1 Echo Power (dBm)',font_size=fs1) ; Plot curtain pnum=3 pow1=reform(pow[*,sub_pt,idx_ovp]) data_var=rotate(pow1,3) ;doesn't work 0,2 ;upside down 1 ;*** cut at surface ;data_var=data_var[*,min(bin_echo_peak_sub[idx]):-1] sfc_echo_idx=min(bin_echo_peak_sub[idx_ovp]) data_var=data_var[*,259-sfc_echo_idx:259] s=size(data_var,/dimensions) nbin_sfc=s[1] ; Bytscal the data result=where(data_var gt -290,count) dmin=-114;min(data_var[result]);-120 dmax=-95;max(data_var[result]);-20 data_image=bytscl(data_var,top=top_color,min=dmin,max=dmax) result=where(data_var eq -299.99,count) if count gt 0 then data_image[result]=253 result=where(data_var eq -9999,count) if count gt 0 then data_image[result]=255 result=where(data_var eq -300,count) if count gt 0 then data_image[result]=254 ; gpm height array height_radar=(findgen(nbin_sfc)*deltah)/1000.0 hidx=where(height_radar lt 10) p1=image(data_image[*,hidx],/current,image_dimensions=[isx,isy],position=pos[pnum,*],rgb_table=mytable) c1=contour(data_var[*,hidx],lat_sub[idx_ovp],height_radar[hidx],ytitle='height (km)',xtitle='latitude',$ xstyle=1,ystyle=1,/nodata,/current,position=pos[pnum,*],font_size=fs1,$ irregular=0,xtickdir=1,ytickdir=1,axis_style=2,xrange=[lat_sub[idx_ovp[0]],lat_sub[idx_ovp[-1]]]) ;c2=colorbar(orientation=1,rgb_table=mycbtable,range=[dmin_radar,dmax_radar],title='dBm',$ ; position=reform(cbpos[pnum,*]),font_size=fs1,tickdir=1,border_on=1) t1=text(pos[pnum,0]+dx,pos[pnum,3]+1*dy,'Sub Satellite point',font_size=fs1) endif ;end of pku ;p0.save,imagedir+'temp_image.png',height=pydim ;stop ;************** ; Image strings ;************** stime=min(julian_day) etime=max(julian_day) 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,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) t3=text(0.1,0.98-4.*dy,'Reflectivity Factor with attenuation correction (dBZ)',font_size=fs1) imagename='ovp.'+datastream+'.dbz.'+$ odate_str+'.'+otime_str2+'.'+forbit if do_gmi eq 'yes' then imagename=imagename+'.'+channel_gmi imagename=imagename+'.ovp.png' print,imagedir+imagename p0.save,imagedir+imagename,height=pydim ; Write out ovp to a file openw,1,imagedir+ovp_filename,/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 ;end of profiles in overpass region endif ;profiles in region endif ; ship in this overpass file endfor ;end of loop through files endif ;files exist end