/*** File libwcs/fitswcs.c *** March 14, 2011 *** By Jessica Mink, jmink@cfa.harvard.edu *** Harvard-Smithsonian Center for Astrophysics *** Copyright (C) 1996-2011 *** Smithsonian Astrophysical Observatory, Cambridge, MA, USA This library is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version. This library is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details. You should have received a copy of the GNU Lesser General Public License along with this library; if not, write to the Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA Correspondence concerning WCSTools should be addressed as follows: Internet email: jmink@cfa.harvard.edu Postal address: Jessica Mink Smithsonian Astrophysical Observatory 60 Garden St. Cambridge, MA 02138 USA * Module: fitswcs.c (FITS file WCS reading and deleting) * Purpose: Read and delete FITS image world coordinate system keywords * * Subroutine: GetWCSFITS (filename, verbose) * Open a FITS or IRAF image file and returns its WCS structure * Subroutine: GetFITShead (filename, verbose) * Open a FITS or IRAF image file and returns a FITS header * Subroutine: DelWCSFITS (header, verbose) * Delete all standard WCS keywords from a FITS header * Subroutine: PrintWCS (header, verbose) * Check the WCS fields and print any that are found if verbose. * Subroutine: SetFITSWCS (header, wcs) * Set FITS WCS keywords from WCS data structure */ #include #include #include #include "fitsfile.h" #include "wcs.h" struct WorldCoor * GetWCSFITS (filename, verbose) char *filename; /* FITS or IRAF file filename */ int verbose; /* Print extra information if nonzero */ { char *header; /* FITS header */ struct WorldCoor *wcs; /* World coordinate system structure */ char *GetFITShead(); char *cwcs; /* Multiple wcs string (name or character) */ /* Read the FITS or IRAF image file header */ header = GetFITShead (filename, verbose); if (header == NULL) return (NULL); /* Set the world coordinate system from the image header */ cwcs = strchr (filename, '%'); if (cwcs != NULL) cwcs++; wcs = wcsinitn (header, cwcs); if (wcs == NULL) { setwcsfile (filename); if (verbose) wcserr (); } free (header); return (wcs); } char * GetFITShead (filename, verbose) char *filename; /* FITS or IRAF file filename */ int verbose; /* Print error messages if nonzero */ { char *header; /* FITS header */ int lhead; /* Maximum number of bytes in FITS header */ char *irafheader; /* IRAF image header */ int nbiraf, nbfits; /* Open IRAF image if .imh extension is present */ if (isiraf (filename)) { if ((irafheader = irafrhead (filename, &nbiraf)) != NULL) { if ((header = iraf2fits (filename, irafheader, nbiraf, &lhead)) == NULL) { if (verbose) fprintf (stderr, "Cannot translate IRAF header %s\n",filename); free (irafheader); irafheader = NULL; return (NULL); } free (irafheader); irafheader = NULL; } else { if (verbose) fprintf (stderr, "Cannot read IRAF header file %s\n", filename); return (NULL); } } else if (istiff (filename) || isgif (filename) || isjpeg (filename)) { if ((header = fitsrtail (filename, &lhead, &nbfits)) == NULL) { if (verbose) fprintf (stderr, "TIFF file %s has no appended header\n", filename); return (NULL); } } /* Open FITS file if .imh extension is not present */ else { if ((header = fitsrhead (filename, &lhead, &nbfits)) == NULL) { if (verbose) /* fprintf (stderr, "Cannot read FITS file %s\n", filename); */ fitserr (); return (NULL); } } return (header); } /* delete all the C* fields. * return 0 if at least one such field is found, else -1. */ int DelWCSFITS (header, verbose) char *header; int verbose; { static char flds[19][8]; char keyword[8]; int i; int n, nfields; double eq; char rastr[32],decstr[32]; n = 0; /* Delete standard WCS fields */ nfields = 19; strcpy (flds[0], "CTYPE1"); strcpy (flds[1], "CTYPE2"); strcpy (flds[2], "CRVAL1"); strcpy (flds[3], "CRVAL2"); strcpy (flds[4], "CDELT1"); strcpy (flds[5], "CDELT2"); strcpy (flds[6], "CRPIX1"); strcpy (flds[7], "CRPIX2"); strcpy (flds[8], "CROTA1"); strcpy (flds[9], "CROTA2"); strcpy (flds[10], "IMWCS"); strcpy (flds[11], "CD1_1"); strcpy (flds[12], "CD1_2"); strcpy (flds[13], "CD2_1"); strcpy (flds[14], "CD2_2"); strcpy (flds[15], "PC1_1"); strcpy (flds[16], "PC1_2"); strcpy (flds[17], "PC2_1"); strcpy (flds[18], "PC2_2"); for (i = 0; i < nfields; i++) { if (hdel (header, flds[i])) { n++; if (verbose) fprintf (stderr,"%s: deleted\n", flds[i]); } } /* Delete projection parameters */ for (i = 0; i < MAXPV; i++) { sprintf (keyword, "PV1_%d", i); if (hdel (header, keyword)) { n++; if (verbose) fprintf (stderr,"%s: deleted\n", keyword); } sprintf (keyword, "PV2_%d", i); if (hdel (header, keyword)) { n++; if (verbose) fprintf (stderr,"%s: deleted\n", keyword); } } /* Delete rotation matrix, if present */ if (ksearch (header, "PC001001")) { int i, j; char keyword[16]; for (i = 1; i < 6; i++) { for (j = 1; j < 6; j++) { sprintf (keyword,"PC%03d%03d", i, j); if (hdel (header, keyword)) { if (verbose) fprintf (stderr,"%s deleted\n", keyword); n++; } } } } if (verbose && n == 0) fprintf (stderr,"DelWCSFITS: No WCS in header\n"); /* Delete RA DEC EPOCH, replacing with saved values, if present */ if (ksearch (header,"WRA")) { hdel (header, "RA"); n++; hchange (header, "WRA","RA"); if (ksearch (header,"WDEC")) { hdel (header, "DEC"); n++; hchange (header, "WDEC", "DEC"); } eq = 0.0; if (!hgetr8 (header, "WEQUINOX", &eq)) { if (!hgetr8 (header, "WEPOCH", &eq)) { if (!hgetr8 (header, "EQUINOX", &eq)) { hgetr8 (header, "EPOCH", &eq); } } } if (ksearch (header,"WEQUINOX")) { hdel (header, "EPOCH"); n++; hchange (header, "WEQUINOX", "EPOCH"); } else if (ksearch (header,"WEPOCH")) { hdel (header, "EPOCH"); n++; hchange (header, "WEPOCH", "EPOCH"); } else if (eq > 0.0) { hputr8 (header, "EPOCH", eq); } if (ksearch (header, "EQUINOX")) { hdel (header, "EQUINOX"); n++; if (verbose) fprintf (stderr,"EQUINOX deleted\n"); } if (hdel (header, "RADECSYS")) { if (verbose) fprintf (stderr,"RADECSYS deleted\n"); n++; } if (hdel (header, "SECPIX1")) { if (verbose) fprintf (stderr,"SECPIX1 deleted\n"); n++; } if (hdel (header, "SECPIX2")) { if (verbose) fprintf (stderr,"SECPIX2 deleted\n"); n++; } if (verbose) { hgets (header,"RA", 31, rastr); hgets (header,"DEC", 31, decstr); eq = 0.0; hgetr8 (header,"EPOCH",&eq); if (eq == 0.0) hgetr8 (header,"EQUINOX",&eq); fprintf (stderr,"DelWCS: Center reset to %s %s %.1f\n", rastr,decstr, eq); } } else if (ksearch (header, "EPOCH") && !ksearch (header, "PLTRAH")) { if (hdel (header,"EQUINOX")) { if (verbose) fprintf (stderr,"EQUINOX: deleted\n"); n++; } else if (verbose) fprintf (stderr,"DelWCS: EPOCH, but not EQUINOX found\n"); } /* Delete IMWCS result keywords, if present */ if (ksearch (header, "WCSMATCH")) { hdel (header, "WCSMATCH"); hdel (header, "WCSRFCAT"); hdel (header, "WCSIMCAT"); hdel (header, "WCSNREF"); hdel (header, "WCSTOL"); hdel (header, "WCSSEP"); } /* Delete SAO polynomial, if present */ if (ksearch (header, "CO1_1")) { int i; char keyword[16]; for (i = 1; i < 13; i++) { sprintf (keyword,"CO1_%d", i); if (hdel (header, keyword)) { if (verbose) fprintf (stderr,"%s deleted\n", keyword); n++; } } for (i = 1; i < 13; i++) { sprintf (keyword,"CO2_%d", i); if (hdel (header, keyword)) { if (verbose) fprintf (stderr,"%s deleted\n", keyword); n++; } } } if (n > 0 && verbose) fprintf (stderr,"%d keywords deleted\n", n); return (n); } /* check the WCS fields and print any that are found if verbose. * return 0 if all are found, else -1. */ int PrintWCS (header, verbose) char *header; /* FITS header */ int verbose; /* 1 to print WCS header keyword values */ { char str[80]; double v; int n, i; char keyword[16]; n = 0; if (hgets (header,"IMWCS",80,str)) { if (verbose) fprintf (stderr,"IMWCS = %s\n", str); n++; } if (hgets (header,"CTYPE1",16,str)) { if (verbose) fprintf (stderr,"CTYPE1 = %s\n", str); n++; } if (hgetr8 (header, "CRVAL1", &v)) { if (verbose) fprintf (stderr,"CRVAL1 = %.8f\n", v); n++; } if (hgetr8 (header, "CRPIX1", &v)) { if (verbose) fprintf (stderr,"CRPIX1 = %.8f\n", v); n++; } if (hgets (header,"CTYPE2",16,str)) { if (verbose) fprintf (stderr,"CTYPE2 = %s\n", str); n++; } if (hgetr8 (header, "CRVAL2", &v)) { if (verbose) fprintf (stderr,"CRVAL2 = %.8f\n", v); n++; } if (hgetr8 (header, "CRPIX2", &v)) { if (verbose) fprintf (stderr,"CRPIX2 = %.8f\n", v); n++; } /* Polynomial plate fit */ if (hgetr8 (header, "CO1_1", &v)) { if (verbose) fprintf (stderr,"CO1_1 = %.8g\n", v); for (i = 1; i < 20; i++) { sprintf (keyword,"CO1_%d",i+1); if (hgetr8 (header, keyword, &v)) { if (verbose) fprintf (stderr,"%s = %.8g\n", keyword, v); n++; } } } if (hgetr8 (header, "CO2_1", &v)) { if (verbose) fprintf (stderr,"CO2_1 = %.8g\n", v); for (i = 1; i < 20; i++) { sprintf (keyword,"CO2_%d",i+1); if (hgetr8 (header, keyword, &v)) { if (verbose) fprintf (stderr,"%s = %.8g\n", keyword, v); n++; } } } /* Plate scale and rotation from CD matrix */ if (hgetr8 (header, "CD1_1", &v)) { if (verbose) fprintf (stderr,"CD1_1 = %.8g\n", v); n++; if (hgetr8 (header, "CD1_2", &v)) { if (verbose) fprintf (stderr,"CD1_2 = %.8g\n", v); n++; } if (hgetr8 (header, "CD2_1", &v)) { if (verbose) fprintf (stderr,"CD2_1 = %.8g\n", v); n++; } if (hgetr8 (header, "CD2_2", &v)) { if (verbose) fprintf (stderr,"CD2_2 = %.8g\n", v); n++; } } /* Plate scale and rotation from CDELTn and CROTAn */ else { if (hgetr8 (header, "CDELT1", &v)) { if (verbose) fprintf (stderr,"CDELT1 = %.8f\n", v); n++; } if (hgetr8 (header, "CROTA1", &v)) { if (verbose) fprintf (stderr,"CROTA1 = %.3f\n", v); n++; } if (hgetr8 (header, "CDELT2", &v)) { if (verbose) fprintf (stderr,"CDELT2 = %.8f\n", v); n++; } if (hgetr8 (header, "CROTA2", &v)) { if (verbose) fprintf (stderr,"CROTA2 = %.3f\n", v); n++; } } return (n > 8 ? 0 : -1); } static char wcsproj[8]="TAN"; /* WCS projection name */ void setwcsproj (type) char *type; { strcpy (wcsproj, type); return; } /* Set FITS C* fields, assuming RA/DEC refers to the reference pixel, CRPIX1/CRPIX2 */ void SetFITSWCS (header, wcs) char *header; /* Image FITS header */ struct WorldCoor *wcs; /* WCS structure */ { double ep; char wcstemp[16]; char *wcsdist; /* Rename old center coordinates */ if (!ksearch (header,"WRA") && ksearch (header,"RA")) hchange (header,"RA","WRA"); if (!ksearch (header,"WDEC") && ksearch (header,"DEC")) hchange (header,"DEC","WDEC"); if (!ksearch (header,"WEQUINOX") && ksearch (header,"EQUINOX")) hchange (header, "EQUINOX", "WEQUINOX"); /* Only change EPOCH if it is used instead of EQUINOX */ else if (!ksearch (header,"WEPOCH") && ksearch (header,"EPOCH")) hchange (header, "EPOCH", "WEPOCH"); /* Set new center coordinates */ if (wcs->xref < 0) wcs->xref = 360.0 + wcs->xref; hputra (header,"RA",wcs->xref); hputdec (header,"DEC",wcs->yref); hputr8 (header, "EQUINOX", wcs->equinox); if (hgetr8 (header, "WEPOCH", &ep)) hputr8 (header, "EPOCH", wcs->equinox); else if (!hgetr8 (header, "EPOCH", &ep)) hputr8 (header, "EPOCH", wcs->equinox); if (wcs->radecsys[0] == 'B' || wcs->radecsys[0] == 'b') hputs (header, "RADECSYS", "FK4"); else if (wcs->radecsys[0] == 'I' || wcs->radecsys[0] == 'i') hputs (header, "RADECSYS", "IRCS"); else if (wcs->radecsys[0] == 'J' || wcs->radecsys[0] == 'j') hputs (header, "RADECSYS", "FK5"); else hputs (header, "RADECSYS", wcs->radecsys); /* Set standard FITS WCS keywords */ wcsdist = getdistcode (wcs); /* FITS WCS distortion code */ strcpy (wcstemp, "RA---"); strcat (wcstemp, wcsproj); if (wcsdist != NULL) strcat (wcstemp, wcsdist); hputs (header, "CTYPE1", wcstemp); strcpy (wcstemp, "DEC--"); strcat (wcstemp, wcsproj); if (wcsdist != NULL) strcat (wcstemp, wcsdist); hputs (header, "CTYPE2", wcstemp); /* Reference pixel in WCS and image coordinates */ hputnr8 (header, "CRVAL1", 9, wcs->xref); hputnr8 (header, "CRVAL2", 9, wcs->yref); hputnr8 (header, "CRPIX1", 4, wcs->xrefpix); hputnr8 (header, "CRPIX2", 4, wcs->yrefpix); /* CD matrix (proposed FITS standard) */ if (wcs->rotmat) { hputnr8 (header, "CD1_1", 12, wcs->cd[0]); hputnr8 (header, "CD1_2", 12, wcs->cd[1]); hputnr8 (header, "CD2_1", 12, wcs->cd[2]); hputnr8 (header, "CD2_2", 12, wcs->cd[3]); hdel (header, "CDELT1"); hdel (header, "CDELT2"); hdel (header, "CROTA1"); hdel (header, "CROTA2"); } /* Scale and rotation (old FITS standard) */ else { hputnr8 (header, "CDELT1", 12, wcs->xinc); hputnr8 (header, "CDELT2", 12, wcs->yinc); hputnr8 (header, "CROTA1", 6, wcs->rot); hputnr8 (header, "CROTA2", 6, wcs->rot); hputnr8 (header, "CD1_1", 12, wcs->cd[0]); hputnr8 (header, "CD1_2", 12, wcs->cd[1]); hputnr8 (header, "CD2_1", 12, wcs->cd[2]); hputnr8 (header, "CD2_2", 12, wcs->cd[3]); /* hdel (header, "CD1_1"); hdel (header, "CD1_2"); hdel (header, "CD2_1"); hdel (header, "CD2_2"); */ } /* Plate scale at reference pixel */ if (-wcs->xinc != wcs->yinc) { if (ksearch (header,"SECPIX")) hdel (header,"SECPIX"); hputnr8 (header, "SECPIX1", 4, -wcs->xinc*3600.0); hputnr8 (header, "SECPIX2", 4, wcs->yinc*3600.0); } else { if (ksearch (header,"SECPIX1")) hdel (header,"SECPIX1"); if (ksearch (header,"SECPIX2")) hdel (header,"SECPIX2"); hputnr8 (header, "SECPIX", 6, wcs->yinc*3600.0); } /* Plate fit coefficients, if present */ if (wcs->ncoeff1 > 0) { char keyword[16]; int i; for (i = 0; i < wcs->ncoeff1; i++) { sprintf (keyword, "CO1_%d",i+1); hputr8 (header, keyword, wcs->x_coeff[i]); } } if (wcs->ncoeff2 > 0) { char keyword[16]; int i; for (i = 0; i < wcs->ncoeff2; i++) { sprintf (keyword, "CO2_%d",i+1); hputr8 (header, keyword, wcs->y_coeff[i]); } } return; } /* May 29 1996 Change name from delWCSFITS to DelWCSFITS * May 31 1996 Print single message if no WCS is found in header * May 31 1996 Use stream I/O instead of standard I/O * Jun 10 1996 Combine imgetwcs.c and imdelwcs.c into fitswcs.c * Jun 17 1996 Delete IMWCS record, too * Jul 16 1996 Update arguments for header-reading subroutines * Aug 6 1996 Fixed small defects after lint * Aug 8 1996 Restore old image center after deleting WCS * Aug 26 1996 Fix subroutine arguments after lint * * Feb 21 1997 Check pointers against NULL explicitly for Linux * Feb 21 1997 Add GetFITShead() subroutine and use it * Mar 20 1997 Remove unused variables * Nov 6 1997 Add PrintWCS() from IMWCS * * Jan 7 1998 Return NULL WCS structure if no FITS header can be read * Feb 18 1998 Move SetFITSWCS() here from imsetwcs.c * Feb 24 1998 Delete CD matrix in DelWCS() * Mar 20 1998 Write CD matrix in SetFITSWCS() * Mar 27 1998 Add plate constants in SetFITSWCS() * Mar 27 1998 Delete plate constants in DelFITSWCS() * Apr 6 1998 Change coefficient keywords from PLTij to COi_j * Apr 7 1998 Change amd_i_coeff to i_coeff * Apr 10 1998 Write out polynomial coefficients correctly * Apr 13 1998 Print polynomial coefficients, if in header * Apr 16 1998 Drop NCOEFF header parameter * Apr 17 1998 Do not write W* keywords if they are already there * May 27 1998 Include fitsio.h instead of fitshead.h * Jun 1 1998 Print error message if WCS cannot be initialized * Jun 11 1998 Change WCSTYPE to WCSPROJ to avoid conflict * Jul 23 1998 In DelWCS, delete specific number of fields * Jul 24 1998 Make irafheader char instead of int * Jul 27 1998 Set irafheader pointer to NULL after use * Aug 6 1998 Change fitsio.h to fitsfile.h * Oct 5 1998 Use isiraf() to determine file type * Oct 28 1998 Delete EQUINOX, RADECSYS, SECPIX1, SECPIX2 from imwcs * * Apr 7 1999 Add file name to error message if WCS error * Jul 8 1999 Write RADECSYS as FK5 or FK4 instead of J2000 or B1950 * Jul 21 1999 Add SECPIX plate scale output to SetFITSWCS() * Oct 21 1999 Fix declarations after lint * * Mar 23 2000 Fix bug in IRAF header error message * * Jan 11 2001 Print all messages to stderr * Jan 31 2001 Add code to extract WCS name or character from filename * Mar 8 2001 Change WCS character separator from : to % in FITS filenames * Jul 11 2001 Add PC matrix to keywords deleted by DelWCS() * * Apr 23 2002 Always write CD matrix in SetFITSWCS() * Jun 19 2002 Add verbose argument to GetWCSFITS() and GetFITShead() * * Feb 10 2003 Print 12 decimal places instead of 9 for CD matrix and CDELT * Oct 23 2003 Add PCi_j to DelWCSFITS() * Nov 3 2003 In SetFITSWCS(), add distortion code if in WCS * Dec 5 2003 Fix bug, delete projection parameters in DelWCSFITS() * * Jul 19 2004 Print error message in verbose mode only * Sep 16 2004 Add 360.0 to negative right ascensions in SetFITSWCS() * Nov 1 2005 Set RADECSYS to ICRS if appropriate * * Feb 23 2006 Add code to read FITS header appended to TIFF file * Jun 1 2006 Fix bug so CD matrix is deleted by DelWCSFITS() * Sep 26 2006 Increase length of rastr and destr from 16 to 32 * * Apr 2 2007 Fix DelWCSFITS() argument description at top of file * Apr 18 2007 Delete WCS result keywords with DelWCSFITS() * Dec 20 2007 Print error message set by fitsrhead() * * Apr 22 2010 Fix parameter deletion errors found by Paul Liptack * Apr 22 2010 Only print deletion confirmation if deletion occurred * Apr 29 2010 Fx loop index test bug i -> j in PC matrix deletion * * Mar 14 2011 Delete up to j=MAXPV PVi_j parameters (for SCAMP from Ed Los) */