search

ratt-ru / meqtrees

A library for implementing radio astronomical Measurement Equations
http://meqtrees.net
10 stars 2 forks source link

FITS reading/writing issues #601

Closed gijzelaerr closed 10 years ago

gijzelaerr commented 10 years ago
at 2007-10-04 00:02:32 Tony Willis reported:

FITS reading/writing issues

gijzelaerr commented 10 years ago

Original comment thread migrated from bugzilla

at 2007-10-04 00:02:32 Tony Willis replied:

Hi Sarod

I'm just copying my e-mail of 10 August here so that things don't get lost forwever!

Subject: RE: question about FITSWriter node

I looked again at the FITS input / output situation. I think you are still making some errors in both your input and output ...

Just to make sure ... I looked at the FITS standard (version of Dec 9, 2005). There the comment to Figure 4.1 states: "Arrays of more than one dimension shall consist of a sequence such that the index along axis 1 varies most rapidly and those along subsequent axes progressively less rapidly." Of course FITS was designed in 1979 and is based on FORTRAN. The python numarrays that we see in the browser, of course, are based on C, and numarray data varies most rapidly along the last axis etc - annoyingly just the opposite of FORTRAN (and FITS!). (Actually I would have thought that the cfitsio package should be handling this stuff in your FITS nodes ...)

To make sure that I could distinguish the last two axes, I modified the veidt_fpa FITS files to have 101 elements in L and 81 in M. So the FITS header looked like:

SIMPLE = T / conforms to FITS standard BITPIX = -64 / array data type NAXIS = 4 / number of array dimensions NAXIS1 = 81 NAXIS2 = 101 NAXIS3 = 1 NAXIS4 = 1 EXTEND = T CTYPE1 = 'M ' CDELT1 = -0.171887338539247 / in degrees for Meq.FITSImage CRPIX1 = 41 / reference pixel (one relative) CRVAL1 = 0.0 / M = 0 at beam peak CUNIT1 = 'deg ' CTYPE2 = 'L ' CDELT2 = 0.171887338539247 / in degrees for Meq.FITSImage CRPIX2 = 51 / reference pixel (one relative) CRVAL2 = 0.0 / L = 0 at beam peak CUNIT2 = 'deg ' CTYPE3 = 'TIME ' CDELT3 = 1 / in sec CRPIX3 = 1 / in pixels (one relative) CRVAL3 = 1.0 / equates to grid point CTYPE4 = 'FREQ ' CDELT4 = 1 / in Hz CRPIX4 = 1 / in pixels (one relative) CRVAL4 = 1.0 / equates to grid point CPLX = 0 / false as data is real CELLS = 1 / true as we want cells END

The Meq.FITSImage node converts this into a 1 x 1 x 81 x 101 numarray with the L axis on the 2D image display clearly smaller than that of the M , but I think you should end up with 1 x 1 x 101 x 81, and the M axis smaller than the L axis.

and the FITSWriter is presently converting your 1 x 1 x 81 x 101 (incorrect) internal numarray to the output

SIMPLE = T / file does conform to FITS standard BITPIX = -64 / number of bits per data pixel NAXIS = 3 / number of data axes NAXIS1 = 1 / length of data axis 1 NAXIS2 = 81 / length of data axis 2 NAXIS3 = 101 / length of data axis 3

etc etc

whereas the output should actually resemble the input header shown above.

Interestingly the 'kvis' program seems remarkably smart about handling incorrectly specified FITS files! Unfortunately some other programs (including our own DRAO fits reader program - which cannot handle the current FITSWriter output, and where I first noticed this problem) fall over.

I have put the FPA data with different sizes in L,M as veidt_fpa_diff_LM.tar.gz in my home directory on lofar9.

Cheers

Tony

On Fri, 10 Aug 2007 yatawatta@astro.rug.nl wrote:

Hi Tony, I will fix this and also your FPA problem (I found a solution). cheers, sarod.

Hi Again

I must have encountered this issue before. I found that I have a python script which converts FITS files from BITPIX = -64 to BITPIX = -32. So this is a very low priority. However I think it would be a useful option for users.

However, I do note one potential problem. The array I'm writing out to FITS has a MeqTrees structure of 1 x 1 x 101 x 101. So I would have thought that NAXIS should be 4 , with NAXIS1 and NAXIS2 both having values of 1. Instead, you generate NAXIS =3, etc ....

SIMPLE = T / file does conform to FITS standard BITPIX = -64 / number of bits per data pixel NAXIS = 3 / number of data axes NAXIS1 = 1 / length of data axis 1 NAXIS2 = 101 / length of data axis 2 NAXIS3 = 101 / length of data axis 3

Cheers

Tony

at 2007-10-16 20:53:19 Tony Willis replied:

I went back and ran the WH/contrib/SBY/demos/demo_FITSImage.py script. Indeed 'kvis' and the browser display give consistent results for L,M (browser displays m==X axis and L == Y axis, so browser display is rotated clockwise 90 deg with respect to kvis display) when I load the image with

ns.image<<Meq.FITSImage(filename="demo.fits",cutoff=1.0,mode=2)

So input seems OK for a square image.

at 2007-10-17 01:05:30 Tony Willis replied:

In fact things seem OK with non-square images. I tried a file with the header:

SIMPLE = T / conforms to FITS standard BITPIX = -64 / array data type NAXIS = 4 / number of array dimensions NAXIS1 = 81 NAXIS2 = 101 NAXIS3 = 1 NAXIS4 = 1 EXTEND = T CTYPE1 = 'L ' CDELT1 = -0.171887338539247 / in degrees for Meq.FITSImage CRPIX1 = 41 / reference pixel (one relative) CRVAL1 = 0.0 / L = 0 at beam peak CUNIT1 = 'deg ' CTYPE2 = 'M ' CDELT2 = 0.171887338539247 / in degrees for Meq.FITSImage CRPIX2 = 51 / reference pixel (one relative) CRVAL2 = 0.0 / M = 0 at beam peak CUNIT2 = 'deg ' CTYPE3 = 'TIME ' CDELT3 = 1 / in sec CRPIX3 = 1 / in pixels (one relative) CRVAL3 = 1.0 / equates to grid point CTYPE4 = 'FREQ ' CDELT4 = 1 / in Hz CRPIX4 = 1 / in pixels (one relative) CRVAL4 = 1.0 / equates to grid point CPLX = 0 / false as data is real CELLS = 1 / true as we want cells END

So I think the input is generally OK, as long as the user understands that L gets translated to the third (internal) axis and M to the fourth.

But the output layout (see my first message in the sequence) definitely has problems.

at 2008-06-24 20:37:13 Tony Willis replied:

Another FITS reading issue ...

I've noticed that when I have a FITS file with 4 frequencies and an L, M image with each frequency, leading to the following header

SIMPLE = T / conforms to FITS standard BITPIX = -64 / array data type NAXIS = 4 / number of array dimensions NAXIS1 = 101 NAXIS2 = 101 NAXIS3 = 1 NAXIS4 = 4 EXTEND = T CTYPE1 = 'L ' CDELT1 = -0.171887338539247 / in degrees for Meq.FITSImage CRPIX1 = 51 / reference pixel (one relative) CRVAL1 = 0.0 / M = 0 at beam peak CUNIT1 = 'deg ' CTYPE2 = 'M ' CDELT2 = 0.171887338539247 / in degrees for Meq.FITSImage CRPIX2 = 51 / reference pixel (one relative) CRVAL2 = 0.0 / L = 0 at beam peak CUNIT2 = 'deg ' CTYPE3 = 'STOKES ' CDELT3 = 1 / in IQUV CRPIX3 = 1 / in pixels (one relative) CRVAL3 = 1.0 / equates to grid point CTYPE4 = 'FREQ ' CDELT4 = 200000000.0 / in MHz CRPIX4 = 1 / in pixels (one relative) CRVAL4 = 1100000000.0 / equates to grid point CPLX = 0 / false as data is real CELLS = 1 / true as we want cells

You read this in and assign the reference frequency to the 2nd frequency - i.e. I get frequencies 900, 1100, 1300 and 1500 MHz. However the FITS reference points are all one-based, due to the ancient FORTRAN heritage, so I think you should be generating the sequence 1100, 1300, 1500, 1700 MHz. (So I also wonder if there is an off-by-one issue with the other axes.)

at 2009-01-08 20:52:19 Tony Willis replied:

Hi Sarod

I think that the FITS reading / writing issues have to be bumped up to critical. I think it extremely important that these issues be resolved before there can be a 'production 'release of MeqTrees. I just tried to write out a node with I,Q,U and V Vells. Although the size of the file on disk seems to be correct, when I try to read the file back in, the first image found only has value 1 for the third dimension. I can get the remaining images by various pyfits incantations, but ...

at 2009-01-12 14:13:18 Sarod Yatawatta replied:

Hi Tony, could you please point me to your scripts so that I could reproduce this. Ill fix it as soon as I get back. cheers, sarod

at 2009-01-14 18:48:58 Sarod Yatawatta replied:

Hi Tony, I just went through the documentation on the wiki page: here's the extract: " FITSWriter

This node is used to save a Result to a FITS file. The _(first)_ VellSet and the Cells will be saved in the FITS file. The initrec of this node requires one parameter,

*

  filename : The name of the file to write to. If this file already exists,

the node will fail. However, you can prepend a ! to your filename so that it will be overwritten by default. For example, instead of giving myfile.fits, you can give !myfile.fits so that the file will be overwritten the next time you run this node.

FITSReader

This node is used to read in a FITS file and convert the values to a Result that include _one_ VellSet and Cells. The initrec of this node requires one parameter,

*

  filename : The name of the file to read.

"""

If you look at the _highlighted_ words, youll see that the writing/reading only works for the first vellset. Writing multiple vellsets to a single FITS file has a fundamental problem : some vellsets can be a constant, or the dimensions might not match up (eg 1D, 2D, 4D)

So what you have encountered is a missing feature. The best way to get around this is to write each vellset to a different FITS file. This will not reproduce an image, but writer/reader was never meant to be used to create FITS images.

cheers, sarod

at 2009-01-15 03:42:07 Tony Willis replied:

Hi Sarod - you should be able to fix the 'additional comment #3 re reading the frequency correctly . In the WH/contrib/RJN directory there is the MG_RJN_UVBrick_FITS.py script. I also checked in a small data cube point_source_cube.fits which has l x m = 8 x 8 (NAXIS1 and NAXIS2). NAXIS3 is Stokes(sorry - in the file is wrongly labeled as 'THIRD COORDINATE' and NAXIS4 should be frequency (this stuff is correctly labeled in the CRVAL4 etc stuff. CRVAL4 is given as 1400000000.0 and CRPIX4 is given as 1. If you execute the script (you may have to edit line 116 to specify the infile directory location properly) and look in the 'corr' bookmarks node, you will see that you have assigned the first frequency point as 1350, and not 1400. This leads me to believe that point 1 in a zero-relative system has value 1400, whereas FITS assumes that point 1 is the first point as FITS was originally defined when FORTRAN was the dominant language.

I'll get busy on an example of the writing out problem ...

at 2009-01-15 05:42:08 Tony Willis replied:

Created an attachment (id=54) script that demonstrates fits writer problems

This script can be used in conjunction with the quadrant.fits image to check whether the FITS writer produces output consistent with the input.

at 2009-01-15 05:43:45 Tony Willis replied:

Created an attachment (id=55) a simple fits file with four non-zero pixels

This fits file can be used to check if the fits writer node is giving sensible output

at 2009-01-15 05:51:15 Tony Willis replied:

Created an attachment (id=56) a comparison of what's currently readin via the FITSImage node vs what gets written out with the FITS writer

This attachment shows (left) the original quadrant.fits image. The right frame shows what one gets when the image gets passed though the fits_writer_test.py script. Amongst other things there is a 90 deg rotation. However we are only lucky enough to get this because kvis seems remarkable sophisticated at trying to produce something useful. Try looking at this with the 'ds9' viewer!

at 2009-01-15 06:05:40 Tony Willis replied:

Hi Sarod

The quadrant.fits file has

NAXIS = 4 / NUMBER OF AXES
NAXIS1 = 8 / # PIXELS PER ROW (E.G. RA)
NAXIS2 = 8 / # ROWS (E.G. DEC)
NAXIS3 = 1
NAXIS4 = 1

which is a valid FITS header and which you correctly convert into an array inside MeqLand with dimensions 1 x 1 x 8 x 8. When one uses the browser to look at the internal record structure of ns.fits in the fits_writer_test.py its still has correctly an array with 1 x 1 x 8 x 8. However the FITS file that is written out has

NAXIS = 3 / number of data axes
NAXIS1 = 1 / length of data axis 1
NAXIS2 = 8 / length of data axis 2
NAXIS3 = 8 / length of data axis 3

So you have lost an axis for starters. The other issue is this damn FORTRAN ordering. In the incoming array NAXIS1 and NAXIS2 are the quantities that are adjacent in memory and you correctly assign them to l, m etc. However, on the output side something is going astray!! You have an NAXIS1 which is the fastest changing index with a dimension of 1! So ds9 interprets the output as a 1x8 vector with 8 somethings (say 8 frequencies) and tries to create a movie with 8 one-dimensional frames!

at 2009-02-11 17:13:33 Oleg Smirnov replied:

Is this working to Tony's satisfaction now? Need to make sure before release.

at 2009-02-11 20:02:21 Tony Willis replied:

The input frequencies are now OK.

I think we have agreed to disagree on whether the output can be considered valid.

In any case, it's easy enough to fiddle the output via PyFits to get the fits file into what I would consider a 'valid' format. So the bug is certainly no longer 'critical' or a blocker for the release.

at 2009-02-12 17:00:45 Oleg Smirnov replied:

In any case, it's easy enough to fiddle the output via PyFits to get the fits file into what I would consider a 'valid' format. So the bug is certainly no longer 'critical' or a blocker for the release.

I have downgraded the bug accordingly. I'll leave it to you (Tony & Sarod) to decide when it should be closed.

at 2009-02-12 23:17:26 Oleg Smirnov replied:

Things seem resolved. Additional discussion here: http://www.astron.nl/meqwiki/MeqImage