Changeset 2459 for branches/hpc33

Timestamp:

04/09/12 18:49:51 (13 years ago)

Author:

Kana Sugimoto

Message:

merged Takeshi's fix to asapgrid and STGrid class form trunk

Location:

branches/hpc33

Files:

• python/asapgrid.py (modified) (10 diffs, 1 prop)
• src/STGrid.cpp (modified) (16 diffs, 1 prop)
• src/STGrid.h (modified) (1 diff, 1 prop)

branches/hpc33/python/asapgrid.py

@@ -105 +105 @@
         numerical value and unit, e.g. '0.5arcmin', or numerical value.
         If those values are specified as numerical value, their units
-        will be assumed to 'arcmin'. If which of those is not specified,
+        will be assumed to 'arcsec'. If which of those is not specified,
         it will be set to the same value as the other. If none of them
         are specified, it will be determined from map extent and number
@@ -130 +130 @@
         """
         if not isinstance( cellx, str ):
-            cellx = '%sarcmin'%(cellx)
+            cellx = '%sarcsec'%(cellx)
         if not isinstance( celly, str ):
-            celly = '%sarcmin'%(celly)
+            celly = '%sarcsec'%(celly)
         self.gridder._defineimage( nx, ny, cellx, celly, center )
 
@@ -145 +145 @@
            'box': 1 pixel
            'sf': 3 pixels
-           'gauss': 3 pixels (width is used as HWHM)
+           'gauss': 1 pixel (width is used as HWHM)
 
         func -- Function type ('box', 'sf', 'gauss').
@@ -265 +265 @@
         idx = 0
         d0 = s.get_direction( 0 ).split()[-1]
-        while ( s.get_direction(self.npol*idx).split()[-1] == d0 ):  
+        while ( s.get_direction(self.npol*idx) is not None \
+                and s.get_direction(self.npol*idx).split()[-1] == d0 ):
             idx += 1
         
@@ -276 +277 @@
         #print self.blc
         #print self.trc
-        incrx = s.get_directionval( self.npol )
-        incry = s.get_directionval( self.nx*self.npol ) 
+        if nrow > 1:
+            incrx = s.get_directionval( self.npol )
+            incry = s.get_directionval( self.nx*self.npol )
+        else:
+            incrx = [0.0,0.0]
+            incry = [0.0,0.0]
         self.cellx = abs( self.blc[0] - incrx[0] )
         self.celly = abs( self.blc[1] - incry[1] )
@@ -287 +292 @@
         else:
             opt = 'pol %s'%(pol)
-        if chan < 0:
+        if type(chan) is list:
+            opt += ', averaged over channel %s-%s'%(chan[0],chan[1])
+        elif chan < 0:
             opt += ', averaged over channel'
         else:
             opt += ', channel %s'%(chan)
-        data = self.getData( chan, pol ) 
+        data = self.getData( chan, pol )
+        data = numpy.fliplr( data )
         title = 'Gridded Image (%s)'%(opt)
         pl.figure(10)
@@ -319 +327 @@
                     #print irow
         # show image
-        extent=[self.blc[0]-0.5*self.cellx,
-                self.trc[0]+0.5*self.cellx,
+        extent=[self.trc[0]+0.5*self.cellx,
+                self.blc[0]-0.5*self.cellx,
                 self.blc[1]-0.5*self.celly,
                 self.trc[1]+0.5*self.celly]
+        deccorr = 1.0/numpy.cos(0.5*(self.blc[1]+self.trc[1]))
         pl.imshow(data,extent=extent,origin='lower',interpolation='nearest')
         pl.colorbar()
         pl.xlabel('R.A. [rad]')
         pl.ylabel('Dec. [rad]')
+        ax = pl.axes()
+        ax.set_aspect(deccorr)
         pl.title( title )
 
@@ -359 +370 @@
 
     def getData( self, chan=-1, pol=-1 ):
-        if chan == -1:
+        if type(chan) == list:
+            spectra = self.__chanAverage(start=chan[0],end=chan[1])
+        elif chan == -1:
             spectra = self.__chanAverage()
         else:
@@ -370 +383 @@
         return retval
 
-    def __chanAverage( self ):
+    def __chanAverage( self, start=-1, end=-1 ):
         s = scantable( self.outfile, average=False )
         nrow = s.nrow() 
@@ -376 +389 @@
         irow = 0
         sp = [0 for i in xrange(self.nchan)]
+        if start < 0:
+            start = 0
+        if end < 0:
+            end = self.nchan
         for i in xrange(nrow/self.npol):
             for ip in xrange(self.npol):
-                sp = s._getspectrum( irow )
+                sp = s._getspectrum( irow )[start:end]
                 spectra[ip,i] = numpy.mean( sp )
                 irow += 1
+            
         return spectra
 

branches/hpc33/src/STGrid.cpp

@@ -1 @@
-#include <iostream>
-#include <fstream>
-#include <cfloat>
-
+//
+// C++ Implementation: STGrid
+//
+// Description:
+//
+//
+// Author: Takeshi Nakazato <takeshi.nakazato@nao.ac.jp>, (C) 2011
+//
+// Copyright: See COPYING file that comes with this distribution
+//
+//
 #include <casa/BasicSL/String.h>
 #include <casa/Arrays/Vector.h>
-#include <casa/Arrays/Matrix.h>
-#include <casa/Arrays/Cube.h>
 #include <casa/Arrays/ArrayMath.h>
-#include <casa/Arrays/ArrayIter.h>
 #include <casa/Quanta/Quantum.h>
 #include <casa/Quanta/QuantumHolder.h>
@@ -16 +20 @@
 #include <tables/Tables/Table.h>
 #include <tables/Tables/TableRecord.h>
+#include <tables/Tables/TableRow.h>
 #include <tables/Tables/ExprNode.h>
 #include <tables/Tables/ScalarColumn.h>
 #include <tables/Tables/ArrayColumn.h>
-#include <tables/Tables/TableIter.h>
+#include <tables/Tables/TableCopy.h>
 
 #include <measures/Measures/MDirection.h>
 
-#include <MathUtils.h>
+#include "MathUtils.h"
+#include <atnf/PKSIO/SrcType.h>
 
 #include "STGrid.h"
@@ -407 +413 @@
   selectData() ;
   t1 = mathutil::gettimeofday_sec() ;
-  os << "selectData: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
+  os << LogIO::DEBUGGING << "selectData: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
 
   setupGrid() ;
@@ -640 +646 @@
   setConvFunc( common.convFunc ) ;
   t1 = mathutil::gettimeofday_sec() ;
-  os << "setConvFunc: elapsed time is " << t1-t0 << " sec." << LogIO::POST ; 
+  os << LogIO::DEBUGGING << "setConvFunc: elapsed time is " << t1-t0 << " sec." << LogIO::POST ; 
 
   // for performance check
@@ -681 +687 @@
     os << "end table " << ifile << LogIO::POST ;   
   }
-  os << "initPol: elapsed time is " << eInitPol << " sec." << LogIO::POST ; 
-  os << "getData: elapsed time is " << eGetData_-eToInt-eGetWeight << " sec." << LogIO::POST ; 
-  os << "toPixel: elapsed time is " << eToPixel_ << " sec." << LogIO::POST ; 
-  os << "ggridsd: elapsed time is " << eGGridSD_ << " sec." << LogIO::POST ; 
-  os << "toInt: elapsed time is " << eToInt << " sec." << LogIO::POST ;
-  os << "getWeight: elapsed time is " << eGetWeight << " sec." << LogIO::POST ;
+  os << LogIO::DEBUGGING << "initPol: elapsed time is " << eInitPol << " sec." << LogIO::POST ; 
+  os << LogIO::DEBUGGING << "getData: elapsed time is " << eGetData_-eToInt-eGetWeight << " sec." << LogIO::POST ; 
+  os << LogIO::DEBUGGING << "toPixel: elapsed time is " << eToPixel_ << " sec." << LogIO::POST ; 
+  os << LogIO::DEBUGGING << "ggridsd: elapsed time is " << eGGridSD_ << " sec." << LogIO::POST ; 
+  os << LogIO::DEBUGGING << "toInt: elapsed time is " << eToInt << " sec." << LogIO::POST ;
+  os << LogIO::DEBUGGING << "getWeight: elapsed time is " << eGetWeight << " sec." << LogIO::POST ;
   
   delete chanMap ;
@@ -790 +796 @@
   setConvFunc( common.convFunc ) ;
   t1 = mathutil::gettimeofday_sec() ;
-  os << "setConvFunc: elapsed time is " << t1-t0 << " sec." << LogIO::POST ; 
+  os << LogIO::DEBUGGING << "setConvFunc: elapsed time is " << t1-t0 << " sec." << LogIO::POST ; 
 
   // for performance check
@@ -831 +837 @@
     os << "end table " << ifile << LogIO::POST ;   
   }
-  os << "initPol: elapsed time is " << eInitPol << " sec." << LogIO::POST ; 
-  os << "getData: elapsed time is " << eGetData_-eToInt-eGetWeight << " sec." << LogIO::POST ; 
-  os << "toPixel: elapsed time is " << eToPixel_ << " sec." << LogIO::POST ; 
-  os << "ggridsd2: elapsed time is " << eGGridSD_ << " sec." << LogIO::POST ; 
-  os << "toInt: elapsed time is " << eToInt << " sec." << LogIO::POST ;
-  os << "getWeight: elapsed time is " << eGetWeight << " sec." << LogIO::POST ;
+  os << LogIO::DEBUGGING << "initPol: elapsed time is " << eInitPol << " sec." << LogIO::POST ; 
+  os << LogIO::DEBUGGING << "getData: elapsed time is " << eGetData_-eToInt-eGetWeight << " sec." << LogIO::POST ; 
+  os << LogIO::DEBUGGING << "toPixel: elapsed time is " << eToPixel_ << " sec." << LogIO::POST ; 
+  os << LogIO::DEBUGGING << "ggridsd2: elapsed time is " << eGGridSD_ << " sec." << LogIO::POST ; 
+  os << LogIO::DEBUGGING << "toInt: elapsed time is " << eToInt << " sec." << LogIO::POST ;
+  os << LogIO::DEBUGGING << "getWeight: elapsed time is " << eGetWeight << " sec." << LogIO::POST ;
   
   delete chanMap ;
@@ -855 +861 @@
               common.clipCMax ) ;
   t1 = mathutil::gettimeofday_sec() ;
-  os << "clipMinMax: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
+  os << LogIO::DEBUGGING << "clipMinMax: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
 //   os << "AFTER CLIPPING" << LogIO::POST ;
 //   os << "gdataArrC=" << common.gdataArrC << LogIO::POST ;
@@ -960 +966 @@
   }
   else 
-    ptab_ = tab_( tab_.col("POLNO") == (uInt)ipol ) ;
+    ptab_ = tab_( tab_.col("POLNO") == pollist_[ipol] ) ;
 
   attach( ptab_ ) ;
@@ -996 +1002 @@
   data_.putStorage( gwgt_p, b2 ) ;
   t1 = mathutil::gettimeofday_sec() ;
-  os << "setData: elapsed time is " << t1-t0 << " sec." << LogIO::POST ; 
+  os << LogIO::DEBUGGING << "setData: elapsed time is " << t1-t0 << " sec." << LogIO::POST ; 
 }
 
@@ -1111 +1117 @@
   }
   cellx_ = qcellx.getValue( "rad" ) ;
+  // DEC correction
+  cellx_ /= cos( center_[1] ) ;
   celly_ = qcelly.getValue( "rad" ) ;
+  //os << "cellx_=" << cellx_ << ", celly_=" << celly_ << ", cos("<<center_(1)<<")=" << cos(center_(1)) << LogIO::POST ;
   if ( nx_ < 0 ) {
     if ( wx == 0.0 ) {
@@ -1156 +1165 @@
   Int ifno = ifno_ ;
   tableList_.resize( nfile_ ) ;
-  if ( ifno == -1 ) {
+  if ( ifno_ == -1 ) {
     Table taborg( infileList_[0] ) ;
     ROScalarColumn<uInt> ifnoCol( taborg, "IFNO" ) ;
-    ifno = ifnoCol( 0 ) ;
+    ifno_ = ifnoCol( 0 ) ;
     os << LogIO::WARN
-       << "IFNO is not given. Using default IFNO: " << ifno << LogIO::POST ;
+       << "IFNO is not given. Using default IFNO: " << ifno_ << LogIO::POST ;
   }
   for ( uInt i = 0 ; i < nfile_ ; i++ ) {
     Table taborg( infileList_[i] ) ;
     TableExprNode node ;
-    if ( isMultiIF( taborg ) ) {
+    if ( ifno != -1 || isMultiIF( taborg ) ) {
       os << "apply selection on IFNO" << LogIO::POST ;
-      node = taborg.col("IFNO") == ifno ;
+      node = taborg.col("IFNO") == ifno_ ;
     }
     if ( scanlist_.size() > 0 ) {
@@ -1180 +1189 @@
       tableList_[i] = taborg( node ) ;
     }
-    os << "tableList_[" << i << "].nrow()=" << tableList_[i].nrow() << LogIO::POST ;
+    os << LogIO::DEBUGGING << "tableList_[" << i << "].nrow()=" << tableList_[i].nrow() << LogIO::POST ;
     if ( tableList_[i].nrow() == 0 ) {
       os << LogIO::SEVERE
-         << "No corresponding rows for given selection: IFNO " << ifno ;
+         << "No corresponding rows for given selection: IFNO " << ifno_ ;
       if ( scanlist_.size() > 0 )
         os << " SCANNO " << scanlist_ ;
@@ -1565 +1574 @@
     // to take into account Gaussian tail
     if ( convSupport_ < 0 )
-      convSupport_ = 12 ; // 3 * 4
+      convSupport_ = 4 ; // 1 * 4
     else {
       convSupport_ = userSupport_ * 4 ;
@@ -1650 +1659 @@
 
   t1 = mathutil::gettimeofday_sec() ;
-  os << "saveData: elapsed time is " << t1-t0 << " sec." << LogIO::POST ; 
+  os << LogIO::DEBUGGING << "saveData: elapsed time is " << t1-t0 << " sec." << LogIO::POST ; 
+
+  fillMainColumns( tab ) ;
 
   return outfile_ ;
@@ -1660 +1671 @@
   t.deepCopy( name, Table::New, False, t.endianFormat(), True ) ;
   tab = Table( name, Table::Update ) ;
-}
-
-}
+  // 2012/02/13 TN
+  // explicitly copy subtables since no rows including subtables are 
+  // copied by Table::deepCopy with noRows=True
+  TableCopy::copySubTables( tab, t ) ;
+}
+
+void STGrid::fillMainColumns( Table &tab ) 
+{
+  // values for fill
+  Table t( infileList_[0], Table::Old ) ;
+  Table tsel = t( t.col( "IFNO" ) == (uInt)ifno_, 1 ) ;
+  ROTableRow row( tsel ) ;
+  row.get( 0 ) ;
+  const TableRecord &rec = row.record() ;
+  uInt freqId = rec.asuInt( "FREQ_ID" ) ;
+  uInt molId = rec.asuInt( "MOLECULE_ID" ) ;
+  uInt tcalId = rec.asuInt( "TCAL_ID" ) ;
+  uInt focusId = rec.asuInt( "FOCUS_ID" ) ;
+  uInt weatherId = rec.asuInt( "WEATHER_ID" ) ;
+  String srcname = rec.asString( "SRCNAME" ) ;
+  String fieldname = rec.asString( "FIELDNAME" ) ;
+  Vector<Float> defaultTsys( 1, 1.0 ) ;
+  // @todo how to set flagtra for gridded spectra?
+  Vector<uChar> flagtra = rec.asArrayuChar( "FLAGTRA" ) ;
+  flagtra = (uChar)0 ;
+  Float opacity = rec.asFloat( "OPACITY" ) ;
+  Double srcvel = rec.asDouble( "SRCVELOCITY" ) ;
+  Vector<Double> srcpm = rec.asArrayDouble( "SRCPROPERMOTION" ) ;
+  Vector<Double> srcdir = rec.asArrayDouble( "SRCDIRECTION" ) ;
+  Vector<Double> scanrate = rec.asArrayDouble( "SCANRATE" ) ;
+  Double time = rec.asDouble( "TIME" ) ;
+  Double interval = rec.asDouble( "INTERVAL" ) ;
+
+  // fill columns
+  Int nrow = tab.nrow() ;
+  ScalarColumn<uInt> scannoCol( tab, "SCANNO" ) ;
+  ScalarColumn<uInt> ifnoCol( tab, "IFNO" ) ;
+  ScalarColumn<uInt> freqIdCol( tab, "FREQ_ID" ) ;
+  ScalarColumn<uInt> molIdCol( tab, "MOLECULE_ID" ) ;
+  ScalarColumn<uInt> tcalidCol( tab, "TCAL_ID" ) ;
+  ScalarColumn<Int> fitidCol( tab, "FIT_ID" ) ;
+  ScalarColumn<uInt> focusidCol( tab, "FOCUS_ID" ) ;
+  ScalarColumn<uInt> weatheridCol( tab, "WEATHER_ID" ) ;
+  ArrayColumn<uChar> flagtraCol( tab, "FLAGTRA" ) ;
+  ScalarColumn<uInt> rflagCol( tab, "FLAGROW" ) ;
+  ArrayColumn<Float> tsysCol( tab, "TSYS" ) ;
+  ScalarColumn<String> srcnameCol( tab, "SRCNAME" ) ;
+  ScalarColumn<String> fieldnameCol( tab, "FIELDNAME" ) ;
+  ScalarColumn<Int> srctypeCol( tab, "SRCTYPE" ) ;
+  ScalarColumn<Float> opacityCol( tab, "OPACITY" ) ;
+  ScalarColumn<Double> srcvelCol( tab, "SRCVELOCITY" ) ;
+  ArrayColumn<Double> srcpmCol( tab, "SRCPROPERMOTION" ) ;
+  ArrayColumn<Double> srcdirCol( tab, "SRCDIRECTION" ) ;
+  ArrayColumn<Double> scanrateCol( tab, "SCANRATE" ) ;
+  ScalarColumn<Double> timeCol( tab, "TIME" ) ;
+  ScalarColumn<Double> intervalCol( tab, "INTERVAL" ) ;
+  for ( Int i = 0 ; i < nrow ; i++ ) {
+    scannoCol.put( i, (uInt)i ) ;
+    ifnoCol.put( i, (uInt)ifno_ ) ;
+    freqIdCol.put( i, freqId ) ;
+    molIdCol.put( i, molId ) ;
+    tcalidCol.put( i, tcalId ) ;
+    fitidCol.put( i, -1 ) ;
+    focusidCol.put( i, focusId ) ;
+    weatheridCol.put( i, weatherId ) ;
+    flagtraCol.put( i, flagtra ) ;
+    rflagCol.put( i, 0 ) ;
+    tsysCol.put( i, defaultTsys ) ;
+    srcnameCol.put( i, srcname ) ;
+    fieldnameCol.put( i, fieldname ) ;
+    srctypeCol.put( i, (Int)SrcType::PSON ) ;
+    opacityCol.put( i, opacity ) ;
+    srcvelCol.put( i, srcvel ) ;
+    srcpmCol.put( i, srcpm ) ;
+    srcdirCol.put( i, srcdir ) ;
+    scanrateCol.put( i, scanrate ) ;
+    timeCol.put( i, time ) ;
+    intervalCol.put( i, interval ) ;
+  }
+}
+
+}

branches/hpc33/src/STGrid.h

@@ -195 +195 @@
   void initTable( uInt idx ) ;
   Bool isMultiIF( Table &tab ) ;
+  void fillMainColumns( Table &tab ) ;
   static bool produceChunk(void *ctx) throw(concurrent::PCException);
   static void consumeChunk(void *ctx) throw(concurrent::PCException);