;
;  This example code illustrates how to access and visualize PO.DAAC
;  SMAP L2B 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 SMAP_L2B_SSS_32221_20210211T162812_R17000_V5.0.h5.idl
;
; Tested under: IDL 8.7.2
; Last updated: 2021-02-16

; Open file.
file_name='SMAP_L2B_SSS_32221_20210211T162812_R17000_V5.0.h5'
file_id=H5F_OPEN(file_name)


; Retrieve data.
datafield_name='/smap_spd'
data_id=H5D_OPEN(file_id, datafield_name)
dataf=H5D_READ(data_id)

; Get '_FillValue' attribute.
fv_id=H5A_OPEN_NAME(data_id, '_FillValue')
fv=H5A_READ(fv_id)
H5A_CLOSE, fv_id

; Get 'valid_min' attribute.
vmin_id=H5A_OPEN_NAME(data_id, 'valid_min')
valid_min=H5A_READ(vmin_id)
H5A_CLOSE, vmin_id

; Get 'valid_max' attribute.
vmax_id=H5A_OPEN_NAME(data_id, 'valid_max')
valid_max=H5A_READ(vmax_id)
H5A_CLOSE, vmax_id

; Get 'long_name' attribute.
title_id=H5A_OPEN_NAME(data_id, 'long_name')
long_name=H5A_READ(title_id)
H5A_CLOSE, title_id

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

; Close dataset.
H5D_CLOSE, data_id

; Retrieve lat/lon.
lat_name='/lat'
lat_id=H5D_OPEN(file_id, lat_name)
lat=H5D_READ(lat_id)

; Get '_FillValue' attribute.
lat_fv_id=H5A_OPEN_NAME(data_id, '_FillValue')
lat_fv=H5A_READ(lat_fv_id)
H5A_CLOSE, lat_fv_id

H5D_CLOSE, lat_id

lon_name='/lon'
lon_id=H5D_OPEN(file_id, lon_name)
lon=H5D_READ(lon_id)

; Get '_FillValue' attribute.
lon_fv_id=H5A_OPEN_NAME(data_id, '_FillValue')
lon_fv=H5A_READ(lon_fv_id)
H5A_CLOSE, lon_fv_id

H5D_CLOSE, lon_id

; Close file.
H5F_CLOSE, file_id

; Process missing value, convert dataf that are equal to missingvaluef to NaN
idx = where(dataf eq fv(0), cnt)
if cnt gt 0 then lat[idx] = !Values.F_NAN
if cnt gt 0 then lon[idx] = !Values.F_NAN
if cnt gt 0 then dataf[idx] = !Values.F_NAN

idx = where(lat eq lat_fv(0), cnt)
if cnt gt 0 then dataf[idx] = !Values.F_NAN
if cnt gt 0 then lat[idx] = !Values.F_NAN
if cnt gt 0 then lon[idx] = !Values.F_NAN

idx = where(lon eq lon_fv(0), cnt)
if cnt gt 0 then lat[idx] = !Values.F_NAN
if cnt gt 0 then lon[idx] = !Values.F_NAN
if cnt gt 0 then dataf[idx] = !Values.F_NAN

; Process valid_range values.
idx=where(dataf LT valid_min(0) OR dataf GT valid_max(0), cnt)
if cnt gt 0 then lat[idx] = !Values.F_NAN
if cnt gt 0 then lon[idx] = !Values.F_NAN
if cnt gt 0 then dataf[idx] = !Values.F_NAN

; Pick only valid points. Otherwise, scatter plot will be wrong.
gidx = WHERE(FINITE(dataf), count)
lats = lat[gidx]
lons = lon[gidx]
datas = dataf[gidx]
datamin = MIN(datas)
datamax = MAX(datas)
m = MAP('Geographic', TITLE=file_name, FONT_SIZE=9, /BUFFER)
ct = COLORTABLE(72, /reverse)
t1 = TEXT(0.35, 0.2, long_name)
c1 = SCATTERPLOT(lons[*], lats[*], OVERPLOT=m, $
                 MAGNITUDE=datas[*], $
                 RGB_TABLE=ct, $
                 POSITION=[0.1, 0.1, 0.83, 0.9], $
                 /SYM_FILLED, SYMBOL='o', SYM_SIZE=0.1)
mc = MAPCONTINENTS()
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