;
; This example code illustrates how to access and visualize GES DISC
; MLS v4 HDF-EOS5 Swath file in IDL. 
;
; If you have any questions, suggestions, comments  on this example, 
; please use the HDF-EOS Forum (http://hdfeos.org/forums).
;
; If you would like to see an  example of any other NASA HDF/HDF-EOS 
; data product that is not listed in the HDF-EOS Comprehensive Examples page
; (http://hdfeos.org/zoo), feel free to contact us at eoshelp@hdfgroup.org
; or post it at the HDF-EOS Forum (http://hdfeos.org/forums).

; Usage:
;
;  $idl MLS-Aura_L2GP-H2O_v04-20-c01_2013d003.he5.idl
;
; Tested under: IDL 8.7.2
; Last updated: 2020-11-09


; Open file.
file_name = 'MLS-Aura_L2GP-H2O_v04-20-c01_2013d003.he5'
file_id = H5F_OPEN(file_name)

; Read H2O dataset.
datafield_name = '/HDFEOS/SWATHS/H2O/Data Fields/L2gpValue'
data_id = H5D_OPEN(file_id,datafield_name)
data = H5D_READ(data_id)

; Get title attribute.
attr_id = H5A_OPEN_NAME(data_id, 'Title')
long_name = H5A_READ(attr_id)
H5A_Close, attr_id

; Get units attribute.
units_id = H5A_OPEN_NAME(data_id, 'Units')
units = H5A_READ(units_id)
H5A_Close, units_id

; Get missing value attribute.
missingvalue_id = H5A_OPEN_NAME(data_id,'MissingValue')
missingvalue = H5A_READ(missingvalue_id)
H5A_Close, missingvalue_id
H5D_Close, data_id

; Open pressure dataset.
pressure_name = '/HDFEOS/SWATHS/H2O/Geolocation Fields/Pressure'
pressure_id = H5D_OPEN(file_id, pressure_name)
pressure = H5D_READ(pressure_id)

; Get units
units_id = H5A_OPEN_NAME(pressure_id, 'Units')
units_p = H5A_READ(units_id)
H5A_Close, units_id
H5D_Close, pressure_id

; Read latitude.
lat_name = '/HDFEOS/SWATHS/H2O/Geolocation Fields/Latitude'
lat_id = H5D_OPEN(file_id, lat_name)
lat = H5D_READ(lat_id)
H5D_Close, lat_id

; Read longitude.
lon_name = '/HDFEOS/SWATHS/H2O/Geolocation Fields/Longitude'
lon_id = H5D_OPEN(file_id, lon_name)
lon = H5D_READ(lon_id)  
H5D_Close, lon_id

; Convert 2D data to 1D.
data = data(0, *)

; Convert data type.
dataf = float(data)
missingvaluef = float(missingvalue(0))

; Process missing value.
idx = where(dataf eq missingvaluef(0), cnt)
if cnt gt 0 then dataf[idx] = !Values.F_NAN

; Generate the plot.
m = MAP('Geographic', TITLE=file_name, FONT_SIZE=9, /BUFFER)
ct = COLORTABLE(72, /reverse)

; Create title label.
t_label = long_name + ' at Pressure = ' + STRING(pressure(0)) + ' ' + units_p
t1 = TEXT(0.35, 0.2, t_label)

; Get max and min value of data for color bar.
datamin = MIN(dataf, /NAN)
datamax = MAX(dataf, /NAN)

; We use SCATTERPLOT because data is 2-d lat/lon swath.
; lon[*]/lat[*]/dataf[*] will make 1-d dataset.
c1 = SCATTERPLOT(lon, lat, OVERPLOT=m, $
                 MAGNITUDE = dataf, $
                 RGB_TABLE=ct, $
                 POSITION=[0.1, 0.1, 0.83, 0.9],$
                 /SYM_FILLED, SYMBOL='o', SYM_SIZE=0.1)
mc = MAPCONTINENTS()

; We need custom colorbar because we use SCATTERPLOT().
; We cannot use TARGET=c1.
cb = COLORBAR(RGB_TABLE=ct, RANGE=[datamin, datamax], /BORDER, ORIENTATION=1, TEXTPOS=1, POSITION=[0.85,0.2,0.87,0.8], TITLE=units)

png = file_name + '.idl.png'
c1.save, png, HEIGHT=600, WIDTH=800
EXIT