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