;------------------------------------------------------------- ;+ ; NAME: ; SAT_TRAJEC2 ; ; PURPOSE: ; overlay flight track of aircraft and few backtrajectories ; and loop through satellite images to detect cloud encounter ; (built for LNOX/PNOX analysis of PEM-Tropics) ; VERSION with 1/2 day time interval and IR images ; ; CATEGORY: ; data analysis ; ; CALLING SEQUENCE: ; SAT_TRAJEC2 [,keywords] ; ; INPUTS: ; FLIGHTDATA --> array with the merged data for the selected flight ; or all flights. Must be accompanied with FLIGHTHEADER ; parameter. If this parameter contains less than 2 elements, ; the file ~mgs/terra/chem1d/INPUT/allnmhc.data is read. ; ; FLIGHTHEADER --> vector with variable names of FLIGHTDATA. Must ; contain the following entries: FLIGHT, LON, LAT, ALTP, PSMB, ; and JDAY ; ; TRAJDATA --> trajetory data for this case. This is a structure ; containing the fields LON, LAT, TIME, ALT, and FLAG as ; returned by READ_TRAJ.PRO. Don't worry about it, simply ; pass a named variable (so it won't have to be read every ; time) and use the /LOADNEW keyword if you move to another day. ; ; KEYWORD PARAMETERS: ; FLIGHTNO --> flight number to be analyzed ; ; LON, LAT, ALT --> geographic location of point to be analyzed ; ; /LOADNEW --> set thiskeyword to force reading of new trajectory ; data. (sure: the program could keep track of day changes, ; but I want to get some sleep as well ;-) ; ; OUTPUTS: ; plots a satellite image and overlays flight TRAJEC of ; aircraft ; ; SUBROUTINES: ; ; REQUIREMENTS: ; uses procedure MAP_IMAGE ; ; NOTES: ; ; EXAMPLE: ; ; MODIFICATION HISTORY: ; mgs, 28 Jan 1998: VERSION 1.00 ; ;- ; Copyright (C) 1998, Martin Schultz, Harvard University ; This software is provided as is without any warranty ; whatsoever. It may be freely used, copied or distributed ; for non-commercial purposes. This copyright notice must be ; kept with any copy of this software. If this software shall ; be used commercially or sold as part of a larger package, ; please contact the author to arrange payment. ; Bugs and comments should be directed to mgs@io.harvard.edu ; with subject "IDL routine sat_trajec" ;------------------------------------------------------------- pro sat_trajec2,flightdata,header,trajdata,flightno=flightno, $ lon=lon,lat=lat,alt=alt,loadnew=loadnew ; program requires FLIGHTNO, LON, LAT, and ALT if (n_elements(flightno) le 0 OR n_elements(LON) le 0 $ OR n_elements(LON) le 0 OR n_elements(LON) le 0) then $ message,'SAT_TRAJEC: requires FLIGHTNO, LON, LAT, and ALT keywords.' ; parse parameters, set defaults if (n_elements(waittime) le 0) then waittime = 2. ; if flightdata is not passed, read it now if(n_elements(flightdata) lt 2) then $ readdata,'~/terra/chem1d/INPUT/allnmhc.data',flightdata,header, $ delim=' ',skp1=1,skp2=3 ; if trajdata is not passed, read it now, else extract from structure if(n_elements(trajdata) le 0 OR keyword_set(loadnew)) then begin read_traj,'/data/pem-t/traj_dc8/fs101d'+ $ string(flightno,format='(i2.2)')+ $ '.pmt',ttime,tlon,tlat,talt,flag=flag trajdata = { time:ttime, lon:tlon, lat:tlat, alt:talt, flag:flag } endif else begin ; check if trajdata is valid structure test = size(trajdata) if (test(test(0)+1) ne 8) then begin print,'trajdata is not a valid structure.' stop endif ttime = trajdata.time tlon = trajdata.lon tlat = trajdata.lat talt = trajdata.alt flag = trajdata.flag endelse ; extract position variables from aircraft data set aflight = flightdata(where(header eq 'FLIGHT'),*) aday = flightdata(where(header eq 'JDAY'),*) alon = flightdata(where(header eq 'LON'),*) alat = flightdata(where(header eq 'LAT'),*) aalt = flightdata(where(header eq 'ALTP'),*) apsmb = flightdata(where(header eq 'PSMB'),*) aind = where(aflight eq flightno) if (aind(0) lt 0) then $ message,'Invalid FLIGHTNO !' ; extract day of flight jday = aday(aind(0)) ; retrieve (mean) pressure level for altitude tmp = aalt(aind) ind = where(abs(alt-tmp) lt 0.15) if (ind(0) lt 0) then begin print,'cannot find pressure for altitude ',alt stop endif meanp = apsmb(aind) meanp = total(meanp(ind))/float(n_elements(ind)) ; look for trajectory starting points near case under investigation tlon0 = tlon(0,*) tlat0 = tlat(0,*) talt0 = talt(0,*) tind = where(abs(lon-tlon0) lt 2. AND abs(lat-tlat0) lt 2. $ AND abs(meanp-talt0) lt 50.) print,' Selected ',n_elements(tind)*(tind(0) ge 0),' trajectories.' if (tind(0) lt 0) then begin print,'Cannot find suitable trajectories.' stop endif ttime = ttime(*,tind) tlon = tlon(*,tind) tlat = tlat(*,tind) talt = talt(*,tind) flag = flag(*,tind) ; load color table grey scale and change first 20 colors myct,0 ; load grey scale color table ; gamma_ct,1.6 ; enhance contrast ; open graphics window of nice size window,0,xsize=900,ysize=850 ; set up animation XInterAnimate,set=[900,850,11],/cycle,/track, $ title='Trajectory cloud chasing' ; loop through satellite images and plot composite image ; use only 1 daily visible image for now for i=0,9 do begin sday = string(jday-(i+1)/2,format='(i3.3)') stime = '1500' if (fix(i/2.) eq i/2.) then stime = '0300' satfile = '~/download/gte/'+sday+'_'+stime+'ful3.jpg' if (not file_exist(satfile)) then begin print,'cannot find satellite image '+satfile goto,nextone endif read_jpeg,satfile,satim ; cut off border ; values for visible image 2100ful1 a = satim if (stime eq '2100') then a = satim(93:1041,23:927) ; values for infrared image 0300ful3 if (stime eq '0300') then begin a = satim(0:744-1-66,24:*) endif ; values for infrared image 1500ful3 if (stime eq '1500') then a = satim(66:*,24:*) p = [ 0.02, 0.02, 0.98, 0.98 ] tvimage,a,position=p,/keep_aspect,ncolors=60,bottom=20 !p.position=p map_set,0,-135.,/satellite,sat_p=[5.5,0.,-.5],/noerase map_continents,color=1 map_grid,color=1,latdel=15,londel=15 ; overlay flight tracks, higlight current flight oplot,alon,alat,color=5 oplot,alon(aind),alat(aind),color=3 ; overlay trajectories, highlight current time points for t=0,n_elements(tind)-1 do begin okind = where(flag(*,t) eq 0) oplot,tlon(okind,t),tlat(okind,t),color=4,thick=2 endfor for t=0,n_elements(tind)-1 do begin dind = where(ttime(*,t) eq i/2.) if (dind(0) ge 0) then $ oplot,tlon(dind,t),tlat(dind,t),color=2, $ psym=sym(1),symsize=0.8 endfor ; wait,waittime nextone: print,'******** FRAME = ',i ; load image in XInterAnimate buffer XInterAnimate,frame=i,window=0 endfor ; image loop ; start animation XInterAnimate return end