PRO GRID_MODIS

COMPILE_OPT IDL2

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

;- 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]
  print, input_file
  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, *]
  help, lat, lon, data

  ;- Get list of grid indices
  lonlat_to_index, lon, lat, resolution, grid_index

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

  ;- Loop over grid cells
  for cell_index = 0, n_elements(grid) - 1 do begin

    ;- 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
      print, cell_index, mean
    endif
   
  endfor

  ;- Display the grid
  map_set, 30, 130, /orthographic, /isotropic, title=input_file
  image = map_image(grid, x0, y0, compress=1)
  loadct, 39, /silent
  tv, bytscl(image, min=0.0, max=5.0), x0, y0
  map_continents

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

endfor

END