pro gpm_overpasses_2019 ;************************************* ; Input ;************************************* ; This is the region covering the entire ship track of the experiment region='investigator_2019' if region eq 'investigator_2019' then begin ulat=-42.0 & llat=-67.0 & llon=137.0 & rlon=153.0 ;big region for 2019 shipdir='/uufs/chpc.utah.edu/common/home/u0079358/gpm/' shipfilename='investigator_ship_track.20190118.210915.20190304.211315.1min.cdf' 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='/uufs/chpc.utah.edu/common/home/mace-group4/field_programs/socrates/ship_track/' shipfilename='socrates_ship_track.20180110.232300.20180221.052900.cdf' endif ; If you put in a gpm orbit, this is the only one that will run ;forbit='022218' ;20180125 ;forbit='022245' ;20180127 ;no cband ;forbit='022263' ;20180128 ;forbit='022368' ;20180204 ;forbit='022370' ;20180204 ;forbit='022398' ;20180206 ;forbit='022429' ;20180208 ;forbit='022432' ;20180208 ;forbit='022445' ;20180209 ;forbit='022448' ;20180209 ;forbit='022460' ;20180210 ;forbit='022463' ;20180210 ;forbit='022479' ;20100211 ;forbit='022494' ;20180212 ;forbit='022509' ;20180213 ;forbit='022522' ;20180214 ;forbit='022525' ;20180214 ;forbit='022567' ;20180217 ;forbit='027903' ;20190126.06 ;forbit='027932' ;20190128.04 forbit='027947' ;20190129.03 ;forbit='027949' ;20190129.05 output_dir='' ;file output directory imagedir='/uufs/chpc.utah.edu/common/home/u0079358/gpm/' ;imagedir='/uufs/chpc.utah.edu/common/home/u0079358/public_html/gpm/overpasses_power/' ;************************************* ; Constants ;************************************* ; directory path_prefix='/uufs/chpc.utah.edu/common/home/' ;chpc ; Choose which radar radar='ku' ;radar='ka' if radar eq 'ku' then begin 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 endif else if radar eq 'ka' then begin radarname='KaPR' radarfreq='35.5 GHz' radarswath='120 km' gpm_dir='/mace-group5/gpm/GPM_L2/GPM_2AKa.05/' ; Sub satellite point sub_pt=12 ;deltah=250 ;m for MS deltah=500 ;m for HS sfc_idx=87 ;surface plane endif ; 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 directory 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 ;****************************************************************************************** ; 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 ;****************************************************************************************** ; 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 then begin file_string='2A.GPM.Ku.V8-20180723.2019*HDF5' ;file_string='2A.GPM.Ku.V8-20180723.*HDF5' endif else begin ;file_string='2A.GPM.Ku.V8-20180723.2019*'+forbit+'*.HDF5' file_string='2A.GPM.Ku.V8-20180723.*'+forbit+'*.HDF5' endelse ; Find matching gpm files filesd=file_search(fdir+file_string,count=num_files) if num_files eq 0 then print,'NO FILES AT ALL' ;*********************************************** ; Loop through each gpm file ;*********************************************** ; If gpm files exist if num_files gt 0 then begin ; Loop through gpm files for i=0,num_files-1 do begin print,filesd[i] ; GPM orbit number 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 ; nray has to be the correct size to continue if nray ge sub_pt*2 then begin ; ; 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 and lon ne -9999.90,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,*]) ;***************** ; 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=where(lat_sub ge llat and lat_sub le ulat and $ colon_sub ge lcolon and colon_sub le rcolon and $ julian_day ge julian_day_ship[0] and $ julian_day le julian_day_ship[-1],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' ; 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 caldat,julian_day[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[-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_sub),/float,value=1e7) for ii=0,n_elements(lon_sub)-1 do begin if lon_sub[ii] gt -9999 and lat_sub[ii] gt -9999 then begin delta_dist[ii]=map_2points(lon_ship_ovp,lat_ship_ovp,$ lon_sub[ii],lat_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[didx] lat_ovp=lat_sub[didx] lon_ovp=lon_sub[didx] colon_ovp=colon_sub[didx] print,'min_dist',min_dist/1000.0,'km' print,'min_dist idx',didx ; Calculate ovp box to be delta_degrees aroud the ship delta_degree=0.5 ;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 ; Check again to see if closest approach is close to ship 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 eq 0 then begin print,'no satellite in ship ovp box' endif if count_ovp gt 1 then begin print,'satellite data in 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_sub ge llat_ovp and lat_sub le ulat_ovp and $ colon_sub ge lcolon_ovp and colon_sub le rcolon_ovp,count_ovp) ; plot range rings around ship location 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 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 ;************************ ; 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 ;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([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]) ;***************************** ; Plot map - Ku or Ka - dbz ;***************************** pnum=4 ; pull out surface var=reform(dbz[sfc_idx,*,*]) ;MS 1 if 1 eq 1 then begin ; 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=4,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,c) 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 print,'start griddata' grid_var=griddata(good_colon,good_lat,good_var,xout=grid_colon,yout=grid_lat,$ /grid,/degrees,/sphere,$ /kriging,min_points=8,sectors=8,empty_sectors=4,$ missing=-9999,triangles=triangles);$ ; Bytscal the data dmin_radar=10;min(grid_var[result]);10 to 35? dmax_radar=35;max(grid_var[result]) 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 ;white result=where(grid_mask eq -8888,count) if count gt 0 then data_image[result]=253 ;gray p0=image(data_image,grid_colon,grid_lat,$ limit=[llat_ovp,lcolon_ovp,ulat_ovp,rcolon_ovp],$ center_longitude=(rcolon_ovp-lcolon_ovp)/2.0,$ /current,dimensions=[pxdim,pydim],grid_units='degrees',$ rgb_table=mytable,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 endif ;try new plotting 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='black',linestyle=2) 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 - Ku or Ka - dbz ;***************************** 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',irregular=0,xtickdir=1,ytickdir=1,$ xstyle=1,ystyle=1,/nodata,/current,position=pos[pnum,*],font_size=fs1,$ 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 map - Ku or Ka - z factor measured ;************************** pnum=6 ; pull out surface var=reform(dbz_measured[sfc_idx-15,*,*]) ; DO NOT TRIANGULATE AGAIN - counterclockwise error ;qhull,newcolon,newlat,triangles,/delaunay,sphere=s ;triangulate,colon,lat,triangles,/degrees,sphere=s ; 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) good_colon=colon[r] good_lat=lat[r] good_var=var[r] ; Triangulate the cloud data triangles=!null qhull,good_colon,good_lat,triangles,/delaunay,sphere=s print,'start griddata' grid_var=griddata(good_colon,good_lat,good_var,xout=grid_colon,yout=grid_lat,$ /grid,/degrees,/sphere,$ ;method='NaturalNeighbor',$ ;/inverse_distance,min_points=10,sectors=8,empty_sectors=3,$ ;works ;/inverse_distance,min_points=15,sectors=8,empty_sectors=3,max_per_sector=3,$ /kriging,min_points=8,sectors=8,empty_sectors=4,$ missing=-9999,triangles=triangles);$ print,'end griddata' ; Bytscal the data dmin_radar=-10;min(grid_var[result]);10 to 35? dmax_radar=50;35;max(grid_var[result]) 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 ;white result=where(grid_mask eq -8888,count) if count gt 0 then data_image[result]=253 ;gray 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,$ 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 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='black',linestyle=2) 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 - Ku or Ka - z factor measured ;***************************** 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) ;****************** ; Get the cband radar ;****************** do_cband_pyart='yes' if do_cband_pyart eq 'yes' then begin read_cband_pyart,julian_day_ship[result_ship[0]],$ ulat_ovp,llat_ovp,llon_ovp,rlon_ovp,dbz_cband,$ lat_cband,lon_cband,alt_cband,julian_day_cband,colon_cband if dbz_cband ne !null then begin pnum=0 sidx=4 dmin_cband=-5 & dmax_cband=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,count) if count gt 0 then data_image[result]=255 ;gray p1=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=(rcolon_ovp-lcolon_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) p2=mapcontinents(color='black') ;p0.mapgrid.box_axes=1 ;p0.mapgrid.box_color='black' ;0.mapgrid.color='black' ;p0.mapgrid.label_color='black' p1.mapgrid.linestyle='dot' ;p0.mapgrid.label_position=0 p8=plot(lon_sub[idx_ovp],lat_sub[idx_ovp],/overplot,color='black',linestyle=2,thick=2) t1=text(pos[pnum,0]-1*dx,pos[pnum,3]+1*dy,string(alt_cband[0,0,sidx]/1000.0,format='(F4.2)')+'km',$ font_size=fs1) for r=0,3 do begin p7=plot(lon_cir[r,*],lat_cir[r,*],/current,/overplot) endfor c2=colorbar(orientation=0,rgb_table=mycbtable,range=[dmin_cband,dmax_cband],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,1]-2*dy,'el='+scan_str,font_size=fs1) t1=text(((pos[pnum,2]-pos[pnum,0])/2)+pos[pnum,0],$ pos[pnum,1]-4*dy,date_time,font_size=fs1) t1=text(((pos[pnum,2]-pos[pnum,0])/2)+pos[pnum,0],$ pos[pnum,1]-6*dy,'CBand',font_size=fs1) t1=text(pos[pnum,0]+1*dx,pos[pnum,3]-1*dy,string(alt_cband[0,0,sidx]/1000.0,format='(F4.2)')+'km',$ font_size=fs1) ;************************** ; Plot cband curtain ;************************** pnum=1 ; 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_sub[idx_ovp[k]]-lat1),idx_lat) closest_lon=min(abs(lon_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,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],$ ;min_value=dmin_cband,max_value=dmax_cband,$ position=pos[pnum,*],rgb_table=mytable) c1=contour(data_var[*,hidx],lon_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_sub[idx_ovp],make_array(count_ovp,value=height_cband[sidx]),$ /overplot,color='gray',thick=3,linestyle=4) ;p10=plot(x,y,overplot=p1,color='gray',thick=3,linestyle=4) ;t1=text(pos[pnum,0]+dx,pos[pnum,3]+2*dy,'Constant Latitude '+string(xlat,format='(F6.1)'),font_size=fs1) ; Colorbar ;cb0=colorbar(rgb_table=mycbtable,range=[dmin_cband,dmax_cband],$ ; orientation=0,position=[0.55,0.97,0.9,0.98],font_size=fs1,$ ; tickdir=1,border_on=1,title='dBZ') p0.save,'test_test2.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 endif ;do the cband plot ;************** ; GMI ;************** do_gmi='no' 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 if n_elements(tb_gmi) gt 5 then begin ;; Calculate regular grid array ;delta=0.05 ;ynum=(ulat-llat)/delta ;grid_lat=findgen(ynum)*delta+llat ;xnum=(rcolon-lcolon)/delta ;grid_colon=findgen(xnum)*delta+lcolon ;grid_lon=grid_colon ;result=where(grid_lon gt 180,count) ;if count gt 0 then grid_lon[result]=grid_lon[result]-360.0 ; 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') ;if count_ship gt 0 then $ ; p1=plot(lon_ship[result_ship],lat_ship[result_ship],/overplot,color='black',thick=3,symbol='o') 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 gmi points exist. 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,$ lat_sub_pow,colon_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) ; 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_pow]) 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 for r=0,3 do begin p7=plot(lon_cir[r,*],lat_cir[r,*],/current,/overplot) endfor 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_pow]) 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_pow]) 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_pow[idx_ovp_pow],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_pow[idx_ovp_pow[0]],lat_sub_pow[idx_ovp_pow[-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.4.cb.png' ;imagedir='/uufs/chpc.utah.edu/common/home/u0079358/gpm/' print,imagedir+imagename p0.save,imagedir+imagename,height=pydim ; Write out ovp to a file ;openw,1,imagedir+'ship_gpm_overpasses_2018.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 ;end of profiles in overpass region endif ;profiles in region endif ; ship in this overpass file endif ;subset is good, nbins is right number endfor ;end of loop through files endif ;files exist end