; ; This program processes a one dimensional variable and writes it to the output ; netcdf file. Determines the number of data points, number of good data ; points,mode, minimum, and maximum value for each time interval ; pro one_dimensional_mode2,dataname common ncdf_file, oldcdfname,twoD,time_dim_name,height_dim_name,$ height_var_name,avgfilename,curdate common old_grid, num_times,num_heights,base_time,time_offset,$ height,site,files,num_files,num_days common new_grid, numtimes, newtimes, numheights, newheights,$ timespacing, heightspacing common bracket, startbst common data, avgdata ; ; Loop through all the matching files in the directory ; j=-1 ;counts number of files containing the variable for i=0,num_files-1 do begin ; ; Don't try to process any gif files that may be in the directory ; if strpos(files[i],'gif') eq -1 and $ strpos(files[i],'png') eq -1 then begin print,i,j,files[i] ; ; Open the netcdf file ; cdf_id=ncdf_open(files[i]) ; ; Get the variable id ; data_id=ncdf_varid(cdf_id,dataname) if data_id ne -1 then begin print,'Processing 1-D ',dataname j=j+1 ; ; Get the variable ; ncdf_varget,cdf_id,data_id,data ; ; If it is the first file ; if j eq 0 then begin combined_data=data datainfo=ncdf_varinq(cdf_id, data_id) ;inquire about the variable numatts=datainfo.natts ;number of attributes assigned to the variable parts=strsplit(files[i],'/',/extract) part=strsplit(parts[5],'.',/extract) origfile=part[0]+'.'+part[1] ; ; Store the attribute names ; dataatr=strarr(numatts) ;string array to hold attribute names for k=0,numatts-1 do begin ; Returns the name of the attribute dataatr[k]=ncdf_attname(cdf_id, data_id, k) endfor ; ; It is not the first file ; endif else combined_data=[combined_data,data] endif ;end of variable exists in file ; ; Close the netcdf file ; ncdf_close,cdf_id endif ;end of not a gif file endfor ;end of looping through the matching files ; ; Continue processing data because some data exists ; if j gt -1 then begin data=combined_data ; ; Remove bad values ; result=where(data ge 1e18, count) if count ne 0 then data[result]=-9999 ; ; Convert cloudtop and cloudbase from msl to agl ; if dataname eq 'RadarFirstTop' or dataname eq 'CloudBaseBestEstimate' then begin if site eq 'sgp' then begin result=where(data ne -9999, count) if count ne 0 then data[result]=data[result]-318. endif if site eq 'nsa' then begin result=where(data ne -9999, count) if count ne 0 then data[result]=data[result]-7. endif if site eq 'twpc1' then begin result=where(data ne -9999, count) if count ne 0 then data[result]=data[result]-3. endif if site eq 'twpc2' then begin result=where(data ne -9999, count) if count ne 0 then data[result]=data[result]-6. endif endif ; ; Give all bad or missing data vales the same bad flag of -9999 ; if dataname eq 'CloudBaseBestEstimate' or dataname eq 'RadarFirstTop' then begin result=where(data le 0, count) if count ne 0 then data[result]=-9999 endif if dataname eq 'liq' then begin result=where(data le 0, count) if count ne 0 then data[result]=-9999 endif ; ; Find the mode value ; possible_points=fltarr(numtimes) ;all possible data points within the time interval, both useable and nonuseable data values data_points=fltarr(numtimes) ;all valid data points within the time interval, only useable data values avg_data=fltarr(numtimes) ;the mode for the time interval hi_ext=fltarr(numtimes) ;the highest data value in the interval lo_ext=fltarr(numtimes) ;the lowest data value in the interval for t=0,numtimes-2 do begin result=where(time_offset ge newtimes[t] and time_offset lt newtimes[t+1], count_all) ;find all possible data within time interval possible_points[t]=count_all ;all possible points within the time interval result1=where(time_offset ge newtimes[t] and time_offset lt newtimes[t+1] and data ne -9999, count1) ;find good data in the time interval data_points[t]=count1 ;possible points with a valid data point if count1 ne 0 then begin sub_data=data[result1] ;good data hist=histogram(sub_data,omax=hi_extreme,omin=lo_extreme) bins=findgen(n_elements(hist))+min(sub_data) hist_max=max(hist,ind) if hist_max eq 1 and count1 gt 1 then avg_data[t]=-9999 else $ avg_data[t]=bins(ind) ;mode value hi_ext[t]=hi_extreme lo_ext[t]=lo_extreme endif else begin avg_data[t]=-9999 hi_ext[t]=-9999 lo_ext[t]=-9999 ;print,avg_data[t],hi_ext[t],lo_ext[0] endelse endfor ; ; Calculate the average for the last data point ; result_all=where(time_offset ge newtimes[numtimes-1] and time_offset lt (newtimes[numtimes-1]+300), count_all) possible_points[numtimes-1]=count_all result=where(time_offset ge newtimes[numtimes-1] and time_offset lt (newtimes[numtimes-1]+300) and data ne -9999, count) data_points[numtimes-1]=count if count ne 0 then begin sub_data=data[result] ;good data hist=histogram(sub_data,omax=hi_extreme,omin=lo_extreme) bins=findgen(n_elements(hist))+min(sub_data) hist_max=max(hist,ind) if hist_max eq 1 and count gt 1 then avg_data[numtimes-1]=-9999 else avg_data[numtimes-1]=bins(ind) ;mode value hi_ext[numtimes-1]=hi_extreme lo_ext[numtimes-1]=lo_extreme endif else begin avg_data[numtimes-1]=-9999 hi_ext[numtimes-1]=-9999 lo_ext[numtimes-1]=-9999 endelse ; ; Plot to check data ; plot1='no' ;flag this yes if you want to see plots to the screen, no if you don't if plot1 eq 'yes' then begin ; ; Find the minimum data value ; result=where(lo_ext eq -9999, count) if count ne 0 then lo_ext[result]=1e8 dminv=min(lo_ext) result=where(lo_ext eq 1e8, count) if count ne 0 then lo_ext[result]=-9999 ; ; Find the data maximum ; dmaxv=max(hi_ext) loadct,5 window,2,title='avg '+dataname plot,newtimes,avg_data,$ psym=6,$;color=255,$ yrange=[dminv,dmaxv] oplot,newtimes,hi_ext,psym=4,color=355 oplot,newtimes,lo_ext,psym=5,color=255 window,3,title='raw '+dataname plot,time_offset,data,$ psym=6,$ yrange=[dminv,dmaxv],$ xrange=[newtimes[0],newtimes[numtimes-1]] endif ;plots to the screen ; ; Open the netcdf file ; avg_id=ncdf_open(avgfilename, /write) ; ; Get the variable dimension ; time_id=ncdf_dimid(avg_id,'time') ; ; Put the netcdf file into define mode ; ncdf_control, avg_id, /redef ; ; Define variable ; possible_points_id=ncdf_vardef(avg_id,dataname+'PossiblePoints',time_id,/short) data_points_id=ncdf_vardef(avg_id,dataname+'UseablePoints',time_id,/short) avg_data_id=ncdf_vardef(avg_id, dataname+'Mode', time_id) hi_ext_id=ncdf_vardef(avg_id, dataname+'MaxExtreme', time_id) lo_ext_id=ncdf_vardef(avg_id, dataname+'MinExtreme', time_id) ; ; Go back to an old file and get the information ; for a variable. Necessary for attribute copying ; for i=0,num_files-1 do begin if strpos(files[i],'gif') eq -1 and $ strpos(files[i],'png') eq -1 then begin cdf_id=ncdf_open(files[i]) data_id=ncdf_varid(cdf_id,dataname) if data_id ne -1 then goto, found_var ncdf_close,cdf_id endif endfor found_var: ;skip out of loop because found a file containing the variable ; ; Define attributes ; for i=0,numatts-1 do begin nowatt=dataatr[i] mytry=ncdf_attcopy(cdf_id, data_id, nowatt, avg_id, avg_data_id) mytry=ncdf_attcopy(cdf_id, data_id, nowatt, avg_id, hi_ext_id) mytry=ncdf_attcopy(cdf_id, data_id, nowatt, avg_id, lo_ext_id) endfor ncdf_close,cdf_id ncdf_attput, avg_id, avg_data_id, 'original data source file ',origfile ncdf_attput, avg_id, possible_points_id, 'long_name','number of possible data points within the time interval' ncdf_attput, avg_id, data_points_id, 'long_name','number of useable points withing the time interval' ; ; Change from define mode to data mode ; ncdf_control, avg_id, /endef ; ; Write variable ; ncdf_varput, avg_id, avg_data_id, avg_data ncdf_varput, avg_id, hi_ext_id, hi_ext ncdf_varput, avg_id, lo_ext_id, lo_ext ncdf_varput, avg_id, possible_points_id, possible_points ncdf_varput, avg_id, data_points_id, data_points ; ; Close the netcdf file ; ncdf_close, avg_id endif ;variable does exist in at least one file return end