EASY READ/EXTRACT HDF DATA IN IDL

The Example IDL Program below, along with the 
associated IDL functions from Ed Santiago, 
show how one can simply and quickly extract 
data from an HDF file(s).

The routine below only requires the names of 
your HDF file(s), Vdata and an (optional) 
list of desired fields.

The single line call that does all the work 
for you is:
   s = read_asc_hdf( filepath, vdataname, lstofflds )

The above line may be incorporated into any IDL 
routine of your own design.  Just remember to 
allow your IDL session to compile the associated 
routines:
   read_asc_hdf.pro
   get_fieldinfo.pro
   get_fieldstruct.pro
   split.pro

The source for these associated routines is
included at the bottom of this file after the 
Example IDL Program.
Or, you can download the code from
http://www.srl.caltech.edu/ACE/ASC/DATA/level2/asc_idl.zip


;EXAMPLE IDL PROGRAM
   
; $Id$
; Author: Glenn Hamell, Ed Santiago 
; NAME: get_HDF_data.pro 
; Copyright (c) 2003, California Institute of Technology. All rights
; reserved. Unauthorized reproduction prohibited.
;
;
; CALLING SEQUENCE:  get_HDF_data
;
; PURPOSE:
;	Quickly & easily access HDF data
;
; INPUT:
;	The HDF file(s) named in the "filepath" variable below.
;
; OUTPUT:
;	A structure of the desired data fields. Optional plots or whatever
;	code modifications you would like to add.	
;
; RETURN:
;	
;
; MODIFICATION HISTORY: 
; ===================== 
; 2003Feb21-Glenn Hamell, Created. Utilizes routines by Ed Santiago.
;	 
;
;################################################
;

PRO get_HDF_data

filepath  = '/home/ull9/glennh/apps/asc/mag_data_1hr_1997.hdf'
vdataname = 'MAG_data_1hr'
lstofflds = 'year,day,hr,fp_doy,Bgse_x,Bgse_y,Bgse_z' ; NOTE: no spaces allowed

s = read_asc_hdf( filepath, vdataname, lstofflds ); GET ONLY FIELDS NAMED 
;s = read_asc_hdf( filepath, vdataname )	  ; GET ALL FIELDS IN vdataname

HELP, s, /st
;PRINT, s

PLOT, s.fp_doy, s.Bgse_x, MIN_VALUE=-999

END
;------------------------------------------------;



;ASSOCIATED IDL 

; = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
;+
; NAME:
;     READ_ASC_HDF
;
; IDENT:
;     $Id: read_asc_hdf.pro,v 1.1 2003/02/14 20:19:48 glennh Exp $
;
; PURPOSE:
;     read in a range of ASC (hdfgen'ed) HDF files
;
; AUTHOR:
;      Ed Santiago
;
; CALLING SEQUENCE:
;     x = READ_ASC_HDF(filelist, basename, fields)
;
; INPUTS:
;     filelist     - array of input files
;     basename     - "hdfgen" name of vdata/sdata, eg, "swepam_dswi"
;     fields       - list of fields to read; if empty all fields are read
;
; OUTPUTS:
;     extracted data is in an array of structures
;
; SIDE EFFECTS:
;
; EXAMPLE:
;
;-
FUNCTION read_asc_hdf, files, basename, fields, Quiet=quiet, Debug=debug

  IF N_Elements(quiet) EQ 0 THEN quiet = 0
  IF N_Elements(debug) EQ 0 THEN debug = 0

  ; Error conditions returned when things don't work out right with one HDF
  ERR_NO_SUCH_VD = 'IDL_M_HDF_VS_SET_FIELDS'
  ERR_NO_SUCH_SD = 'IDL_M_HDF_SD_SELECT'

  nfiles = N_Elements(files)

  FOR i=0,nfiles-1 DO BEGIN
      IF NOT quiet THEN BEGIN
          status = string(format='("-> ",A)', files[i])
          IF nfiles GT 1 THEN $
               status = status + string(format='(" (",I0," of ",I0,")")',$
					i+1, nfiles)
          print, status
      ENDIF

      ; Initialize the HDF thingies
      fid   = HDF_Open(files[i])
      IF fid LE 0 THEN BEGIN
          MESSAGE, 'HDF_Open(' + files[i] + ') failed!', /Continue
          RETURN, 0
      ENDIF

      v     = HDF_VD_Find(fid,basename)
      IF v LE 0 THEN BEGIN
        MESSAGE, 'HDF_Find('+basename+') in '+files[i]+' failed!', /Continue
        HDF_Close, fid
        RETURN, 0
      ENDIF

      vdata = HDF_VD_Attach(fid, v)
      IF vdata LE 0 THEN BEGIN
        HDF_Close, fid
        MESSAGE, 'HDF_VD_Attach('+files[i]+') failed!'
      ENDIF

      ; First time around, figure out what the data types are
      IF N_Elements(fieldinfo) EQ 0 THEN BEGIN
          fieldinfo   = get_fieldinfo(files[i], basename, vdata, fields)
	  fieldstruct = get_fieldstruct(fieldinfo)
      ENDIF

      ; read the file contents, and extract fields
      HDF_VD_Get, vdata, count=nr
      IF nr LE 0 THEN GOTO, skipThisFile

      today = REPLICATE(fieldstruct, nr)

      ; Loop over each field, reading it in using the appropriate HDF interface
      FOR j=0,N_Elements(fieldinfo)-1 DO BEGIN
          ; If field is a multimensional array, get it from the SD interface
          IF fieldinfo[j].ndims GT 1 THEN BEGIN
              sd = -1
	      Catch, err_sts
	      IF err_sts EQ 0 THEN BEGIN
                  sd  = HDF_SD_Start(files[i])
	          sdi = HDF_SD_NameToIndex(sd, basename+'_'+fieldinfo[j].name)
	          sds = HDF_SD_Select(sd, sdi)
	          IF debug THEN print,'Reading SD: '+fieldinfo[j].name+'..'
	          HDF_SD_GetData, sds, data
                  IF debug THEN print,'Copying...'
                  today.(j+1) = temporary(data)
	          IF debug THEN print,'done!'
                  HDF_SD_EndAccess, sds
                  HDF_SD_End, sd
	          sd = -1
	      ENDIF ELSE IF !Error_State.Name EQ ERR_NO_SUCH_SD THEN BEGIN
	          print,'['+files[i]+': no such SD: ' + fieldinfo[j].name + ']'
	      ENDIF ELSE BEGIN
help,!error_state,/st
	          Catch, /Cancel
	          print, 'hdf_sd_getdata(): err=',err_sts
		  HDF_Close, fid
	          Message, !ERR_STRING
	      ENDELSE
	      Catch, /Cancel
	      IF sd GE 0 THEN HDF_SD_End, sd
          ENDIF ELSE BEGIN
	      ; Field is scalar or vector; read it using the VDATA interface
	      Catch, err_sts
	      IF err_sts EQ 0 THEN BEGIN
                  HDF_VD_Seek, vdata, 0
                  status = HDF_VD_Read(vdata, data, fields=fieldinfo[j].name)
	          IF status NE nr THEN MESSAGE, 'mismatch'
                  today.(j+1) = temporary(data)
	      ENDIF ELSE IF !Error_State.Name EQ ERR_NO_SUCH_VD THEN BEGIN
	          print,'['+files[i]+': no such VD: ' + fieldinfo[j].name + ']'
	      ENDIF ELSE BEGIN
	          Catch, /Cancel
	          print, 'hdf_vd_read(): err=',err_sts
	          Message, !ERR_STRING
	      ENDELSE
	      Catch, /Cancel
          ENDELSE
      ENDFOR

      ; Done reading all fields.  Make a list of all structure elements.
      IF N_Elements(matchlist) EQ 0 THEN BEGIN
          matchlist = today
      ENDIF ELSE BEGIN
          matchlist = [ matchlist, today ]
      ENDELSE

      ; done.
skipThisFile:
      HDF_VD_Detach, vdata
      HDF_Close, fid
  ENDFOR

  IF N_Elements(matchlist) EQ 0 THEN BEGIN
      IF NOT quiet THEN MESSAGE, 'Could not find any data in this range',/Cont
      return, 0
  ENDIF

  RETURN, matchlist
END

; = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
;+
; NAME:
;     GET_FIELDINFO
;
; IDENT:
;     $Id: get_fieldinfo_hdf.pro,v 1.1 2003/02/14 20:17:45 glennh Exp $
;
; PURPOSE:
;     determine the dimensions and types of HDF file field thingies
;
; AUTHOR:
;     Ed Santiago
;
; CALLING SEQUENCE:
;     x = GET_FIELDINFO( filename, basename, vdata, fields )
;
; INPUTS:
;     filename     - name of the HDF input file
;     basename     - name of the ASC thingy (eg, "swepam_i")
;     vdata        - HDF vdata thingy
;     fields       - comma-separated string of fields to be extracted
;
; OUTPUTS:
;     an array of structures
;
; SIDE EFFECTS:
;
; EXAMPLE:
;      fid   = HDF_Open(files[i])
;      v     = HDF_VD_Find(fid,'swepam_i')
;      vdata = HDF_VD_Attach(fid, v)
;      fields = 'foo,bar'
;
;      fieldinfo   = get_fieldinfo(files[i], 'swepam_i', vdata, fields)
;
;-
FUNCTION get_fieldinfo, filename, basename, vdata, fields
  ; define the "fieldinfo" structure, containing field name/type/dimensions
  fieldinfo = { fieldinfo, name:'', type:'', ndims:0L, dims:LonArr(8) }
  fieldlist = [ fieldinfo ]

  ; Yuk.  IDL is case-insensitive, but HDFs are not.  The fields are
  ; stored in case-sensitive manner, so if we are called with "bmag"
  ; and the field is actually called "Bmag", we will fail.
  ;
  ; To avoid this, we always read in all the VD and SD field names,
  ; then do a case-insensitive comparison of the fields passed in.
  ; This allows us to find out the *real* HDF name of the field.

  ; obtain the VDATA fields
  HDF_VD_Get, vdata, fields=fields_all

  ; obtain the SDATA fields
  sd = HDF_SD_Start(filename)
  HDF_SD_FileInfo, sd, datasets, attributes
  FOR i=0,datasets-1 DO BEGIN
      sds = HDF_SD_Select(sd,i)
      HDF_SD_GetInfo, sds, name=nnn
      IF STRMID(nnn,0,STRLEN(basename)+1) EQ basename + '_' THEN BEGIN
          fields_all = fields_all + ',' + STRMID(nnn, STRLEN(basename)+1, 999)
      ENDIF
      HDF_SD_EndAccess, sds
  ENDFOR
  HDF_SD_End, sd


  ; If no fields were given, assume user wants them all
  IF N_Elements(fields) EQ 0 THEN fields = fields_all

  ;
  ; Split the list of fields into an array.  For each element, do some 
  ; HDF info-get stuff to figure out what it is.
  ;
  fieldarr     = split(',', fields)
  fieldarr_all = split(',', fields_all)
  FOR i = 0, N_Elements(fieldarr)-1 DO BEGIN
      ; Find out the real (HDF case-sensitive) name of the field
      tmp = where(StrUpCase(fieldarr_all) EQ StrUpCase(fieldarr[i]), c)
      IF c NE 1 THEN BEGIN
        MESSAGE, 'no such VD or SD: ' + fieldarr[i], /Continue
        GOTO, skipThisField
      ENDIF
      fieldinfo.name  = fieldarr_all[tmp]

      ; HDF routines choke on errors, but we want to proceed.
      Catch, err_sts

      IF err_sts EQ 0 THEN BEGIN
          ;
          ; First, see if it's a VDATA.  If it is, life is simple.  If
          ; it isn't, the following line will choke, and control will 
          ; pass on to the ELSE below.
          ;
          HDF_VD_GetInfo, vdata, fieldinfo.name, Type=t, Order=o
          IF o EQ 1 THEN fieldinfo.ndims = 0 ELSE fieldinfo.ndims = 1
          fieldinfo.dims[0] = o
      ENDIF ELSE IF !Error_State.Name EQ 'IDL_M_HDF_VS_FEXIST' THEN BEGIN
	  ;
	  ; Sigh, the HDF_VD_GetInfo call failed.  That means that the
	  ; field in question is not a VDATA.  Now try to see if it's 
	  ; an SDATA.  Since this HDF routine chokes also, we need a
	  ; new CATCH.
	  ;
          Catch, err_sts
          IF err_sts EQ 0 THEN BEGIN
	      ;
	      ; Initialize the SDATA interface, and try to get info.  If
	      ; it works, we set the dimensions and type.
	      ;
              sd  = HDF_SD_Start(filename)
              sdi = HDF_SD_NameToIndex(sd, basename + '_' + fieldinfo.name)
              sds = HDF_SD_Select(sd, sdi)

              HDF_SD_GetInfo, sds, Type=t, Dims=o
              fieldinfo.ndims = N_Elements(o) - 1
              fieldinfo.dims  = o
	      HDF_SD_EndAccess, sds
	      HDF_SD_End, sd
          ENDIF ELSE IF !Error_State.Name EQ 'IDL_M_HDF_SD_SELECT' THEN BEGIN
	      ; Nope, could not find a VDATA or SDATA by that name.
              print, '['+filename+': no such VD or SD: ' + fieldarr[i] + ']'
          ENDIF ELSE BEGIN
	      ; Weird unexpected error
              print, 'hdf_sd failed, err=',err_sts
              stop
          ENDELSE
      ENDIF ELSE BEGIN
	  ; Weird unexpected error
	  Catch, /Cancel
          print, 'hdf_vd_info failed with status ', err_sts, fieldarr[i]
	  Message, !ERR_STRING, /NoName
      ENDELSE
      Catch, /Cancel

      IF err_sts EQ 0 THEN BEGIN
          fieldinfo.type  = t
          fieldlist = [ fieldlist, fieldinfo ]
      ENDIF
skipThisField:
  ENDFOR	; i=0,nfields

  RETURN, fieldlist[1:*]
END

; = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
;+
; NAME:
;     GET_FIELDSTRUCT
;
; IDENT:
;     $Id: get_fieldstruct_hdf.pro,v 1.1 2003/02/14 20:18:57 glennh Exp $
;
; PURPOSE:
;     Given a fieldinfo struct (see GET_FIELDINFO.PRO), build a data
;     structure of the right type.
;
; AUTHOR:
;     Ed Santiago
;
; CALLING SEQUENCE:
;     data = GET_FIELDSTRUCT( fieldinfo )
;
; INPUTS:
;     fieldinfo    - data structure returned by GET_FIELDINFO()
;
; OUTPUTS:
;     a structure template
;
; SIDE EFFECTS:
;
; EXAMPLE:
;
;-
FUNCTION get_fieldstruct, fieldinfo
  mystruct = create_struct( 'dnum', 0D )

  FOR i=0, N_Elements(fieldinfo)-1 DO BEGIN
    type  = fieldinfo[i].type		; String describing the type
    ndims = fieldinfo[i].ndims		; Number of array dims (eg, 1, 2, 3)
    dims  = fieldinfo[i].dims		; The dimensions themselves (eg, [9x3])

    ; Sigh.  Rather than give us something *USEFUL*, such as its internal
    ; data type number (float=4,int=2,etc), IDL gives us a string.  It's
    ; up to us to convert that to a real data type.  We do so in a gross
    ; CASE statement, converting each string descriptor to a variable
    ; of the correct IDL type.
    ;
    ; To make things worse, with the introduction of unsigned types
    ; in 5.2, IDL has added the 'UINT' and 'ULONG' strings.  Piece
    ; o' cake to add those, along with their respective definitions
    ; as "0U" and "0UL", right?  Well, no.  IDL <5.2 gives a compile-
    ; time error with those, so we can't do that if we want to remain
    ; compatible with 5.1.  The solution is to define the type as a
    ; string, then EXECUTE() it.  This way we don't get the compilation
    ; errors, nor will EXECUTE() fail, since UINT/ULONG don't exist in <5.2.
    CASE type OF 
      'BYTE':   t = '0B'
      'INT':    t = '0'
      'UINT':	t = '0U'
      'LONG':   t = '0L'
      'ULONG':  t = '0UL'
      'FLOAT':  t = '0.0'
      'DOUBLE': t = '0.0D'
    ENDCASE

    r = execute('xxx = ' + t)

    ; If field is an array, make it so.
    IF ndims NE 0 THEN $
        xxx = Make_Array(size=[ndims, dims[0:ndims-1], (size(xxx))[1], 0])

    ; Append this new field to the Master Data Structure
    mystruct = create_struct( mystruct, fieldinfo[i].name, xxx)
  ENDFOR

  RETURN, mystruct
END

; = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
;+
; NAME:
;     SPLIT
;
; IDENT:
;     $Id: split_hdf.pro,v 1.1 2003/02/14 20:20:07 glennh Exp $
;
; PURPOSE:
;     Just like perl.  Splits a string into an array of strings.
;
; AUTHOR:
;     Ed Santiago
;
; CALLING SEQUENCE:
;     stringarr = split(delimiter, string)
;
; INPUTS:
;     delimiter  -- character to split on
;     string     -- string to split
;
; OUTPUTS:
;     an array of strings
;
; SIDE EFFECTS:
;
; EXAMPLE:
;     IDL> x=split(',', 'this,is,a,test')
;     IDL> print, N_Elements(x), x
;           4 this is a test
;
;-


FUNCTION split, delimiter, string	; Just like perl
  arr     = [ 'x' ] ; sigh

  len     = strlen(string)
  lastpos = 0
  WHILE lastpos LT len DO BEGIN
      pos     = STRPOS(string, delimiter, lastpos)
      IF pos EQ -1 THEN pos = len
      arr     = [ arr, STRMID(string, lastpos, pos-lastpos) ]

      ; Collapse multiple spaces into one
      IF delimiter EQ ' ' THEN WHILE StrMid(string,pos+1,1) EQ ' ' DO pos=pos+1

      lastpos = pos+1
  ENDWHILE

  ; Always guaranteed at least one hit, unless string is null
  RETURN, arr[1:*]
END