plot_modis.pro

PRO PLOT_MODIS, LAT, LON

COMPILE_OPT IDL2

;- Crear angulos grillados iguales
;- (la esquina inferior izquierda de la grilla esta a 180W, 90S, en efecto, -180.0, -90.0
resolution = 1.0
ncol = long(360.0 / resolution)
nrow = long(180.0 / resolution)
grid = fltarr(ncol, nrow)

;- Obtener la lista de archivos grillados
list = file_search('G:/AOD MODIS/AQUA/grid/*.grd')
help, list

;- Crear matrices de salida
data = fltarr(n_elements(list))
time = fltarr(n_elements(list))
hour = intarr(n_elements(list))
date = strarr(n_elements(list))

;- Calcular los indices de la celda requerida
lonlat_to_index, lon, lat, resolution, index
help, index

;- Bucle sobre los archivos grillados
for i = 0, n_elements(list) - 1 do begin

  ;- Leer un archivo
  input_file = list[i]
  openr, lun, input_file, /get_lun
  readu, lun, grid
  free_lun, lun

  ;- Obtener data de la celda deseada
  data[i] = grid[index]

  ;- Obtener el dia y año del nombre del archivo
  year = long(strmid(basename(input_file), 0, 4))
  day_of_year = long(strmid(basename(input_file), 4, 3))

  ;- Obtener el mes y el dia del mes
  ydn2md, year, day_of_year, month, day

  ;- Obtener dia Juliano
  time[i] = julday(month, day, year)
  hour[i] = strmid(basename(input_file),8,4)
  date [i] = string(day)+'-'+string(month)+'-'+string(year)
endfor

help, data
help, time

;- Establecer los valores faltantes a NaN (Not a Number)
loc = where(data lt 0.0, count)
if (count gt 0) then data[loc] = !VALUES.F_NAN

;- Calcular la serie de tiempo, suavizada sobre 16 dias
smoothed_data =smooth(data, 16, /nan, /edge_truncate)
print, smoothed_data
;- Cargar la tabla de escala de grises
loadct, 0

;- Plotear la data
result = label_date(date_format=['%M %Y'])
plot, time, smoothed_data, $
  xrange=[min(time), max(time)], $
  xstyle=1, $
  psym=1, $
  color=0, background=255, $
  xtickunits='months', $
  xtickformat='label_date', $
  ytitle='Aerosol Optical Depth', $
  title='Serie de tiempo MODIS en' + string(lat, lon, format='(2f8.2)')

;- Guardar la data a plotear
output_file = "G:/AOD MODIS/AQUA/" + string(fix(lat)) + "_" + string(fix(lon)) + ".csv"
;openw, lun, output_file, /get_lun
write_csv, output_file, date,hour,smoothed_data

END