IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 24190


Ignore:
Timestamp:
May 14, 2009, 4:46:29 PM (17 years ago)
Author:
Paul Price
Message:

Updating scripts to provide detections to MOPS following developments in diff tables and CMF file format.

Location:
trunk/ippScripts/scripts
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/ippScripts/scripts/ipp_mops_translate.pl

    r21201 r24190  
    1515
    1616use constant EXTNAME => 'MOPS_TRANSIENT_DETECTIONS'; # Extension name for output table
    17 use constant POSITION_ERROR => 0.2 / 3600; # Error in positions, degrees
    1817use constant ZERO_POINT => 25;  # Magnitude zero point
    19 use constant MAG_ERROR => 0.1;  # Error in magnitudes
    2018use constant OBSERVATORY_CODE => 566; # IAU Observatory Code
    2119use constant FAKE_LIMITING_MAG => 23.0; # Fake limiting magnitude to report
    … …  
    2523     $extname,                  # Name of extension containing photometry
    2624     $skycell,                  # Skycell file with WCS
    27      $diff_image_id,            # difference image id
    2825     $output,                   # Name of output file
    2926     $save_temps,               # Save temporary files?
    … …  
    3330           'input=s'    => \$input,
    3431           'extname=s'  => \$extname,
    35            'skycell=s'  => \$skycell,
    3632           'output=s'   => \$output,
    37            'diff_image_id=s' => \$diff_image_id,
    3833           'save-temps' => \$save_temps,
    3934) or pod2usage( 2 );
    4035
    4136pod2usage( -msg => "Unknown option: @ARGV", -exitval => 2 ) if @ARGV;
    42 pod2usage( -msg => "Required options: --input --extname --skycell --diff_image_id --output",
     37pod2usage( -msg => "Required options: --input --extname --output",
    4338           -exitval => 3)
    4439    unless defined $input
    4540    and defined $extname
    46     and defined $skycell
    47     and defined $diff_image_id,
    4841    and defined $output;
    4942
    … …  
    6457# Header translation table
    6558my $headers = {
     59    # IPP name     => MOPS name, type, comment
    6660    'FPA.RA'       => { name => 'RA',       type => TSTRING, comment => 'Right ascension of boresight'},
    6761    'FPA.DEC'      => { name => 'DEC',      type => TSTRING, comment => 'Declination of boresight' },
    68     'FPA.OBS'      => { name => 'FPA_ID',   type => TSTRING, comment => 'Exposure identifier' },
    6962    'FPA.FILTER'   => { name => 'FILTER',   type => TSTRING, comment => 'Filter name' },
    7063    'EXPTIME'      => { name => 'EXPTIME',  type => TDOUBLE, comment => 'Exposure time' },
    … …  
    7265    'FPA.ALT'      => { name => 'TEL_ALT',  type => TDOUBLE, comment => 'Telescope altitude' },
    7366    'FPA.AZ'       => { name => 'TEL_AZ',   type => TDOUBLE, comment => 'Telescope azimuth' },
    74             };
    75 
    76 
    77 # Read the skycell header for the WCS
    78 my $skycellResolved = $ipprc->file_resolve( $skycell );
    79 my $status = 0;
    80 my $skyfile_fits = Astro::FITS::CFITSIO::open_file( $skycellResolved, READONLY, $status );
    81     die("failed to open skycell file: $skycellResolved: $status") if $status;
    82 
    83 my $wcsheader;
    84 ($wcsheader, $status) = Astro::FITS::CFITSIO::fits_read_header( $skyfile_fits );
    85 die("Unable to read skycell header: $status") if $status;
    86 
    87 # Get the useful header keywords
    88 my $naxis1 = $$wcsheader{'NAXIS1'};
    89 my $naxis2;
    90 if ($naxis1) {
    91     $naxis2 = $$wcsheader{'NAXIS2'} or die("Can't find NAXIS2");
    92 } else {
    93     # if the skyfile is compressed then the WCS won't be in the primary header
    94     my $hdutype;
    95     $skyfile_fits->movrel_hdu(1, $hdutype, $status);
    96     die("Unable to movrel_hdu: $status") if $status;
    97 
    98     ($wcsheader, $status) = Astro::FITS::CFITSIO::fits_read_header( $skyfile_fits );
    99     die("Unable to read extension header: $status") if $status;
    100     my $xtension = $$wcsheader{'XTENSION'} or die("Can't find XTENSION");
    101     die("XTENSION found: $xtension") if $xtension ne "'BINTABLE'";
    102     $naxis1 = $$wcsheader{'ZNAXIS1'} or die("Can't find ZNAXIS1");
    103     $naxis2 = $$wcsheader{'ZNAXIS2'} or die("Can't find ZNAXIS2");
    104 }
    105 
    106 my $ctype1 = $$wcsheader{'CTYPE1'} or die("Can't find CTYPE1");
    107 my $ctype2 = $$wcsheader{'CTYPE2'} or die("Can't find CTYPE2");
    108 my $cdelt1 = $$wcsheader{'CDELT1'} or die("Can't find CDELT1");
    109 my $cdelt2 = $$wcsheader{'CDELT2'} or die("Can't find CDELT2");
    110 my $crval1 = deg2rad($$wcsheader{'CRVAL1'}) or die("Can't find CRVAL1");
    111 my $crval2 = deg2rad($$wcsheader{'CRVAL2'}) or die("Can't find CRVAL2");
    112 my $crpix1 = $$wcsheader{'CRPIX1'} or die("Can't find CRPIX1");
    113 my $crpix2 = $$wcsheader{'CRPIX2'} or die("Can't find CRPIX2");
    114 my $pc11 = $$wcsheader{'PC001001'} or die("Can't find PC001001");
    115 my $pc12 = $$wcsheader{'PC001002'} or die("Can't find PC001002");
    116 my $pc21 = $$wcsheader{'PC002001'} or die("Can't find PC002001");
    117 my $pc22 = $$wcsheader{'PC002002'} or die("Can't find PC002002");
    118 my $crota1 = $$wcsheader{'CROTA1'};
    119 my $crota2 = $$wcsheader{'CROTA2'};
    120 
    121 die("Unexpected projection: $ctype1 and $ctype2.") unless $ctype1 =~ /^\'RA---TAN\s*\'$/ and
    122     $ctype2 =~ /^\'DEC--TAN\s*\'$/;
    123 die("Can't determine size of skycell ($naxis1,$naxis2)") unless $naxis1 > 0 and $naxis2 > 0;
    124 die("Can't determine scale of skycell ($cdelt1,$cdelt2)") if not defined $cdelt1 or $cdelt1 == 0 or
    125     not defined $cdelt2 or $cdelt2 == 0;
    126 die("We don't know how to handle rotations ($crota1,$crota2)") if defined $crota1 or defined $crota2;
    127 
    128 
    129 # Read the input table
    130 my $inputResolved = $ipprc->file_resolve($input);
     67    'IMAGEID'      => { name => 'DIFFIMID', type => TINT,    comment => 'Difference image identifier'},
     68    'FPA.OBS'      => { name => 'FPA_ID',   type => TSTRING, comment => 'Exposure identifier' },
     69};
     70
     71
     72# Read header
     73my $inputResolved = $ipprc->file_resolve($input); # Resolved filename
     74my $status = 0;                 # CFITSIO status
    13175my $inFits = Astro::FITS::CFITSIO::open_file( $inputResolved, READONLY, $status ); # FITS file handle
     76die("failed to open input file: $inputResolved: $status") if $status;
    13277my $inHeader = $inFits->read_header(); # Header for input
    13378check_fitsio($status);
     79
     80# Plate scales
     81my $cdelt1 = $$inHeader{'CDELT1'} or die("Can't find CDELT1");
     82my $cdelt2 = $$inHeader{'CDELT2'} or die("Can't find CDELT2");
     83
     84# Read table data
    13485$inFits->movnam_hdu(BINARY_TBL, $extname, 0, $status) and check_fitsio($status);
    13586my $numRows;                    # Number of rows in table
    13687$inFits->get_num_rows($numRows, $status) and check_fitsio($status);
    13788
    138 my ($xCol, $yCol, $magCol, $ensCol);  # Column numbers for x and y
    139 $inFits->get_colnum(0, 'X_PSF', $xCol, $status) and check_fitsio($status);
    140 $inFits->get_colnum(0, 'Y_PSF', $yCol, $status) and check_fitsio($status);
    141 $inFits->get_colnum(0, 'PSF_INST_MAG', $magCol, $status) and check_fitsio($status);
    142 $inFits->get_colnum(0, 'EXT_NSIGMA', $ensCol, $status) and check_fitsio($status);
    143 
    144 my ($xType, $yType, $magType, $ensType);  # Column types
    145 $inFits->get_coltype($xCol, $xType, undef, undef, $status) and check_fitsio($status);
    146 $inFits->get_coltype($yCol, $yType, undef, undef, $status) and check_fitsio($status);
    147 $inFits->get_coltype($magCol, $magType, undef, undef, $status) and check_fitsio($status);
    148 $inFits->get_coltype($ensCol, $ensType, undef, undef, $status) and check_fitsio($status);
    149 
    150 my ($x, $y, $mag, $ens);              # Sources arrays
    151 $inFits->read_col($xType, $xCol, 1, 1, $numRows, 0, $x, undef, $status) and check_fitsio($status);
    152 $inFits->read_col($yType, $yCol, 1, 1, $numRows, 0, $y, undef, $status) and check_fitsio($status);
    153 $inFits->read_col($magType, $magCol, 1, 1, $numRows, 0, $mag, undef, $status) and check_fitsio($status);
    154 $inFits->read_col($ensType, $ensCol, 1, 1, $numRows, 0, $ens, undef, $status) and check_fitsio($status);
     89my $ra = column($inFits, 'RA_PSF', $numRows); # Right Ascension, degrees
     90my $dec = column($inFits, 'DEC_PSF', $numRows); # Declination, degrees
     91my $mag = column($inFits, 'PSF_INST_MAG', $numRows); # Magnitude
     92my $magErr = column($inFits, 'PSF_INST_MAG_SIG', $numRows); # Magnitude error
     93my $ens = column($inFits, 'EXT_NSIGMA', $numRows); # Significance of extension
     94my $xErr = column($inFits, 'X_PSF_SIG', $numRows); # Error in x position, pixels
     95my $yErr = column($inFits, 'Y_PSF_SIG', $numRows); # Error in y position, pixels
     96my $scale = column($inFits, 'PLTSCALE', $numRows); # Plate scale, arcsec/pixel
     97my $angle = column($inFits, 'POSANGLE', $numRows); # Position angle, degrees
    15598
    15699$inFits->close_file( $status );
     100
     101my ($raErr, $decErr);           # Error in ra, dec
     102for (my $i = 0; $i < $numRows; $i++) {
     103    my $cosPA = cos($$angle[$i]);
     104    my $sinPA = sin($$angle[$i]);
     105
     106    # XXX Not sure about the transformation here --- check
     107    $$raErr[$i] = $$scale[$i] * ($cosPA * $xErr + $sinPA * $yErr) / 3600;
     108    $$decErr[$i] = $$scale[$i] * ($sinPA * $xErr - $cosPA * $yErr) / 3600;
     109}
     110
    157111
    158112# Parse the list of columns
    … …  
    169123# Convert the input data into the output formats
    170124for (my $i = 0; $i < $numRows; $i++) {
    171     my ($ra, $dec) = pixels_to_sky($$x[$i], $$y[$i]); # Sky coordinates
    172 
    173125    push @{$colData{'RA_DEG'}}, $ra;
    174126    push @{$colData{'DEC_DEG'}}, $dec;
    175     push @{$colData{'RA_SIG'}}, POSITION_ERROR;
    176     push @{$colData{'DEC_SIG'}}, POSITION_ERROR;
     127    push @{$colData{'RA_SIG'}}, $raErr;
     128    push @{$colData{'DEC_SIG'}}, $decErr;
    177129    push @{$colData{'FLUX'}}, $$mag[$i] + ZERO_POINT;
    178     push @{$colData{'FLUX_SIG'}}, MAG_ERROR;
     130    push @{$colData{'FLUX_SIG'}}, $magErr;
    179131    push @{$colData{'ANG'}}, 0.0;
    180132    push @{$colData{'ANG_SIG'}}, 0.0;
    … …  
    211163}
    212164
    213 # set DIFFIMID
    214 # XXX TODO: we could put the diff_image_id in the .cmf file and get it from there
    215 $outFits->write_key( TLONGLONG, 'DIFFIMID', $diff_image_id, 'Difference image identifier', $status);
    216 check_fitsio( $status );
    217 
    218165# Adjust the time from start to mid-point
    219166{
    … …  
    260207
    261208
     209# Read a column specified by name
     210sub column
     211{
     212    my $fits = shift;           # FITS file
     213    my $name = shift;           # Name of column
     214    my $rows = shift;           # Number of rows
     215    my $status = 0;             # Status of CFITSIO
     216
     217    my $num;                    # Column number
     218    $fits->get_colnum(0, $name, $num, $status) and check_fitsio($status);
     219    my $type;                   # Type of column
     220    $inFits->get_coltype($num, $type, undef, undef, $status) and check_fitsio($status);
     221    my $data;                   # Array with data
     222    $inFits->read_col($type, $num, 1, 1, $rows, 0, $data, undef, $status) and check_fitsio($status);
     223
     224    return $data;
     225}
     226
     227
     228
    262229
    263230# From Astro::FITS::CFITSIO demo
    … …  
    273240}
    274241
    275 # Convert pixel coordinates to sky coordinates
    276 sub pixels_to_sky
    277 {
    278     my ($x, $y) = @_;           # Coordinates
    279 
    280     # Pixel coordinate relative to reference
    281     $x -= $crpix1;
    282     $y -= $crpix2;
    283 
    284     # Coordinates on tangent plane
    285     my $xi = $pc11 * ($x) + $pc12 * ($y);
    286     my $eta = $pc21 * ($x) + $pc22 * ($y);
    287     $xi *= $cdelt1;
    288     $eta *= $cdelt2;
    289 
    290     # Coordinates on rotated celestial sphere
    291     my $phi = atan2($eta,$xi) + pi/2;
    292     my $theta = atan(180 / pi / sqrt($xi**2 + $eta**2));
    293 
    294     # Coordinates on celestial sphere
    295     my $ra = $crval1 + atan2(cos($theta) * sin($phi),
    296                              sin($theta) * cos($crval2) + cos($theta) * sin($crval2) * cos($phi));
    297     my $dec = asin(sin($theta) * sin($crval2) - cos($theta) * cos($crval2) * cos($phi));
    298 
    299     $ra = rad2deg($ra);
    300     $dec = rad2deg($dec);
    301 
    302     $ra += 360 if $ra < 0;
    303     $dec += 360 if $dec < 0;
    304     $ra -= 360 if $ra >= 360;
    305     $dec -= 360 if $dec >= 360;
    306 
    307     return ($ra, $dec);
    308 }
    309 
    310 
    311242__END__
    312243

  • trunk/ippScripts/scripts/ipp_serial_mops.pl

    r19717 r24190  
    7171    die "Unable to connect to database: $DBI::errstr";
    7272
    73 my $sql = "SELECT exp_id, diff_id, skycell_id, warp_id, diffSkyfile.path_base AS diff_path_base, warpSkyfile.path_base AS warp_path_base FROM diffSkyfile JOIN diffRun using(diff_id) JOIN diffInputSkyfile USING(diff_id,tess_id,skycell_id) JOIN warpRun USING(warp_id,tess_id) JOIN warpSkyfile USING(warp_id,skycell_id,tess_id) JOIN fakeRun using(fake_id) JOIN camRun USING(cam_id) JOIN chipRun USING(chip_id) JOIN rawExp USING(exp_id) WHERE diffSkyfile.fault = 0 AND rawExp.camera = \'$camera\'";
    74 $sql .= " AND diffRun.label = \'$label\'" if defined $label;
    75 $sql .= ';';
     73my $where_label = defined $label ? "AND label = '$label'" : ""; # WHERE for label
    7674
    77 my $rows = $db->selectall_arrayref( $sql, { Slice => {} } ) or die "Unable to execute SQL: $DBI::errstr";
    78 $db->disconnect;
     75my $sql = "
     76-- Get a list of exposures on which magic may be performed
     77SELECT
     78    exp_id,
     79    MAX(diffWarps.diff_id) AS diff_id,
     80    -- The following trick pulls out the 'inverse' value for the maximum diff_id
     81    CONVERT(SUBSTRING_INDEX(GROUP_CONCAT(diffWarps.inverse ORDER BY diffWarps.diff_id), ',', 1), UNSIGNED) AS inverse
     82FROM (
     83    -- Forward diffs
     84    SELECT
     85        diffRun.diff_id,
     86        warp1 AS warp_id,
     87        0 AS inverse
     88    FROM diffRun
     89    JOIN diffInputSkyfile USING(diff_id)
     90    WHERE diffInputSkyfile.warp1 IS NOT NULL
     91        AND diffRun.exposure = 1
     92        $where_label
     93    UNION
     94    -- Backward diffs
     95    SELECT
     96        diffRun.diff_id,
     97        warp2 AS warp_id,
     98        1 AS inverse
     99    FROM diffRun
     100    JOIN diffInputSkyfile USING(diff_id)
     101    WHERE diffInputSkyfile.warp2 IS NOT NULL
     102        AND diffRun.exposure = 1
     103        AND diffRun.bothways = 1
     104        $where_label
     105    ) AS diffWarps
     106JOIN warpRun USING(warp_id)
     107JOIN fakeRun USING(fake_id)
     108JOIN camRun USING(cam_id)
     109JOIN chipRun USING(chip_id)
     110WHERE diffRun.state = 'full'
     111    AND rawExp.camera = '$camera'
     112GROUP BY exp_id;";
    79113
    80 print "Selected " . scalar @$rows . " rows.\n";
     114my $diffs = $db->selectall_arrayref( $sql, { Slice => {} } ) or die "Unable to execute SQL: $DBI::errstr";
     115
     116print "Selected " . scalar @$diffs . " rows.\n";
    81117
    82118$ipprc->outroot_prepare( $outroot );
    … …  
    84120my ($dsFile, $dsName) = tempfile( "$outrootResolved.dslist.XXXX", UNLINK => !$save_temps);
    85121
    86 foreach my $row ( @$rows ) {
    87     my $exp_id = $row->{exp_id};
    88     my $diff_id = $row->{diff_id};
    89     my $skycell_id = $row->{skycell_id};
    90     my $warp_id = $row->{warp_id};
    91     my $diff_base = $row->{diff_path_base};
    92     my $warp_base = $row->{warp_path_base};
     122foreach my $diff ( @$diffs ) {
     123    my $exp_id = $diff->{exp_id};
     124    my $diff_id = $diff->{diff_id};
     125    my $inverse = $diff->{inverse};
    93126
    94     my $input = $ipprc->filename("PSPHOT.OUT.CMF.MEF", $diff_base);
    95     die "Can't find $input\n" unless $ipprc->file_exists($input);
    96     $input = $ipprc->file_resolve($input);
     127    my $sql = "SELECT * FROM diffSkyfile WHERE diff_id = $diff_id;";
     128    my $skycells = $db->selectall_arrayref( $sql, { Slice => {} } ) or die "Unable to execute SQL: $DBI::errstr";
    97129
    98     my $skycell = $ipprc->filename("SKYCELL.TEMPLATE", $warp_base);
    99     die "Can't find $skycell\n" unless $ipprc->file_exists($skycell);
    100     $skycell = $ipprc->file_resolve($skycell);
     130    foreach my $skycell ( @$skycells ) {
     131        my $skycell_id = $diff->{skycell_id};
     132        my $path_base = $diff->{path_base};
    101133
    102     my $output = caturi( $outroot, "mops.$exp_id.$skycell_id.$diff_id.fits" );
    103     $output = $ipprc->file_resolve($output);
     134        my $sources = $inverse ? "PPSUB.INVERSE.SOURCES" : "PPSUB.OUTPUT.SOURCES";
    104135
    105     unless ($no_op) {
    106         $ipprc->file_prepare($output);
    107         my $command = "$mops --input $input --extname SkyChip.psf --skycell $skycell --output $output";
    108         my $success = run( command => $command, verbose => $verbose );
    109         die "Couldn't translate $input\n" unless $success;
     136        my $input = $ipprc->filename($sources, $path_base);
     137        die "Can't find $input\n" unless $ipprc->file_exists($input);
     138        $input = $ipprc->file_resolve($input);
     139
     140        my $output = caturi( $outroot, "mops.$exp_id.$skycell_id.$diff_id.fits" );
     141        $output = $ipprc->file_resolve($output);
     142
     143        unless ($no_op) {
     144            $ipprc->file_prepare($output);
     145            my $command = "$mops --input $input --extname SkyChip.psf --output $output";
     146            my $success = run( command => $command, verbose => $verbose );
     147            die "Couldn't translate $input\n" unless $success;
     148        }
     149
     150        # format: filename|filesize|md5sum|filetype|
     151        # note: since we omit filesize and md5sum, dsreg will calculate them
     152        print $dsFile "${output}|||ipp-mops|\n";
    110153    }
    111 
    112     # format: filename|filesize|md5sum|filetype|
    113     # note: since we omit filesize and md5sum, dsreg will calculate them
    114     print $dsFile "${output}|||ipp-mops|\n";
    115154}
    116155close $dsFile;
     156$db->disconnect;
    117157
    118158# Register new files with the data store
Note: See TracChangeset for help on using the changeset viewer.