MATLAB routines to read FITS

[Stephen Walton]

Lewis C. Roberts wrote:

> Hi,
>
> Does anyone have a Matlab routine to read in FITS files that
> they would be willing to share?

Here's one which was posted on comp.soft-sys.matlab a while back.  It is
short enough to put here, I hope, and be found useful.  If  Netscape
wordwraps it in some weird way, let me know and I'll e-mail a clean
copy.

function output_image=fitsread(filename,n_hdr_areas)
% FITSREAD reads a very simple FITS image into a Matlab variable.
%
%  data=fitsread(filename)
%  data=fitsread(filename,n_hdr_areas)
%
%  The second form explicitly indicates the number of ascii header
%  areas present in the file. Explicitly telling FITSREAD
%  how many header areas are present greatly speeds up
%  the file access, since in this case the function does not need to
%  parse each header line to determine where the header ends
%  and the data begins.
%
%  The following syntax is also acceptible:
%
%    data=fitsread(filename,'hst')
%    data=fitsread(filename,'stsci')
%
%  where the 'hst' and 'stsci'  keywords indicate FITSREAD
%  should use the number of HDUs (=6) in some typical FITS files
%  originating from the Space Telescope Science Institute.
%
%
%  Useful FITS Definitions:
%
%  "card" = 80 byte line of ASCII data in a file. This
%           contains keyword/value pairs and/or comments.
%
%  "Header area" or "HDU" = a set of 36 "cards". Note
%          that FITS files must have an integer number of header
%          areas (typically between 1 and 6, depending on how much
%          header information is stored in the file.)
%
%
% Known deficiences:
%
%  1. does not use bscale or bzero -- if your data needs to be
%     rescaled you can add code to look for these keywords
%
%  2. does not read anything more complicated than a simple
%     image, ie. no binary tables or data cubes.


% Version 2.1
% R. Abraham, Institute of Astronomy, Cambridge University
%
% Changes
% Oct. 98 -- added code to detect 80X86 and Alpha machines and to read
%            data into these with the big endian switch


%The first few cards of the first header/data unit must give information

%in a pre-defined order (eg. the data format, number of axes,
%size etc) but after that the header keywords can come in any
%order. The end of the last card giving information is flagged by
%the keyword END, after which blank lines are added to pad the
%last header so it also contains 36 cards. After the last card
%the data begins, in a format specified by the BITPIX keyword.
%The dimensions of the data are specified by the NAXIS1 and
%NAXIS2 keywords.
%
%Reference:   NASA/Science Office of Standards and Technology
%           "Definition of the Flexible Image Transport System"
%                            NOST 100-1.0
%
%           This and other FITS documents are available on-line at:
%           http://www.gsfc.nasa.gov/astro/fits/basics_info.html

%Set flag indicating unknown number of HDUs
if nargin<2
   n_hdr_areas='unknown';
end;

if isstr(n_hdr_areas)
   n_hdr_areas=upper(n_hdr_areas);
end;

%Allow user to specify a few keywords to indicate number of header areas

if strmatch(n_hdr_areas,'HST') | strmatch(n_hdr_areas,'STSCI') | ...
   strmatch(n_hdr_areas,'MDS')
   n_hdr_areas=6;
end
if strmatch(n_hdr_areas,'FREI')
   n_hdr_areas=2;
end

% try to figure out if we need to swap bytes. This is
% imperfect as I don't know the endian-ness of each
% architecture, so I'm only looking for ones I know for
% sure.
friend = computer;
if strmatch(friend,'PCWIN')
   bswap = 'b';
elseif strmatch(friend,'LNX86')
   bswap = 'b';
elseif strmatch(friend,'ALPHA')
   bswap = 'b';
else
   bswap = 'l';
end

%Open the file
fid=-1;
if ~isstr(filename)
   filename=setstr(filename);
end;
if (isempty(findstr(filename,'.'))==1)
   filename=[filename,'.fits'];
end
[file,message] = fopen(filename,'r','l');
if file == -1
   error(message);
end

%First five cards must contain specific information
[d,simple,d]=parse_card(setstr(fread(file,80,'uchar')'));
[d,bitpix,d]=parse_card(setstr(fread(file,80,'uchar')'));
[d,naxis,d]=parse_card(setstr(fread(file,80,'uchar')'));
[d,naxis1,d]=parse_card(setstr(fread(file,80,'uchar')'));
[d,naxis2,d]=parse_card(setstr(fread(file,80,'uchar')'));

if strmatch(n_hdr_areas,'UNKNOWN')
   %Keep reading cards until one turns up with the keyword 'END'.
   n_card=5;
   keyword='NIL';
   while(~strcmp(deblank(upper(keyword)),'END'))
      n_card=n_card+1;
      card=setstr(fread(file,80,'uchar')');
      [keyword]=parse_card(card);
   end;
   %Go past the blank lines of padding before the start of the data
   if (rem(n_card,36) ~= 0)
      n_blanks = 36 - rem(n_card,36);
      dummy=fread(file,n_blanks*80,'uchar');
   end;
else
   dummy=fread(file,((n_hdr_areas*36)-5)*80,'uchar');
end;

%Read the data.
if bitpix==-64
   X=fread(file,naxis1*naxis2,'float32',bswap);
elseif bitpix==-32
   X=fread(file,naxis1*naxis2,'float',bswap);
elseif bitpix==8
   X=fread(file,naxis1*naxis2,'uint8',bswap);
elseif bitpix==16
   X=fread(file,naxis1*naxis2,'short',bswap);
elseif bitpix==32
   X=fread(file,naxis1*naxis2,'long',bswap);
else
   error('data type specified by BITPIX keyword is not -64, -32, 8, 16,
or 32');
end;

%Clean up and output data
fclose(file);

output_image=reshape(X,naxis1,naxis2);



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [keyword,value,comment] = parse_card(s)
%Parses a FITS header card.
%Reference:
%                NASA/Science Office of Standards and Technology
%           "Definition of the Flexible Image Transport System (FITS)"
%                        NOST 100-1.0    June 19, 1993

%Set defaults
keyword=[];value=[];comment=[];

%Get keyword in bytes 1 - 8
keyword=s(1:8);
if nargout==1
   return;
end

%If keyword is blank then the line is a comment
if strmatch(keyword,'       ')
    keyword=[];
 value=[];
 comment=deblank(s(11:80));
 return;
end;


%Keyword is non-blank. Check if there is a corresponding value/comment.
%If not then the only possibilities are that bytes 11:80 are a comment
%or that they are blank
if ~strmatch(s(9:10),'= ')
    keyword=deblank(keyword);
 value=[];
 comment=deblank(s(11:80));
 return;
end;

%Card is a standard keyword/value/comment structure. Break the
value/comment
%string (bytes 11 - 80) into separate strings by tokenizing on "/"
character.
%Remove the leading and trailing blanks on the value and the trailing
blanks
%on the comment.

keyword=deblank(keyword);
[value,comment]=strtok(s(11:80),'/');
comment=deblank(comment);
value=fliplr(deblank(fliplr(deblank(value))));

%Now figure out whether to output the value as a string or as a number.
%The FITS standard requires all values that are strings to be in single
%quotes like this: 'foo bar', so I can simply look for occurences of a
%single quote to flag a string. However, logical variables can take the
%values T or F without having any single quotes, so I'll have to look
%out for those also.

%Test for logical. Return logical as a string.
if strmatch(upper(value),'T') | strmatch(upper(value),'F')
 return;
end;

%Test for string. Return string unconverted.
if length(findstr('''',value)) ~= 0
 return;
end;

%Only thing left is a number. Convert string to number.
value=str2num(value);