QuantitiesByAst.hxx

Go to the documentation of this file.
00001 
00002 
00003 
00004 //-----------------------------------------------------------------------------
00005 // SVN: $Id: QuantitiesByAst.hxx 6830 2010-04-16 14:33:45Z EckhardSutorius $
00006 //-----------------------------------------------------------------------------
00007 #ifndef QUANTITIES_HXX
00008 #define QUANTITIES_HXX
00009 
00010 // We will make use of AST for operations involving WCS
00011 extern "C"
00012 {
00013   #include "ast.h"
00014 }
00015 
00016 #include <cmath>
00017 #include <string>
00018 
00019 #include "Constants.h"
00020 #include "Converter.h"
00021 #include "FitsFile.h"
00022 #include "Logger.h"
00023 #include "DataMethod.hxx"
00024 #include "DataOps.h"
00025 #include "StringOps.h"
00026 //-----------------------------------------------------------------------------
00034 template<typename DataType>
00035 class QuantitiesByAst : public DataMethod<DataType>
00036 {
00037 public:
00039   QuantitiesByAst(const FitsFile& aFile)
00040     : mFile(aFile)
00041     { DataMethod<DataType>::methodName = "QuantitiesByAst"; }
00042 
00044   virtual ~QuantitiesByAst() { }
00045 
00048   void setSource(IntMap& tbls, IntMap& cols, StringMap units)
00049     { extnumColno = cols[0]; }
00050 
00052   void setTarget(IntMap& cols, StringMap units)
00053   {
00054     raColno     = cols[0];
00055     decColno    = cols[1];
00056     posangColno = cols[2];
00057     extentColno = cols[3];
00058     xsizeColno  = cols[4];
00059     ysizeColno  = cols[5];
00060     minRaColno  = cols[6];
00061     minDecColno = cols[7];
00062     maxRaColno  = cols[8];
00063     maxDecColno = cols[9];
00064 
00065     // AST will calculate angles in radians, so we need the conversion
00066     // factors for the required units in the schema
00067     Converter conv;
00068     raFact     = conv.getFactor("RADIANS", units[0]);
00069     decFact    = conv.getFactor("RADIANS", units[1]);
00070     posangFact = conv.getFactor("RADIANS", units[2]);
00071     extentFact = conv.getFactor("RADIANS", units[3]);
00072     xsizeFact  = conv.getFactor("RADIANS", units[4]);
00073     ysizeFact  = conv.getFactor("RADIANS", units[5]);
00074     minRaFact  = conv.getFactor("RADIANS", units[6]);
00075     minDecFact = conv.getFactor("RADIANS", units[7]);
00076     maxRaFact  = conv.getFactor("RADIANS", units[8]);
00077     maxDecFact = conv.getFactor("RADIANS", units[9]);
00078   }
00079 
00081   void doit(TableData<DataType>& data, int row1, int row2)
00082   {
00083     Logger log;
00084 
00085     std::string filename = mFile.getFileName();
00086     int numHDUs = mFile.getNumHdus();
00087     double crval1[numHDUs],crval2[numHDUs];
00088     int npixx[numHDUs], npixy[numHDUs];
00089 
00090     // Use these arrays for AST input
00091     double xpix[5], ypix[5], xsky[5], ysky[5];
00092 
00093     // Set up AST
00094     astBegin;
00095 
00096     // Create a new AST fits channel
00097     AstFitsChan *fitschan = astFitsChan(NULL, NULL, "%.*s", 0, "NULL");
00098     for (int rowNo = row1; rowNo < row2; ++rowNo)
00099     {
00100 
00101       // Get the FITS extension number for the data in this row. Note
00102       // that the extension number is assumed to be in the fitsio
00103       // convention
00104       Numeric value;
00105       data.value(extnumColno, rowNo, value);
00106       int hduNo = (int)value;
00107 
00108       int extNum = hduNo - 2;
00109 
00110 
00111       AstFrameSet *wcsinfo = 0;
00112 
00113       // Move to this extension on the FITS file
00114       mFile.movabsHdu(hduNo);
00115 
00116       // Put all cards in the header of this extension into the AST
00117       // fits channel
00118       int ncards = mFile.getNumCards();
00119       for (int n = 0; n < ncards; ++n) {
00120           char card[80];
00121           mFile.readCard(n, card);
00122           astPutFits(fitschan, card, 0);
00123       }
00124       // Make AST wcsinfo object
00125       astClear( fitschan, "Card" );
00126 
00127       mFile.readKey(FitsFile::CRVAL1, crval1[extNum]);
00128       mFile.readKey(FitsFile::CRVAL2, crval2[extNum]);
00129 
00130       if (crval1[extNum] == 0 && crval2[extNum] == 0) {
00131         if (filename.find("dark") == string::npos) {
00132           log.addError("Could not read WCS info in QBA for " + filename);
00133         }
00134       } else {
00135           wcsinfo = (AstFrameSet*)astRead(fitschan);
00136 
00137         // Centre of image in pixel coordinates
00138 
00139         if (mFile.hasKey("ZIMAGE"))
00140         {
00141           mFile.readKey(FitsFile::ZNAXIS1, npixx[extNum]);
00142           mFile.readKey(FitsFile::ZNAXIS2, npixy[extNum]);
00143         }
00144         else
00145         {
00146           mFile.readKey(FitsFile::NAXIS1, npixx[extNum]);
00147           mFile.readKey(FitsFile::NAXIS2, npixy[extNum]);
00148         }
00149       }
00150 
00151       if (crval1[extNum] != 0 || crval2[extNum] != 0)
00152       {
00153         // Not sure how pixel coordinates are defined. What is (x,y) of the
00154         // centre of the first pixel?
00155         xpix[0] = (double)npixx[extNum] / 2;
00156         ypix[0] = (double)npixy[extNum] / 2;
00157 
00158         // Pixel coordinates of the four corners
00159         xpix[1] = 0;
00160         ypix[1] = 0;
00161         xpix[2] = npixx[extNum];
00162         ypix[2] = 0;
00163         xpix[3] = 0;
00164         ypix[3] = npixy[extNum];
00165         xpix[4] = npixx[extNum];
00166         ypix[4] = npixy[extNum];
00167 
00168         // Get the sky coordinates for the array of pixel coordinates. Note
00169         // these will be in radians. Angles can also be negative
00170         astTran2(wcsinfo, 5, xpix, ypix, 1, xsky, ysky);
00171 
00172         // Make sure all RAs are in the range 0 to 2*pi
00173         for (int n = 0; n < 5; ++n)
00174         {
00175           if (xsky[n] < 0)
00176           {
00177             xsky[n] += 2*M_PI;
00178           }
00179         }
00180         // Get the central ra and dec in units as required in the schema
00181         Numeric raCen  = raFact * xsky[0];
00182         Numeric decCen = decFact * ysky[0];
00183 
00184         // Calculate max and min of RA and DEC from xsky, ysky
00185         // first find max and min values
00186         double xskyMin = 9999.999;
00187         double xskyMax = -9999.999;
00188         double yskyMin = 9999.999;
00189         double yskyMax = -9999.999;
00190         for (int n = 1; n < 5; ++n)
00191         {
00192           if (xsky[n] < xskyMin) xskyMin = xsky[n];
00193           if (xsky[n] > xskyMax) xskyMax = xsky[n];
00194           if (ysky[n] < yskyMin) yskyMin = ysky[n];
00195           if (ysky[n] > yskyMax) yskyMax = ysky[n];
00196         }
00197         // then project RA in 0 to 3*pi for fields crossing RA=0
00198         double h = 0.0;
00199         if (xskyMax - xskyMin > M_PI)
00200         {
00201           h = xskyMax;
00202           xskyMax = xskyMin + 2*M_PI;
00203           xskyMin = h;
00204         }
00205         Numeric raMin  = minRaFact * xskyMin;
00206         Numeric decMin = minDecFact * yskyMin;
00207         Numeric raMax  = maxRaFact * xskyMax;
00208         Numeric decMax = maxDecFact * yskyMax;
00209 
00210         // Find angular separation between the centre and each of the
00211         // four corners and find the maximum of these
00212         double maxsep = 0;
00213         for (int n = 0; n < 5; ++n)
00214         {
00215           double ang = DataOps::calcAngsep(xsky[0], ysky[0],
00216                                            xsky[n], ysky[n]);
00217           if (ang > maxsep) maxsep = ang;
00218         }
00219         // Take image extent as maximum separation and put in required units
00220         Numeric imageExtent = maxsep * extentFact;
00221 
00222         // Find the posAngle
00223         // Define N-S line by incrementing the ysky coordinate (usually dec)
00224         // by a small amount
00225         xsky[1] = xsky[0];
00226         ysky[1] = ysky[0] + 0.01;
00227         // Find equiavelent pixel coordinates for the N-S point
00228         astTran2(wcsinfo, 5, xsky, ysky, 0, xpix, ypix);
00229         // Determine angle between N-S line and positive x-axis with origin
00230         // at the centre of the image
00231         Numeric posAng = posangFact * atan2(ypix[1] - ypix[0],
00232                                             xpix[1] - xpix[0]);
00233 
00234         // X and Y pixel size is calculated by incrementing the
00235         // respective pixel sizes at the image centre
00236         double centre[2]={xpix[0],ypix[0]};
00237         double convfact[2]={xsizeFact,ysizeFact};
00238 
00239         Numeric xPixSize,yPixSize;
00240 
00241         DataOps::calcPixSize(wcsinfo,centre,convfact,xPixSize,yPixSize);
00242 
00243         data.assign(raColno, rowNo, raCen);
00244         data.assign(decColno, rowNo, decCen);
00245         data.assign(minRaColno, rowNo, raMin);
00246         data.assign(minDecColno, rowNo, decMin);
00247         data.assign(maxRaColno, rowNo, raMax);
00248         data.assign(maxDecColno, rowNo, decMax);
00249         data.assign(extentColno, rowNo, imageExtent);
00250         data.assign(posangColno, rowNo, posAng);
00251         data.assign(xsizeColno, rowNo, xPixSize);
00252         data.assign(ysizeColno, rowNo, yPixSize);
00253       }
00254     }
00255     // Finshed with AST
00256     astEnd;
00257   }
00258   // No data from raw data included
00259   void doit2(TableData<DataType>& rawData, TableData<DataType>& data, int row1, int row2) { }
00260 
00261 private:
00262   FitsFile mFile;
00263   Logger log;
00264 
00265   int extnumColno, raColno, decColno, posangColno;
00266   int extentColno, xsizeColno, ysizeColno;
00267   int minRaColno,maxRaColno,minDecColno,maxDecColno;
00268   double raFact, decFact, posangFact, extentFact, xsizeFact, ysizeFact;
00269   double minRaFact, maxRaFact, minDecFact, maxDecFact;
00270 };
00271 //-----------------------------------------------------------------------------
00272 #endif
00273 //-----------------------------------------------------------------------------
00274 // Change log:
00275 //
00276 //  7-Jun-2005, ETWS: Original version.
00277 // 22-Mar-2007, RSC:  Robustified against dodgy WCS values.
00278 //  2-May-2007,  RSC: Reverted to original design of allowing dodgy WCS values.
00279 //  7-Apr-2008, ETWS: Upgraded to use new detection table layout.
Generated on Mon Oct 4 10:38:34 2010 for WfcamSrc by  doxygen 1.6.3