Context Navigation

MSWriter.cpp@ 1976

Visit:

Last change on this file since 1976 was 1975, checked in by Takeshi Nakazato, 14 years ago

New Development: Yes

JIRA Issue: Yes CAS-2718

Ready for Test: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes/No

Module(s): Module Names change impacts.

Description: Describe your changes here...

Writer improvement.

File size: 31.5 KB

Rev	Line
[1974]	1	//
	2	// C++ Interface: MSWriter
	3	//
	4	// Description:
	5	//
	6	// This class is specific writer for MS format
	7	//
	8	// Takeshi Nakazato <takeshi.nakazato@nao.ac.jp>, (C) 2010
	9	//
	10	// Copyright: See COPYING file that comes with this distribution
	11	//
	12	//
	13
	14	#include <casa/OS/File.h>
	15	#include <casa/OS/RegularFile.h>
	16	#include <casa/OS/Directory.h>
	17	#include <casa/OS/SymLink.h>
	18	#include <casa/BasicSL/String.h>
	19
[1975]	20	#include <tables/Tables/ExprNode.h>
[1974]	21	#include <tables/Tables/TableDesc.h>
	22	#include <tables/Tables/SetupNewTab.h>
	23	#include <tables/Tables/TableIter.h>
	24	#include <tables/Tables/RefRows.h>
	25
	26	#include <ms/MeasurementSets/MeasurementSet.h>
	27	#include <ms/MeasurementSets/MSColumns.h>
	28	#include <ms/MeasurementSets/MSPolIndex.h>
	29	#include <ms/MeasurementSets/MSDataDescIndex.h>
[1975]	30	#include <ms/MeasurementSets/MSSourceIndex.h>
[1974]	31
	32	#include "MSWriter.h"
	33	#include "STHeader.h"
	34	#include "STFrequencies.h"
[1975]	35	#include "STMolecules.h"
[1974]	36
	37	using namespace casa ;
	38
	39	namespace asap
	40	{
	41
	42	MSWriter::MSWriter(CountedPtr<Scantable> stable)
	43	: table_(stable)
	44	{
	45	os_ = LogIO() ;
	46	os_.origin( LogOrigin( "MSWriter", "MSWriter()", WHERE ) ) ;
	47	os_ << "MSWriter::MSWriter()" << LogIO::POST ;
	48
	49	// initialize writer
	50	init() ;
	51	}
	52
	53	MSWriter::~MSWriter()
	54	{
	55	os_.origin( LogOrigin( "MSWriter", "~MSWriter()", WHERE ) ) ;
	56	os_ << "MSWriter::~MSWriter()" << LogIO::POST ;
	57	}
	58
	59	bool MSWriter::write(const string& filename, const Record& rec)
	60	{
	61	os_.origin( LogOrigin( "MSWriter", "write()", WHERE ) ) ;
	62	os_ << "MSWriter::write()" << LogIO::POST ;
	63
	64	filename_ = filename ;
	65
	66	// parsing MS options
	67	Bool overwrite = False ;
	68	if ( rec.isDefined( "ms" ) ) {
	69	Record msrec = rec.asRecord( "ms" ) ;
	70	if ( msrec.isDefined( "overwrite" ) ) {
	71	overwrite = msrec.asBool( "overwrite" ) ;
	72	}
	73	}
	74
	75	os_ << "Parsing MS options" << endl ;
	76	os_ << " overwrite = " << overwrite << LogIO::POST ;
	77
	78	File file( filename_ ) ;
	79	if ( file.exists() ) {
	80	if ( overwrite ) {
	81	os_ << filename_ << " exists. Overwrite existing data... " << LogIO::POST ;
	82	if ( file.isRegular() ) RegularFile(file).remove() ;
	83	else if ( file.isDirectory() ) Directory(file).removeRecursive() ;
	84	else SymLink(file).remove() ;
	85	}
	86	else {
	87	os_ << LogIO::SEVERE << "ERROR: " << filename_ << " exists..." << LogIO::POST ;
	88	return False ;
	89	}
	90	}
	91
	92	// set up MS
	93	setupMS() ;
	94
	95	// subtables
	96	// OBSERVATION
	97	fillObservation() ;
	98
	99	// ANTENNA
	100	fillAntenna() ;
	101
[1975]	102	// PROCESSOR
	103	fillProcessor() ;
	104
	105	// SOURCE
	106	fillSource() ;
	107
[1974]	108	// MAIN
	109	// Iterate over several ids
	110	Vector<uInt> processedFreqId( 0 ) ;
	111	Int defaultFieldId = 0 ;
	112	//
	113	// ITERATION: FIELDNAME
	114	//
	115	Int added0 = 0 ;
	116	Int current0 = mstable_->nrow() ;
	117	TableIterator iter0( table_->table(), "FIELDNAME" ) ;
	118	while( !iter0.pastEnd() ) {
	119	Table t0( iter0.table() ) ;
[1975]	120	ROScalarColumn<String> sharedStrCol( t0, "FIELDNAME" ) ;
	121	String fieldName = sharedStrCol( 0 ) ;
	122	sharedStrCol.attach( t0, "SRCNAME" ) ;
	123	String srcName = sharedStrCol( 0 ) ;
	124	ROScalarColumn<Double> timeCol( t0, "TIME" ) ;
	125	Double minTime = min( timeCol.getColumn() ) ;
	126	ROArrayColumn<Double> scanRateCol( t0, "SCANRATE" ) ;
	127	Vector<Double> scanRate = scanRateCol.getColumn()[0] ;
[1974]	128	String::size_type pos = fieldName.find( "__" ) ;
	129	Int fieldId = -1 ;
	130	if ( pos != String::npos ) {
	131	os_ << "fieldName.substr( pos+2 )=" << fieldName.substr( pos+2 ) << LogIO::POST ;
	132	fieldId = String::toInt( fieldName.substr( pos+2 ) ) ;
	133	fieldName = fieldName.substr( 0, pos ) ;
	134	}
	135	else {
	136	os_ << "use default field id" << LogIO::POST ;
	137	fieldId = defaultFieldId ;
	138	defaultFieldId++ ;
	139	}
	140	os_ << "fieldId" << fieldId << ": " << fieldName << LogIO::POST ;
	141	//
	142	// ITERATION: BEAMNO
	143	//
	144	Int added1 = 0 ;
	145	Int current1 = mstable_->nrow() ;
	146	TableIterator iter1( t0, "BEAMNO" ) ;
	147	while( !iter1.pastEnd() ) {
	148	Table t1( iter1.table() ) ;
	149	ROScalarColumn<uInt> beamNoCol( t1, "BEAMNO" ) ;
[1975]	150	uInt beamNo = beamNoCol( 0 ) ;
[1974]	151	os_ << "beamNo = " << beamNo << LogIO::POST ;
	152	//
	153	// ITERATION: SCANNO
	154	//
	155	Int added2 = 0 ;
	156	Int current2 = mstable_->nrow() ;
	157	TableIterator iter2( t1, "SCANNO" ) ;
	158	while( !iter2.pastEnd() ) {
	159	Table t2( iter2.table() ) ;
	160	ROScalarColumn<uInt> scanNoCol( t2, "SCANNO" ) ;
[1975]	161	uInt scanNo = scanNoCol( 0 ) ;
[1974]	162	os_ << "scanNo = " << scanNo << LogIO::POST ;
	163	//
	164	// ITERATION: CYCLENO
	165	//
	166	Int added3 = 0 ;
	167	Int current3 = mstable_->nrow() ;
	168	TableIterator iter3( t2, "CYCLENO" ) ;
	169	while( !iter3.pastEnd() ) {
	170	Table t3( iter3.table() ) ;
	171	//
	172	// ITERATION: IFNO
	173	//
	174	Int added4 = 0 ;
	175	Int current4 = mstable_->nrow() ;
	176	TableIterator iter4( t3, "IFNO" ) ;
	177	while( !iter4.pastEnd() ) {
	178	Table t4( iter4.table() ) ;
	179	ROScalarColumn<uInt> ifNoCol( t4, "IFNO" ) ;
[1975]	180	uInt ifNo = ifNoCol( 0 ) ;
[1974]	181	os_ << "ifNo = " << ifNo << LogIO::POST ;
	182	ROScalarColumn<uInt> freqIdCol( t4, "FREQ_ID" ) ;
[1975]	183	uInt freqId = freqIdCol( 0 ) ;
[1974]	184	os_ << "freqId = " << freqId << LogIO::POST ;
	185	//
	186	// ITERATION: TIME
	187	//
	188	Int added5 = 0 ;
	189	Int current5 = mstable_->nrow() ;
	190	TableIterator iter5( t4, "TIME" ) ;
	191	while( !iter5.pastEnd() ) {
	192	Table t5( iter5.table().sort("POLNO") ) ;
	193	Int prevnr = mstable_->nrow() ;
	194	Int nrow = t5.nrow() ;
	195	os_ << "nrow = " << nrow << LogIO::POST ;
	196
	197	// add row
	198	mstable_->addRow( 1, True ) ;
	199
	200	Vector<Int> polnos( nrow ) ;
	201	indgen( polnos, 0 ) ;
	202	Int polid = addPolarization( polnos ) ;
	203	os_ << "polid = " << polid << LogIO::POST ;
	204	//
	205	// LOOP: POLNO
	206	//
	207	for ( Int ipol = 0 ; ipol < nrow ; ipol++ ) {
	208	}
	209
	210	// TIME and TIME_CENTROID
	211	ROScalarMeasColumn<MEpoch> timeCol( t5, "TIME" ) ;
	212	ScalarMeasColumn<MEpoch> msTimeCol( *mstable_, "TIME" ) ;
[1975]	213	msTimeCol.put( prevnr, timeCol( 0 ) ) ;
[1974]	214	msTimeCol.attach( *mstable_, "TIME_CENTROID" ) ;
[1975]	215	msTimeCol.put( prevnr, timeCol( 0 ) ) ;
[1974]	216
	217	// INTERVAL and EXPOSURE
	218	ROScalarColumn<Double> intervalCol( t5, "INTERVAL" ) ;
	219	ScalarColumn<Double> msIntervalCol( *mstable_, "INTERVAL" ) ;
[1975]	220	msIntervalCol.put( prevnr, intervalCol( 0 ) ) ;
[1974]	221	msIntervalCol.attach( *mstable_, "EXPOSURE" ) ;
[1975]	222	msIntervalCol.put( prevnr, intervalCol( 0 ) ) ;
[1974]	223
	224	// add DATA_DESCRIPTION row
	225	Int ddid = addDataDescription( polid, ifNo ) ;
	226	os_ << "ddid = " << ddid << LogIO::POST ;
	227	ScalarColumn<Int> ddIdCol( *mstable_, "DATA_DESC_ID" ) ;
	228	ddIdCol.put( prevnr, ddid ) ;
	229
	230	added5 += 1 ;
	231	os_ << "added5 = " << added5 << " current5 = " << current5 << LogIO::POST ;
	232	iter5.next() ;
	233	}
	234
	235	// add SPECTRAL_WINDOW row
	236	if ( allNE( processedFreqId, freqId ) ) {
	237	uInt vsize = processedFreqId.size() ;
	238	processedFreqId.resize( vsize+1, True ) ;
	239	processedFreqId[vsize] = freqId ;
	240	addSpectralWindow( ifNo, freqId ) ;
	241	}
	242
	243	added4 += added5 ;
	244	os_ << "added4 = " << added4 << " current4 = " << current4 << LogIO::POST ;
	245	iter4.next() ;
	246	}
	247	added3 += added4 ;
	248	os_ << "added3 = " << added3 << " current3 = " << current3 << LogIO::POST ;
	249	iter3.next() ;
	250	}
	251
	252	// SCAN_NUMBER
	253	RefRows rows3( current3, current3+added3-1 ) ;
	254	Vector<Int> scanNum( added3, scanNo ) ;
	255	ScalarColumn<Int> scanNumCol( *mstable_, "SCAN_NUMBER" ) ;
	256	scanNumCol.putColumnCells( rows3, scanNum ) ;
	257
	258	added2 += added3 ;
	259	os_ << "added2 = " << added2 << " current2 = " << current2 << LogIO::POST ;
	260	iter2.next() ;
	261	}
	262
	263	// FEED1 and FEED2
	264	RefRows rows2( current2, current2+added2-1 ) ;
	265	Vector<Int> feedId( added2, beamNo ) ;
	266	ScalarColumn<Int> feedCol( *mstable_, "FEED1" ) ;
	267	feedCol.putColumnCells( rows2, feedId ) ;
	268	feedCol.attach( *mstable_, "FEED2" ) ;
	269	feedCol.putColumnCells( rows2, feedId ) ;
	270
	271	// add FEED row
	272	addFeed( beamNo ) ;
	273
	274	added1 += added2 ;
	275	os_ << "added1 = " << added1 << " current1 = " << current1 << LogIO::POST ;
	276	iter1.next() ;
	277	}
	278
	279	// FIELD_ID
	280	RefRows rows1( current1, current1+added1-1 ) ;
	281	Vector<Int> fieldIds( added1, fieldId ) ;
	282	ScalarColumn<Int> fieldIdCol( *mstable_, "FIELD_ID" ) ;
	283	fieldIdCol.putColumnCells( rows1, fieldIds ) ;
	284
[1975]	285	// add FIELD row
	286	addField( fieldId, fieldName, srcName, minTime, scanRate ) ;
	287
[1974]	288	added0 += added1 ;
	289	os_ << "added0 = " << added0 << " current0 = " << current0 << LogIO::POST ;
	290	iter0.next() ;
	291	}
	292
	293	// OBSERVATION_ID is always 0
	294	ScalarColumn<Int> sharedIntCol( *mstable_, "OBSERVATION_ID" ) ;
	295	Vector<Int> sharedIntArr( added0, 0 ) ;
	296	sharedIntCol.putColumn( sharedIntArr ) ;
	297
	298	// ANTENNA1 and ANTENNA2 are always 0
	299	sharedIntArr = 0 ;
	300	sharedIntCol.attach( *mstable_, "ANTENNA1" ) ;
	301	sharedIntCol.putColumn( sharedIntArr ) ;
	302	sharedIntCol.attach( *mstable_, "ANTENNA2" ) ;
	303	sharedIntCol.putColumn( sharedIntArr ) ;
	304
	305	// ARRAY_ID is tentatively set to 0
	306	sharedIntArr = 0 ;
	307	sharedIntCol.attach( *mstable_, "ARRAY_ID" ) ;
	308	sharedIntCol.putColumn( sharedIntArr ) ;
	309
[1975]	310	// PROCESSOR_ID is tentatively set to 0
	311	sharedIntArr = 0 ;
	312	sharedIntCol.attach( *mstable_, "PROCESSOR_ID" ) ;
	313	sharedIntCol.putColumn( sharedIntArr ) ;
	314
[1974]	315	return True ;
	316	}
	317
	318	void MSWriter::init()
	319	{
	320	// os_.origin( LogOrigin( "MSWriter", "init()", WHERE ) ) ;
	321	os_ << "MSWriter::init()" << LogIO::POST ;
	322
	323	// access to scantable
	324	header_ = table_->getHeader() ;
	325
	326	// FLOAT_DATA? or DATA?
	327	if ( header_.npol > 2 ) {
	328	isFloatData_ = False ;
	329	isData_ = True ;
	330	}
	331	else {
	332	isFloatData_ = True ;
	333	isData_ = False ;
	334	}
	335
	336	// polarization type
	337	polType_ = header_.poltype ;
	338	if ( polType_ == "" )
	339	polType_ = "stokes" ;
	340	else if ( polType_.find( "linear" ) != String::npos )
	341	polType_ = "linear" ;
	342	else if ( polType_.find( "circular" ) != String::npos )
	343	polType_ = "circular" ;
	344	else if ( polType_.find( "stokes" ) != String::npos )
	345	polType_ = "stokes" ;
	346	else if ( polType_.find( "linpol" ) != String::npos )
	347	polType_ = "linpol" ;
	348	else
	349	polType_ = "notype" ;
	350
	351	}
	352
	353	void MSWriter::setupMS()
	354	{
	355	// os_.origin( LogOrigin( "MSWriter", "setupMS()", WHERE ) ) ;
	356	os_ << "MSWriter::setupMS()" << LogIO::POST ;
	357
	358	TableDesc msDesc = MeasurementSet::requiredTableDesc() ;
	359
	360	// any additional columns?
	361	// TODO: add FLOAT_DATA column if npol == 2 otherwise add DATA column
	362
	363	SetupNewTable newtab( filename_, msDesc, Table::New ) ;
	364
	365	mstable_ = new MeasurementSet( newtab ) ;
	366
	367	// create subtables
	368	TableDesc antennaDesc = MSAntenna::requiredTableDesc() ;
	369	SetupNewTable antennaTab( mstable_->antennaTableName(), antennaDesc, Table::New ) ;
	370	mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::ANTENNA ), Table( antennaTab ) ) ;
	371
	372	TableDesc dataDescDesc = MSDataDescription::requiredTableDesc() ;
	373	SetupNewTable dataDescTab( mstable_->dataDescriptionTableName(), dataDescDesc, Table::New ) ;
	374	mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::DATA_DESCRIPTION ), Table( dataDescTab ) ) ;
	375
	376	TableDesc dopplerDesc = MSDoppler::requiredTableDesc() ;
	377	SetupNewTable dopplerTab( mstable_->dopplerTableName(), dopplerDesc, Table::New ) ;
	378	mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::DOPPLER ), Table( dopplerTab ) ) ;
	379
	380	TableDesc feedDesc = MSFeed::requiredTableDesc() ;
	381	SetupNewTable feedTab( mstable_->feedTableName(), feedDesc, Table::New ) ;
	382	mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::FEED ), Table( feedTab ) ) ;
	383
	384	TableDesc fieldDesc = MSField::requiredTableDesc() ;
	385	SetupNewTable fieldTab( mstable_->fieldTableName(), fieldDesc, Table::New ) ;
	386	mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::FIELD ), Table( fieldTab ) ) ;
	387
	388	TableDesc flagCmdDesc = MSFlagCmd::requiredTableDesc() ;
	389	SetupNewTable flagCmdTab( mstable_->flagCmdTableName(), flagCmdDesc, Table::New ) ;
	390	mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::FLAG_CMD ), Table( flagCmdTab ) ) ;
	391
	392	TableDesc freqOffsetDesc = MSFreqOffset::requiredTableDesc() ;
	393	SetupNewTable freqOffsetTab( mstable_->freqOffsetTableName(), freqOffsetDesc, Table::New ) ;
	394	mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::FREQ_OFFSET ), Table( freqOffsetTab ) ) ;
	395
	396	TableDesc historyDesc = MSHistory::requiredTableDesc() ;
	397	SetupNewTable historyTab( mstable_->historyTableName(), historyDesc, Table::New ) ;
	398	mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::HISTORY ), Table( historyTab ) ) ;
	399
	400	TableDesc observationDesc = MSObservation::requiredTableDesc() ;
	401	SetupNewTable observationTab( mstable_->observationTableName(), observationDesc, Table::New ) ;
	402	mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::OBSERVATION ), Table( observationTab ) ) ;
	403
	404	TableDesc pointingDesc = MSPointing::requiredTableDesc() ;
	405	SetupNewTable pointingTab( mstable_->pointingTableName(), pointingDesc, Table::New ) ;
	406	mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::POINTING ), Table( pointingTab ) ) ;
	407
	408	TableDesc polarizationDesc = MSPolarization::requiredTableDesc() ;
	409	SetupNewTable polarizationTab( mstable_->polarizationTableName(), polarizationDesc, Table::New ) ;
	410	mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::POLARIZATION ), Table( polarizationTab ) ) ;
	411
	412	TableDesc processorDesc = MSProcessor::requiredTableDesc() ;
	413	SetupNewTable processorTab( mstable_->processorTableName(), processorDesc, Table::New ) ;
	414	mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::PROCESSOR ), Table( processorTab ) ) ;
	415
	416	TableDesc sourceDesc = MSSource::requiredTableDesc() ;
[1975]	417	MSSource::addColumnToDesc( sourceDesc, MSSourceEnums::TRANSITION, 1 ) ;
	418	MSSource::addColumnToDesc( sourceDesc, MSSourceEnums::REST_FREQUENCY, 1 ) ;
	419	MSSource::addColumnToDesc( sourceDesc, MSSourceEnums::SYSVEL, 1 ) ;
[1974]	420	SetupNewTable sourceTab( mstable_->sourceTableName(), sourceDesc, Table::New ) ;
	421	mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::SOURCE ), Table( sourceTab ) ) ;
	422
	423	TableDesc spwDesc = MSSpectralWindow::requiredTableDesc() ;
	424	SetupNewTable spwTab( mstable_->spectralWindowTableName(), spwDesc, Table::New ) ;
	425	mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::SPECTRAL_WINDOW ), Table( spwTab ) ) ;
	426
	427	TableDesc stateDesc = MSState::requiredTableDesc() ;
	428	SetupNewTable stateTab( mstable_->stateTableName(), stateDesc, Table::New ) ;
	429	mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::STATE ), Table( stateTab ) ) ;
	430
	431	// TODO: add TCAL_SPECTRUM and TSYS_SPECTRUM if necessary
	432	TableDesc sysCalDesc = MSSysCal::requiredTableDesc() ;
	433	SetupNewTable sysCalTab( mstable_->sysCalTableName(), sysCalDesc, Table::New ) ;
	434	mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::SYSCAL ), Table( sysCalTab ) ) ;
	435
	436	TableDesc weatherDesc = MSWeather::requiredTableDesc() ;
	437	SetupNewTable weatherTab( mstable_->weatherTableName(), weatherDesc, Table::New ) ;
	438	mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::WEATHER ), Table( weatherTab ) ) ;
	439
	440	mstable_->initRefs() ;
	441
	442	}
	443
	444	void MSWriter::fillObservation()
	445	{
	446	os_ << "set up OBSERVATION subtable" << LogIO::POST ;
	447	// only 1 row
	448	mstable_->observation().addRow( 1, True ) ;
	449	MSObservationColumns msObsCols( mstable_->observation() ) ;
	450	msObsCols.observer().put( 0, header_.observer ) ;
	451	// tentatively put antennaname (from ANTENNA subtable)
	452	String hAntennaName = header_.antennaname ;
	453	String::size_type pos = hAntennaName.find( "//" ) ;
	454	String telescopeName ;
	455	if ( pos != String::npos ) {
	456	telescopeName = hAntennaName.substr( 0, pos ) ;
	457	}
	458	else {
	459	pos = hAntennaName.find( "@" ) ;
	460	telescopeName = hAntennaName.substr( 0, pos ) ;
	461	}
	462	os_ << "telescopeName = " << telescopeName << LogIO::POST ;
	463	msObsCols.telescopeName().put( 0, telescopeName ) ;
	464	msObsCols.project().put( 0, header_.project ) ;
	465	//ScalarMeasColumn<MEpoch> timeCol( table_->table().sort("TIME"), "TIME" ) ;
	466	Table sortedtable = table_->table().sort("TIME") ;
	467	ScalarMeasColumn<MEpoch> timeCol( sortedtable, "TIME" ) ;
	468	Vector<MEpoch> trange( 2 ) ;
	469	trange[0] = timeCol( 0 ) ;
	470	trange[1] = timeCol( table_->nrow()-1 ) ;
	471	msObsCols.timeRangeMeas().put( 0, trange ) ;
	472	}
	473
	474	void MSWriter::fillAntenna()
	475	{
	476	os_ << "set up ANTENNA subtable" << LogIO::POST ;
	477	// only 1 row
	478	mstable_->antenna().addRow( 1, True ) ;
	479	MSAntennaColumns msAntCols( mstable_->antenna() ) ;
	480
	481	String hAntennaName = header_.antennaname ;
	482	String::size_type pos = hAntennaName.find( "//" ) ;
	483	String antennaName ;
	484	String stationName ;
	485	if ( pos != String::npos ) {
	486	hAntennaName = hAntennaName.substr( pos+2 ) ;
	487	}
	488	pos = hAntennaName.find( "@" ) ;
	489	if ( pos != String::npos ) {
	490	antennaName = hAntennaName.substr( 0, pos ) ;
	491	stationName = hAntennaName.substr( pos+1 ) ;
	492	}
	493	else {
	494	antennaName = hAntennaName ;
	495	stationName = hAntennaName ;
	496	}
	497	os_ << "antennaName = " << antennaName << LogIO::POST ;
	498	os_ << "stationName = " << stationName << LogIO::POST ;
	499
	500	msAntCols.name().put( 0, antennaName ) ;
	501	msAntCols.station().put( 0, stationName ) ;
	502
	503	os_ << "antennaPosition = " << header_.antennaposition << LogIO::POST ;
	504
	505	msAntCols.position().put( 0, header_.antennaposition ) ;
	506	}
	507
[1975]	508	void MSWriter::fillProcessor()
	509	{
	510	os_ << "set up PROCESSOR subtable" << LogIO::POST ;
	511
	512	// only add empty 1 row
	513	MSProcessor msProc = mstable_->processor() ;
	514	msProc.addRow( 1, True ) ;
	515	}
	516
	517	void MSWriter::fillSource()
	518	{
	519	os_ << "set up SOURCE subtable" << LogIO::POST ;
	520
	521
	522	// access to MS SOURCE subtable
	523	MSSource msSrc = mstable_->source() ;
	524
	525	// access to MOLECULE subtable
	526	STMolecules stm = table_->molecules() ;
	527
	528	Int srcId = 0 ;
	529	Vector<Double> restFreq ;
	530	Vector<String> molName ;
	531	Vector<String> fMolName ;
	532	//
	533	// ITERATION: SRCNAME
	534	//
	535	Int added0 = 0 ;
	536	Int current0 = msSrc.nrow() ;
	537	TableIterator iter0( table_->table(), "SRCNAME" ) ;
	538	while( !iter0.pastEnd() ) {
	539	Table t0( iter0.table() ) ;
	540
	541	// get necessary information
	542	ROScalarColumn<String> srcNameCol( t0, "SRCNAME" ) ;
	543	String srcName = srcNameCol( 0 ) ;
	544	ROArrayColumn<Double> sharedDArrRCol( t0, "SRCPROPERMOTION" ) ;
	545	Vector<Double> srcPM = sharedDArrRCol( 0 ) ;
	546	sharedDArrRCol.attach( t0, "SRCDIRECTION" ) ;
	547	Vector<Double> srcDir = sharedDArrRCol( 0 ) ;
	548	ROScalarColumn<Double> srcVelCol( t0, "SRCVELOCITY" ) ;
	549	Double srcVel = srcVelCol( 0 ) ;
	550
	551	//
	552	// ITERATION: MOLECULE_ID
	553	//
	554	Int added1 = 0 ;
	555	Int current1 = msSrc.nrow() ;
	556	TableIterator iter1( t0, "MOLECULE_ID" ) ;
	557	while( !iter1.pastEnd() ) {
	558	Table t1( iter1.table() ) ;
	559
	560	// get necessary information
	561	ROScalarColumn<uInt> molIdCol( t1, "MOLECULE_ID" ) ;
	562	uInt molId = molIdCol( 0 ) ;
	563	stm.getEntry( restFreq, molName, fMolName, molId ) ;
	564
	565	//
	566	// ITERATION: IFNO
	567	//
	568	Int added2 = 0 ;
	569	Int current2 = msSrc.nrow() ;
	570	TableIterator iter2( t1, "IFNO" ) ;
	571	while( !iter2.pastEnd() ) {
	572	Table t2( iter2.table() ) ;
	573	uInt prevnr = msSrc.nrow() ;
	574
	575	// get necessary information
	576	ROScalarColumn<uInt> ifNoCol( t2, "IFNO" ) ;
	577	uInt ifno = ifNoCol( 0 ) ; // IFNO = SPECTRAL_WINDOW_ID
	578	MEpoch midTimeMeas ;
	579	Double interval ;
	580	getValidTimeRange( midTimeMeas, interval, t2 ) ;
	581
	582	// add row
	583	msSrc.addRow( 1, True ) ;
	584
	585	// fill SPECTRAL_WINDOW_ID
	586	ScalarColumn<Int> sSpwIdCol( msSrc, "SPECTRAL_WINDOW_ID" ) ;
	587	sSpwIdCol.put( prevnr, ifno ) ;
	588
	589	// fill TIME, INTERVAL
	590	ScalarMeasColumn<MEpoch> sTimeMeasCol( msSrc, "TIME" ) ;
	591	sTimeMeasCol.put( prevnr, midTimeMeas ) ;
	592	ScalarColumn<Double> sIntervalCol( msSrc, "INTERVAL" ) ;
	593	sIntervalCol.put( prevnr, interval ) ;
	594
	595	added2++ ;
	596	iter2.next() ;
	597	}
	598
	599	// fill NUM_LINES, REST_FREQUENCY, TRANSITION, SYSVEL
	600	RefRows rows2( current2, current2+added2-1 ) ;
	601	uInt numFreq = restFreq.size() ;
	602	Vector<Int> numLines( added2, numFreq ) ;
	603	ScalarColumn<Int> numLinesCol( msSrc, "NUM_LINES" ) ;
	604	numLinesCol.putColumnCells( rows2, numLines ) ;
	605	if ( numFreq != 0 ) {
	606	Array<Double> rf( IPosition( 2, numFreq, added2 ) ) ;
	607	Array<String> trans( IPosition( 2, numFreq, added2 ) ) ;
	608	Array<Double> srcVelArr( IPosition( 2, numFreq, added2 ), srcVel ) ;
	609	for ( uInt ifreq = 0 ; ifreq < numFreq ; ifreq++ ) {
	610	Slicer slice( Slice(ifreq),Slice(0,added2,1) ) ;
	611	rf( slice ) = restFreq[ifreq] ;
	612	String transStr = fMolName[ifreq] ;
	613	if ( transStr.size() == 0 ) {
	614	transStr = molName[ifreq] ;
	615	}
	616	trans( slice ) = transStr ;
	617	}
	618	ArrayColumn<Double> sharedDArrCol( msSrc, "REST_FREQUENCY" ) ;
	619	sharedDArrCol.putColumnCells( rows2, rf ) ;
	620	sharedDArrCol.attach( msSrc, "SYSVEL" ) ;
	621	sharedDArrCol.putColumnCells( rows2, srcVelArr ) ;
	622	ArrayColumn<String> transCol( msSrc, "TRANSITION" ) ;
	623	transCol.putColumnCells( rows2, trans ) ;
	624	}
	625
	626	added1 += added2 ;
	627	iter1.next() ;
	628	}
	629
	630	// fill NAME, SOURCE_ID
	631	RefRows rows1( current1, current1+added1-1 ) ;
	632	Vector<String> nameArr( added1, srcName ) ;
	633	Vector<Int> srcIdArr( added1, srcId ) ;
	634	ScalarColumn<String> sNameCol( msSrc, "NAME" ) ;
	635	ScalarColumn<Int> srcIdCol( msSrc, "SOURCE_ID" ) ;
	636	sNameCol.putColumnCells( rows1, nameArr ) ;
	637	srcIdCol.putColumnCells( rows1, srcIdArr ) ;
	638
	639	// fill DIRECTION, PROPER_MOTION
	640	Array<Double> srcPMArr( IPosition( 2, 2, added1 ) ) ;
	641	Array<Double> srcDirArr( IPosition( 2, 2, added1 ) ) ;
	642	for ( uInt i = 0 ; i < 2 ; i++ ) {
	643	Slicer slice( Slice(i),Slice(0,added1,1) ) ;
	644	srcPMArr( slice ) = srcPM[i] ;
	645	srcDirArr( slice ) = srcDir[i] ;
	646	}
	647	ArrayColumn<Double> sharedDArrCol( msSrc, "DIRECTION" ) ;
	648	sharedDArrCol.putColumnCells( rows1, srcDirArr ) ;
	649	sharedDArrCol.attach( msSrc, "PROPER_MOTION" ) ;
	650	sharedDArrCol.putColumnCells( rows1, srcPMArr ) ;
	651
	652	// fill TIME, INTERVAL
	653
	654	// increment srcId if SRCNAME changed
	655	srcId++ ;
	656
	657	added0 += added1 ;
	658	iter0.next() ;
	659	}
	660	}
	661
[1974]	662	void MSWriter::addFeed( Int id )
	663	{
	664	os_ << "set up FEED subtable" << LogIO::POST ;
	665
	666	// add row
	667	MSFeed msFeed = mstable_->feed() ;
	668	msFeed.addRow( 1, True ) ;
	669	Int nrow = msFeed.nrow() ;
	670
	671	MSFeedColumns msFeedCols( mstable_->feed() ) ;
	672
	673	msFeedCols.feedId().put( nrow-1, id ) ;
	674	msFeedCols.antennaId().put( nrow-1, 0 ) ;
	675	}
	676
	677	void MSWriter::addSpectralWindow( Int spwid, Int freqid )
	678	{
	679	os_ << "set up SPECTRAL_WINDOW subtable" << LogIO::POST ;
	680
	681	// add row
	682	MSSpectralWindow msSpw = mstable_->spectralWindow() ;
	683	while( (Int)msSpw.nrow() <= spwid ) {
	684	msSpw.addRow( 1, True ) ;
	685	}
	686
	687	MSSpWindowColumns msSpwCols( msSpw ) ;
	688
	689	STFrequencies stf = table_->frequencies() ;
	690
	691	// MEAS_FREQ_REF
	692	msSpwCols.measFreqRef().put( spwid, stf.getFrame( True ) ) ;
	693
	694	Double refpix ;
	695	Double refval ;
	696	Double inc ;
	697	stf.getEntry( refpix, refval, inc, (uInt)freqid ) ;
	698
	699	// NUM_CHAN
	700	Int nchan = refpix * 2 ;
	701	msSpwCols.numChan().put( spwid, nchan ) ;
	702
	703	// TOTAL_BANDWIDTH
	704	Double bw = nchan * inc ;
	705	msSpwCols.totalBandwidth().put( spwid, bw ) ;
	706
	707	// REF_FREQUENCY
	708	Double refFreq = refval - refpix * inc ;
	709	msSpwCols.refFrequency().put( spwid, refFreq ) ;
	710
	711	// NET_SIDEBAND
	712	// tentative: USB->0, LSB->1
	713	Int netSideband = 0 ;
	714	if ( inc < 0 )
	715	netSideband = 1 ;
	716	msSpwCols.netSideband().put( spwid, netSideband ) ;
	717
	718	// RESOLUTION, CHAN_WIDTH, EFFECTIVE_BW
	719	Vector<Double> sharedDoubleArr( nchan, inc ) ;
	720	msSpwCols.resolution().put( spwid, sharedDoubleArr ) ;
	721	msSpwCols.chanWidth().put( spwid, sharedDoubleArr ) ;
	722	msSpwCols.effectiveBW().put( spwid, sharedDoubleArr ) ;
	723
	724	// CHAN_FREQ
	725	indgen( sharedDoubleArr, refFreq, inc ) ;
	726	msSpwCols.chanFreq().put( spwid, sharedDoubleArr ) ;
	727	}
	728
[1975]	729	void MSWriter::addField( Int fid, String fieldname, String srcname, Double t, Vector<Double> rate )
	730	{
	731	os_ << "set up FIELD subtable" << LogIO::POST ;
	732
	733	MSField msField = mstable_->field() ;
	734	while( (Int)msField.nrow() <= fid ) {
	735	msField.addRow( 1, True ) ;
	736	}
	737	MSFieldColumns msFieldCols( msField ) ;
	738
	739	// Access to SOURCE table
	740	MSSource msSrc = mstable_->source() ;
	741
	742	// fill target row
	743	msFieldCols.name().put( fid, fieldname ) ;
	744	msFieldCols.time().put( fid, t ) ;
	745	Int numPoly = 0 ;
	746	if ( anyNE( rate, 0.0 ) )
	747	numPoly = 1 ;
	748	msFieldCols.numPoly().put( fid, numPoly ) ;
	749	MSSourceIndex msSrcIdx( msSrc ) ;
	750	Int srcId = -1 ;
	751	Vector<Int> srcIdArr = msSrcIdx.matchSourceName( srcname ) ;
	752	if ( srcIdArr.size() != 0 ) {
	753	srcId = srcIdArr[0] ;
	754	MSSource msSrcSel = msSrc( msSrc.col("SOURCE_ID") == srcId ) ;
	755	ROMSSourceColumns msSrcCols( msSrcSel ) ;
	756	Vector<Double> srcDir = msSrcCols.direction()( 0 ) ;
	757	Array<Double> srcDirA( IPosition( 2, 2, 1+numPoly ) ) ;
	758	os_ << "srcDirA = " << srcDirA << LogIO::POST ;
	759	os_ << "sliced srcDirA = " << srcDirA( Slicer( Slice(0,2,1), Slice(0) ) ) << LogIO::POST ;
	760	srcDirA( Slicer( Slice(0,2,1), Slice(0) ) ) = srcDir.reform( IPosition(2,2,1) ) ;
	761	os_ << "srcDirA = " << srcDirA << LogIO::POST ;
	762	if ( numPoly != 0 )
	763	srcDirA( Slicer( Slice(0,2,1), Slice(1) ) ) = rate ;
	764	msFieldCols.phaseDir().put( fid, srcDirA ) ;
	765	msFieldCols.referenceDir().put( fid, srcDirA ) ;
	766	msFieldCols.delayDir().put( fid, srcDirA ) ;
	767	}
	768	msFieldCols.sourceId().put( fid, srcId ) ;
	769	}
	770
[1974]	771	Int MSWriter::addPolarization( Vector<Int> polnos )
	772	{
[1975]	773	os_ << "set up POLARIZATION subtable" << LogIO::POST ;
[1974]	774
	775	os_ << "polnos = " << polnos << LogIO::POST ;
	776	MSPolarization msPol = mstable_->polarization() ;
	777	uInt nrow = msPol.nrow() ;
	778
	779	Vector<Int> corrType = toCorrType( polnos ) ;
	780	os_ << "corrType = " << corrType << LogIO::POST ;
	781
	782	MSPolarizationIndex msPolIdx( msPol ) ;
	783	Vector<Int> polids = msPolIdx.matchCorrType( corrType ) ;
	784	os_ << "polids = " << polids << LogIO::POST ;
	785
	786	Int polid = -1 ;
	787
	788	if ( polids.size() == 0 ) {
	789	// add row
	790	msPol.addRow( 1, True ) ;
	791	polid = (Int)nrow ;
	792
	793	MSPolarizationColumns msPolCols( msPol ) ;
	794
	795	// CORR_TYPE
	796	msPolCols.corrType().put( nrow, corrType ) ;
	797
	798	// NUM_CORR
	799	uInt npol = corrType.size() ;
	800	msPolCols.numCorr().put( nrow, npol ) ;
	801
	802	// CORR_PRODUCT
	803	Matrix<Int> corrProd( 2, npol, -1 ) ;
	804	if ( npol == 1 ) {
	805	corrProd = 0 ;
	806	}
	807	else if ( npol == 2 ) {
[1975]	808	corrProd.column( 0 ) = 0 ;
	809	corrProd.column( 1 ) = 1 ;
[1974]	810	}
	811	else {
[1975]	812	corrProd.column( 0 ) = 0 ;
	813	corrProd.column( 3 ) = 1 ;
	814	corrProd( 0,1 ) = 0 ;
	815	corrProd( 1,1 ) = 1 ;
	816	corrProd( 0,2 ) = 1 ;
	817	corrProd( 1,2 ) = 0 ;
[1974]	818	}
	819	msPolCols.corrProduct().put( nrow, corrProd ) ;
	820
	821	}
	822	else {
	823	polid = polids[0] ;
	824	}
	825
	826	return polid ;
	827	}
	828
[1975]	829	Int MSWriter::addDataDescription( Int polid, Int spwid )
	830	{
	831	os_ << "set up SPECTRAL_WINDOW subtable" << LogIO::POST ;
	832
	833	MSDataDescription msDataDesc = mstable_->dataDescription() ;
	834	uInt nrow = msDataDesc.nrow() ;
	835
	836	MSDataDescIndex msDataDescIdx( msDataDesc ) ;
	837
	838	Vector<Int> ddids = msDataDescIdx.matchSpwIdAndPolznId( spwid, polid ) ;
	839	os_ << "ddids = " << ddids << LogIO::POST ;
	840
	841	Int ddid = -1 ;
	842	if ( ddids.size() == 0 ) {
	843	msDataDesc.addRow( 1, True ) ;
	844	MSDataDescColumns msDataDescCols( msDataDesc ) ;
	845	msDataDescCols.polarizationId().put( nrow, polid ) ;
	846	msDataDescCols.spectralWindowId().put( nrow, spwid ) ;
	847	ddid = (Int)nrow ;
	848	}
	849	else {
	850	ddid = ddids[0] ;
	851	}
	852
	853	return ddid ;
	854	}
	855
[1974]	856	Vector<Int> MSWriter::toCorrType( Vector<Int> polnos )
	857	{
	858	uInt npol = polnos.size() ;
	859	Vector<Int> corrType( npol, Stokes::Undefined ) ;
	860
	861	if ( npol == 4 ) {
	862	if ( polType_ == "linear" ) {
	863	for ( uInt ipol = 0 ; ipol < npol ; ipol++ ) {
	864	if ( polnos[ipol] == 0 )
	865	corrType[ipol] = Stokes::XX ;
	866	else if ( polnos[ipol] == 1 )
	867	corrType[ipol] = Stokes::XY ;
	868	else if ( polnos[ipol] == 2 )
	869	corrType[ipol] = Stokes::YX ;
	870	else if ( polnos[ipol] == 3 )
	871	corrType[ipol] = Stokes::YY ;
	872	}
	873	}
	874	else if ( polType_ == "circular" ) {
	875	for ( uInt ipol = 0 ; ipol < npol ; ipol++ ) {
	876	if ( polnos[ipol] == 0 )
	877	corrType[ipol] = Stokes::RR ;
	878	else if ( polnos[ipol] == 1 )
	879	corrType[ipol] = Stokes::RL ;
	880	else if ( polnos[ipol] == 2 )
	881	corrType[ipol] = Stokes::LR ;
	882	else if ( polnos[ipol] == 3 )
	883	corrType[ipol] = Stokes::LL ;
	884	}
	885	}
	886	else if ( polType_ == "stokes" ) {
	887	for ( uInt ipol = 0 ; ipol < npol ; ipol++ ) {
	888	if ( polnos[ipol] == 0 )
	889	corrType[ipol] = Stokes::I ;
	890	else if ( polnos[ipol] == 1 )
	891	corrType[ipol] = Stokes::Q ;
	892	else if ( polnos[ipol] == 2 )
	893	corrType[ipol] = Stokes::U ;
	894	else if ( polnos[ipol] == 3 )
	895	corrType[ipol] = Stokes::V ;
	896	}
	897	}
	898	}
	899	else if ( npol == 2 ) {
	900	if ( polType_ == "linear" ) {
	901	for ( uInt ipol = 0 ; ipol < npol ; ipol++ ) {
	902	if ( polnos[ipol] == 0 )
	903	corrType[ipol] = Stokes::XX ;
	904	else if ( polnos[ipol] == 1 )
	905	corrType[ipol] = Stokes::YY ;
	906	}
	907	}
	908	else if ( polType_ == "circular" ) {
	909	for ( uInt ipol = 0 ; ipol < npol ; ipol++ ) {
	910	if ( polnos[ipol] == 0 )
	911	corrType[ipol] = Stokes::RR ;
	912	else if ( polnos[ipol] == 1 )
	913	corrType[ipol] = Stokes::LL ;
	914	}
	915	}
	916	else if ( polType_ == "stokes" ) {
	917	for ( uInt ipol = 0 ; ipol < npol ; ipol++ ) {
	918	if ( polnos[ipol] == 0 )
	919	corrType[ipol] = Stokes::I ;
	920	else if ( polnos[ipol] == 1 )
	921	corrType[ipol] = Stokes::V ;
	922	}
	923	}
	924	else if ( polType_ == "linpol" ) {
	925	for ( uInt ipol = 0 ; ipol < npol ; ipol++ ) {
	926	if ( polnos[ipol] == 1 )
	927	corrType[ipol] = Stokes::Plinear ;
	928	else if ( polnos[ipol] == 2 )
	929	corrType[ipol] = Stokes::Pangle ;
	930	}
	931	}
	932	}
	933	else if ( npol == 1 ) {
	934	if ( polType_ == "linear" )
	935	corrType[0] = Stokes::XX ;
	936	else if ( polType_ == "circular" )
	937	corrType[0] = Stokes::RR ;
	938	else if ( polType_ == "stokes" )
	939	corrType[0] = Stokes::I ;
	940	}
	941
	942	return corrType ;
	943	}
	944
[1975]	945	void MSWriter::getValidTimeRange( MEpoch &me, Double &interval, Table &tab )
[1974]	946	{
[1975]	947	ROScalarMeasColumn<MEpoch> timeMeasCol( tab, "TIME" ) ;
	948	ROScalarColumn<Double> timeCol( tab, "TIME" ) ;
	949	String refStr = timeMeasCol( 0 ).getRefString() ;
	950	Vector<Double> timeArr = timeCol.getColumn() ;
	951	MEpoch::Types meType ;
	952	MEpoch::getType( meType, refStr ) ;
	953	Unit tUnit = timeMeasCol.measDesc().getUnits()( 0 ) ;
	954	Double minTime ;
	955	Double maxTime ;
	956	minMax( minTime, maxTime, timeArr ) ;
	957	Double midTime = 0.5 * ( minTime + maxTime ) ;
	958	me = MEpoch( Quantity( midTime, tUnit ), meType ) ;
	959	interval = Quantity( maxTime-minTime, tUnit ).getValue( "s" ) ;
	960	}
[1974]	961
	962	}

Note: See TracBrowser for help on using the repository browser.

Context Navigation

source: trunk/src/MSWriter.cpp@ 1976

Download in other formats: