Visualization_Using_Ozone_Data_of_HDF5_Format.pro

pro Visualization_Using_Aerosol_Data_of_HDF5_Format

    cd, 'C:\SYSTEM\PROG\IDL\HDF5'

    dim = [1934, 1544]
    
    hdf4_file='INPUT/coms_cn_geos_lonlat.hdf'
    hdf5_file = "INPUT/coms_mi_le2_aod_cn_201402010000.h5"
    
    file_id = h5f_open(hdf5_file)
    data_id = h5d_open(file_id, '/Product/Aerosol_Optical_Depth')
    data = h5d_read(data_id)
    h5f_close, file_id
    data = float(data)
  
    valid_min = 0.000000
    valid_max = 5.00000
    fill = 32767
    offset = 0.000000
    scale = 0.00488759
    
    loc1 = where( ((data lt valid_min[0]) or (data gt valid_max[0]) or (data eq fill[0])), nloc1, complement=loc2, ncomplement=nloc2)
    if (nloc1 gt 0) then data[loc1] = !values.f_nan
    if (nloc2 gt 0) then data[loc2] = scale[0] * (data[loc2] - offset[0])
    val2D = data
  
    sd_id=hdf_sd_start(hdf4_file, /read)
    sds_index = hdf_sd_nametoindex(sd_id, 'Lat')
    sds_id = hdf_sd_select(sd_id, sds_index)
    hdf_sd_getdata, sds_id, data
    data = float(data)
  
    loc1 = where( (data eq -999.0), nloc1, complement=loc2, ncomplement=nloc2)
    if (nloc1 gt 0) then data[loc1] = !values.f_nan
    if (nloc2 gt 0) then data[loc2] = data[loc2]
    lat2D = data
  
    sd_id=hdf_sd_start(hdf4_file, /read)
    sds_index = hdf_sd_nametoindex(sd_id, 'Lon')
    sds_id = hdf_sd_select(sd_id, sds_index)
    hdf_sd_getdata, sds_id, data
    data = float(data)
  
    loc1 = where( (data eq -999.0), nloc1, complement=loc2, ncomplement=nloc2)
    if (nloc1 gt 0) then data[loc1] = !values.f_nan
    if (nloc2 gt 0) then data[loc2] = data[loc2]
    lon2D = data
    
    hdf_sd_endaccess, sds_id
    hdf_sd_end, sd_id
  
    dims = size(var, /dimensions)

    print, min(lon2D, /nan), mean(lon2D, /nan), max(lon2D, /nan)
    print, min(lat2D, /nan), mean(lat2D, /nan), max(lat2D, /nan)
    print, min(val2D, /nan), mean(val2D, /nan), max(val2D, /nan)
    
    help, lon2D, lat2D, val2D
  
  
    ROI = [-20, 60, 79.9, 200] 
    cnLat = 0
    cnLon = 128.2
  
    dtDateYmd = STRMID(hdf5_file, 25, 12)
    psName = "Img01_" + dtDateYmd
    mainName = "HDF5 Aerosol Optical Depth : " + dtDateYmd

    fnMakePsPlot, lon2D, lat2D, val2D, dim, cnLat, cnLon, ROI, psName, mainName, 0, 2, 0.1, 33

end

;====================================================================
; Subroutine : Function Make Postscript to Png Graph
;====================================================================
pro fnMakePsPlot, lon2D, lat2D, val2D, dim, cnLat, cnLon, ROI, psName, mainName, zmin, zmax, interval, color_table

  set_plot, "ps"

  device, filename = psName + ".ps", decomposed = 0, bits = 8, /color, xsize = 16, ysize = 11, /inches, font_size = 13, /Helvetica

  !p.font = 0 & !p.charsize = 2.0 & !p.charthick = 1.6  & !p.multi = [0,1,1] & !p.background = 255

  xvert = [0.1, 0.1, -0.1, -0.1 ,0.1]
  yvert = [-0.1, 0.1, 0.1, -0.1, -0.1]

  usersym, xvert, yvert, /fill

  start_color = 0 & end_color = 255 & colorn = end_color - start_color + 1

  cgloadct, color_table

  latmin = ROI(0) & latmax= ROI(1) & lonmin = ROI(2) & lonmax= ROI(3)

  cgMAP_SET, cnLat, cnLon, /CYLINDRICAL $
    , limit = [latmin, lonmin, latmax, lonmax] $
    , position = [0.05, 0.0, 0.85, 0.90], /noborder, /ISOTROPIC

  for i = 0L, dim[0] - 1 do begin
    for j = 0L, dim[1] - 1 do begin
      
      nLon = lon2D[i, j]
      nLat = lat2D[i, j]
      nVal = val2D[i, j]

      if (lonmin gt nLon or nLon gt lonmax or finite(nLon, /nan) eq 1) then continue
      if (latmin gt nLat or nLat gt latmax or finite(nLat, /nan) eq 1) then continue
      if (finite(nVal, /nan) eq 1) then continue

      plots, nLon, nLat, psym = 8, symsize = 2, color = BYTSCL(nVal, zmin, zmax)

    endfor
  endfor

  cgColorbar, NColors = 255, BOTTOM = 1 $
    , Divisions = 5, COLOR = "black" $
    , Position = [0.1, 0.90, 0.8, 0.94] $
    , Range = [zmin, zmax] $
    , CHARSIZE = 1.5,CHARthick = 2 $
    , VERTICAL = 1, RIGHT = 1

;  cgLoadct, 2

;  contour, val2D, lon, lat, /overplot, C_THICK = 2 $
;    , LEVELS = zmin + findgen(10) * interval $
;    , C_LABELS = zmin + findgen(10) * interval $
;    , C_CHARSIZE = 1.5, C_CHARTHICK = 2 $
;    , C_COLORS = BYTSCL(zmin + findgen(10) * interval, zmin, zmax)

  cgMAP_CONTINENTS, /coast, /countries, COLOR = "black"

  lats = [-90:90:15]
  lons = [-180:360:20]

  lats_names = strarr(n_elements(lats))
  lons_names = strarr(n_elements(lons))

  for i = 0L, n_elements(lats) - 1 do begin
    if (lats[i] gt 0) then begin
      lats_names[i] = textoidl(string(lats[i]) + "\circN")
    endif else if (lats[i] eq 0) then begin
      lats_names[i] = textoidl(string(lats[i]) + "\circ")
    endif else begin
      lats_names[i] = textoidl(string(-lats[i]) + "\circS")
    endelse
  endfor

  for i = 0L,n_elements(lons) - 1  do begin
    if (lons[i] gt 0) then begin
      lons_names[i] = textoidl(string(lons[i]) + "\circE")
    endif else if (lons[i] eq 0) then begin
      lons_names[i] = textoidl(string(lons[i]) + "\circ")
    endif else begin
      lons_names[i] = textoidl(string(-lons[i]) + "\circW")
    endelse
  endfor

  cgMap_grid, color = "black", charsize = 1.8, thick = 2.5, lats = lats,  latnames = lats_names, lons = lons, lonnames = lons_names $
             , linestyle = 1, bthick = 300, /box_axes; , /no_grid

  cgText, 0.45, 0.95, mainName, /normal, charsize = 2, CHARTHICK = 2, color = "black", alignment = 0.5

  device, /close_file

  com = "convert -flatten -background white " + psName + ".ps" + " " + "./FIG/" + file_basename(psName, ".ps") + ".png"
  spawn, /hide, com

  file_delete, ps_name + ".ps"

  return

end