;
; This program processes a two dimensional variable and writes it to the output
; netcdf file.  Determines the nearest neighbor value for stratus retrieval.
;

pro two_dimensional_nearest_neighbor,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, avg_data

;
;  Loop through all the matching files in the directory
;
m=-1  ;counts number of cdf files without the variable
p=-1  ;counts number of cdf files with variable
j=-1  ;counts number of cdf files
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  
    j=j+1
    print,i,j,files[i]
;
;  Open the netcdf file
;
    cdf_id=ncdf_open(files[i])
;
;  Get the dimension id's
;
    time_did=ncdf_dimid(cdf_id,time_dim_name)
    height_did=ncdf_dimid(cdf_id,height_dim_name)
;
;  Get the dimensions
;
    ncdf_diminq,cdf_id,time_did,charstring,ntimes
    ncdf_diminq,cdf_id,height_did,charstring,nheights
;
;  Create array for nodata value
;
    dummy=make_array(ntimes,nheights,/float,value=-8888)
;
;  Get the variable id
;
    data_id=ncdf_varid(cdf_id,dataname)
;
;  Get the variable if it exists
;
    if data_id ne -1 then begin
      print,'Processing 2-D ',dataname
      ncdf_varget, cdf_id, data_id, data
      data=transpose(data)
      p=p+1
    endif else begin
      print, 'Variable not in this file'
      m=m+1
    endelse
;
;  If it is the first real variable then get the attributes
;
    if p eq 0 and data_id ne -1 then begin
      print, 'first variable, getting the attributes'
        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
    endif
;
;  Combine the data
;
    if p eq 0 and m lt 0 and data_id ne -1 then begin
      print,'first file and variable exists' 
      combined_data=data
    endif else if p ge 0 and data_id ne -1 then begin
      print,'adding on existing variable'
      combined_data=[combined_data,data]
    endif else if p lt 0 and m eq 0 and data_id eq -1 then begin
      print,'first file and variable doesnt exist'
      combined_data=dummy
    endif else if m ge 0 and data_id eq -1 then begin
      print,'adding on dummy variable'
      combined_data=[combined_data,dummy]
    endif
;
;  Close the netcdf file
;
    ncdf_close,cdf_id
  endif  ;end of not a gif file
endfor  ;end of looping though the matching files

;
;  Continue processing data because some data exists
;

if p gt -1 then begin 

;
;  Manipulate the big data array
;

data=combined_data
combined_data=0  ;release memory

;********************************
;  Nearest neighbor interpolation
;********************************

;
;  Create the array of closest time subscripts
;

close_time_sub=fltarr(numtimes)
for t=0,numtimes-1 do begin
  timespaces=ABS(time_offset-newtimes[t])
  closest=min(timespaces,closest_sub)
  close_time_sub[t]=closest_sub
endfor

;
;  Create the array of closest height subscripts
;

close_height_sub=fltarr(numheights)
for h=0,numheights-1 do begin
  heightspaces=abs(height-newheights[h])
  closest=min(heightspaces,closest_sub)
  close_height_sub[h]=closest_sub
endfor

;
;  Create the new data array
;

avg_data=make_array(numtimes,numheights, /float, value=-9999)

for h=0,numheights-1 do begin
for t=0,numtimes-1 do begin

  ;if nearest data point is more than 1 min away, then give -8888 data value
  ;if nearest data point is more than 50 meters away, then give -8888 data value
  deltat=abs(time_offset[close_time_sub[t]]-newtimes[t])
  deltah=abs(height[close_height_sub[h]]-newheights[h])
  if (deltat gt 60) or (deltah gt 50) then begin
    avg_data[t,h]=-8888 
  endif else begin
    avg_data[t,h]=data[close_time_sub[t],close_height_sub[h]]
  endelse
 
endfor
endfor

;**************************************
;  Add variable to netcdf file
;**************************************

;
;  Open the netcdf file
;

avg_id=ncdf_open(avgfilename, /write)

;
;  Get the variable dimensions
;

time_id=ncdf_dimid(avg_id,'time')
height2_did=ncdf_dimid(avg_id,'height')
datadim=intarr(2)
datadim[1]=time_id
datadim[0]=height2_did

;
;  Put the netcdf file into define mode
;

ncdf_control, avg_id, /redef

;
;  Define variable
;

avg_data_id=ncdf_vardef(avg_id, dataname, datadim)

;
;  Define attributes
;

;
;  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
for i=0,numatts-1 do begin
  nowatt=dataatr[i]
  mytry=ncdf_attcopy(cdf_id, data_id, nowatt, avg_id, avg_data_id)
endfor
ncdf_close,cdf_id  ;close the source file that attributes were copied from

ncdf_attput, avg_id, avg_data_id, 'missing ','-8888'
ncdf_attput, avg_id, avg_data_id, 'original data source file ',origfile

;
;  Change from define mode to data mode
;

ncdf_control, avg_id, /endef

;
;  Write variable
;

avg_data=transpose(avg_data)  ;transpose for writing to netcdf
ncdf_varput, avg_id, avg_data_id, avg_data
avg_data=transpose(avg_data)  ;retrun to original form for plotting

;
;  Close the netcdf file
;

ncdf_close, avg_id

;
;  Plot to check data
;

scatterp='no'  ;flag this yes if you want plots to the screen, no if you don't
if scatterp eq 'yes' then begin

;
;  Scatter plot to check the data
;

if dataname eq 'ReflectivityBestEstimate' then begin
	miv=-5000
	mav=max(avg_data)
endif else begin
	mav=max(avg_data)
	miv=0
endelse
window,2,title='avg '+dataname
surface,avg_data,newtimes,newheights,min_value=miv,max_value=mav,/horizontal

window,3,title='raw '+dataname

n1=0
while time_offset[n1] lt newtimes[0] do  n1=n1+1
n2=n1
while time_offset[n2] lt newtimes[numtimes-1] do begin
	n2=n2+1
        if n2 eq num_times then begin
		n2=n2-1
		goto, enough
	endif
endwhile
enough:
sub_data=data[n1:n2,*]
sub_time_offset=time_offset[n1:n2]
zoomsize=size(data)
if zoomsize[1] ge 2 then surface,sub_data,sub_time_offset,height,min_value=miv,max_value=mav,/horizontal

endif  ;scatter plots

;
;  image plots to further check the data
;

imagep='no'  ;flag this yes if you want image plots to the screen, no if you don't
if imagep eq 'yes' then begin

;
;  Find the original data range and set the null values to 1e6
;

dmax=max(sub_data)  ;max data value
dmin=min(sub_data)  ;null value
result=where(sub_data eq dmin,count)
if count ne 0 then sub_data[result]=1e6
dmin=min(sub_data)  ;min data value
;print, 'data min', dmin, 'data max', dmax

image=bytscl(sub_data,max=dmax,min=dmin)

window,3,xsize=600,ysize=600
loadct,5
contour, sub_data, sub_time_offset, height, /nodata,$
	xstyle=1,ystyle=1
	px=!x.window*!d.x_vsize
	py=!y.window*!d.y_vsize
	sx=px[1]-px[0]+1
	sy=py[1]-py[0]+1
tv,congrid(image,sx,sy),px[0],py[0]

;
;  Regridded data range is commented out so to be scaled
;  according to the original data range.
;

;dmin=min(newdata)  ;null value
;result=where(newdata eq dmin,count)
;if count ne 0 then newdata[result]=1e6
;dmin=min(newdata)  ;min data value
;print, 'newdata min', dmin, 'newdata max', dmax

;
;  Scale the data vales and insert then into image
;

image=bytscl(newdata,max=dmax,min=dmin)

;
;  Plot the new regridded data
;

window,4,xsize=600,ysize=600
loadct,5
contour, newdata, newtimes, newheights, /nodata,$
	xstyle=1,ystyle=1
px=!x.window*!d.x_vsize
py=!y.window*!d.y_vsize
sx=px[1]-px[0]+1
sy=py[1]-py[0]+1
tv, congrid(newimage,sx,sy),px[0],py[0]
endif  ;image plots

endif  ;variable does exist in at least one file

return
end