vircam_getstds.c

00001 /* $Id: vircam_getstds.c,v 1.21 2008/05/06 08:40:10 jim Exp $
00002  *
00003  * This file is part of the VIRCAM Pipeline
00004  * Copyright (C) 2005 Cambridge Astronomy Survey Unit
00005  *
00006  * This program is free software; you can redistribute it and/or modify
00007  * it under the terms of the GNU General Public License as published by
00008  * the Free Software Foundation; either version 2 of the License, or
00009  * (at your option) any later version.
00010  *
00011  * This program is distributed in the hope that it will be useful,
00012  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00013  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00014  * GNU General Public License for more details.
00015  *
00016  * You should have received a copy of the GNU General Public License
00017  * along with this program; if not, write to the Free Software
00018  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00019  */
00020 
00021 /*
00022  * $Author: jim $
00023  * $Date: 2008/05/06 08:40:10 $
00024  * $Revision: 1.21 $
00025  * $Name:  $
00026  */
00027 
00028 /* Includes */
00029 
00030 #ifdef HAVE_CONFIG_H
00031 #include <config.h>
00032 #endif
00033 
00034 #include <sys/stat.h>
00035 #include <sys/types.h>
00036 #include <unistd.h>
00037 #include <string.h>
00038 #include <cpl.h>
00039 
00040 #include <math.h>
00041 
00042 #include "vircam_mods.h"
00043 #include "vircam_utils.h"
00044 #include "vircam_wcsutils.h"
00045 
00046 #define CACHEDIR "catcache"
00047 #define CACHEIND "catcache/index"
00048 #define SZBUF 1024
00049 
00050 static cpl_table *vircam_2mass_extract(char *path, float ramin, float ramax,
00051                                        float decmin, float decmax);
00052 static cpl_table *check_cache(char *catname, float ra1_im, float ra2_im, 
00053                               float dec1_im, float dec2_im);
00054 static void addto_cache(cpl_table *stds, char *catname, float ramin, 
00055                         float ramax, float decmin, float decmax);
00056 
00059 /*---------------------------------------------------------------------------*/
00108 /*---------------------------------------------------------------------------*/
00109 
00110 extern int vircam_getstds(cpl_propertylist *plist, int cache, char *path,
00111                           char *catname, cpl_table **stds, int *status) {
00112     const char *fctid = "vircam_getstds";
00113     double x1,x2,y1,y2,r,d,xx,yy;
00114     float ramin,ramax,decmin,decmax,*ra,*dec,*x,*y;
00115     cpl_wcs *wcs;
00116     int n,i;
00117     cpl_propertylist *p;
00118 
00119     /* Inherited status */
00120 
00121     *stds = NULL;
00122     if (*status != VIR_OK)
00123         return(*status);
00124 
00125     /* Get the coverage of the input WCS */
00126 
00127     (void)vircam_coverage(plist,0,&x1,&x2,&y1,&y2,status);
00128     ramin = (float)x1;
00129     ramax = (float)x2;
00130     decmin = (float)y1;
00131     decmax = (float)y2;
00132 
00133     /* If using the cache, then check if for a previously used table */
00134 
00135     *stds = NULL;
00136     if (cache)
00137         *stds = check_cache(catname,ramin,ramax,decmin,decmax);
00138 
00139     /* If there was nothing in the cache (or you didn't use it) then search
00140        the standard 2mass catalogues */
00141 
00142     if (*stds == NULL) {        
00143 
00144         /* Read the standards from the catalogue */
00145 
00146         *stds = vircam_2mass_extract(path,ramin,ramax,decmin,decmax);
00147         if (*stds == NULL) {
00148             cpl_msg_error(fctid,"Unable to extract data in %s\n",path);
00149             FATAL_ERROR
00150         }
00151         
00152         /* Add this table to the cache if you want to */
00153 
00154         if (cache)
00155             addto_cache(*stds,catname,ramin,ramax,decmin,decmax);
00156     }
00157 
00158     /* If there are no rows in the table, this may be a cause for concern.
00159        So add the columns into the table and then set a warning return 
00160        status */
00161 
00162     n = cpl_table_get_nrow(*stds);
00163     if (n == 0) {
00164         cpl_table_new_column(*stds,"xpredict",CPL_TYPE_FLOAT);
00165         cpl_table_new_column(*stds,"ypredict",CPL_TYPE_FLOAT);
00166         WARN_RETURN
00167     }
00168 
00169     /* Now fill the coordinates in */
00170 
00171     wcs = cpl_wcs_new_from_propertylist((const cpl_propertylist *)plist);
00172     ra = cpl_table_get_data_float(*stds,"RA");
00173     dec = cpl_table_get_data_float(*stds,"Dec");
00174     x = cpl_malloc(n*sizeof(*x));
00175     y = cpl_malloc(n*sizeof(*y));
00176     for (i = 0; i < n; i++) {
00177         r = (double)ra[i];
00178         d = (double)dec[i];
00179         vircam_radectoxy(wcs,r,d,&xx,&yy);
00180         x[i] = (float)xx;
00181         y[i] = (float)yy;
00182     }
00183     cpl_wcs_delete(wcs);
00184     
00185     /* Add the predicted x,y coordinates columns to the table */
00186 
00187     cpl_table_wrap_float(*stds,x,"xpredict");
00188     cpl_table_set_column_unit(*stds,"xpredict","pixels");
00189     cpl_table_wrap_float(*stds,y,"ypredict");
00190     cpl_table_set_column_unit(*stds,"ypredict","pixels");
00191 
00192     /* Finally sort this by ypredict */
00193 
00194     p = cpl_propertylist_new();
00195     cpl_propertylist_append_bool(p,"ypredict",0);
00196     cpl_table_sort(*stds,p);
00197     cpl_propertylist_delete(p);
00198 
00199     /* Get out of here */
00200 
00201     GOOD_STATUS
00202 }
00203 
00204 
00205 /*---------------------------------------------------------------------------*/
00233 /*---------------------------------------------------------------------------*/
00234 
00235 static cpl_table *vircam_2mass_extract(char *path, float ramin1, float ramax1, 
00236                                        float decmin, float decmax) {
00237     cpl_table *t,*s,*o;
00238     int i,nrows,start,finish,first_index,last_index,irow,init,j;
00239     int first_index_ra,last_index_ra,wrap,iwrap;
00240     float dectest,ratest,ramin,ramax;
00241     char fullname[SZBUF];
00242     cpl_array *a;
00243     char *deccol[] = {"Dec"};
00244     cpl_propertylist *p;
00245 
00246     /* Create an output table */
00247 
00248     o = cpl_table_new(0);
00249     init = 1;
00250 
00251     /* Create a cpl array */
00252 
00253     a = cpl_array_wrap_string(deccol,1);
00254 
00255     /* Is there a wrap around problem? */
00256 
00257     wrap = (ramin1 < 0.0 && ramax1 > 0.0) ? 2 : 1;
00258 
00259     /* Loop for each query. If there is a wrap around problem then we need 2
00260        queries. If not, then we only need 1 */
00261 
00262     for (iwrap = 0; iwrap < wrap; iwrap++) {
00263         if (wrap == 2) {
00264             if (iwrap == 0) {
00265                 ramin = ramin1 + 360.0;
00266                 ramax = 360.0;
00267             } else {
00268                 ramin = 0.000001;
00269                 ramax = ramax1;
00270             }
00271         } else {
00272             ramin = ramin1;
00273             ramax = ramax1;
00274         }
00275 
00276         /* Find out where in the index to look */
00277 
00278         first_index_ra = (int)ramin;
00279         last_index_ra = min((int)ramax,359);
00280         
00281         /* Look at the min and max RA and decide which files need to be 
00282            opened. */
00283 
00284         for (i = first_index_ra; i <= last_index_ra; i++) {
00285 
00286             /* Ok, we've found one that needs opening. Read the file with 
00287                the relevant CPL call */
00288 
00289             (void)snprintf(fullname,SZBUF,"%s/npsc%03d.fits",path,i);
00290 
00291             /* Read the propertylist so that you know how many rows there
00292                are in the table */
00293 
00294             p = cpl_propertylist_load(fullname,1);
00295             if (p == NULL) {
00296                 freetable(o);
00297                 cpl_array_unwrap(a);
00298                 return(NULL);
00299             }
00300             nrows = cpl_propertylist_get_int(p,"NAXIS2");
00301             cpl_propertylist_delete(p);
00302 
00303             /* Load various rows until you find the Dec range that you 
00304                have specified. First the minimum Dec */
00305 
00306             start = 0;
00307             finish = nrows;
00308             first_index = nrows/2;
00309             while (finish - start >= 2) {
00310                 t = cpl_table_load_window(fullname,1,0,a,first_index,1);
00311                 dectest = cpl_table_get_float(t,"Dec",0,NULL);
00312                 cpl_table_delete(t);
00313                 if (dectest < decmin) {
00314                     start = first_index;
00315                     first_index = (first_index + finish)/2;
00316                 } else {
00317                     finish = first_index;
00318                     first_index = (first_index + start)/2;
00319                 }
00320             }
00321 
00322             /* Load various rows until you find the Dec range that you 
00323                have specified. Now the maximum Dec */
00324 
00325             start = first_index;
00326             finish = nrows;
00327             last_index = start + (finish - start)/2;
00328             while (finish - start >= 2) {
00329                 t = cpl_table_load_window(fullname,1,0,a,last_index,1);
00330                 dectest = cpl_table_get_float(t,"Dec",0,NULL);
00331                 cpl_table_delete(t);
00332                 if (dectest < decmax) {
00333                     start = last_index;
00334                     last_index = (last_index + finish)/2;
00335                 } else {
00336                     finish = last_index;
00337                     last_index = (last_index + start)/2;
00338                 }
00339             }    
00340             if (last_index < first_index)
00341                 last_index = first_index;
00342 
00343             /* Ok now now load all the rows in the relevant dec limits */
00344 
00345             nrows = last_index - first_index + 1;
00346             if ((t = cpl_table_load_window(fullname,1,0,NULL,first_index,
00347                                            nrows)) == NULL) {
00348                 freetable(o);
00349                 cpl_array_unwrap(a);
00350                 return(NULL);
00351             }
00352             cpl_table_unselect_all(t);
00353 
00354             /* Right, we now know what range of rows to search. Go through 
00355                these and pick the ones that are in the correct range of RA.
00356                If a row qualifies, then 'select' it. */
00357 
00358             for (j = 0; j < nrows; j++) {
00359                 ratest = cpl_table_get_float(t,"RA",j,NULL);
00360                 if (cpl_error_get_code() != CPL_ERROR_NONE) {
00361                     cpl_table_delete(t);
00362                     cpl_array_unwrap(a);
00363                     freetable(o);
00364                     return(NULL);
00365                 }
00366                 if (ratest >= ramin && ratest <= ramax)
00367                     cpl_table_select_row(t,j);
00368             }
00369 
00370             /* Extract the rows that have been selected now and append them
00371                onto the output table */
00372 
00373             s = cpl_table_extract_selected(t);
00374             if (init == 1) {
00375                 cpl_table_copy_structure(o,t);
00376                 init = 0;
00377             }
00378             irow = cpl_table_get_nrow(o) + 1;
00379             cpl_table_insert(o,s,irow);
00380 
00381             /* Tidy up */
00382 
00383             cpl_table_delete(t);
00384             cpl_table_delete(s);
00385         }
00386     }
00387 
00388     /* Ok, now just return the table and get out of here */
00389 
00390     cpl_array_unwrap(a);
00391     return(o);
00392 }
00393         
00394 /*---------------------------------------------------------------------------*/
00422 /*---------------------------------------------------------------------------*/
00423 
00424 static cpl_table *check_cache(char *catname, float ra1_im, float ra2_im, 
00425                               float dec1_im, float dec2_im) {
00426     int wrap1,wrap2;
00427     FILE *fd;
00428     char fname[BUFSIZ],catname2[SZBUF],cat_cache[SZBUF];
00429     float best,ra1_cat,ra2_cat,dec1_cat,dec2_cat,d1,d2,fra,fdec,ftot;
00430     cpl_table *out_cat;
00431 
00432     /* Open the index file.  NB the path and file name are hardcoded */
00433 
00434     fd = fopen(CACHEIND,"r");
00435     if (fd == NULL)
00436         return(NULL);
00437 
00438     /* Check to see if there is wrap around in the coordinates */
00439 
00440     wrap1 = (ra1_im < 0.0 ? 1 : 0);
00441 
00442     /* Now see if you have any matching entries */
00443 
00444     best = 0.0;
00445     while (fscanf(fd,"%s %s %g %g %g %g",fname,catname2,&ra1_cat,&ra2_cat,
00446                   &dec1_cat,&dec2_cat) != EOF) {
00447         wrap2 = (ra1_cat < 0.0 ? 1 : 0);
00448         if (wrap1 != wrap2)
00449             continue;
00450         if (strcmp(catname,catname2))
00451             continue;
00452         
00453         /* Check to see if there is at least some overlap */
00454 
00455         if (!(((ra1_im >= ra1_cat && ra1_im <= ra2_cat) ||
00456              (ra2_im >= ra1_cat && ra2_im <= ra2_cat)) &&
00457             ((dec1_im >= dec1_cat && dec1_im <= dec2_cat) ||
00458              (dec2_im >= dec1_cat && dec2_im <= dec2_cat))))
00459             continue;
00460 
00461         /* Work out exactly how much there is in each coordinate */
00462 
00463         d1 = max(0.0,ra1_cat-ra1_im);
00464         d2 = max(0.0,ra2_im-ra2_cat);
00465         fra = 1.0 - (d1 + d2)/(ra2_im - ra1_im);
00466         d1 = max(0.0,dec1_cat-dec1_im);
00467         d2 = max(0.0,dec2_im-dec2_cat);
00468         fdec = 1.0 - (d1 + d2)/(dec2_im - dec1_im);
00469         ftot = fra*fdec;
00470 
00471         /* Keep track of which is the best one */
00472 
00473         if (ftot > best) {
00474             snprintf(cat_cache,SZBUF,"%s/%s",CACHEDIR,fname);
00475             best = ftot;
00476         }
00477     }
00478     fclose(fd);
00479 
00480     /* Return a bad status if there isn't sufficient overlap */
00481 
00482     if (best < 0.9)
00483         return(NULL);
00484 
00485     /* If there is enough overlap, then try and read the FITS table. If it
00486        reads successfully, then return the table pointer */
00487 
00488     out_cat = cpl_table_load(cat_cache,1,0);
00489     return(out_cat);
00490 }
00491     
00492 /*---------------------------------------------------------------------------*/
00518 /*---------------------------------------------------------------------------*/
00519 
00520 static void addto_cache(cpl_table *stds, char *catname, float ramin, 
00521                         float ramax, float decmin, float decmax) {
00522     FILE *fd;
00523     char newname[SZBUF];
00524     int i;
00525 
00526     /* Check to see if the cache directory exists.  If it doesn't, then create
00527        it. */
00528 
00529     if (access(CACHEDIR,0) != 0)
00530         mkdir(CACHEDIR,0755);
00531 
00532     /* Open the index file with 'append' access */
00533 
00534     fd = fopen(CACHEIND,"a");
00535 
00536     /* Check the files in the directory to get a number that isn't already
00537        being used */
00538 
00539     i = 0;
00540     while (i >= 0) {
00541         i++;
00542         snprintf(newname,SZBUF,"%s/cch_%08d",CACHEDIR,i);
00543         if (access(newname,F_OK) != 0) 
00544             break;
00545     }
00546 
00547     /* Now write the current entry and make a copy of the file into the 
00548        directory */
00549 
00550     snprintf(newname,SZBUF,"%s/cch_%08d",CACHEDIR,i);
00551     cpl_table_save(stds,NULL,NULL,newname,CPL_IO_DEFAULT);
00552     snprintf(newname,SZBUF,"cch_%08d",i);
00553     if (cpl_error_get_code() == CPL_ERROR_NONE)
00554         fprintf(fd, "%s %s %g %g %g %g\n",newname,catname,ramin-0.0005,
00555                 ramax+0.0005,decmin-0.0005,decmax+0.0005);
00556     fclose(fd);
00557 }
00558 
00562 /* 
00563 
00564 $Log: vircam_getstds.c,v $
00565 Revision 1.21  2008/05/06 08:40:10  jim
00566 Modified to use cpl_wcs interface
00567 
00568 Revision 1.20  2007/10/25 17:34:00  jim
00569 Modified to remove lint warnings
00570 
00571 Revision 1.19  2007/10/19 09:25:09  jim
00572 Fixed problems with missing includes
00573 
00574 Revision 1.18  2007/10/19 06:55:06  jim
00575 Modifications made to use new method for directing the recipes to the
00576 standard catalogues using the sof
00577 
00578 Revision 1.17  2007/10/15 12:50:53  jim
00579 Modified to read new version of 2mass catalogue files
00580 
00581 Revision 1.16  2007/03/29 12:19:39  jim
00582 Little changes to improve documentation
00583 
00584 Revision 1.15  2007/03/01 12:42:41  jim
00585 Modified slightly after code checking
00586 
00587 Revision 1.14  2007/02/25 06:34:20  jim
00588 Plugged memory leak
00589 
00590 Revision 1.13  2007/01/17 23:54:00  jim
00591 Plugged some memory leaks
00592 
00593 Revision 1.12  2006/12/18 17:14:17  jim
00594 Another tidy of some error messages
00595 
00596 Revision 1.11  2006/12/18 16:41:47  jim
00597 Added bit to check for the existence of the index file rather than asking cpl
00598 to do that check with loading the table
00599 
00600 Revision 1.10  2006/12/18 12:51:20  jim
00601 Tightened up some of the error reporting
00602 
00603 Revision 1.9  2006/12/15 12:03:27  jim
00604 Fixed bug in 2mass_extract where index was offset by -1
00605 
00606 Revision 1.8  2006/11/27 12:07:22  jim
00607 Argument list now contains catname. Modified caching routines to take this
00608 into account.
00609 
00610 Revision 1.7  2006/08/11 12:45:41  jim
00611 Modified to use wcslib
00612 
00613 Revision 1.6  2006/07/04 09:19:05  jim
00614 replaced all sprintf statements with snprintf
00615 
00616 Revision 1.5  2006/06/08 14:50:09  jim
00617 Initialised output table to NULL and fixed a problem in check_cache that
00618 resets the cpl error code
00619 
00620 Revision 1.4  2006/03/23 21:18:47  jim
00621 Minor changes mainly to comment headers
00622 
00623 Revision 1.3  2006/03/22 13:58:31  jim
00624 Cosmetic fixes to keep lint happy
00625 
00626 Revision 1.2  2006/01/23 16:05:33  jim
00627 tidied up error handling
00628 
00629 Revision 1.1  2006/01/23 10:31:56  jim
00630 New file
00631 
00632 
00633 */

Generated on Wed Apr 10 04:01:55 2013 for VIRCAM Pipeline by  doxygen 1.5.1