vircam_wcsutils.c

/* $Id: vircam_wcsutils.c,v 1.28 2008/07/10 13:05:53 jim Exp $
 *
 * This file is part of the VIRCAM Pipeline
 * Copyright (C) 2005 Cambridge Astronomy Survey Unit
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program 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 General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 */

/*
 * $Author: jim $
 * $Date: 2008/07/10 13:05:53 $
 * $Revision: 1.28 $
 * $Name:  $
 */

/* Includes */

#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <cpl.h>
#include <math.h>

#include "vircam_wcsutils.h"
#include "vircam_pfits.h"
#include "vircam_utils.h"

#define WCS_FATAL_ERR(_p) {cpl_msg_error(fctid,"Unable to find keyword %s",_p); cpl_error_reset(); return(wcs);}

/* WCS keywords that should be removed from FITS tables*/

#define NNOTABKEYS 6
static const char *notabkeys[NNOTABKEYS] = {"^CRVAL[1-2]*$","^CRPIX[1-2]*",
                                            "^CD[1-2]*_[1-2]*","^CDELT[1-2]*",
                                            "^CTYPE[1-2]*","^PV2_[1-5]*"};

/* Define a list of VIRCAM WCS header keywords */

#define NWCSKEYS 15
static const char *vir_wcskeys[NWCSKEYS] = {"CTYPE1","CTYPE2","CRVAL1",
                                            "CRVAL2","CRPIX1","CRPIX2","CD1_1",
                                            "CD1_2","CD2_1","CD2_2","PV2_1",
                                            "PV2_2","PV2_3","PV2_4","PV2_5"};

#define SZKEY 9

/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

extern void vircam_xytoradec(cpl_wcs *wcs, double x, double y, double *ra,
                             double *dec) {
    double *xy,*radec;
    cpl_matrix *from,*to;
    cpl_array *status;
    int stat;

    /* Load up the information */
    
    from = cpl_matrix_new(1,2);
    xy = cpl_matrix_get_data(from);
    xy[0] = x;
    xy[1] = y;

    /* Call the conversion routine */

    cpl_wcs_convert(wcs,from,&to,&status,CPL_WCS_PHYS2WORLD);

    /* Pass it back now */

    radec = cpl_matrix_get_data(to);
    *ra = radec[0];
    *dec = radec[1];

    /* Tidy and exit */

    cpl_matrix_delete(from);
    cpl_matrix_delete(to);
    cpl_array_delete(status);
    return;
}

/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

extern void vircam_radectoxy(cpl_wcs *wcs, double ra, double dec, double *x,
                             double *y) {
    double *xy,*radec;
    cpl_matrix *from,*to;
    cpl_array *status;
    int stat;

    /* Load up the information */
    
    from = cpl_matrix_new(1,2);
    radec = cpl_matrix_get_data(from);
    radec[0] = ra;
    radec[1] = dec;

    /* Call the conversion routine */

    cpl_wcs_convert(wcs,from,&to,&status,CPL_WCS_WORLD2PHYS);

    /* Pass it back now */

    xy = cpl_matrix_get_data(to);
    *x = xy[0];
    *y = xy[1];

    /* Tidy and exit */

    cpl_matrix_delete(from);
    cpl_matrix_delete(to);
    cpl_array_delete(status);
}

/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

extern void vircam_radectoxieta(cpl_wcs *wcs, double ra, double dec, 
                                double *xi, double *eta) {

    double *xy,*radec;
    cpl_matrix *from,*to;
    cpl_array *status;
    int stat;

    /* Load up the information */
    
    from = cpl_matrix_new(1,2);
    radec = cpl_matrix_get_data(from);
    radec[0] = ra;
    radec[1] = dec;

    /* Call the conversion routine */

    cpl_wcs_convert(wcs,from,&to,&status,CPL_WCS_WORLD2STD);

    /* Pass it back now */

    xy = cpl_matrix_get_data(to);
    *xi = xy[0]/DEGRAD;
    *eta = xy[1]/DEGRAD;

    /* Tidy and exit */

    cpl_matrix_delete(from);
    cpl_matrix_delete(to);
    cpl_array_delete(status);
}

/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

extern int vircam_coverage(cpl_propertylist *plist, int fudge, double *ra1, 
                           double *ra2, double *dec1, double *dec2, 
                           int *status) {

    cpl_wcs *wcs;
    double ra,dec,x,y,dra,ddec,boxfudge,min_4q,max_1q;
    int first_quad,fourth_quad,*naxes;
    long i,j;
    cpl_array *a;

    /* Initialise these in case of failure later*/

    *ra1 = 0.0;
    *ra2 = 0.0;
    *dec1 = 0.0;
    *dec2 = 0.0;
    if (*status != VIR_OK)
        return(*status);

    /* Grab the WCS info from the property list */

    wcs = cpl_wcs_new_from_propertylist(plist);
    if (wcs == NULL) 
        FATAL_ERROR

    /* Get the size of the data array */

    a = cpl_wcs_get_image_dims(wcs);
    naxes = cpl_array_get_data_int(a);

    /* Find the RA and Dec limits of the image */

    *ra1 = 370.0;
    *ra2 = -370.0;
    *dec1 = 95.0;
    *dec2 = -95.0;
    first_quad = 0;
    fourth_quad = 0;
    min_4q = 370.0;
    max_1q = 0.0;
    for (j = 1; j < naxes[1]; j += 10) {
        j = min(j,naxes[1]);
        y = (double)j;
        for (i = 1; i < naxes[0]; i += 10) {
            i = min(i,naxes[0]);
            x = (double)i;
            vircam_xytoradec(wcs,x,y,&ra,&dec);
            if (ra >= 0.0 && ra <= 90.0) {
                first_quad = 1;
                max_1q = max(ra,max_1q);
            } else if (ra >= 270.0 && ra <= 360.0) {
                fourth_quad = 1;
                min_4q = min((ra-360.0),min_4q);
            }
            *ra1 = min(*ra1,ra);
            *ra2 = max(*ra2,ra);
            *dec1 = min(*dec1,dec);
            *dec2 = max(*dec2,dec);
        }
    }
    cpl_wcs_delete(wcs);

    /* Now have a look to see if you had RA values in both the first and
       fourth quadrants.  If you have, then make the minimum RA a negative
       value.  This will be the signal to the caller that you have the
       wraparound... */

    if (first_quad && fourth_quad) {
        *ra1 = min_4q;
        *ra2 = max_1q;
    }

    /* Pad out search a bit */

    if (fudge) {
        boxfudge = 0.01*(float)fudge;
        dra = 0.5*boxfudge*(*ra2 - *ra1);
        *ra1 -= dra;
        *ra2 += dra;
        ddec = 0.5*boxfudge*(*dec2 - *dec1);
        *dec1 -= ddec;
        *dec2 += ddec;
    }

    /* Exit */

    GOOD_STATUS
}

/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

extern int vircam_crpixshift(cpl_propertylist *p, double scalefac, 
                             double sh[]) {
    int i;
    double val;
    char key[SZKEY];
    cpl_type type;
    const char *fctid = "vircam_crpixshift";

    /* Make sure that the scale factor isn't zero */

    if (scalefac == 0.0) {
        cpl_msg_error(fctid,"Scale factor is zero!");
        return(VIR_FATAL);
    }

    /* Loop through both axes. Shift and scale the values of crpix */

    for (i = 1; i <= 2; i++) {
        snprintf(key,SZKEY,"CRPIX%d",i);

        /* First make sure the property exists. It damn well better exist! */

        if (cpl_propertylist_has(p,key) == 0) {
            cpl_msg_error(fctid,"Header is missing WCS key %s",key);
            return(VIR_FATAL);
        }

        /* Now find the type... */

        type = cpl_propertylist_get_type(p,key);

        /* Switch for the type. If it's neither float nor double, then 
           signal an error as this is clearly nonsense. */

        switch (type) {
        case CPL_TYPE_FLOAT:
            val = (double)cpl_propertylist_get_float(p,key);
            break;
        case CPL_TYPE_DOUBLE:
            val = cpl_propertylist_get_double(p,key);
            break;
        default:
            cpl_msg_error(fctid,"Header has WCS key %s as non-floating point!",
                          key);
            return(VIR_FATAL);
        }

        /* Right, now do the actual work */

        val = (val - sh[i-1])/scalefac - 1.0;

        /* And update the property now */

        switch (type) {
        case CPL_TYPE_FLOAT:
            cpl_propertylist_update_float(p,key,(float)val);
            break;
        case CPL_TYPE_DOUBLE:
            cpl_propertylist_update_double(p,key,val);
            break;
        default:
            break;
        }
    }
    return(VIR_OK);
}

/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

extern int vircam_rescalecd(cpl_propertylist *p, double scalefac) {
    int i,j;
    cpl_type type;
    char key[SZKEY];
    const char *fctid = "vircam_rescalecd";
    double val;

    /* Make sure that the scale factor isn't zero */

    if (scalefac == 0.0) {
        cpl_msg_error(fctid,"Scale factor is zero!");
        return(VIR_FATAL);
    }

    /* Loop through both axes. Shift and scale the values of cd */

    for (i = 1; i <= 2; i++) {
        for (j = 1; j <= 2; j++) {
            snprintf(key,SZKEY,"CD%d_%d",i,j);

            /* First make sure the property exists. It damn well better exist! */

            if (cpl_propertylist_has(p,key) == 0) {
                cpl_msg_error(fctid,"Header is missing WCS key %s",key);
                return(VIR_FATAL);
            }

            /* Now find the type... */

            type = cpl_propertylist_get_type(p,key);

            /* Switch for the type. If it's neither float nor double, then 
               signal an error as this is clearly nonsense. */

            switch (type) {
            case CPL_TYPE_FLOAT:
                val = (double)cpl_propertylist_get_float(p,key);
                break;
            case CPL_TYPE_DOUBLE:
                val = cpl_propertylist_get_double(p,key);
                break;
            default:
                cpl_msg_error(fctid,"Header has WCS key %s as non-floating point!",
                              key);
                return(VIR_FATAL);
            }

            /* Right, now do the actual work */

            val *= scalefac;

            /* And update the property now */

            switch (type) {
            case CPL_TYPE_FLOAT:
                cpl_propertylist_update_float(p,key,(float)val);
                break;
            case CPL_TYPE_DOUBLE:
                cpl_propertylist_update_double(p,key,val);
                break;
            default:
                break;
            }
        }
    }
    return(VIR_OK);
}
    
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

extern int vircam_diffxywcs(cpl_wcs *wcs, cpl_wcs *wcsref, float *xoff,
                            float *yoff, int *status) {
    double xc,yc,ra,dec,xnew,ynew;
    cpl_array *a;
    int *dims;
    const char *fctid = "vircam_diffxywcs";

    /* Inherited status */

    *xoff = 0.0;
    *yoff = 0.0;
    if (*status != VIR_OK) 
        return(*status);

    /* Check for rubbish input */

    if (wcs == NULL || wcsref == NULL) {
        cpl_msg_error(fctid,"NULL wcs information");
        FATAL_ERROR
    }

    /* Work out the ra and dec of the central pixel in the reference image */

    a = cpl_wcs_get_image_dims(wcsref);
    dims = cpl_array_get_data_int(a);
    xc = 0.5*(double)dims[0];
    yc = 0.5*(double)dims[1];
    vircam_xytoradec(wcsref,xc,yc,&ra,&dec);

    /* Ok, calculate where in x,y space these equatorial coordinates fall on
       the programme image */

    vircam_radectoxy(wcs,ra,dec,&xnew,&ynew);

    /* Right, now return the offsets */

    *xoff = (float)(xc - xnew);
    *yoff = (float)(yc - ynew);
    GOOD_STATUS
}

/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

extern int vircam_removewcs(cpl_propertylist *p, int *status) {
    int i;
    const char *fctid = "vircam_removewcs";

    /* Inherited status */

    if (*status != VIR_OK)
        return(*status);
    if (p == NULL) {
        cpl_msg_error(fctid,"Propertylist passed is NULL\nProgramming error");
        FATAL_ERROR
    }

    /* Loop through all the template keywords and remove them */

    for (i = 0; i < NNOTABKEYS; i++) 
        cpl_propertylist_erase_regexp(p,notabkeys[i],0);

    GOOD_STATUS
}

/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

extern int vircam_tabwcs(cpl_propertylist *p, int xcol, int ycol, 
                         int *status) {
    int i;
    char key[9],key2[9],*fctid="vircam_tabwcs";

    /* Inherited status */

    if (*status != VIR_OK)
        return(*status);
    if (p == NULL) {
        cpl_msg_error(fctid,"Propertylist passed is NULL\nProgramming error");
        FATAL_ERROR
    }

    /* If either of the columns are -1, then just get rid of the image WCS 
       and get out of here */

    if (xcol == -1 || ycol == -1) {
        vircam_removewcs(p,status);
        GOOD_STATUS
    }

    /* Go through the standard VIRCAM WCS header keywords one by one
       and translate them.  Start with CTYPE */

    (void)snprintf(key,8,"TCTYP%d",xcol);
    vircam_rename_property(p,"CTYPE1",key);
    (void)snprintf(key,8,"TCTYP%d",ycol);
    vircam_rename_property(p,"CTYPE2",key);


    /* Now CRVAL */

    (void)snprintf(key,8,"TCRVL%d",xcol);
    vircam_rename_property(p,"CRVAL1",key);
    (void)snprintf(key,8,"TCRVL%d",ycol);
    vircam_rename_property(p,"CRVAL2",key);

    /* Now CRPIX */

    (void)snprintf(key,8,"TCRPX%d",xcol);
    vircam_rename_property(p,"CRPIX1",key);
    (void)snprintf(key,8,"TCRPX%d",ycol);
    vircam_rename_property(p,"CRPIX2",key);

    /* Now PV matrix */

    for (i = 1; i <= 5; i++) {
        (void)snprintf(key2,8,"PV2_%d",i);
        (void)snprintf(key,8,"TV%d_%d",ycol,i);
        if (cpl_propertylist_has(p,key2))
            vircam_rename_property(p,key2,key);
    }

    /* Now the CD matrix */

    (void)snprintf(key,8,"TC%d_%d",xcol,xcol);
    vircam_rename_property(p,"CD1_1",key);
    (void)snprintf(key,8,"TC%d_%d",xcol,ycol);
    vircam_rename_property(p,"CD1_2",key);
    (void)snprintf(key,8,"TC%d_%d",ycol,xcol);
    vircam_rename_property(p,"CD2_1",key);
    (void)snprintf(key,8,"TC%d_%d",ycol,ycol);
    vircam_rename_property(p,"CD2_2",key);

    /* Now get out of here */

    GOOD_STATUS
    
}


/*

$Log: vircam_wcsutils.c,v $
Revision 1.28  2008/07/10 13:05:53  jim
Modified to use v4.2 version of cpl_wcs

Revision 1.27  2008/06/20 11:14:01  jim
Fixed dodgy call to cpl_wcs_get_image_dims

Revision 1.26  2008/05/06 08:40:10  jim
Modified to use cpl_wcs interface

Revision 1.25  2007/10/25 17:34:01  jim
Modified to remove lint warnings

Revision 1.24  2007/10/15 12:50:28  jim
Modified for compatibility with cpl_4.0

Revision 1.23  2007/05/14 15:25:28  jim
Fixed formatting of wcs floating point parameters in vircam_wcsopen

Revision 1.22  2007/05/03 11:15:33  jim
Fixed little problem with table wcs

Revision 1.21  2007/05/02 12:53:11  jim
typo fixes in docs

Revision 1.20  2007/05/02 09:16:12  jim
Added vircam_tabwcs

Revision 1.19  2007/03/01 12:42:42  jim
Modified slightly after code checking

Revision 1.18  2007/01/17 23:54:01  jim
Plugged some memory leaks

Revision 1.17  2007/01/08 19:12:49  jim
Fixed stupid calloc

Revision 1.16  2006/10/05 09:23:00  jim
Small modifications to a couple of cpl calls to bring them into line with
cpl v3.0

Revision 1.15  2006/08/11 12:45:41  jim
Modified to use wcslib

Revision 1.14  2006/07/04 09:19:06  jim
replaced all sprintf statements with snprintf

Revision 1.13  2006/07/03 09:33:19  jim
Fixed a few things to keep the compiler happy

Revision 1.12  2006/06/14 14:13:58  jim
fixed minor doc problem

Revision 1.11  2006/06/13 14:09:54  jim
Added vircam_removewcs

Revision 1.10  2006/05/24 13:36:15  jim
Added vircam_diffxywcs

Revision 1.9  2006/04/20 11:30:29  jim
Added naxis1 and naxis2 to the structure

Revision 1.8  2006/03/22 13:58:32  jim
Cosmetic fixes to keep lint happy

Revision 1.7  2006/02/18 11:48:55  jim
*** empty log message ***

Revision 1.6  2006/01/23 10:30:50  jim
Mainly documentation mods

Revision 1.5  2005/12/14 22:17:33  jim
Updated docs

Revision 1.4  2005/11/29 11:37:38  jim
fixed bug in error macro

Revision 1.3  2005/11/25 15:33:22  jim
Some code fixes to keep splint happy

Revision 1.2  2005/11/03 13:28:51  jim
All sorts of changes to tighten up error handling

Revision 1.1.1.1  2005/08/05 08:29:09  jim
Initial import


*/

---