On Fri, 10 Aug 2007, David Berry wrote:
> On 09/08/07, Peter W. Draper <[log in to unmask]> wrote:
>> On Thu, 9 Aug 2007, David Berry wrote:
>>
>>>
>>> Not with you. So when you come to write out the new FrameSet to the
>>> FitsChan, how does the presence of PROJP keywords in the FitsChan
>>> cause the output FITS header to be written in FITS-PC?
>>
>> It's not the write that's the problem, that's fine, it's the read when the
>> FITS file is re-opened.
>>
>> Perhaps the title should be:
>>
>> Re: FITS-WCS written, but read back as FITS-PC.
>>
>>> Do you set Encoding before writing the WCS FrameSet out to the FitsChan,
>>> or do you use the default?
>>
>> I set it to FITS-WCS.
>>
>> So (condensed headers from the files), initially I have:
>>
>> CTYPE1 = 'RA---ZPN' / Algorithm type for axis 1
>> CTYPE2 = 'DEC--ZPN' / Algorithm type for axis 2
>> CRPIX1 = 3.431523900000000E+03 / Dither offset Y
>> CRPIX2 = -4.111654500000001E+03 / Dither offset Y
>> CRVAL1 = 2.8992249E+02 / [deg] Right ascension at the reference pixel
>> CRVAL2 = 1.5713925E+01 / [deg] Declination at the reference pixel
>> CRUNIT1 = 'deg ' / Unit of right ascension co-ordinates
>> CRUNIT2 = 'deg ' / Unit of declination co-ordinates
>> CD1_1 = 2.0097707E-07 / Transformation matrix element
>> CD1_2 = 5.5739441E-05 / Transformation matrix element
>> CD2_1 = -5.5726214E-05 / Transformation matrix element
>> CD2_2 = 2.7165242E-07 / Transformation matrix element
>> PV2_1 = 1. / Pol.coeff. for pixel -> celestial coord
>> PV2_2 = 0.000000E+00 / Pol.coeff. for pixel -> celestial coord
>> PV2_3 = -50.
>> PROJP1 = 1. /
>> PROJP3 = -50. /
>>
>> which reads as a FITS-WCS. That WCS is then changed in GAIA to refine the
>> positions, and re-written, as FITS-WCS (after a destructive read on the
>> above) and becomes:
>>
>> CTYPE1 = 'DEC--TAN' / Type of co-ordinate on axis 1
>> CTYPE2 = 'RA---TAN' / Type of co-ordinate on axis 2
>> CRPIX1 = 150.0 / Reference pixel on axis 1
>> CRPIX2 = 150.0 / Reference pixel on axis 2
>> CRVAL1 = 15.8980643331 / Value at ref. pixel on axis 1
>> CRVAL2 = 290.169112711 / Value at ref. pixel on axis 2
>> CRUNIT1 = 'deg ' / Unit of right ascension co-ordinates
>> CRUNIT2 = 'deg ' / Unit of declination co-ordinates
>> PROJP1 = 1.0
>> PROJP3 = -50.0
>> NUMZPT = 6642 / Number of standards used
>> CDELT1 = -5.59098742690742E-5/ Pixel size on axis 1
>> CDELT2 = 5.59121728822906E-5 / Pixel size on axis 2
>> PC1_1 = 1.0 / Transformation matrix element
>> PC1_2 = -0.00442691692661737 / Transformation matrix element
>> PC2_1 = 0.00212609392345515 / Transformation matrix element
>> PC2_2 = 1.0 / Transformation matrix element
>
>
> Surely this is a bug in AST???? In the text above the headers you say
> that you re-write is as FITS-WCS, by which I presume you mean that you
> set Encoding explicitly to FITS-WCS before calling astWrite.
That's right.
> If so, then why does the above list of headers contain FITS-PC keywords
> rather than FITS-WCS? You're asking for FITS-WCS but getting FITS-PC...
> I could understand this behaviour if in fact you were not setting an
> explicit Encoding before calling astWrite, but if you *are* setting
> Encoding, then it looks like a bug to me.
I thought having PC or CD values was now FITS-WCS. "Paper I" says that.
>> so the PROJP1/3 keywords survive. Now when this latter header is read AST
>> says ha, that's FITS-PC, not the FITS-WCS that was written, and uses the
>> PROJP1/3 values and screws up.
>>
>> So what I'd like is to clean the first set of headers of all WCS related
>> cards, and then write the FITS-WCS encoding, which would get rid of the
>> PROJP1/3 cards, plus any other potential conflicts.
>
> Yes, that can be done but I also want to make sure there are no other
> bugs. Also, requiring an application to make an extra call to ensure
> they get the correct encoding is not such a good solution as fixing
> the problem in AST so that you get the correct encoding all the time.
OK, if that's what this turns out to be I agree.
> The other issue is that doing a "spring clean" on other peoples headers
> sounds a bit out of the FITS style.
Still think that's wrong. When I write an updated WCS to the FITS channel
all the other partial or even full encodings become invalid, so what
better than a spring clean first.
Compile and link the attached program and run it on the attached headers:
./astchange before.head after.head FITS-WCS
before.head is:
CTYPE1 = 'RA---ZPN' / Algorithm type for axis 1
CTYPE2 = 'DEC--ZPN' / Algorithm type for axis 2
CRPIX1 = 3.431523900000000E+03 / Dither offset Y
CRPIX2 = -4.111654500000001E+03 / Dither offset Y
CRVAL1 = 2.8992249E+02 / [deg] Right ascension at the reference
pixel
CRVAL2 = 1.5713925E+01 / [deg] Declination at the reference pixel
CRUNIT1 = 'deg ' / Unit of right ascension co-ordinates
CRUNIT2 = 'deg ' / Unit of declination co-ordinates
CD1_1 = 2.0097707E-07 / Transformation matrix element
CD1_2 = 5.5739441E-05 / Transformation matrix element
CD2_1 = -5.5726214E-05 / Transformation matrix element
CD2_2 = 2.7165242E-07 / Transformation matrix element
PV2_1 = 1. / Pol.coeff. for pixel -> celestial coord
PV2_2 = 0.000000E+00 / Pol.coeff. for pixel -> celestial coord
PV2_3 = -50.
PROJP1 = 1. /
PROJP3 = -50. /
and after running it I get after.head:
CTYPE1 = 'DEC--ZPN' / Type of co-ordinate on axis 1
CTYPE2 = 'RA---ZPN' / Type of co-ordinate on axis 2
CRPIX1 = 3431.5239 / Dither offset Y
CRPIX2 = -4111.6545 / Dither offset Y
CRVAL1 = 15.713925 / Value at ref. pixel on axis 1
CRVAL2 = 289.92249 / Value at ref. pixel on axis 2
CRUNIT1 = 'deg ' / Unit of right ascension co-ordinates
CRUNIT2 = 'deg ' / Unit of declination co-ordinates
PROJP1 = 1.0
PROJP3 = -50.0
WCSAXES = 2 / Number of WCS axes
CDELT1 = -5.57268761174625E-5/ Pixel size on axis 1
CDELT2 = 5.57398033255864E-5 / Pixel size on axis 2
PC1_1 = 1.0 / Transformation matrix element
PC1_2 = -0.0048747110716872 / Transformation matrix element
PC2_1 = 0.00360562933500328 / Transformation matrix element
PC2_2 = 1.0 / Transformation matrix element
PV1_1 = 1.0 / Projection parameter
PV1_2 = 0.0 / Projection parameter
PV1_3 = -50.0 / Projection parameter
RADESYS = 'ICRS ' / Reference frame for RA/DEC values
So CD -> PC and encoding is FITS-WCS.
Cheers,
Peter.
#include <stdio.h>
#include "ast.h"
int main( int argc, char *argv[] )
{
FILE *infile;
FILE *outfile;
char *encod;
char line[1025];
AstFitsChan *fitsChan = NULL;
AstFrameSet *wcsinfo;
/* Check command-line */
if ( argc < 4 ) {
fprintf( stderr, "Usage: astchange infile outfile encoding\n" );
return;
}
infile = fopen( argv[1], "r" );
outfile = fopen( argv[2], "w" );
encod = argv[3];
astBegin;
fitsChan = astFitsChan( NULL, NULL, "" );
while( fgets( line, 1024, infile ) != NULL ) {
astPutFits( fitsChan, line, 0 );
}
astClear( fitsChan, "Card" );
printf( "Initial Encoding = %s\n", astGetC( fitsChan, "Encoding" ) );
wcsinfo = astRead( fitsChan );
astSet( fitsChan, "Card=1,Encoding=%s", encod );
printf( "Final Encoding = %s\n", astGetC( fitsChan, "Encoding" ) );
astWrite( fitsChan, wcsinfo );
astClear( fitsChan, "Card" );
while ( astFindFits( fitsChan, "%f", line, 1 ) ) {
fprintf( outfile, "%s\n", line);
}
close( infile );
close( outfile );
astEnd;
}
CTYPE1 = 'RA---ZPN' / Algorithm type for axis 1
CTYPE2 = 'DEC--ZPN' / Algorithm type for axis 2
CRPIX1 = 3.431523900000000E+03 / Dither offset Y
CRPIX2 = -4.111654500000001E+03 / Dither offset Y
CRVAL1 = 2.8992249E+02 / [deg] Right ascension at the reference pixel
CRVAL2 = 1.5713925E+01 / [deg] Declination at the reference pixel
CRUNIT1 = 'deg ' / Unit of right ascension co-ordinates
CRUNIT2 = 'deg ' / Unit of declination co-ordinates
CD1_1 = 2.0097707E-07 / Transformation matrix element
CD1_2 = 5.5739441E-05 / Transformation matrix element
CD2_1 = -5.5726214E-05 / Transformation matrix element
CD2_2 = 2.7165242E-07 / Transformation matrix element
PV2_1 = 1. / Pol.coeff. for pixel -> celestial coord
PV2_2 = 0.000000E+00 / Pol.coeff. for pixel -> celestial coord
PV2_3 = -50.
PROJP1 = 1. /
PROJP3 = -50. /
|