Changeset 1458 for branches/alma


Ignore:
Timestamp:
12/18/08 13:15:57 (16 years ago)
Author:
Takeshi Nakazato
Message:

New Development: No

JIRA Issue: Yes CAS-1043

Ready to Release: Yes

Interface Changes: Yes/No

What Interface Changed: Please list interface changes

Test Programs: Execute following commands:


asap_init()
s=sd.scantable(filename,average=False)

where filename should be NRO data.

Put in Release Notes: Yes/No

Description: Modified STFiller class is able to import

NRO OTF raw data. STFiller automatically
identifies if a type of input file is NRO
raw data or not.


Location:
branches/alma/src
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • branches/alma/src/STFiller.cpp

    r1446 r1458  
    2525#include <tables/Tables/TableRow.h>
    2626
     27#include <measures/Measures/MDirection.h>
     28#include <measures/Measures/MeasConvert.h>
     29
    2730#include <atnf/PKSIO/PKSreader.h>
    2831#include <casa/System/ProgressMeter.h>
    29 
     32#include <atnf/PKSIO/NROReader.h>
     33
     34#include <time.h>
    3035
    3136#include "STDefs.h"
    … …  
    4247  reader_(0),
    4348  header_(0),
    44   table_(0)
     49  table_(0),
     50  nreader_(0)
    4551{
    4652}
    … …  
    4955  reader_(0),
    5056  header_(0),
    51   table_(stbl)
     57  table_(stbl),
     58  nreader_(0)
    5259{
    5360}
    … …  
    5663  reader_(0),
    5764  header_(0),
    58   table_(0)
     65  table_(0),
     66  nreader_(0)
    5967{
    6068  open(filename, whichIF, whichBeam);
    … …  
    8896  Vector<Bool> beams, ifs;
    8997  Vector<uInt> nchans,npols;
     98
     99  //
     100  // if isNRO_ is true, try NROReader
     101  //
     102  // 2008/11/11 Takeshi Nakazato
     103  isNRO_ = fileCheck() ;
     104  if ( isNRO_ ) {
     105    if ( (nreader_ = getNROReader( inName, format )) == 0 ) {
     106      throw(AipsError("Creation of NROReader failed")) ;
     107    }
     108    else {
     109      openNRO( whichIF, whichBeam ) ;
     110      return ;
     111    }
     112  }
     113  //
     114
    90115  if ( (reader_ = getPKSreader(inName, 0, 0, format, beams, ifs,
    91116                              nchans, npols, haveXPol_,haveBase, haveSpectra
    … …  
    231256{
    232257  delete reader_;reader_=0;
     258  delete nreader_;nreader_=0;
    233259  delete header_;header_=0;
    234260  table_ = 0;
    … …  
    238264{
    239265  int status = 0;
     266
     267  //
     268  // for NRO data
     269  //
     270  // 2008/11/12 Takeshi Nakazato
     271  if ( isNRO_ ) {
     272    status = readNRO() ;
     273    return status ;
     274  }
     275  //
    240276
    241277  Int    beamNo, IFno, refBeam, scanNo, cycleNo;
    … …  
    400436}
    401437
     438/**
     439 * For NRO data
     440 *
     441 * 2008/11/11 Takeshi Nakazato
     442 **/
     443void STFiller::openNRO( int whichIF, int whichBeam )
     444{
     445  // open file
     446  // DEBUG
     447  time_t t0 ;
     448  time( &t0 ) ;
     449  tm *ttm = localtime( &t0 ) ;
     450 
     451  cout << "STFiller::openNRO()  Start time = " << t0
     452       << " ("
     453       << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday
     454       << " "
     455       << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec
     456       << ")" << endl ;
     457  //
     458  if ( nreader_->open() != 0 ) {
     459    throw( AipsError( "Error while opening file "+filename_ ) ) ;
     460  }
     461
     462  isNRO_ = true ;
     463
     464  // store data into  NROHeader and NRODataset
     465  if ( nreader_->readHeader() != 0 ) {
     466    throw( AipsError( "Error while reading file "+filename_ ) ) ;
     467  }
     468
     469  // get header data
     470  NROHeader *nheader =  nreader_->getHeader() ;
     471
     472  // fill STHeader
     473  header_ = new STHeader() ;
     474
     475  header_->nchan = nreader_->getSpectrum( 0 ).size() ;
     476  header_->npol = nreader_->getPolarizationNum() ;
     477  header_->observer = nheader->getOBSVR() ;
     478  header_->project = nheader->getPROJ() ;
     479  header_->obstype = nheader->getSWMOD() ;
     480  header_->antennaname = nheader->getSITE() ;
     481  // TODO: should be investigated antenna position since there are
     482  //       no corresponding information in the header
     483  // 2008/11/13 Takeshi Nakazato
     484  //
     485  // INFO: tentative antenna posiiton is obtained for NRO 45m from ITRF website
     486  // 2008/11/26 Takeshi Nakazato
     487  vector<double> pos = nreader_->getAntennaPosition() ;
     488  header_->antennaposition = pos ;
     489  char *eq = nheader->getEPOCH() ;
     490  if ( strncmp( eq, "B1950", 5 ) == 0 )
     491    header_->equinox = 1950.0 ;
     492  else if ( strncmp( eq, "J2000", 5 ) == 0 )
     493    header_->equinox = 2000.0 ;
     494  header_->freqref = nheader->getVREF() ;
     495  header_->reffreq = nheader->getF0CAL()[0] ;
     496  header_->bandwidth = nheader->getBEBW()[0] ;
     497  header_->utc = nreader_->getStartTime() ;
     498  header_->fluxunit = "K" ;
     499  header_->epoch = "UTC" ; 
     500  char *poltp = nheader->getPOLTP()[0] ;
     501  if ( strcmp( poltp, "" ) == 0 )
     502    poltp = "None" ;
     503  header_->poltype = poltp ;
     504  // DEBUG
     505  cout << "STFiller::openNRO()  poltype = " << header_->poltype << endl ;
     506  //
     507
     508  vector<Bool> ifs = nreader_->getIFs() ;
     509  ifOffset_ = 0;
     510  nIF_ = ifs.size() ;
     511  if ( whichIF >= 0 ) {
     512    if ( whichIF >= 0 && whichIF < nIF_ ) {
     513      for ( int i = 0 ; i < nIF_ ; i++ )
     514        ifs[i] = False ;
     515      ifs[whichIF] = True ;
     516      header_->nif = 1;
     517      nIF_ = 1;
     518      ifOffset_ = whichIF;
     519    } else {
     520      delete reader_;
     521      reader_ = 0;
     522      delete header_;
     523      header_ = 0;
     524      throw(AipsError("Illegal IF selection"));
     525    }
     526  }
     527
     528  beamOffset_ = 0;
     529  vector<Bool> beams = nreader_->getBeams() ;
     530  nBeam_ = beams.size() ;
     531  if (whichBeam>=0) {
     532    if (whichBeam>=0 && whichBeam<nBeam_) {
     533      for ( int i = 0 ; i < nBeam_ ; i++ )
     534        beams[i] = False ;
     535      beams[whichBeam] = True;
     536      header_->nbeam = 1;
     537      nBeam_ = 1;
     538      beamOffset_ = whichBeam;
     539    } else {
     540      delete reader_;
     541      reader_ = 0;
     542      delete header_;
     543      header_ = 0;
     544      throw(AipsError("Illegal Beam selection"));
     545    }
     546  }
     547  header_->nbeam = nBeam_ ;
     548  header_->nif = nIF_ ;
     549
     550  // set header
     551  table_->setHeader( *header_ ) ;
     552
     553  // DEBUG
     554  time_t t1 ;
     555  time( &t1 ) ;
     556  ttm = localtime( &t1 ) ;
     557  cout << "STFiller::openNRO()  End time = " << t1
     558       << " ("
     559       << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday
     560       << " "
     561       << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec
     562       << ")" << endl ;
     563  cout << "STFiller::openNRO()  Elapsed time = " << t1 - t0 << " sec" << endl ;
     564  //
     565
     566  return ;
     567}
     568
     569int STFiller::readNRO()
     570{
     571  //return 0 ;
     572
     573  // DEBUG
     574  time_t t0 ;
     575  time( &t0 ) ;
     576  tm *ttm = localtime( &t0 ) ;
     577  cout << "STFiller::readNRO()  Start time = " << t0
     578       << " ("
     579       << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday
     580       << " "
     581       << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec
     582       << ")" << endl ;
     583  //
     584
     585  // get header
     586  //vector<NRODataset *> data = nreader_->getData() ;
     587  NROHeader *h = nreader_->getHeader() ;
     588
     589  // fill row
     590  uInt id ;
     591  uInt imax = nreader_->getRowNum() ;
     592  //uInt imax = 30 ;
     593  uInt i = 0 ;
     594  int count = 0 ;
     595  for ( i = 0 ; i < imax ; i++ ) {
     596    NRODataset *d = nreader_->getData( i ) ;
     597
     598    //char *scanType = d->getSCANTP() ;
     599    char *scanType = d->SCANTP ;
     600    Int srcType = -1 ;
     601    if ( strncmp( scanType, "ON", 2 ) == 0 ) {
     602      srcType = 0 ;
     603    }
     604    else if ( strncmp( scanType, "OFF", 3 ) == 0 ) {
     605      //cout << "OFF srcType: " << i << endl ;
     606      srcType = 1 ;
     607    }
     608    else if ( strncmp( scanType, "ZERO", 4 ) == 0 ) {
     609      //cout << "ZERO srcType: " << i << endl ;
     610      srcType = 2 ;
     611    }
     612    else {
     613      //cout << "Undefined srcType: " << i << endl ;
     614      srcType = 3 ;
     615    }
     616 
     617    // if srcType is 2 (ZERO scan), ignore scan
     618    if ( srcType != 2 && srcType != -1 && srcType != 3 ) {
     619      TableRow row( table_->table() ) ;
     620      TableRecord& rec = row.record();
     621
     622      RecordFieldPtr<Int> srctCol(rec, "SRCTYPE") ;
     623      *srctCol = srcType ;     
     624      RecordFieldPtr<Array<Double> > srateCol(rec, "SCANRATE");
     625      Array<Double> srcarr( IPosition( 1, 2 ) ) ;
     626      srcarr = 0.0 ;
     627      *srateCol = srcarr ;
     628      RecordFieldPtr<Array<Double> > spmCol(rec, "SRCPROPERMOTION") ;
     629      *spmCol = srcarr ;
     630      RecordFieldPtr<Array<Double> > sdirCol(rec, "SRCDIRECTION") ;
     631      *sdirCol = nreader_->getSourceDirection() ;
     632      RecordFieldPtr<Double> svelCol(rec, "SRCVELOCITY") ;
     633      *svelCol = h->getURVEL() ;   // [m/s]
     634      RecordFieldPtr<Int> fitCol(rec, "FIT_ID") ;
     635      *fitCol = -1 ;
     636      RecordFieldPtr<uInt> scannoCol( rec, "SCANNO" ) ;
     637      //*scannoCol = (uInt)(d->getISCAN()) ;
     638      *scannoCol = (uInt)(d->ISCAN) ;
     639      RecordFieldPtr<uInt> cyclenoCol(rec, "CYCLENO");
     640      RecordFieldPtr<Double> mjdCol( rec, "TIME" ) ;
     641      *mjdCol = Double( nreader_->getStartIntTime( i ) ) ;
     642      RecordFieldPtr<Double> intervalCol( rec, "INTERVAL" ) ;
     643      *intervalCol = Double( h->getIPTIM() ) ;
     644      RecordFieldPtr<String> srcnCol(rec, "SRCNAME") ;
     645      *srcnCol = String( h->getOBJ() ) ;
     646      RecordFieldPtr<String> fieldnCol(rec, "FIELDNAME");
     647      *fieldnCol = String( h->getOBJ() ) ;
     648      // BEAMNO is 0-base
     649      RecordFieldPtr<uInt> beamCol(rec, "BEAMNO") ;
     650      string arryt = string( d->ARRYT ) ;
     651      string sbeamno = arryt.substr( 1, arryt.size()-1 ) ;
     652      uInt ibeamno = atoi( sbeamno.c_str() ) ;
     653      *beamCol = ibeamno - 1 ;
     654      RecordFieldPtr<Int> rbCol(rec, "REFBEAMNO") ;
     655      RecordFieldPtr<uInt> ifCol(rec, "IFNO") ;
     656      RecordFieldPtr<uInt> mfreqidCol(rec, "FREQ_ID") ;
     657      vector<double> fqs = nreader_->getFrequencies( i ) ;
     658      id = table_->frequencies().addEntry( Double( fqs[0] ),
     659                                           Double( fqs[1] ),
     660                                           Double( fqs[2] ) ) ;
     661      *mfreqidCol = id ;
     662      *ifCol = id ;
     663      RecordFieldPtr<uInt> molidCol(rec, "MOLECULE_ID") ;
     664      Vector<Double> restfreq( IPosition( 1, 1 ) ) ;
     665      restfreq( 0 ) = d->FREQ0 ;
     666      id = table_->molecules().addEntry( restfreq ) ;
     667      *molidCol = id ;
     668      RecordFieldPtr<uInt> mcalidCol(rec, "TCAL_ID") ;
     669      //
     670      // No Tcal information in the data
     671      //
     672      // 2008/11/20 Takeshi Nakazato
     673      //
     674      *mcalidCol = 0 ;
     675      RecordFieldPtr<uInt> mweatheridCol(rec, "WEATHER_ID") ;
     676
     677      id = table_->weather().addEntry( Float( d->TEMP ),
     678                                       Float( d->PATM ),
     679                                       Float( d->PH2O ),
     680                                       Float( d->VWIND ),
     681                                       Float( d->DWIND ) ) ;
     682
     683      *mweatheridCol = id ;         
     684      RecordFieldPtr<uInt> mfocusidCol(rec, "FOCUS_ID") ;
     685      RecordFieldPtr< Array<Double> > dirCol(rec, "DIRECTION") ;
     686      *dirCol = nreader_->getDirection( i ) ;
     687      RecordFieldPtr<Float> azCol(rec, "AZIMUTH") ;
     688      *azCol = d->RAZ ;
     689      RecordFieldPtr<Float> elCol(rec, "ELEVATION") ;
     690      *elCol = d->REL ;
     691      RecordFieldPtr<Float> parCol(rec, "PARANGLE") ;
     692      RecordFieldPtr< Array<Float> > specCol(rec, "SPECTRA") ;
     693      vector<double> spec = nreader_->getSpectrum( i ) ;
     694      Array<Float> sp( IPosition( 1, spec.size() ) ) ;
     695      int index = 0 ;
     696      for ( Array<Float>::iterator itr = sp.begin() ; itr != sp.end() ; itr++ ) {
     697        *itr = spec[index++] ;
     698      }
     699      *specCol = sp ;
     700      RecordFieldPtr< Array<uChar> > flagCol(rec, "FLAGTRA") ;
     701      Array<uChar> flag( sp.shape() ) ;
     702      flag.set( 0 ) ;
     703      *flagCol = flag ;
     704      RecordFieldPtr< uInt > polnoCol(rec, "POLNO") ;
     705      *polnoCol = 0 ;
     706      RecordFieldPtr< Array<Float> > tsysCol(rec, "TSYS") ;
     707      Array<Float> tsys( IPosition( 1, 1 ), d->TSYS ) ;
     708      *tsysCol = tsys ;
     709
     710      table_->table().addRow() ;
     711      row.put(table_->table().nrow()-1, rec) ;
     712    }
     713    else {
     714      count++ ;
     715    }
     716  }
     717
     718  // DEBUG
     719  time_t t1 ;
     720  time( &t1 ) ;
     721  ttm = localtime( &t1 ) ;
     722  cout << "STFiller::readNRO()  Processed " << i << " rows" << endl ;
     723  cout << "STFiller::readNRO()  Added " << i - count << " rows (ignored "
     724       << count << " \"ZERO\" scans)" << endl ;
     725  cout << "STFiller::readNRO()  End time = " << t1
     726       << " ("
     727       << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday
     728       << " "
     729       << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec
     730       << ")" << endl ;
     731  cout << "STFiller::readNRO()  Elapsed time = " << t1 - t0 << " sec" << endl ;
     732  //
     733
     734  return 0 ;
     735}
     736
     737Bool STFiller::fileCheck()
     738{
     739  // if filename_ is directory, return false
     740  File inFile( filename_ ) ;
     741  if ( inFile.isDirectory() )
     742    return false ;
     743 
     744  // if beginning of header data is "RW", return true
     745  // otherwise, return false ;
     746  FILE *fp = fopen( filename_.c_str(), "r" ) ;
     747  char buf[9] ;
     748  fread( buf, 4, 1, fp ) ;
     749  buf[4] = '\0' ;
     750  if ( ( strncmp( buf, "RW", 2 ) == 0 ) )
     751    return true ;
     752
     753  return false ;
     754}
     755
    402756}//namespace asap

  • branches/alma/src/STFiller.h

    r1388 r1458  
    2727
    2828class PKSreader;
     29class NROReader;
    2930
    3031namespace asap {
    … …  
    6162    */
    6263  explicit STFiller( const std::string& filename, int whichIF=-1,
    63                   int whichBeam=-1 );
     64                     int whichBeam=-1 );
    6465
    6566  /**
    … …  
    9394  casa::CountedPtr<Scantable> getTable() const { return table_;}
    9495
     96  /**
     97   * For NRO data
     98   *
     99   * 2008/11/11 Takeshi Nakazato
     100   *
     101   * openNRO  : NRO version of open(), which performs to open file and
     102   *            read header data.
     103   * 
     104   * readNRO  : NRO version of read(), which performs to read scan
     105   *            records.
     106   *
     107   * fileCheck: Identify a type (NRO data or not) of filename_.
     108   **/
     109  void openNRO( int whichIF=-1, int whichBeam=-1 ) ;
     110  int readNRO() ;
     111  casa::Bool fileCheck() ;
     112
    95113private:
    96114
    … …  
    102120  casa::uInt ifOffset_, beamOffset_;
    103121  casa::Vector<casa::Bool> haveXPol_;
     122  NROReader *nreader_ ;
     123  casa::Bool isNRO_ ;
    104124};
    105125
Note: See TracChangeset for help on using the changeset viewer.