4

I am trying to read FITS file created by HEALPix software. The file contains a spherical map of a scalar field.

There is the online support on how to use the file in IDL, but not for Mathematica.

I understand FITS is partly binary format, but I expect it is still standardised and contains readable description on how the data is stored.

I tried reading Elements from the file in this way

Import["test.fits", #] & /@ 
  Import["test.fits", "Elements"]

but the output is quite short and mainly empty lists, especially for

Import["test.fits", "Data"]
Import["test.fits", "TableData"]

The size of my file is 12.6 MB. Here is a link to the file.

Is there another way to read the content of FITS file, besides these general 'Elements'?

Do I need to know something about the structure of my FITS file in order to be able to read it with Mathematica? Could it be the problem is that Mathematica does not read all types of FITS files?

Vladimir
  • 1,503
  • 1
  • 12
  • 28
  • 1
    For a binary file, structure is in general needed in order to parse and read the file correctly. Otherwise Mathematica does not know if the bit it reads is a part of an integer or a real number, or where the number starts. You can have a look at BinaryRead and BinaryReadList once you know the structure of the file. – Yi Wang Sep 01 '14 at 10:14
  • Sure, Yi, only that this is a FITS file format – Vladimir Sep 01 '14 at 14:09
  • 1
    Can you post the file? – mfvonh Sep 01 '14 at 15:22
  • FITS is supported, http://reference.wolfram.com/language/ref/format/FITS.html . By the way the"fits" tag may not be related to this question. – rhermans Sep 03 '14 at 22:35
  • 1
    @mfvonh: Unfortunatelly, uploading never finishes. It seems 13MB is too much. I'll try to make a smaller example file. However, I discovered that neither IDL can open the FITS file from HEALPix. It seems these guys modified the format a little bit and have even created a specific function in IDL to work with their format. The only hope know is to use a header of the FITS file and understand the structure of the rest of the file which is binary. – Vladimir Sep 03 '14 at 22:42
  • @rhermans: Thank you for the tag correction! – Vladimir Sep 03 '14 at 22:42
  • @Vladimir Just in case you have not considered this, people often use direct links from services like dropbox to share larger files here. – mfvonh Sep 04 '14 at 00:27
  • @mfvonh: Thank you. Here is the FITS file, which makes the trouble. – Vladimir Sep 04 '14 at 10:10
  • related http://mathematica.stackexchange.com/questions/43791/exporting-data-extensions-in-fits-file/44460#44460 – george2079 Sep 04 '14 at 13:47

1 Answers1

3

Your file seems to be corrupt, there is one extra byte in the header section. Open in a binary clean text editor, search for "uK^2" and delete exactly one space following the 2.

The file can then be read by this: (adapted from link in comment)

 f = OpenRead["test.fits", BinaryFormat -> True];
 parsehead[hh_] := (metadat[#[[1]]] = 
    #[[2]]; #) & /@ ({StringReplace[#[[1]], " " -> ""],
    If[Length[#] > 1, StringReplace[#[[2]], " " -> ""], ""],
    If[Length[#] > 2, #[[3 ;;]], ""]} &@ 
   StringSplit[ # , {"=", "/"}] & /@ hh)
 getheaders[f_] := 
 parsehead[First@Last@
       Reap[ While[ (Sow[
            StringJoin@BinaryRead[f, Table["Character8", {80}]]] != 
                StringJoin["END", Table[" ", {77}]])] ]]
 hh = getheaders[f];
 Print["lines read = ", nh = Length[hh]]
 blkpad = BinaryRead[f, Table["Byte", {2880 - Mod[nh 80, 2880]}]];
 Union[blkpad] == {32} (* verfy we just read blank *)
 Print["gettingheaders"]
 hh = getheaders[f];
 Print["lines read = ", nh = Length[hh]]
 hh // MatrixForm
 BinaryRead[f, Table["Byte", {2880 - Mod[nh 80, 2880]}]];

(data size/type hard coded based on looking at the header data..)

 data = BinaryRead[f, Table["UnsignedInteger8", {3072}, {4096}], ByteOrdering -> 1];
 blkpad = BinaryRead[f, Table["Byte", {2880 - Mod[ (3072 4096 ) , 2880]}], 
       ByteOrdering -> 1];
 Union[blkpad] == {32}
 StreamPosition[f] == FileByteCount[First@f] (* verfy entire file read *)
 Close[f]

The resulting image is pretty much snow, but then I don't know what its supposed to be.

 Image[data/2^8]

enter image description here

the headers:

enter image description here enter image description here

george2079
  • 38,913
  • 1
  • 43
  • 110
  • Thank you, George. This is extensive code. The image is a simulated CMB map, so it should not look like this at all. It has colours and I think the file should contain 5 (or 7) columns of 8bit numbers that can be plotted as colours. As I said, it is probably to difficult to get this, since HEALPix is probably using a modified FITS format (which is open source), and hence it is not standard any more... – Vladimir Sep 05 '14 at 00:01
  • 1
    the raw data is a 3072x4096 byte array. Once you have that you can sort out how to partition it up. There is also a ton of header information (stored in hh in the above code) - i didn't study it to see if it describes the data structure. – george2079 Sep 05 '14 at 12:05