Exporting data extensions in fits file

Question

Context

In Astronomy the de facto standard for images and data is FITS. I would like to export heterogenous data sets into a single fits file.

Attempt

I am interested in exporting into a fits file some data, say.

 dat = Table[i, {i, 5}]

(* {1,2,3,4,5} *)

and additionally some metadata, say a range of $\nu$s:

 rnus = ("nu" <> ToString[#] -> 0.5 + # & /@ ( Range[5]))

(* {nu1->1.5,nu2->2.5,nu3->3.5,nu4->4.5,nu5->5.5} *)

Following closely the documentation I can write:

 Export["image.fits", {"Data" -> dat, 
 "Metadata" -> Join[rnus, {"Object" -> "something"}]}, "Rules"]

(* image.fits *)

and read back accordingly both the data:

 Import["image.fits", "RawData"][[1]]

(* {1,2,3,4,5} *)

and the metadata:

  "NU1" /. Import["image.fits", "Metadata"]
  "OBJECT" /. Import["image.fits", "Metadata"]

(* 1.5 ) ( Something *)

BUT, what I would like to do is to write different (large) tables of different sizes into a single fits file. So using a single keyword for each element is unrealistic.

If the table were of the same size I could join them into a single list, as in

   Export["image.fits", {"Data" -> {dat,dat}}, "Rules"]     

but in my case the tables are not of the same size.

In principle the documentation says "For FITS files that contain multiple images or data extensions, the above elements are taken to be lists of the respective expressions."

Question

How can I write multiple images or data extensions into a unique fits file? Or in other words, what does the documentation mean?

Update

As requested, here is an example of fits with an extension containing {1,2,3,4,5} and {1,2,3,4,5,6} :

https://dl.dropboxusercontent.com/u/659996/test.fits

or it can be produced in python as

python
import pyfits
import numpy
data=numpy.array([1.,2.,3.,4.,5.])
pyfits.writeto("test.fits",data)
data=numpy.array([1.,2.,3.,4.,5.,6.])
pyfits.append("test.fits",data)
exit

Answer 1

This is a bit of a shot in the dark but this general approach may work with a bit of tweaking. (note another edit fix..)

 dat = Table[i, {i, 5}]
 rnus = ("nu" <> ToString[#] -> 0.5 + # & /@ (Range[5]))
 BinaryWrite[g = OpenWrite["test.fits", BinaryFormat -> True],
          ExportString[{"Data" -> dat, "Metadata" -> Join[rnus,
             {"Object" -> "something", "NEXTEND" -> 1}]}, {"FITS", "Rules"}]]
 BinaryWrite[g,                      
     StringReplace[
         ExportString[
             {"Data" -> dat, 
                "Metadata" -> Join[rnus, {"Object" -> "something"}]}, 
                  {"FITS",  "Rules"}], 
               "SIMPLE  =                    T" -> 
               "XTENSION= 'IMAGE   '          "]]
  Close[g];

Edit -- note that StringReplace must preserve the string length.

Here's what my reader shows (changing the data type to integer) .. note mathematica has added a bunch of extra metadata compared to the example.

Approach 2

This is a more low level approach, forgoing Export

 stringpad[s_, n_] := StringJoin[s, Table[" ", {n - StringLength@s}]];
 writefitshead[f_][hdat_] := (
    BinaryWrite[f, stringpad[StringJoin[stringpad[#[[1]], 8],
      If[Length[#] > 1, StringJoin["= ", #[[2]], " "], ""],
      If[Length[#] > 2, StringJoin["/", #[[3]]], ""]], 80]] & /@ hdat;
    BinaryWrite[f, stringpad["", 80]] & /@ Range[ 36 - Length[hdat]];)
 bytelen["Real64"] = 8;
 bytelen["Integer16"] = 2;
 bytelen["Integer8"] = 1;
 writefitsdat[f_, type_][dat_] := ( 
    BinaryWrite[f, dat, type, ByteOrdering -> 1] ;
    BinaryWrite[f, 
       Table[0, {i, 2880 - Mod[Length[dat]  bytelen[type], 2880]}], "Integer8"] );
 f = OpenWrite["mytest.fits", BinaryFormat -> True];
 writefitshead[f]@{
    {"SIMPLE", "T", "conforms to FITS standard "},
    {"BITPIX", "-64", "array data type"},
    {"NAXIS", "1", " number of array dimensions"},
    {"NAXIS1", "5"},
    {"EXTEND", "T"},
    {"END"}};
 writefitsdat[f, "Real64"]@Table[i, {i, 5}]; (*note the type here has to match BITPIX*)
 writefitshead[f]@{
    {"XTENSION", "'IMAGE   '", "Image extension "},
    {"BITPIX", "-64 ", " array data type  "},
    {"NAXIS", "1", "number of array dimensions"},
    {"NAXIS1", "6"},
    {"PCOUNT", "0", "number of parameters"},
    {"GCOUNT", "1"},
    {"END"}};
 writefitsdat[f, "Real64"]@Table[i, {i, 6}];
 Close[f];

Answer 2

This reads your example: a crude start anyway, i'm not parsing the headers just showing them as strings.

 f = OpenRead["test.fits", BinaryFormat -> True];
 First@Last@
     Reap[ While[ (Sow[
          StringJoin@BinaryRead[f, Table["Character8", {80}]]] != 
           StringJoin["END", Table[" ", {77}]])] ]
 Print["lines read = ", nh = Length[%]]
 BinaryRead[f, Table["Byte", {2880 - Mod[nh 80, 2880]}]];

manually set the data type and length: it looks like "BITPIX = -64" denotes 64 bit floats. Obviously you need to parse the headers to automate this. On the next line the first arg to Mod is 5 values * 8 bytes

 BinaryRead[f, Table["Real64", {5}], ByteOrdering -> 1]
 BinaryRead[f, Table["Byte", {2880 - Mod[5 8, 2880]}]]; 

 First@Last@
    Reap[ While[ (Sow[
       StringJoin@BinaryRead[f, Table["Character8", {80}]]] != 
            StringJoin["END", Table[" ", {77}]])] ]
 Print["lines read = ", nh = Length[%]]
 BinaryRead[f, Table["Byte", {2880 - Mod[nh 80, 2880]}]];

again manually set the length/type

 BinaryRead[f, Table["Real64", {6}], ByteOrdering -> 1]
 Close[f]

Answer 3

I put here a function built from @george2079's answer as it may be of interest to other people.

Clear[FitsExport];
FitsExport[fname_, dats_, descriptions__] :=
Module[{g},
  BinaryWrite[g = OpenWrite[fname, BinaryFormat -> True], 
  ExportString[{"Data" -> dats[[1]], 
  "Metadata" -> {"Field" -> descriptions[[1]], 
   "NEXTEND" -> 1}}, {"FITS", "Rules"}]];
  If[Length[dats] > 1,
  Do[
   BinaryWrite[g, StringReplace[
   ExportString[{"Data" -> dats[[i]], 
     "Metadata" -> {"Field" -> descriptions[[i]]}}, {"FITS", "Rules"}], 
   "SIMPLE  =                    T" -> 
    "XTENSION= 'IMAGE   '          "]],
    {i, 2, Length[dats]}];, ""];
  Close[g];
  ]

  FitsExport[fname_, dats_] :=FitsExport[fname, dats, Table[" ", {Length[dats]}]];

And as an example:

FitsExport["test.fits", 
  {Range[5], Table[i + j + 0.5, {i, 2}, {j, 3}], Range[3] + 0.5},
  {"a", "b", "c"}
 ]

Note that mathematica's Import actually works on such file:

   u = Import["test.fits", "RawData"];

(* {{1, 2, 3, 4, 5}, {{2.5, 3.5, 4.5}, {3.5, 4.5, 5.5}}, {1.5, 2.5, 3.5}} *)

Note finally that one can use fits files to import multiple tensors and their name and use mathematica to do the assignment on the fly as follows.

 Options[FitsImport] = {SetVariables -> False};
  SetVariables::usage = "SetVariables is an option for FitsImport;
  if set the variables names are set from the fits file using the field FIELD;
  default value is False";

 FitsImport::usage = 
 "dat=FitsImport[fname]; read into the extented fits file fname;
  FitsImport[fname,SetVariables\[Rule] True] set variables from \
  extended fits file fname, using as variable name the content of 
  the field FIELD; WARNING it will override existing names";

  FitsImport[fname_, opts : OptionsPattern[]] :=
  Module[{var, head,dat},
  head = Import[fname, "Metadata"];       
  var = Select[Flatten["FIELD" /. head], StringQ];
  var // TableForm // Print;
  dat = Import[fname, "RawData"];
  If[OptionValue[SetVariables],
    Clear /@ var;
   Do[Set[Evaluate@ ToExpression[var[[i]]], dat[[i]]], {i, Length[var]}]; ,
    dat]
   ]

So that

   FitsImport["test.fits", SetVariables -> True]

(* a, b, c *)

    a 

(* {1,2,3,4,5} *)

Clearly for the assignment to be operational the descriptors of the field must be possible variable names. So AGivenName is ok, but A Given Name is not for instance. Once again, when SetVariables is set to True the corresponding variables are cleared first, so beware of side effects.

Answer 4

Update

I use Mr.Wizard dynamicPartition (code included below).

dat1 = Range@5;
dat2 = Range@6;
metadata = {
   "MetaInfo" -> "MyData",
   "MyValue" -> ToString[Dimensions /@ {dat1, dat2}]
   };
Export["foo.fits", {"Data" -> Flatten@{dat1, dat2}, "Metadata" -> metadata}, "Rules"];

Example 1:

data = Import["1Example.fits", "RawData"]
(*
  {{1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 6}}
*)

dims = "MYVALUE" /. Import["1Example.fits", "Metadata"] // ToExpression
(*
  {{5}, {6}}
*)

dynamicPartition[First@data, Times @@@ dims]
(*
  {{1, 2, 3, 4, 5}, {1, 2, 3, 4, 5, 6}}
*)

Example 2:

m1 = Import["http://exampledata.wolfram.com/messier61.fits.gz", "RawData"];
m2 = Import[
   "http://www.atnf.csiro.au/people/mcalabre/data/WCS/1904-66_ARC.fits.gz", "RawData"];

metadata = {
   "MetaInfo" -> "MyData" ,
   "DataSets" -> 2, 
   "MyValue" -> ToString[Dimensions /@ First /@ {m1, m2}]
   };
Export["foo.fits", {"Data" -> Flatten@{m1, m2}, "Metadata" -> metadata}, "Rules"];

data = Import["foo.fits", "RawData"];
dims = "MYVALUE" /. Import["foo.fits", "Metadata"] // ToExpression;

{n1, n2} = MapThread[
   ArrayReshape, {dynamicPartition[First@data, Times @@@ dims], dims}];

{n1} == m1
{n2} === m2
(* 
  True
  True
*)

---

Mr.Wizard's code, included for convenience:

dPcore[L_, p : {q___, _}] := 
 Inner[L[[# ;; #2]] &, {0, q} + 1, p, Head@L]

dPcore[L_, p_, All] := dPcore[L, p]~Append~Drop[L, Last@p]

dPcore[L_, p_, n__] := dPcore[L, p]~Join~Partition[L~Drop~Last@p, n]

dynamicPartition[L_, p : {__Integer}, x___] := 
 dPcore[L, Accumulate@p, x] /; ! Negative@Min@p && Length@L >= Tr@p