grid_modis.pro

PRO GRID_MODIS

COMPILE_OPT IDL2

;- Obtener la lista de todos los archivos de entrada
list = file_search('G:/AOD MODIS/AQUA/out/MYD04_L2.*.dat')
help, list

;- Crear una grilla de igual ángulo
;- (la esquina inferior izquierda es a 180W, 90S; en efecto, -180.0 y -90.0)
resolution = 1.0
ncol = long(360.0 / resolution)
nrow = long(180.0 / resolution)
grid = fltarr(ncol, nrow)

;- Bucle sobre los archivos de entrada
for index = 0, n_elements(list) - 1 do begin

  ;- Abrir cada archivo de entrada y obtener un numero de filas
  input_file = list[index]
  print, input_file
  openr, lun, input_file, /get_lun
  result = fstat(lun)
  nrows = result.size / (4 * 3)

  ;- Leer la tabla de datos
  table = fltarr(3, nrows)
  readu, lun, table
  free_lun, lun

  ;- Extraer, lat, lon y data
  lat = table[0, *]
  lon = table[1, *]
  data = table[2, *]
  help, lat, lon, data

  ;- Obtener la lista de indices de grillado
  lonlat_to_index, lon, lat, resolution, grid_index

  ;- Resetear la matriz de grillado
  grid[*] = -999.9

  ;- Bucle sobre las celdas de grillado
  for cell_index = 0, n_elements(grid) - 1 do begin

    ;- Encontrar los data points dentro de esta celda de grillado
    loc = where(grid_index eq cell_index, count)

    ;- Guardar el promedio para esta celda
    if (count gt 0) then begin
      mean = total(data[loc]) / count
      grid[cell_index] = mean
      print, cell_index, mean
    endif

  endfor

  ;- Mostrar el grillado
  map_set, -10, -90, /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

  ;- Guargar el grillado de salida
  date_name = strmid(basename(input_file), 10, 12)
  output_file =  'G:/AOD MODIS/AQUA/grid/' + date_name + '.grd'
  print, 'Escribiendo ' + output_file
  openw, lun, output_file, /get_lun
  writeu, lun, grid
  free_lun, lun

endfor

END