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