;+ ; NAME: ; ReadEarthFile ; PURPOSE: ; Read ea* file(s) ; CALLING SEQUENCE: FUNCTION ReadEarthFile,file_in,t1,f1,tt_shift=tt_shift, $ tt_range=tt_range, tt_nfile=tt_nfile, $ data_filter=data_filter, data_path=data_path ; INPUTS: ; file scalar; type: string ; (if an array is specified only the first element is used) ; name of input file (may include wildcard) ; OPTIONAL INPUT PARAMETERS: ; These are used if a single record is to be extracted from a group ; of ea* files (see PROCEDURE): ; ; tt_shift=tt_shift array[1]; time difference structure; default: 0 hours ; offset from forecast time in days (see PROCEDURE) ; internally tt_shift is rounded to the nearest quarter day. ; tt_range=tt_range array[1] or array[2]; type: time structure or float ; if array[1] only ea* files prior to tt_range are used ; if array[2] only ea* files between tt_range[0] and ; tt_range[1] will be used ; tt_range can be specified as a time structure or as a ; modified Carrington variable. ; tt_nfile=tt_nfile scalar; type: integer ; max. # files from which a single record is read ; only the most recent tt_nfile files are used. ; OUTPUTS: ; status scalar; type: byte ; 0B: error; t1 and f1 should not be used ; 1B: file(s) succesfully read. ; t1 array[n]; type: time structure ; times read from file ; f1 array[n,4]; type: float ; velocity, density, radial magnetic field, ; tangential magnetic field at times t1 ; INCLUDE: @compile_opt.pro ; On error, return to caller ; CALLS: ; SetFileSpec, GetFileSpec, sfloat, flt_read, tagarray ; TimeSet, TimeGet, TimeUnit, TimeString, Carrington, BadValue, TimeOp ; PROCEDURE: ; The tomographic data at the location of earth are stored in files with ; name ea_1234.678, where 1234.678 indicates a modified Carrington variable. ; ; Each file contains 161 records. The 'forecast time' is the middle record ; (record 80). Records are separated by time steps of 6 hours (0.25 days). ; ; If none of the keywords tt_range, tt_shift or tt_nfile is specified then a ; single ea* file is read into arrays t1 and f1. ; ; If either one of the keywords is specified then from each ea* file the ; record corresponding to the forecast time, plus tt_shift is extracted. ; If tt_range is specified then only files prior to tt_range[0] are used ; (if tt_range has only one element) or all files between tt_range[0] and ; tt_range[1] are used (if tt_range has two elements). ; ; After tt_range has been applied a group of ea* files remain. If tt_nfile is ; specified then only the most recent tt_nfile ea* files are accessed. ; If tt_nfile is not specified then all ea* files are accessed. ; MODIFICATION HISTORY: ; NOV-1999, Paul Hick (UCSD) ; OCT-2000, Paul Hick (UCSD); added option to read single record from a group of ; ea* files, using keywords tt_shift and tt_nfile. ;- single_rec = n_elements(tt_range) ne 0 or n_elements(tt_shift) ne 0 or n_elements(tt_nfile) ne 0 if n_elements(data_path) eq 0 then $ data_path = filepath(root=getenv('NAGOYA'),subdir='slow','raw') if n_elements(data_filter) eq 0 then $ data_filter = 'ea*.*' bad = [-9999.999,-9999.999, -999.999, -999.999, -999.999, -999.999] ; Never modify the input argument file_in !!!!! if single_rec then begin ; Read single record ; Set up the list of files to be accessed. if n_elements(file_in) eq 0 then begin path = data_path filter = data_filter endif else begin SetFileSpec, file_in[0] path = GetFileSpec(upto='dir' ) if path eq '' then path = data_path filter = GetFileSpec(from='name') if filter eq '' then filter = data_filter endelse file = file_search( filepath(root=path, filter), count=count ) if file[0] eq '' then begin ; No ea* files found message, /info, 'no '+filter+' file(s) found in '+path return, 0 endif file = file[ sort(file) ] ; Make sure the files are in chronological order ntt = n_elements(tt_range) if ntt ne 0 then begin SetFileSpec, file name = GetFileSpec(from='name') tmp = sfloat(name[0], fmt=format, /xfora) times = dblarr(n_elements(name)) reads, name, format=format, times times = Carrington(times) trange = Carrington(tt_range,/get_time) ; Convert to time structure case ntt of 1: begin tmp = TimeOp(/subtract, times, trange[0], sign=sign) sign = where(sign lt 0, count) if count eq 0 then begin message, /info, 'no '+filter+' file(s) found in '+path+ $ ' before '+TimeString(trange[0], /ymd) return, 0 endif end 2: begin tmp = TimeOp(/subtract,times, trange[0], sign=sign1) tmp = TimeOp(/subtract,times, trange[1], sign=sign2) sign = where(sign1 gt 0 and sign2 lt 0, count) if count eq 0 then begin message, /info, 'no '+filter+' file(s) found in '+path+ $ ' between '+TimeString(trange[0], /ymd)+' and '+ $ TimeString(trange[1], /ymd) return, 0 endif end endcase file = file[sign] endif if n_elements(tt_nfile) ne 0 then begin ; If # files specified if count gt tt_nfile then begin file = file[count-tt_nfile:count-1]; Keep only last tt_nfile files count = tt_nfile endif endif endif else begin ; Read a single ea* file if n_elements(file_in) eq 0 then begin path = data_path filter = data_filter endif else begin SetFileSpec, file_in[0] ; Only directory specified if GetFileSpec(from='name',upto='name') eq '' then begin path = file_in[0] filter = data_filter endif else begin file = file_search(file_in[0], count=tmp) case tmp of 0: begin message, /info, 'no file(s) found: '+file_in[0] return, 0 end 1: file = file[0] else: begin ; More than one file, probably wildcard path = GetFileSpec(upto='dir' ) filter = GetFileSpec(from='name') ; Filter contains wildcard end endcase endelse endelse if n_elements(path) ne 0 then begin file = (dialog_pickfile(path=path,filter=filter))[0] if file eq '' then begin message, /info, 'no '+filter+' file selected from '+path return, 0 endif endif endelse ; Check the first file for the column values for density, velocity, etc. status = flt_read(file[0], in, comments=comments, /silent, delay=1) if not status or n_elements(comments) lt 2 then return, 0 header = tagarray(strlowcase(comments[1]),' ',/nonull) columns = [where(header eq 'velocity'), where(header eq 'density'), where(header eq 'brad' ), where(header eq 'btang')] if columns[0] eq -1 or columns[1] eq -1 then return, 0 ntime = (size(in))[2] ; # times in each ea* file ; If the column 'rot' (Carrington rotation) is present the file was created with ; the ipsg2 program, which sets the density at Earth to one. unit_density = (where(header eq 'rot'))[0] ne -1 ; Older files may not have any magnetic field data present if columns[2] eq -1 or columns[3] eq -1 then begin message, /info, 'no B data: '+file[0] columns = columns[0:1] ; Pick up velocity and density only endif ncolumn = n_elements(columns) ; # columns with valid data if single_rec then begin if n_elements(tt_shift) eq 0 then tt_shift = TimeSet( /diff, hour=0) ; Extract record itime. (ntime-1)/2 is the forecast time (the middle record in the file). ; tt_shift is used to determine an offset from the middle record assuming that the time step is 6 hour. itime = (ntime-1)/2+round( TimeGet(/diff, tt_shift, TimeUnit(/day))/0.25) itime = itime[0] if itime lt 0 or itime ge ntime then begin message, /info, 'specified time offset not inside time range' return, 0 endif f1 = replicate(1,count)#bad ; Initialize with bad values everywhere t1 = replicate(!TheTime, n_elements(f1[*,0])) for ifile=0,count-1 do begin ; Loop over all files if not flt_read(file[ifile], in, /silent, delay=1) then return, 0 ; Extract data from record itime t1[ifile ] = TimeSet(yr=in[0,itime], doy=in[1,itime], hour=in[2,itime]) f1[ifile,0:ncolumn-1] = in[columns,itime] endfor message, /info, 'processed'+strcompress(count)+' '+filter+' files' endif else begin t1 = reform(TimeSet(yr=in[0,*], doy=in[1,*], hour=in[2,*])) f1 = replicate(1,ntime)#bad ; Initialize with bad values everywhere f1[*,0:ncolumn-1] = transpose(in[columns,*]) endelse ; Scale the densities to a nominal density at Earth of 5. if unit_density then f1[*,0] = 5*f1[*,0] ; Replace bad data values by !values.f_nan. for i=0,n_elements(bad)-1 do begin tmp = where(f1[*,i] eq bad[i]) if tmp[0] ne -1 then f1[tmp,i] = BadValue(bad[i]) endfor if n_elements(file_in) eq 0 then file_in = file return, 1 & end