function READFITS, filename, header, heap, NOSCALE = noscale, $
SILENT = silent, EXTEN_NO = exten_no, NUMROW = numrow, $
POINTLUN = pointlun, STARTROW = startrow, $
NaNvalue = NaNvalue, NSLICE = nslice
;+
; NAME:
; READFITS
; PURPOSE:
; Read a FITS file into IDL data and header variables.
; EXPLANATION:
; Under Unix, READFITS() can also read gzip or Unix compressed FITS files.
;
; CALLING SEQUENCE:
; Result = READFITS( Filename,[ Header, heap, /NOSCALE, EXTEN_NO=,
; NSLICE=, /SILENT , NaNVALUE =, STARTROW =, NUMROW = ] )
;
; INPUTS:
; FILENAME = Scalar string containing the name of the FITS file
; (including extension) to be read. If the filename has
; a *.gz extension, it will be treated as a gzip compressed
; file. If it has a .Z extension, it will be treated as a
; Unix compressed file.
;
; OUTPUTS:
; Result = FITS data array constructed from designated record.
; If the specified file was not found, then Result = -1
;
; OPTIONAL OUTPUT:
; Header = String array containing the header from the FITS file.
; heap = For extensions, the optional heap area following the main
; data array (e.g. for variable length binary extensions).
;
; OPTIONAL INPUT KEYWORDS:
;
; EXTEN_NO - scalar integer specify the FITS extension to read. For
; example, specify EXTEN = 1 or /EXTEN to read the first
; FITS extension. Extensions are read using recursive
; calls to READFITS.
;
; NaNVALUE - This scalar is only needed on architectures (such as VMS)
; that do not recognize the IEEE "not a number" (NaN) convention.
; It specifies the value to translate any IEEE "not a number"
; values in the FITS data array. In addition, if the data is
; stored as integer (BITPIX = 16 or 32), and BSCALE is present,
; then NaNValue gives the values to pixels assigned with the
; BLANK keyword.
;
; NOSCALE - If present and non-zero, then the ouput data will not be
; scaled using the optional BSCALE and BZERO keywords in the
; FITS header. Default is to scale.
;
; NSLICE - An integer scalar specifying which N-1 dimensional slice of a
; N-dimensional array to read. For example, if the primary
; image of a file 'wfpc.fits' contains a 800 x 800 x 4 array,
; then
;
; IDL> im = readfits('wfpc.fits',h, nslice=2)
; is equivalent to
; IDL> im = readfits('wfpc.fits',h)
; IDL> im = im(*,*,2)
; but the use of the NSLICE keyword is much more efficient.
;
; NUMROW - This keyword only applies when reading a FITS extension.
; If specifies the number of rows (scalar integer) of the
; extension table to read. Useful when one does not want to
; read the entire table.
;
; POINT_LUN - Position (in bytes) in the FITS file at which to start
; reading. Useful if READFITS is called by another procedure
; which needs to directly read a FITS extension. Should
; always be a multiple of 2880.
;
; SILENT - Normally, READFITS will display the size the array at the
; terminal. The SILENT keyword will suppress this
;
; STARTROW - This keyword only applies when reading a FITS extension
; It specifies the row (scalar integer) of the extension table at
; which to begin reading. Useful when one does not want to read
; the entire table.
;
; EXAMPLE:
; Read a FITS file TEST.FITS into an IDL image array, IM and FITS
; header array, H. Do not scale the data with BSCALE and BZERO.
;
; IDL> im = READFITS( 'TEST.FITS', h, /NOSCALE)
;
; If the file contain a FITS extension, it could be read with
;
; IDL> tab = READFITS( 'TEST.FITS', htab, /EXTEN )
;
; The function TBGET() can be used for further processing of a binary
; table, and FTGET() for an ASCII table.
; To read only rows 100-149 of the FITS extension,
;
; IDL> tab = READFITS( 'TEST.FITS', htab, /EXTEN,
; STARTR=100, NUMR = 50 )
;
; To read in a file that has been compressed:
;
; IDL> tab = READFITS('test.fits.gz',h)
;
; ERROR HANDLING:
; If an error is encountered reading the FITS file, then
; (1) the system variable !ERROR is set (via the MESSAGE facility)
; (2) the error message is displayed (unless /SILENT is set),
; and the message is also stored in !ERR_STRING
; (3) READFITS returns with a value of -1
; RESTRICTIONS:
; (1) Cannot handle random group FITS
;
; NOTES:
; (1) If data is stored as integer (BITPIX = 16 or 32), and BSCALE
; and/or BZERO keywords are present, then the output array is scaled to
; floating point (unless /NOSCALE is present) using the values of BSCALE
; and BZERO. In the header, the values of BSCALE and BZERO are then
; reset to 1. and 0., while the original values are written into the
; new keywords O_BSCALE and O_BZERO. If the BLANK keyword was
; present, then any input integer values equal to BLANK in the input
; integer image are scaled to NaN (or the value of the NaNValue
; keyword) after the scaling to floating point.
;
; (2) The procedure FXREAD can be used as an alternative to READFITS.
; FXREAD has the option of reading an arbitary subsection of the
; primary FITS data.
;
; (3) The use of the NSLICE keyword is incompatible with the NUMROW
; or STARTROW keywords.
; PROCEDURES USED:
; Functions: SXPAR(), WHERENAN()
; Procedures: IEEE_TO_HOST, SXADDPAR, FDECOMP
;
; MODIFICATION HISTORY:
; MODIFIED, Wayne Landsman October, 1991
; Added call to TEMPORARY function to speed processing Feb-92
; Added STARTROW and NUMROW keywords for FITS tables Jul-92
; Work under "windows" R. Isaacman Jan-93
; Check for SIMPLE keyword in first 8 characters Feb-93
; Removed EOF function for DECNET access Aug-93
; Work under "alpha" Sep-93
; Null array processing fixed: quotes in a message
; properly nested, return added. Affected case when
; readfits called from another procedure. R.S.Hill Jul-94
; Correct size of variable length binary tables W.Landsman Dec-94
; To read in compressed files on Unix systems. J. Bloch Jan-95
; Check that file is a multiple of 2880 bytes Aug-95
; Added FINDFILE check for file existence K.Feggans Oct-95
; Consistent Error Handling W. Landsman Nov-95
; Handle gzip image extensions W. Landsman Apr-96
; Fixed bug reading 1-d data introduced Apr-96 W. Landsman Jun-96
; Don't use FINDFILE (too slow), & check for Blank values WBL Oct-96
; !VALUES wasn't compatible with IDL V3.6 WBL Jan-97
; Added ability to read Unix compressed (.Z) files WBL Jan-97
; Changed a FIX to LONG to handle very large tables WBL Apr-97
; Force use of /bin/sh shell with gzip WBL Apr-97
; Recognize BSCALE, BZERO in IMAGE extensions WBL Jun-97
; Added NSLICE keyword WBL Jul-97
; Added ability to read heap area after extensions WBL Aug-97
; Suppress *all* nonfatal messages with /SILENT WBL Dec-97
; Converted to IDL V5.0 WBL Dec-1997
; Fix NaN assignment for int data C. Gehman/JPL Mar-98
; Fix bug with NaNvalue = 0.0 C. Gehman/JPL Mar-98
;-
On_error,2 ;Return to user
; Check for filename input
if N_params() LT 1 then begin
print,'Syntax - im = READFITS( filename, [ h, heap, /NOSCALE, /SILENT,
print,' NaNValue = ,EXTEN_NO =, STARTROW = , NUMROW='
print,' NSLICE = ]
return, -1
endif
;
; Determine if file exists.
;
silent = keyword_set( SILENT )
openr, lun, filename, ERROR=error,/get_lun
if error EQ 0 then free_lun,lun else begin
message,/con,NoPrint=Silent,' ERROR - Unable to locate file ' + filename
return, -1
end
;
; Determine if the input file is compressed with gzip by the extension
;
fdecomp,filename,disk,dir,name,ext,ver
if (ext EQ "gz") or (ext EQ "Z") then gzip = 1 else gzip = 0
;
;
if not keyword_set( EXTEN_NO ) then exten_no = 0
; Open file and read header information
if gzip then begin
if (not silent) and (exten_no EQ 0) then $
message,"Input file compressed with gzip",/inform
spawn,"gzip -l "+filename,unit=unit
tmp=" "
readf,unit,tmp,format="(A80)"
compsize=0L
filesize=0L
dum = " "
origpcnt = " "
readf,unit,compsize,filesize,origpcnt
free_lun,unit
; Alternate means of obtaining filesize for 'compress'-ed files
if filesize eq -1 then begin
spawn,"zcat "+filename+" | wc -c ",unit=unit
readf,unit,filesize
free_lun,unit
fcompr = (filesize-compsize)/float(filesize)
origpcnt=string(fcompr*100,dir+name,format='(f6.1,"% ",a)')
endif
if (not silent) and (exten_no EQ 0) then begin
message,"Compressed size:"+string(compsize)+" bytes",/inform
message,"Uncompressed size:"+string(filesize)+" bytes",/inform
message,"% Compress/Original Filename: "+origpcnt,/inform
endif
spawn,"gzip -cd "+filename,unit=unit,/sh
file=fstat(unit)
endif else begin
openr, unit, filename, /GET_LUN, /BLOCK
file = fstat(unit)
flen = file.size/2880.
if (long(flen) NE flen) then $
message,'WARNING - File size of ' + strupcase(filename) + $
' is not a multiple of 2880 bytes',/CONT,NOPRINT=silent
endelse
if keyword_set( POINTLUN) then begin
if gzip then begin
tmp=bytarr(pointlun)
readu,unit,tmp
endif else begin
point_lun, unit, pointlun
endelse
endif else pointlun = 0
if gzip then nbytesleft = filesize - pointlun else $
nbytesleft = file.size - pointlun
hdr = bytarr( 80, 36, /NOZERO )
if nbytesleft LT 2880 then begin
free_lun, unit
message,/CON,NoPrint=Silent, $
'ERROR - EOF encountered while reading FITS header'
return, -1
endif
readu, unit, hdr
nbytesleft = nbytesleft - 2880
header = string( hdr > 32b )
if ( pointlun EQ 0 ) then $
if strmid( header[0], 0, 8) NE 'SIMPLE ' then begin
message,/CON,NoPrint=Silent, $
'ERROR - Header does not contain required SIMPLE keyword'
free_lun, unit
return, -1
endif
endline = where( strmid(header,0,8) EQ 'END ', Nend )
if Nend GT 0 then header = header[ 0:endline[0] ]
while Nend EQ 0 do begin
if nbytesleft LT 2880 then begin
message,/CON,NoPrint=Silent, $
'ERROR - EOF encountered while reading FITS header'
free_lun, unit
return, -1
endif
readu, unit, hdr
nbytesleft = nbytesleft - 2880
hdr1 = string( hdr > 32b )
endline = where( strmid(hdr1,0,8) EQ 'END ', Nend )
if Nend GT 0 then hdr1 = hdr1[ 0:endline[0] ]
header = [ header, hdr1 ]
endwhile
; Get parameter values
Naxis = sxpar( header, 'NAXIS' )
bitpix = sxpar( header, 'BITPIX' )
if !ERR EQ -1 then begin
message,/CON,NoPrint=Silent, $
'ERROR - FITS header missing required BITPIX keyword'
free_lun, unit
return, -1
endif
gcount = sxpar( header, 'GCOUNT') > 1
pcount = sxpar( header, 'PCOUNT')
case BITPIX of
8: IDL_type = 1 ; Byte
16: IDL_type = 2 ; Integer*2
32: IDL_type = 3 ; Integer*4
-32: IDL_type = 4 ; Real*4
-64: IDL_type = 5 ; Real*8
else: begin
message,/CON,NoPrint=Silent, $
'ERROR - Illegal value of BITPIX (= ' + $
strtrim(bitpix,2) + ') in FITS header'
free_lun,unit
return, -1
end
endcase
; Check for dummy extension header
if Naxis GT 0 then begin
Nax = sxpar( header, 'NAXIS*' ) ;Read NAXES
ndata = nax[0]
if naxis GT 1 then for i = 2, naxis do ndata = ndata*nax[i-1]
endif else ndata = 0
nbytes = (abs(bitpix)/8) * gcount * (pcount + ndata)
if pointlun EQ 0 then begin
extend = sxpar( header, 'EXTEND')
if !ERR EQ -1 then extend = 0
if not ( SILENT) then begin
if (exten_no GT 0) and (not EXTEND) then message,NoPrint=Silent, $
'ERROR - EXTEND keyword not found in primary header',/CON
endif
endif
if keyword_set( EXTEN_NO ) then begin
nrec = long(( nbytes +2879)/ 2880)
if gzip then pointlun = filesize - nbytesleft else $
point_lun, -unit, pointlun ;Current position
pointlun = pointlun + nrec*2880l ;Next FITS extension
free_lun, unit
im = READFITS( filename, header, heap, POINTLUN = pointlun, $
SILENT = silent, NUMROW = numrow, $
EXTEN = exten_no - 1, STARTROW = startrow )
return, im
endif
if nbytes EQ 0 then begin
if not SILENT then message, $
"FITS header has NAXIS or NAXISi = 0, no data array read",/CON
free_lun, unit
return,-1
endif
; Check for FITS extensions, GROUPS
groups = sxpar( header, 'GROUPS' )
if groups then MESSAGE,'WARNING - FITS file contains random GROUPS', /CON
; If an extension, did user specify row to start reading, or number of rows
; to read?
if not keyword_set(STARTROW) then startrow = 0
if naxis GE 2 then nrow = nax[1] else nrow = ndata
if not keyword_set(NUMROW) then numrow = nrow
if pointlun GT 0 then begin
xtension = strtrim( sxpar( header, 'XTENSION' , Count = N_ext),2)
if N_ext EQ 0 then message, /CON, NoPRINT = Silent, $
'ERROR - Header missing XTENSION keyword'
endif
if (pointlun GT 0) and ((startrow NE 0) or (numrow NE nrow)) then begin
nax[1] = nax[1] - startrow
nax[1] = nax[1] < numrow
sxaddpar, header, 'NAXIS2', nax[1]
if gzip then pointlun = filesize - nbytesleft else $
point_lun, -unit, pointlun ;Current position
pointlun = pointlun + startrow*nax[0] ;Next FITS extension
if gzip then begin
if startrow GT 0 then begin
tmp=bytarr(startrow*nax[0])
readu,unit,tmp
endif
endif else point_lun, unit, pointlun
endif else if keyword_set(NSLICE) then begin
lastdim = nax[naxis-1]
if nslice GE lastdim then message,/CON, NoPRINT = Silent, $
'ERROR - Value of NSLICE must be less than ' + strtrim(lastdim,2)
nax = nax[0:naxis-2]
sxdelpar,header,'NAXIS' + strtrim(naxis,2)
naxis = naxis-1
sxaddpar,header,'NAXIS',naxis
ndata = ndata/lastdim
if gzip then currpoint = filesize - nbytesleft else $
point_lun, -unit, currpoint ;Current position
currpoint = currpoint + nslice*ndata*abs(bitpix/8)
if gzip then begin
tmp = make_array( DIM = nax, TYPE = IDL_type, /NOZERO)
if nslice GT 0 then for i=0,nslice-1 do readu,unit,tmp
endif else point_lun, unit, currpoint
endif
if not (SILENT) then begin ;Print size of array being read
if pointlun GT 0 then message, $
'Reading FITS extension of type ' + xtension, /INF
snax = strtrim(NAX,2)
st = snax[0]
if Naxis GT 1 then for I=1,NAXIS-1 do st = st + ' by '+SNAX[I] $
else st = st + ' element'
st = 'Now reading ' + st + ' array'
if (pointlun GT 0) and (pcount GT 0) then st = st + ' + heap area'
message,/INF,st
endif
; Read Data in a single I/O call
data = make_array( DIM = nax, TYPE = IDL_type, /NOZERO)
if nbytesleft LT N_elements(data) then begin
message,/CON,NoPRINT=Silent, $
'ERROR - End of file encountered while reading data array'
free_lun,unit
return,-1
endif
readu, unit, data
if (pointlun GT 0) and (pcount GT 0) then begin
theap = sxpar(header,'THEAP')
skip = theap - N_elements(data)
if skip GT 0 then begin
temp = bytarr(skip,/nozero)
readu, unit, skip
endif
heap = bytarr(pcount*gcount*abs(bitpix)/8)
readu, unit, heap
endif
free_lun, unit
; If necessary, replace NaN values, and convert to host byte ordering
check_NaN = (bitpix LT 0) and ( N_elements(NaNvalue) GT 0 )
if check_NaN then NaNpts = whereNaN( data, Count)
ieee_to_host, data
if check_NaN then $
if ( Count GT 0 ) then data[ NaNpts] = NaNvalue
; Scale data unless it is an extension, or /NOSCALE is set
; Use "TEMPORARY" function to speed processing.
do_scale = not keyword_set( NOSCALE )
if (do_scale and ( PointLun GT 0)) then do_scale = xtension EQ 'IMAGE'
if do_scale then begin
Nblank = 0
if bitpix GT 0 then begin
blank = sxpar( header, 'BLANK')
if !ERR NE -1 then $
blankval = where( data EQ blank, Nblank)
endif
bscale = float( sxpar( header, 'BSCALE' ))
if !ERR NE -1 then $
if ( Bscale NE 1. ) then begin
data = temporary(data) * Bscale
sxaddpar, header, 'BSCALE', 1.
sxaddpar, header, 'O_BSCALE', Bscale,' Original BSCALE Value'
endif
bzero = float( sxpar ( header, 'BZERO' ) )
if !ERR NE -1 then $
if (Bzero NE 0) then begin
data = temporary( data ) + Bzero
sxaddpar, header, 'BZERO', 0.
sxaddpar, header, 'O_BZERO', Bzero,' Original BZERO Value'
endif
if Nblank GT 0 then begin
if keyword_set(NaNValue) then $
data[blankval] = NaNvalue else $
data[blankval] = float( [127b,192b,0b,0b],0,1);NaN default
endif
endif
; Return array
return, data
end
Return to the shapelets web page or the code help menu.
| Last modified on 2nd Mar 2009 by Richard Massey. |