Changeset 25027 for branches/pap/ippScripts/scripts/ipp_mops_translate.pl

Timestamp:

Aug 7, 2009, 4:08:25 PM (17 years ago)

Author:

Paul Price

Message:

Merging trunk (r25026) to get up-to-date on old branch.

Location:

branches/pap

Files:

• . (modified) (1 prop)
• ippScripts/scripts/ipp_mops_translate.pl (modified) (9 diffs)

branches/pap/ippScripts/scripts/ipp_mops_translate.pl

@@ -15 +15 @@
 
 use constant EXTNAME => 'MOPS_TRANSIENT_DETECTIONS'; # Extension name for output table
-use constant POSITION_ERROR => 0.2 / 3600; # Error in positions, degrees
 use constant ZERO_POINT => 25;  # Magnitude zero point
-use constant MAG_ERROR => 0.1;  # Error in magnitudes
 use constant OBSERVATORY_CODE => 566; # IAU Observatory Code
 use constant FAKE_LIMITING_MAG => 23.0; # Fake limiting magnitude to report
@@ -25 +23 @@
      $extname,                  # Name of extension containing photometry
      $skycell,                  # Skycell file with WCS
-     $diff_image_id,            # difference image id
      $output,                   # Name of output file
      $save_temps,               # Save temporary files?
@@ -33 +30 @@
            'input=s'    => \$input,
            'extname=s'  => \$extname,
-           'skycell=s'  => \$skycell,
            'output=s'   => \$output,
-           'diff_image_id=s' => \$diff_image_id,
            'save-temps' => \$save_temps,
 ) or pod2usage( 2 );
 
 pod2usage( -msg => "Unknown option: @ARGV", -exitval => 2 ) if @ARGV;
-pod2usage( -msg => "Required options: --input --extname --skycell --diff_image_id --output",
+pod2usage( -msg => "Required options: --input --extname --output",
            -exitval => 3)
     unless defined $input
     and defined $extname
-    and defined $skycell
-    and defined $diff_image_id,
     and defined $output;
 
@@ -64 +57 @@
 # Header translation table
 my $headers = {
+    # IPP name     => MOPS name, type, comment
     'FPA.RA'       => { name => 'RA',       type => TSTRING, comment => 'Right ascension of boresight'},
     'FPA.DEC'      => { name => 'DEC',      type => TSTRING, comment => 'Declination of boresight' },
-    'FPA.OBS'      => { name => 'FPA_ID',   type => TSTRING, comment => 'Exposure identifier' },
     'FPA.FILTER'   => { name => 'FILTER',   type => TSTRING, comment => 'Filter name' },
     'EXPTIME'      => { name => 'EXPTIME',  type => TDOUBLE, comment => 'Exposure time' },
@@ -72 +65 @@
     'FPA.ALT'      => { name => 'TEL_ALT',  type => TDOUBLE, comment => 'Telescope altitude' },
     'FPA.AZ'       => { name => 'TEL_AZ',   type => TDOUBLE, comment => 'Telescope azimuth' },
-            };
-
-
-# Read the skycell header for the WCS
-my $skycellResolved = $ipprc->file_resolve( $skycell );
-my $status = 0;
-my $skyfile_fits = Astro::FITS::CFITSIO::open_file( $skycellResolved, READONLY, $status ); 
-    die("failed to open skycell file: $skycellResolved: $status") if $status;
-
-my $wcsheader;
-($wcsheader, $status) = Astro::FITS::CFITSIO::fits_read_header( $skyfile_fits );
-die("Unable to read skycell header: $status") if $status;
-
-# Get the useful header keywords
-my $naxis1 = $$wcsheader{'NAXIS1'};
-my $naxis2;
-if ($naxis1) {
-    $naxis2 = $$wcsheader{'NAXIS2'} or die("Can't find NAXIS2");
-} else {
-    # if the skyfile is compressed then the WCS won't be in the primary header
-    my $hdutype;
-    $skyfile_fits->movrel_hdu(1, $hdutype, $status);
-    die("Unable to movrel_hdu: $status") if $status;
-
-    ($wcsheader, $status) = Astro::FITS::CFITSIO::fits_read_header( $skyfile_fits );
-    die("Unable to read extension header: $status") if $status;
-    my $xtension = $$wcsheader{'XTENSION'} or die("Can't find XTENSION");
-    die("XTENSION found: $xtension") if $xtension ne "'BINTABLE'";
-    $naxis1 = $$wcsheader{'ZNAXIS1'} or die("Can't find ZNAXIS1");
-    $naxis2 = $$wcsheader{'ZNAXIS2'} or die("Can't find ZNAXIS2");
-}
-
-my $ctype1 = $$wcsheader{'CTYPE1'} or die("Can't find CTYPE1");
-my $ctype2 = $$wcsheader{'CTYPE2'} or die("Can't find CTYPE2");
-my $cdelt1 = $$wcsheader{'CDELT1'} or die("Can't find CDELT1");
-my $cdelt2 = $$wcsheader{'CDELT2'} or die("Can't find CDELT2");
-my $crval1 = deg2rad($$wcsheader{'CRVAL1'}) or die("Can't find CRVAL1");
-my $crval2 = deg2rad($$wcsheader{'CRVAL2'}) or die("Can't find CRVAL2");
-my $crpix1 = $$wcsheader{'CRPIX1'} or die("Can't find CRPIX1");
-my $crpix2 = $$wcsheader{'CRPIX2'} or die("Can't find CRPIX2");
-my $pc11 = $$wcsheader{'PC001001'} or die("Can't find PC001001");
-my $pc12 = $$wcsheader{'PC001002'} or die("Can't find PC001002");
-my $pc21 = $$wcsheader{'PC002001'} or die("Can't find PC002001");
-my $pc22 = $$wcsheader{'PC002002'} or die("Can't find PC002002");
-my $crota1 = $$wcsheader{'CROTA1'};
-my $crota2 = $$wcsheader{'CROTA2'};
-
-die("Unexpected projection: $ctype1 and $ctype2.") unless $ctype1 =~ /^\'RA---TAN\s*\'$/ and
-    $ctype2 =~ /^\'DEC--TAN\s*\'$/;
-die("Can't determine size of skycell ($naxis1,$naxis2)") unless $naxis1 > 0 and $naxis2 > 0;
-die("Can't determine scale of skycell ($cdelt1,$cdelt2)") if not defined $cdelt1 or $cdelt1 == 0 or
-    not defined $cdelt2 or $cdelt2 == 0;
-die("We don't know how to handle rotations ($crota1,$crota2)") if defined $crota1 or defined $crota2;
-
-
-# Read the input table
-my $inputResolved = $ipprc->file_resolve($input);
+    'IMAGEID'      => { name => 'DIFFIMID', type => TINT,    comment => 'Difference image identifier'},
+    'FPA.OBS'      => { name => 'FPA_ID',   type => TSTRING, comment => 'Exposure identifier' },
+};
+
+
+# Read header
+my $inputResolved = $ipprc->file_resolve($input); # Resolved filename
+my $status = 0;                 # CFITSIO status
 my $inFits = Astro::FITS::CFITSIO::open_file( $inputResolved, READONLY, $status ); # FITS file handle
+die("failed to open input file: $inputResolved: $status") if $status;
 my $inHeader = $inFits->read_header(); # Header for input
 check_fitsio($status);
+
+# Read table data
 $inFits->movnam_hdu(BINARY_TBL, $extname, 0, $status) and check_fitsio($status);
 my $numRows;                    # Number of rows in table
 $inFits->get_num_rows($numRows, $status) and check_fitsio($status);
 
-my ($xCol, $yCol, $magCol, $ensCol);  # Column numbers for x and y
-$inFits->get_colnum(0, 'X_PSF', $xCol, $status) and check_fitsio($status);
-$inFits->get_colnum(0, 'Y_PSF', $yCol, $status) and check_fitsio($status);
-$inFits->get_colnum(0, 'PSF_INST_MAG', $magCol, $status) and check_fitsio($status);
-$inFits->get_colnum(0, 'EXT_NSIGMA', $ensCol, $status) and check_fitsio($status);
-
-my ($xType, $yType, $magType, $ensType);  # Column types
-$inFits->get_coltype($xCol, $xType, undef, undef, $status) and check_fitsio($status);
-$inFits->get_coltype($yCol, $yType, undef, undef, $status) and check_fitsio($status);
-$inFits->get_coltype($magCol, $magType, undef, undef, $status) and check_fitsio($status);
-$inFits->get_coltype($ensCol, $ensType, undef, undef, $status) and check_fitsio($status);
-
-my ($x, $y, $mag, $ens);              # Sources arrays
-$inFits->read_col($xType, $xCol, 1, 1, $numRows, 0, $x, undef, $status) and check_fitsio($status);
-$inFits->read_col($yType, $yCol, 1, 1, $numRows, 0, $y, undef, $status) and check_fitsio($status);
-$inFits->read_col($magType, $magCol, 1, 1, $numRows, 0, $mag, undef, $status) and check_fitsio($status);
-$inFits->read_col($ensType, $ensCol, 1, 1, $numRows, 0, $ens, undef, $status) and check_fitsio($status);
+my $ra = column($inFits, 'RA_PSF', $numRows); # Right Ascension, degrees
+my $dec = column($inFits, 'DEC_PSF', $numRows); # Declination, degrees
+my $mag = column($inFits, 'PSF_INST_MAG', $numRows); # Magnitude
+my $magErr = column($inFits, 'PSF_INST_MAG_SIG', $numRows); # Magnitude error
+my $ens = column($inFits, 'EXT_NSIGMA', $numRows); # Significance of extension
+my $xErr = column($inFits, 'X_PSF_SIG', $numRows); # Error in x position, pixels
+my $yErr = column($inFits, 'Y_PSF_SIG', $numRows); # Error in y position, pixels
+my $scale = column($inFits, 'PLTSCALE', $numRows); # Plate scale, arcsec/pixel
+my $angle = column($inFits, 'POSANGLE', $numRows); # Position angle, degrees
 
 $inFits->close_file( $status );
+
+my ($raErr, $decErr);           # Error in ra, dec
+for (my $i = 0; $i < $numRows; $i++) {
+    my $cosPA = cos($$angle[$i]);
+    my $sinPA = sin($$angle[$i]);
+
+    # XXX Not sure about the transformation here --- check
+    $$raErr[$i] = $$scale[$i] * ($cosPA * $xErr + $sinPA * $yErr) / 3600;
+    $$decErr[$i] = $$scale[$i] * ($sinPA * $xErr - $cosPA * $yErr) / 3600;
+}
+
+# Plate scales
+#my $cdelt1 = $$inHeader{'CDELT1'} or die("Can't find CDELT1");
+#my $cdelt2 = $$inHeader{'CDELT2'} or die("Can't find CDELT2");
+### XXX WCS wasn't being set in inverse diffs, but it's available elsewhere
+my $cdelt1 = $$scale[0] / 3600;
+my $cdelt2 = $$scale[1] / 3600;
 
 # Parse the list of columns
@@ -169 +125 @@
 # Convert the input data into the output formats
 for (my $i = 0; $i < $numRows; $i++) {
-    my ($ra, $dec) = pixels_to_sky($$x[$i], $$y[$i]); # Sky coordinates
-
     push @{$colData{'RA_DEG'}}, $ra;
     push @{$colData{'DEC_DEG'}}, $dec;
-    push @{$colData{'RA_SIG'}}, POSITION_ERROR;
-    push @{$colData{'DEC_SIG'}}, POSITION_ERROR;
+    push @{$colData{'RA_SIG'}}, $raErr;
+    push @{$colData{'DEC_SIG'}}, $decErr;
     push @{$colData{'FLUX'}}, $$mag[$i] + ZERO_POINT;
-    push @{$colData{'FLUX_SIG'}}, MAG_ERROR;
+    push @{$colData{'FLUX_SIG'}}, $magErr;
     push @{$colData{'ANG'}}, 0.0;
     push @{$colData{'ANG_SIG'}}, 0.0;
@@ -211 +165 @@
 }
 
-# set DIFFIMID
-# XXX TODO: we could put the diff_image_id in the .cmf file and get it from there
-$outFits->write_key( TLONGLONG, 'DIFFIMID', $diff_image_id, 'Difference image identifier', $status);
-check_fitsio( $status );
-
 # Adjust the time from start to mid-point
 {
@@ -260 +209 @@
 
 
+# Read a column specified by name
+sub column
+{
+    my $fits = shift;           # FITS file
+    my $name = shift;           # Name of column
+    my $rows = shift;           # Number of rows
+    my $status = 0;             # Status of CFITSIO
+
+    my $num;                    # Column number
+    $fits->get_colnum(0, $name, $num, $status) and check_fitsio($status);
+    my $type;                   # Type of column
+    $inFits->get_coltype($num, $type, undef, undef, $status) and check_fitsio($status);
+    my $data;                   # Array with data
+    $inFits->read_col($type, $num, 1, 1, $rows, 0, $data, undef, $status) and check_fitsio($status);
+
+    return $data;
+}
+
+
+
 
 # From Astro::FITS::CFITSIO demo
@@ -273 +242 @@
 }
 
-# Convert pixel coordinates to sky coordinates
-sub pixels_to_sky
-{
-    my ($x, $y) = @_;           # Coordinates
-
-    # Pixel coordinate relative to reference
-    $x -= $crpix1;
-    $y -= $crpix2;
-
-    # Coordinates on tangent plane
-    my $xi = $pc11 * ($x) + $pc12 * ($y);
-    my $eta = $pc21 * ($x) + $pc22 * ($y);
-    $xi *= $cdelt1;
-    $eta *= $cdelt2;
-
-    # Coordinates on rotated celestial sphere
-    my $phi = atan2($eta,$xi) + pi/2;
-    my $theta = atan(180 / pi / sqrt($xi**2 + $eta**2));
-
-    # Coordinates on celestial sphere
-    my $ra = $crval1 + atan2(cos($theta) * sin($phi),
-                             sin($theta) * cos($crval2) + cos($theta) * sin($crval2) * cos($phi));
-    my $dec = asin(sin($theta) * sin($crval2) - cos($theta) * cos($crval2) * cos($phi));
-
-    $ra = rad2deg($ra);
-    $dec = rad2deg($dec);
-
-    $ra += 360 if $ra < 0;
-    $dec += 360 if $dec < 0;
-    $ra -= 360 if $ra >= 360;
-    $dec -= 360 if $dec >= 360;
-
-    return ($ra, $dec);
-}
-
-
 __END__