Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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
share|improve this question
    
I'm looking forward to see the answers :) I have a question, do you know how to import multidimensional fits files into MMA? I mean, one file with multiple images for example. I was never able to do this. –  Kuba Mar 19 at 8:03
    
@Kuba No but it would have been my next question. Let's try and trigger some interest from WRI... –  chris Mar 19 at 8:35
    
Wolfram Community is a good place for this but I'm not particiating there, one addiction is enough ;) –  Kuba Mar 19 at 8:39
    
my wag you will need to directly Write the file rather than Export fits.gsfc.nasa.gov/fits_primer.html –  george2079 Mar 19 at 11:57
add comment

4 Answers 4

up vote 5 down vote accepted
+200

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.

enter image description here

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];
share|improve this answer
    
Thank you for your contribution. I am not yet able to access the second table but will try again later. –  chris Mar 21 at 7:39
    
what are you using to try to read it? This is a tough problem to work without access to an external application that reads the format. –  george2079 Mar 21 at 12:55
    
based on the example I changed the XTENSION type to IMAGE -- see if that does it. –  george2079 Mar 21 at 15:13
    
This seems to do the trick; with ` dat2 = Table[i, {i, 6}]` and "Data" -> dat2 in the second write. Pity its not directly available in Export. –  chris Mar 21 at 16:02
add comment

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.

share|improve this answer
add comment

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"]

"foo.fits"

data = Import["1Example.fits", "RawData"];
dims = "MYVALUE" /. Import["1Example.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
share|improve this answer
    
Thank you for your attempt. It seems like a solution, as long as the length of dat1 and dat2 are comparable; but if dat2, say , is much shorter it is a bit of a waste? –  chris Mar 21 at 7:37
    
Yes, it's a waste of space. Do you have an example of a FITS file where the sizes are different? Say one available online? I couldn't find one, but I don't deal with FITS regularly. –  Michael E2 Mar 21 at 10:30
    
I ll get you one. –  chris Mar 21 at 10:33
    
@Chris The update does not waste space, but, like my other one, it must be a nonstandard solution. –  Michael E2 Mar 21 at 12:58
    
Thanks again. It does deal with the original challenge. Indeed it is not standard which would make reading in the data in other language a bit cumbersome. –  chris Mar 21 at 13:35
show 2 more comments

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]

enter image description here

share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.