PRO GRID_MODIS

COMPILE_OPT IDL2

;- Get list of input files
list = file_search('a1.*.dat')

;- Create equal angle grid
;- (lower left corner of grid is at 180W, 90S, i.e. -180.0, -90.0)
resolution = 1.0
ncol = long(360.0 / resolution)
nrow = long(180.0 / resolution)
grid = fltarr(ncol, nrow)

;- Loop over input files
for index = 0, n_elements(list) - 1 do begin

  ;- Open input file and get number of rows
  input_file = list[index]
  openr, lun, input_file, /get_lun
  result = fstat(lun)
  nrows = result.size / (4 * 3)

  ;- Read data table
  table = fltarr(3, nrows)
  readu, lun, table
  free_lun, lun

  ;- Extract, lat, lon, data values
  lat = table[0, *]
  lon = table[1, *]
  data = table[2, *]

  ;- Get a list of grid indices for every lat/lon in the input data
  lonlat_to_index, lon, lat, resolution, grid_index

  ;- Get a list of unique grid indices in the input data
  sort_idx = sort(grid_index)
  uniq_idx = uniq(grid_index[sort_idx])
  unique_grid_index = grid_index[uniq_idx] 

  ;- Reset the grid array
  grid[*] = -999.9

  ;- Loop over the grid cells which are present in the data
  for i = 0, n_elements(unique_grid_index) - 1 do begin

    ;- Get this cell index
    cell_index = unique_grid_index[i]
        
    ;- Find the data points within this grid cell
    loc = where(grid_index eq cell_index, count)

    ;- Save the mean for this cell
    if (count gt 0) then begin
      mean = total(data[loc]) / count
      grid[cell_index] = mean
    endif
   
  endfor

  ;- Save the output grid
  output_file = strmid(input_file, 0, 8) + '.grd'
  openw, lun, output_file, /get_lun
  writeu, lun, grid
  free_lun, lun
  print, 'Created ', output_file
  
endfor

END