PRO C_bd_VAD_s_ang ;-------Last changed on March 19, 2009 by Shelby Frisch ;----Much of the data reading part was made possible by Irina Djalaova Common Block1,a0,a1,b1,vrr,theta; COMMON BLOCK2,flag_thresh,rnsc_win_width COMMON BLOCK3,err,nmpts COMMON Block4,flag COMMON Block5,ht,outfile_1 COMMON Block7,count COMMON Block8,el COMMON Block9,dzzz,dzzthreshold ;=========================================================================================== ;This program uses the NASA software for reading the Ron Browm C-BAND radar files ; It calls a RANSAC type routine to clean up the radial velocity data by randomly taking ;radial velocity data from the 3- 120 degree sectors, one point from each sector. It finds the coef. ;for representing the radial velocity as vr=a0+a1*sin(theta)+b1*cos(theta), that is a0,a1,b1 ;It uses this fit to make a window of vr+delta_vr and Vr-delta_vr. It accepts only data in this window ; for a least squares fit of another representation vr'=a0'+a1'*sin(theta)+b1'*cos(theta) ; We compute a metric r^2=(vr(windowed)-vr')^2 and use the window for the data that yeilds a minimum of ;r^2. This should remove most errors caused by birds, etc in the data. A swarm of insects that cover ; a large fraction of the sweep might cause large errors however, this should be evident by doing a ;VAD display. ;--------------flag_thresh--------------- ;the data flag is for known bad data the flag_thresh should be less than this flag_thresh=-20. ;---------------rnsc_win_width------------------ ; This should be some fraction of nyquist, rnsc_win_width=1.5 ;---------------num_val -------------------- ;the number of non thresholded values to be used in each 120 degree sector num_val=10.0 ;----------------err------------------- ;This is the rms difference between the fitted curve and the actual ransac selected data ;======================================================== ;-----Input and output files ;---must set up the path for the input files to inputfile statement, outfile path for the ;---u and v components, etc, and outfile_2 for some plots of the vad radial velocity with ; the raw data and random consensus selected data. ;files= FILE_SEARCH('G:\Work Dan\TEXAQS_2006\working folder\vad\Input\August 8\*.*',$ ;/ISSUE_ACCESS_ERROR, Count=numFiles) ;======== INPUT FILE PATH============================ ;print,'input the input file path' ;inputfile='' ;READ,inputfile files= FILE_SEARCH('/cruises/TexAQS_2006/RHB/radar/cband/Raw/bad_files/*.*',/ISSUE_ACCESS_ERROR, Count=numFiles) ;======================================================== ;==========OUTPUT DATA FILE PATH============== ;print,'input the output file path' outfil='/cruises/TexAQS_2006/RHB/radar/cband/Raw/bad_files/' ;READ,outfil ;================================================ ;=================OUTPUT PLOT PATH============== ;print,'input the dataplot file path' plotfile='/cruises/TexAQS_2006/RHB/radar/cband/Raw/bad_files/' ;READ,plotfile ;========================================= print,'number of files=',numfiles ;===============SET up output arrays====================== Ux=FLTARR(100) ; East wind component m/s Uy=FLTARR(100) ; North wind component m/s height=FLTARR(100) ; this is the height of the range gate in meters elevation=FLTARR(30) ; radar elevation angle of the vad sweep in degrees Speed=FLTARR(100) ; this is the wind speed in m/s Theta=FLTARR(100) ; this is the direction of the wind speed from the North clockwise in degrees error=FLTARR(100) ; This is a measure of the lack of fit of the model to the data in m/s ct=FLTARR(100); This is the number of points used at each height in the calculations ;===================================================================== ;================Set up the file reading ==================== FOR kkk=0,numfiles-1 DO BEGIN filename=files(kkk) print,'FILENAME=',filename radar=rsl_anyformat_to_radar(filename) ;====================Write time and date information for checking as the program runs ============= PRINT,'TIME (MONTH DAY YEAR HOUR MINUTE SECOND): ' PRINT,radar.h.month,radar.h.day,radar.h.year,radar.h.hour,radar.h.minute,radar.h.sec strmon=STRING(radar.h.month) sstrmon=STRCOMPRESS(strmon,/REMOVE_ALL) strday=STRING(radar.h.day) sstrday=STRCOMPRESS(strday,/REMOVE_ALL) stryr=STRING(radar.h.year) sstryr=STRCOMPRESS(stryr,/REMOVE_ALL) strhr=STRING(radar.h.hour) sstrhr=STRCOMPRESS(strhr,/REMOVE_ALL) strmin=STRING(radar.h.minute) sstrmin=STRCOMPRESS(strmin,/REMOVE_ALL) strsec=STRING(radar.h.sec) sstrsec=STRCOMPRESS(strsec,/REMOVE_ALL) PRINT,'RADAR POSITION: ' PRINT,radar.h.latd+radar.h.latm/60.+radar.h.lats/3600.,$ radar.h.lond+radar.h.lonm/60.+radar.h.lons/3600. PRINT,'CALIBRATION CONSTS: ',radar.volume[0].h.calibr_const PRINT,'MISSING DATA FLAG: ',radar.volume[0].h.no_data_flag nvolumes=radar.h.nvolumes PRINT,'FIELD MIN MAX' FOR ivolume=0,nvolumes-1 DO BEGIN maxvar=MAX(radar.volume[ivolume].sweep.ray.range) data=radar.volume[ivolume].sweep.ray.range index=WHERE(data NE radar.volume[0].h.no_data_flag,count) IF count GT 0 THEN minvar=MIN(data[index]) ELSE minvar=radar.volume[0].h.no_data_flag PRINT,radar.volume[ivolume].h.field_type,minvar,maxvar IF radar.volume[ivolume].h.field_type EQ 'DZ' THEN DZvolume=ivolume IF radar.volume[ivolume].h.field_type EQ 'VR' THEN VRvolume=ivolume IF radar.volume[ivolume].h.field_type EQ 'SW' THEN SWvolume=ivolume ENDFOR ;=================Number of sweeps in the volume, number of rays in every sweep and number of bins in every ray nsweeps=radar.volume[VRvolume].h.nsweeps nrays =radar.volume[VRvolume].sweep[0].h.nrays nbins =radar.volume[VRvolume].sweep[0].ray[0].h.nbins PRINT,'nsweeps=',nsweeps,' nrays=',nrays,' nbins=',nbins ;Meters to the first gate and gate_range PRINT,'first_gate=',radar.volume[VRvolume].sweep[0].ray[0].h.range_bin1,'(Meters)',$ ' gate_size=',radar.volume[VRvolume].sweep[0].ray[0].h.gate_size,'(Meters)' ;==============================Assume that we are extracting the data of the first 10 Km maxrange=10000 ;in Meters ngates=(maxrange-radar.volume[VRvolume].sweep[0].ray[0].h.range_bin1)/radar.volume[VRvolume].sweep[0].ray[0].h.gate_size PRINT,'maxrange=',maxrange,'(Meters) ngates for output=',ngates ;=======================Output Arrays============================== ;TimeARRAY time_array=FLTARR(6,nrays,nsweeps) ;0-month,1-day,2-year,3-hour,4-minute,5-second ;Naval Variables ARRAY navy_array=FLTARR(10,nrays,nsweeps) ;0-pitch,1-roll,2-heading,3-pitch_rate,4-roll_rate,5-heading_rate,6-vel_east, ; 7-vel_north,8-vel_up ,9- vr correction ;DataARRAY data_array=FLTARR(nvolumes,ngates,nrays,nsweeps) ; ============================Output files outfile=outfil+sstrmon+'_'+sstrday+'_'+sstryr+'_'$ +sstrhr+'_'+sstrmin+'_'+sstrsec+'.dat' outfile_1=plotfile+sstrmon+'_'+sstrday+'_'+sstryr+'_'+sstrhr+'_'+sstrmin+'_'+sstrsec ;================print output file information for checking ================ print,'outfile=',outfile print,'plotfile=',outfile_1 ;=====================open output U and V data file =========================== OpenW,iu,outfile,/get_lun FOR isweep=0,nsweeps-1 DO BEGIN PRINT,'isweep=',isweep,' sweep_elevation=',radar.volume[VRvolume].sweep[isweep].ray[0].h.fix_angle FOR iray=0,nrays-1 DO BEGIN ray1=radar.volume[VRvolume].sweep[isweep].ray[iray] time_array[0,iray,isweep]=ray1.h.month time_array[1,iray,isweep]=ray1.h.day time_array[2,iray,isweep]=ray1.h.year time_array[3,iray,isweep]=ray1.h.hour time_array[4,iray,isweep]=ray1.h.minute time_array[5,iray,isweep]=ray1.h.sec navy_array[0,iray,isweep]=ray1.h.pitch navy_array[1,iray,isweep]=ray1.h.roll navy_array[2,iray,isweep]=ray1.h.heading navy_array[3,iray,isweep]=ray1.h.pitch_rate navy_array[4,iray,isweep]=ray1.h.roll_rate navy_array[5,iray,isweep]=ray1.h.heading_rate navy_array[6,iray,isweep]=ray1.h.vel_east navy_array[7,iray,isweep]=ray1.h.vel_north navy_array[8,iray,isweep]=ray1.h.vel_up navy_array[9,iray,isweep]=ray1.h.rvc FOR igate=0,ngates-1 DO BEGIN FOR ivolume=0,nvolumes-1 DO $ data_array[ivolume,igate,iray,isweep]=radar.volume[ivolume].sweep[isweep].ray[iray].range[igate] ENDFOR ENDFOR ENDFOR ;Meters to the first gate and gate_range PRINT,'first_gate=',radar.volume[VRvolume].sweep[0].ray[0].h.range_bin1,'(Meters)',$ ' gate_size=',radar.volume[VRvolume].sweep[0].ray[0].h.gate_size,'(Meters)' first_gate=radar.volume[VRvolume].sweep[0].ray[0].h.range_bin1 gate_size=radar.volume[VRvolume].sweep[0].ray[0].h.gate_size ;Assume that we are extracting the data of the first 10 Km maxrange=10000 ;in Meters ngates=(maxrange-radar.volume[VRvolume].sweep[0].ray[0].h.range_bin1)/radar.volume[VRvolume].sweep[0].ray[0].h.gate_size PRINT,'maxrange=',maxrange,'(Meters) ngates for output=',ngates ;----------Set up variable arrays vr=FLTARR(nrays,ngates) az=FINDGEN(nrays)*360/nrays vrr=FLTARR(nrays) vr_cor=FLTARR(nrays) ; get radial velocity correction vr_cor=navy_array[9,*] elevation=radar.volume[0].sweep[0].ray[0].h.elev elevation=elevation[0:nsweeps-1]+1.0 ;elevation[0]=1.0 ; get Vr in the sweep nearest 10 degree elevation vol_Vr=rsl_get_volume(radar,'vr') Vr_sweep=RSL_GET_SWEEP(radar.volume[VRvolume],elevation(0)) ;--reinitialze height-----to eliminate carryover to next file height[*]=0 ;-------Open the output file and print the header---------------------- PrintF,iu,'(MONTH DAY YEAR HOUR MINUTE SECOND): ' PrintF,iu,radar.h.month,radar.h.day,radar.h.year,radar.h.hour,radar.h.minute,radar.h.sec PrintF,iu,'RADAR POSITION: ' PrintF,iu,radar.h.latd+radar.h.latm/60.+radar.h.lats/3600.,$ radar.h.lond+radar.h.lonm/60.+radar.h.lons/3600. Maxheight=MAX(height) print,'elevation angle=',elevation PrintF,iu,'elev_deg','height','speed','theta','corrUx','corrUy',$ 'error','total_used',format='(8a15)' FOR k=7,ngates-1 DO BEGIN ;-----------Note first gate is smallest k---------------------- FOR j=0,nrays-1 DO BEGIN FOR j=0,nrays-1 DO BEGIN VR_ray=(RSL_GET_RAY_FROM_SWEEP(Vr_sweep,j)).range ; note j is the azimuth angle in degrees, ray=rsl_get_ray(vol_Vr,1.,j) header=ray.h az(j)=header.azimuth*0.017453293 vr(j,k)=VR_ray[k] +vr_cor(j) ; note k is the range gate number ENDFOR ENDFOR ;-------set up vrr array for call to ransac2--------- vrr=vr[*,k] Height(k)=first_gate+gate_size*(k+1)*SIN(0.017453*elevation(0)) ;-----------clean up the fit with ransac2------------------ el=elevation[0] ht=Height(k) ; -=============================--Clean up data call - Ransac4 ==================================================== ; RANSAC4 ;======================================================================================== ct(k)=count ; number of points used error(k)=err IF flag EQ 1 THEN BEGIN Ux(k)=a1/COS(0.017453*elevation[0])+navy_array[6,0] Uy(k)=b1/COS(0.017453*elevation[0])+navy_array[7,0] speed(k)=SQRT(ux(k)^2.0+uy(k)^2.0) theta(k)=180/!PI*ATAN(ux(k)/uy(k)) IF theta(k) LT 0 THEN theta(k)=theta(k)+360 Height(k)=first_gate+gate_size*(k+1)*SIN(0.017453*elevation(0)) ;;---------print the values to the output file-------------- PrintF,iu,elevation(0),Height(k),speed(k),theta(k),Ux(k),$ Uy(k),error(k),ct(k),format='(8f15.5)' print,'kkk=',kkk,'k=',k,'height=',height[k] ENDIF ENDFOR CLOSE,1 ENDFOR END