IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 33133


Ignore:
Timestamp:
Jan 22, 2012, 6:32:39 AM (14 years ago)
Author:
eugene
Message:

make setphot a bit more generic for loading ubercal-style files (read metadata describing correction from the file itself)

Location:
branches/eam_branches/ipp-20111122/Ohana/src/uniphot
Files:
8 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20111122/Ohana/src/uniphot/include/uniphot.h

    r33114 r33133  
    6767int          VERBOSE;
    6868int          UBERCAL; // load the supplied ubercal zero point fits table (with flat-field corrections)
     69int          NO_METADATA; // the supplied ubercal data has no descriptive metadata
    6970int          NLOOP;
    7071int          TimeSelect;
    … …  
    7778time_t       TSTOP;
    7879PhotCode    *photcode;
    79 
    80 // hard-wired values which describe the ubercal analysis
    81 # define NFILTER 5
    82 # define NSEASON 4
    83 # define NCHIP_X 8
    84 # define NCHIP_Y 8
    85 # define NCELL_X 2
    86 # define NCELL_Y 2
    87 
    88 // these are initialized by initialize_setphot (in initialize.c)
    89 e_time tstart    [NSEASON];
    90 e_time tstop     [NSEASON];
    91 char   filter    [NFILTER];
    9280
    9381enum {black, white, red, orange, yellow, green, blue, indigo, violet};
    … …  
    150138
    151139ZptTable *load_zpt_ubercal (char *filename, int *nzpts, FlatCorrectionTable *flatcorrTable);
     140ZptTable *load_zpt_ubercal_nometadata (char *filename, int *nzpts, FlatCorrectionTable *flatcorrTable);
    152141int match_flatcorr_to_images (Image *image, off_t Nimage, FlatCorrectionTable *flatcorrTable);
    153142

  • branches/eam_branches/ipp-20111122/Ohana/src/uniphot/src/args.c

    r33114 r33133  
    8888  }
    8989
     90  NO_METADATA = FALSE;
     91  if ((N = get_argument (argc, argv, "-no-metadata"))) {
     92    remove_argument (N, &argc, argv);
     93    NO_METADATA = TRUE;
     94  }
     95
    9096  UPDATE = FALSE;
    9197  if ((N = get_argument (argc, argv, "-update"))) {

  • branches/eam_branches/ipp-20111122/Ohana/src/uniphot/src/initialize.c

    r33114 r33133  
    4242    exit (1);
    4343  }
    44 
    45   // we have hard-coded the MJD season ranges from Eddie in uniphot.h
    46 
    47   double tstart_mjd[] = {55000.0, 55296.0, 55327.0, 55662.0};
    48   double tstop_mjd[]  = {55296.0, 55327.0, 55662.0, 60000.0};
    49   filter[0] = 'g';
    50   filter[1] = 'r';
    51   filter[2] = 'i';
    52   filter[3] = 'z';
    53   filter[4] = 'y';
    54 
    55   int i;
    56   for (i = 0; i < NSEASON; i++) {
    57     tstart[i] = ohana_mjd_to_sec(tstart_mjd[i]);
    58     tstop[i]  = ohana_mjd_to_sec(tstop_mjd[i]);
    59   }
    6044}
    6145

  • branches/eam_branches/ipp-20111122/Ohana/src/uniphot/src/load_zpt_table.c

    r33114 r33133  
    4242}
    4343
    44 ZptTable *load_zpt_ubercal (char *filename, int *nzpts, FlatCorrectionTable *flatcorrTable) {
     44/* Slightly more generic loader than the ubercal version.  Still assumes Nfilters x Nseasons
     45   the input file must contain the following:
     46   PHU Header : NFILTER, NSEASON, NCHIP_X, NCHIP_Y, NCELL_X, NCELL_Y, TS0_nnnn (season nnnn start mjd), TS1_nnnn (season nnnn end mjd)
     47   
     48   NSEASON * 2 extensions with
     49   TABLE : mjd, zpt; Header: FILTER
     50   IMAGE : 1D array of flat offsets
     51*/
     52
     53ZptTable *load_zpt_ubercal(char *filename, int *nzpts, FlatCorrectionTable *flatcorrTable) {
    4554
    4655  int i, nfilter, nseason, ix, iy, ixc, iyc, Ncol, seq;
    4756  off_t Nrow;
    48   char type[16];
     57  char type[16], filter[80];
    4958  int Nzpts, NZPTS;
    5059  ZptTable *zpts;
    … …  
    5463  FTable ftable;
    5564
     65  // parameters describing the flat-field correction
     66  int NSEASON;
     67  int NFILTER;
     68  int NCHIP_X;
     69  int NCHIP_Y;
     70  int NCELL_X;
     71  int NCELL_Y;
     72  int CHIP_DX;
     73  int CHIP_DY;
     74
    5675  FILE *f;
    5776
    … …  
    7897  NZPTS = 0;
    7998  ALLOCATE (zpts, ZptTable, NZPTS);
     99
     100  // this function would be better if we read the list of filters, seasons, and the dimensions from the header
     101  // for current testing, make fake smfs that correspond to specific chips, filter, and mjd ranges?
     102
     103  if (!gfits_scan (&header, "NSEASON", "%d", 1, &NSEASON)) {
     104    fprintf (stderr, "cannot find NSEASON in header of %s\n", filename);
     105    fclose (f);
     106    return NULL;
     107  }
     108  if (!gfits_scan (&header, "NFILTER", "%d", 1, &NFILTER)) {
     109    fprintf (stderr, "cannot find NFILTER in header of %s\n", filename);
     110    fclose (f);
     111    return NULL;
     112  }
     113  if (!gfits_scan (&header, "NCHIP_X", "%d", 1, &NCHIP_X)) {
     114    fprintf (stderr, "cannot find NCHIP_X in header of %s\n", filename);
     115    fclose (f);
     116    return NULL;
     117  }
     118  if (!gfits_scan (&header, "NCHIP_Y", "%d", 1, &NCHIP_Y)) {
     119    fprintf (stderr, "cannot find NCHIP_Y in header of %s\n", filename);
     120    fclose (f);
     121    return NULL;
     122  }
     123  if (!gfits_scan (&header, "NCELL_X", "%d", 1, &NCELL_X)) {
     124    fprintf (stderr, "cannot find NCELL_X in header of %s\n", filename);
     125    fclose (f);
     126    return NULL;
     127  }
     128  if (!gfits_scan (&header, "NCELL_Y", "%d", 1, &NCELL_Y)) {
     129    fprintf (stderr, "cannot find NCELL_Y in header of %s\n", filename);
     130    fclose (f);
     131    return NULL;
     132  }
     133  if (!gfits_scan (&header, "CHIP_DX", "%d", 1, &CHIP_DX)) {
     134    fprintf (stderr, "cannot find CHIP_DX in header of %s\n", filename);
     135    fclose (f);
     136    return NULL;
     137  }
     138  if (!gfits_scan (&header, "CHIP_DY", "%d", 1, &CHIP_DY)) {
     139    fprintf (stderr, "cannot find CHIP_DY in header of %s\n", filename);
     140    fclose (f);
     141    return NULL;
     142  }
     143
     144  flatcorrTable->Nseason = NSEASON;
     145  ALLOCATE (flatcorrTable->tstart, e_time, NSEASON);
     146  ALLOCATE (flatcorrTable->tstop,  e_time, NSEASON);
     147
     148  for (i = 0; i < NSEASON; i++) {
     149    double mjdvalue;
     150    char name[9];
     151    snprintf (name, 9, "TS0_%04d", i);
     152    if (!gfits_scan (&header, name, "%lf", 1, &mjdvalue)) {
     153      fprintf (stderr, "cannot find %s in header of %s\n", name, filename);
     154      fclose (f);
     155      return NULL;
     156    }
     157    flatcorrTable->tstart[i] = ohana_mjd_to_sec(mjdvalue);
     158
     159    snprintf (name, 9, "TS1_%04d", i);
     160    if (!gfits_scan (&header, name, "%lf", 1, &mjdvalue)) {
     161      fprintf (stderr, "cannot find %s in header of %s\n", name, filename);
     162      fclose (f);
     163      return NULL;
     164    }
     165    flatcorrTable->tstop[i]  = ohana_mjd_to_sec(mjdvalue);
     166  }
    80167
    81168  // we have 5 filters, and 4 flat-field correction sets for each
    … …  
    115202      zpts[i+Nzpts].time = ohana_mjd_to_sec(mjd[i]);
    116203      zpts[i+Nzpts].zpt = zp[i];
    117       zpts[i+Nzpts].zpt_err = zperr[i];
     204      // zpts[i+Nzpts].zpt_err = zperr[i];
    118205    }
    119206    Nzpts += Nrow;
    120207
    121     // the image contains the flat-field corrections
     208    // the image contains the flat-field corrections for a specific filter
    122209
    123210    // load data for this header
    124211    if (!gfits_load_header (f, &header)) return (NULL);
    125212
     213    if (!gfits_scan (&theader, "FILTER", "%s", 1, filter)) {
     214      fprintf (stderr, "cannot find FILTER in header of %s\n", filename);
     215      fclose (f);
     216      return NULL;
     217    }
     218
    126219    // read the fits table bytes
    127220    if (!gfits_fread_matrix (f, &matrix, &header)) return (NULL);
     221    assert (header.bitpix == -64); // hardwired as a double
     222    double *offset = (double *) matrix.buffer;
     223
     224    for (nseason = 0; nseason < NSEASON; nseason++) { // seasons
     225      for (iy = 0; iy < NCHIP_Y; iy++) { // y-chip
     226        for (ix = 0; ix < NCHIP_X; ix++) { // x-chip
     227
     228          // photcode name
     229          char photname[64];
     230          snprintf (photname, 64, "GPC1.%s.XY%d%d", filter, ix, iy);
     231          // note that the XY00, XY07, etc, chips will have photcode values of 0
     232
     233          flatcorrTable->image[Nimage].photcode = GetPhotcodeCodebyName(photname);
     234          flatcorrTable->image[Nimage].Nx = NCELL_X;
     235          flatcorrTable->image[Nimage].Ny = NCELL_Y;
     236          flatcorrTable->image[Nimage].ID = corrID;
     237          flatcorrTable->image[Nimage].DX = CHIP_DX;
     238          flatcorrTable->image[Nimage].DY = CHIP_DY;
     239          flatcorrTable->image[Nimage].tstart = flatcorrTable->tstart[nseason];
     240          flatcorrTable->image[Nimage].tstop  = flatcorrTable->tstop[nseason];
     241         
     242          // This enforces a 180 chip rotation for XY3n - XY7n & is only known to be valid for GPC1 (the XYnn names as well)
     243          for (iyc = 0; iyc < NCELL_Y; iyc++) {
     244            for (ixc = 0; ixc < NCELL_X; ixc++) {
     245              if (ix > 3) {
     246                seq = nfilter*NSEASON*NCHIP_X*NCHIP_Y*NCELL_X*NCELL_Y + nseason*NCHIP_X*NCHIP_Y*NCELL_X*NCELL_Y + (iy*NCELL_Y + NCELL_Y - 1 - iyc)*NCHIP_X*NCELL_X + (ix*NCELL_X + NCELL_X - 1 - ixc);
     247              } else {
     248                seq = nfilter*NSEASON*NCHIP_X*NCHIP_Y*NCELL_X*NCELL_Y + nseason*NCHIP_X*NCHIP_Y*NCELL_X*NCELL_Y + (iy*NCELL_Y + iyc)*NCHIP_X*NCELL_X + (ix*NCELL_X + ixc);
     249              }
     250              assert (!flatcorrTable->corr[seq].ID);
     251              flatcorrTable->corr[seq].x = ixc;
     252              flatcorrTable->corr[seq].y = iyc;
     253              flatcorrTable->corr[seq].offset = offset[seq];
     254              flatcorrTable->corr[seq].ID = corrID;
     255            }
     256          }
     257          corrID ++;
     258          Nimage ++;
     259        }
     260      }
     261    }         
     262  }
     263
     264  /*** convert from corr,image format to offsets ***/
     265  FlatCorrectionInternal (flatcorrTable);
     266
     267  fprintf (stderr, "loaded %d zero points\n", Nzpts);
     268
     269  *nzpts = Nzpts;
     270  return zpts;
     271}
     272
     273ZptTable *load_zpt_ubercal_nometadata (char *filename, int *nzpts, FlatCorrectionTable *flatcorrTable) {
     274
     275  int i, nfilter, nseason, ix, iy, ixc, iyc, Ncol, seq;
     276  off_t Nrow;
     277  char type[16];
     278  int Nzpts, NZPTS;
     279  ZptTable *zpts;
     280  Header header;
     281  Header theader;
     282  Matrix matrix;
     283  FTable ftable;
     284
     285  // hard-wired values which describe the ubercal analysis
     286  int NFILTER = 5;
     287  int NSEASON = 4;
     288  int NCHIP_X = 8;
     289  int NCHIP_Y = 8;
     290  int NCELL_X = 2;
     291  int NCELL_Y = 2;
     292  int CHIP_DX = 4880;
     293  int CHIP_DY = 4864;
     294
     295  flatcorrTable->Nseason = NSEASON;
     296  ALLOCATE (flatcorrTable->tstart, e_time, NSEASON);
     297  ALLOCATE (flatcorrTable->tstop,  e_time, NSEASON);
     298
     299  char filter[] = {'g', 'r', 'i', 'z', 'y'};
     300
     301  double tstart_mjd[] = {55000.0, 55296.0, 55327.0, 55662.0};
     302  double tstop_mjd[]  = {55296.0, 55327.0, 55662.0, 60000.0};
     303
     304  for (i = 0; i < NSEASON; i++) {
     305    flatcorrTable->tstart[i] = ohana_mjd_to_sec(tstart_mjd[i]);
     306    flatcorrTable->tstop[i]  = ohana_mjd_to_sec(tstop_mjd[i]);
     307  }
     308
     309  FILE *f = fopen (filename, "r");
     310  if (!f) {
     311    fprintf (stderr, "ERROR: cannot open zpt table file %s\n", filename);
     312    exit (1);
     313  }
     314
     315  /* load in PHU segment (ignore) */
     316  if (!gfits_fread_header (f, &header)) {
     317    if (VERBOSE) fprintf (stderr, "can't read Flat Correction header\n");
     318    fclose (f);
     319    return (NULL);
     320  }
     321  if (!gfits_fread_matrix (f, &matrix, &header)) {
     322    if (VERBOSE) fprintf (stderr, "can't read Flat Correction matrix\n");
     323    gfits_free_header (&header);
     324    fclose (f);
     325    return (NULL);
     326  }
     327
     328  Nzpts = 0;
     329  NZPTS = 0;
     330  ALLOCATE (zpts, ZptTable, NZPTS);
     331
     332  // this function would be better if we read the list of filters, seasons, and the dimensions from the header
     333  // for current testing, make fake smfs that correspond to specific chips, filter, and mjd ranges?
     334
     335  // we have 5 filters, and 4 flat-field correction sets for each
     336  flatcorrTable->Ncorr = NFILTER*NSEASON*NCHIP_X*NCHIP_Y*NCELL_X*NCELL_Y;
     337  flatcorrTable->Nimage = NFILTER*NSEASON*NCHIP_X*NCHIP_Y;
     338
     339  ALLOCATE (flatcorrTable->corr, FlatCorrection, flatcorrTable->Ncorr);
     340  ALLOCATE (flatcorrTable->image, FlatCorrectionImage, flatcorrTable->Nimage);
     341  memset (flatcorrTable->corr, 0, flatcorrTable->Ncorr*sizeof(FlatCorrection));
     342
     343  int corrID = 1;
     344  int Nimage = 0;
     345  ftable.header = &theader;
     346  for (nfilter = 0; nfilter < NFILTER; nfilter++) {
     347    // load data for this header
     348    if (!gfits_load_header (f, &theader)) return (NULL);
     349
     350    // read the fits table bytes
     351    if (!gfits_fread_ftable_data (f, &ftable, FALSE)) return (NULL);
     352 
     353    // skip over remaining bytes in data segment
     354    fseeko (f, ftable.datasize - ftable.validsize, SEEK_CUR);
     355
     356    // need to create and assign to flat-field correction
     357    double *mjd = gfits_get_bintable_column_data (&theader, &ftable, "mjd_obs", type, &Nrow, &Ncol);
     358    assert (!strcmp(type, "double"));
     359
     360    double *zp = gfits_get_bintable_column_data (&theader, &ftable, "zp", type, &Nrow, &Ncol);
     361    assert (!strcmp(type, "double"));
     362     
     363    // float *zperr = gfits_get_bintable_column_data (&theader, &ftable, "resid", type, &Nrow, &Ncol);
     364    // assert (!strcmp(type, "float"));
     365     
     366    NZPTS += Nrow;
     367    REALLOCATE (zpts, ZptTable, NZPTS);
     368    for (i = 0; i < Nrow; i++) {
     369      zpts[i+Nzpts].time = ohana_mjd_to_sec(mjd[i]);
     370      zpts[i+Nzpts].zpt = zp[i];
     371      // zpts[i+Nzpts].zpt_err = zperr[i];
     372    }
     373    Nzpts += Nrow;
     374
     375    // the image contains the flat-field corrections
     376
     377    // load data for this header
     378    if (!gfits_load_header (f, &header)) return (NULL);
     379
     380    // read the fits table bytes
     381    if (!gfits_fread_matrix (f, &matrix, &header)) return (NULL);
     382    assert (header.bitpix == -64); // hardwired as a double
     383    double *offset = (double *) matrix.buffer;
    128384
    129385    for (nseason = 0; nseason < NSEASON; nseason++) { // seasons
    … …  
    140396          flatcorrTable->image[Nimage].Ny = NCELL_Y;
    141397          flatcorrTable->image[Nimage].ID = corrID;
    142           flatcorrTable->image[Nimage].DX = 4880;
    143           flatcorrTable->image[Nimage].DY = 4864;
    144           flatcorrTable->image[Nimage].tstart = tstart[nseason];
    145           flatcorrTable->image[Nimage].tstop  = tstop[nseason];
     398          flatcorrTable->image[Nimage].DX = CHIP_DX;
     399          flatcorrTable->image[Nimage].DY = CHIP_DY;
     400          flatcorrTable->image[Nimage].tstart = flatcorrTable->tstart[nseason];
     401          flatcorrTable->image[Nimage].tstop  = flatcorrTable->tstop[nseason];
    146402         
    147403          for (iyc = 0; iyc < NCELL_Y; iyc++) {
    … …  
    155411              flatcorrTable->corr[seq].x = ixc;
    156412              flatcorrTable->corr[seq].y = iyc;
    157               flatcorrTable->corr[seq].offset = matrix.buffer[seq];
     413              flatcorrTable->corr[seq].offset = offset[seq];
    158414              flatcorrTable->corr[seq].ID = corrID;
    159415            }

  • branches/eam_branches/ipp-20111122/Ohana/src/uniphot/src/match_flatcorr_to_images.c

    r33114 r33133  
    1818  maxCode ++; // we want the outer bound, not the last value
    1919
    20   ALLOCATE (index, short *, NSEASON);
    21   for (i = 0; i < NSEASON; i++) {
     20  ALLOCATE (index, short *, flatcorrTable->Nseason);
     21  for (i = 0; i < flatcorrTable->Nseason; i++) {
    2222      ALLOCATE (index[i], short, maxCode);
    2323      for (j = 0; j < maxCode; j++) {
    2424        index[i][j] = -1;
    2525      }
    26   }
    27 
    28   for (j = 0; j < NSEASON; j++) {
    29     assert (tstart[j]);
    30     assert (tstop[j]);
    3126  }
    3227
    … …  
    3732
    3833    // which season?
    39     for (j = 0; j < NSEASON; j++) {
    40       if (flatcorrTable->image[i].tstart == tstart[j]) {
     34    for (j = 0; j < flatcorrTable->Nseason; j++) {
     35      if (flatcorrTable->image[i].tstart == flatcorrTable->tstart[j]) {
    4136        assert (index[j][flatcorrTable->image[i].photcode] == -1);
    4237        index[j][flatcorrTable->image[i].photcode] = i;
    … …  
    4540    }
    4641  }
    47 
     42 
    4843  for (i = 0; i < Nimage; i++) {
    4944    if (!image[i].photcode) continue; // skip PHU images
    5045
    5146    int found = FALSE;
    52     for (j = 0; !found && (j < NSEASON); j++) {
    53       if (image[i].tzero < tstart[j]) continue;
    54       if (image[i].tzero > tstop[j]) continue;
     47    for (j = 0; !found && (j < flatcorrTable->Nseason); j++) {
     48      if (image[i].tzero < flatcorrTable->tstart[j]) continue;
     49      if (image[i].tzero > flatcorrTable->tstop[j]) continue;
    5550
    5651      int seq = index[j][image[i].photcode];

  • branches/eam_branches/ipp-20111122/Ohana/src/uniphot/src/match_zpts_to_images.c

    r33114 r33133  
    7272    image[Ni].dMcal = zpts[Nz].zpt_err;
    7373    image[Ni].flags &= ~ID_IMAGE_PHOTOM_NOCAL; // clear the NOCAL flag
     74    // image[Ni].flags |=  ID_IMAGE_PHOTOM_EXTERN; XXX do we want some flag like this?
     75    if (UBERCAL) {
     76      image[Ni].flags |=  ID_IMAGE_PHOTOM_UBERCAL;
     77    }
    7478    zpts[Nz].found = TRUE;
    7579    NImatch ++;

  • branches/eam_branches/ipp-20111122/Ohana/src/uniphot/src/setphot.c

    r33118 r33133  
    1919
    2020  if (UBERCAL) {
    21     zpts = load_zpt_ubercal (argv[1], &Nzpts, &flatcorrTable);
     21    if (NO_METADATA) {
     22      // the simple files from Eddie have no internal metadata describing the corrections,
     23      // so they must be manually encoded
     24      zpts = load_zpt_ubercal_nometadata (argv[1], &Nzpts, &flatcorrTable);
     25    } else {
     26      zpts = load_zpt_ubercal (argv[1], &Nzpts, &flatcorrTable);
     27    }
     28
    2229    char flatcorrfile[64];
    2330    snprintf (flatcorrfile, 64, "%s/flatcorr.fits", CATDIR);

  • branches/eam_branches/ipp-20111122/Ohana/src/uniphot/src/update_catalog_setphot.c

    r33104 r33133  
    9696      catalog[0].measure[m].Mcal = Mcal + Mcal_offset;
    9797      catalog[0].measure[m].dMcal = dMcal;
     98
     99      if (UBERCAL) {
     100        catalog[0].measure[m].dbFlags |=  ID_MEAS_PHOTOM_UBERCAL;
     101      }
     102
    98103      found ++;
    99104    }
Note: See TracChangeset for help on using the changeset viewer.