diff --git a/CEP/MS/include/MS/BaselineSelect.h b/CEP/MS/include/MS/BaselineSelect.h index 80daa0578712324f07b0451ef34a4a6329b7a42b..540ab928f70f1cd81a05cad41cb1045d653770cb 100644 --- a/CEP/MS/include/MS/BaselineSelect.h +++ b/CEP/MS/include/MS/BaselineSelect.h @@ -29,11 +29,18 @@ //# Includes #include <Common/lofar_string.h> +#include <ms/MSSel/MSSelectionErrorHandler.h> +#include <ostream> + //# Forward Declarations namespace casa { template<class T> class Matrix; + template<class T> class Vector; + class Table; + class TableExprNode; + class MPosition; } namespace LOFAR @@ -51,10 +58,48 @@ class BaselineSelect public: // Parse the MSSelection baseline string and create a Matrix telling // which baselines are selected. + // Possible messages from the parser are written to the ostream. static casa::Matrix<bool> convert (const string& msName, - const string& baselineSelection); + const string& baselineSelection, + std::ostream&); + + // Parse the MSSelection baseline string and create a Matrix telling + // which baselines are selected. + // The input is a vector of station names and positions. + // Possible messages from the parser are written to the ostream. + static casa::Matrix<bool> convert (const casa::Vector<casa::String>& names, + const std::vector<casa::MPosition>& pos, + const casa::Vector<casa::Int>& ant1, + const casa::Vector<casa::Int>& ant2, + const string& baselineSelection, + std::ostream&); + +private: + static casa::Matrix<bool> convert (casa::Table& anttab, + casa::TableExprNode& a1, + casa::TableExprNode& a2, + const string& baselineSelection, + std::ostream& os); + }; + + +// This class handles an error from the Casacore's MSAntennaParse. +// It adds the message to the message list of the parent BaselineSelect. +class BaselineSelectErrorHandler : public casa::MSSelectionErrorHandler +{ +public: + BaselineSelectErrorHandler (std::ostream& os) + : itsStream (os) + {} + virtual ~BaselineSelectErrorHandler(); + virtual void reportError (const char *token, const casa::String message); +private: + std::ostream& itsStream; +}; + + // @} } // end namespace diff --git a/CEP/MS/include/MS/MSCreate.h b/CEP/MS/include/MS/MSCreate.h index e91a1f665e5b18b1a6168ca7f08ebd864aeec1b0..f35ce3dda2a2a84b4b7a93b323155aa0fa47b402 100644 --- a/CEP/MS/include/MS/MSCreate.h +++ b/CEP/MS/include/MS/MSCreate.h @@ -84,6 +84,10 @@ public: // Destructor ~MSCreate(); + // Close all subtables to limit the nr of open files. + // This can be done once all subtables have been written. + void closeSubTables(); + // Add the extra columns needed for lwimager for every column not existing. // These are CORRECTED_DATA, MODEL_DATA, and IMAGING_WEIGHT. // Furthermore it sets the CHANNEL_SELECTION keyword for casa::VisSet. diff --git a/CEP/MS/src/BaselineSelect.cc b/CEP/MS/src/BaselineSelect.cc index 8015a4de6833abf2149d2ada2c0f5c70da7dd16f..123085347349b3b7387646dddd6837dc8bed9207 100644 --- a/CEP/MS/src/BaselineSelect.cc +++ b/CEP/MS/src/BaselineSelect.cc @@ -22,23 +22,38 @@ //# Includes #include <MS/BaselineSelect.h> +#include <Common/LofarLogger.h> + #include <ms/MeasurementSets/MeasurementSet.h> +#include <ms/MeasurementSets/MSAntenna.h> +#include <ms/MeasurementSets/MSAntennaColumns.h> #if defined(casacore) #include <ms/MSSel/MSSelection.h> +#include <ms/MSSel/MSAntennaParse.h> +#include <ms/MSSel/MSAntennaGram.h> #else #include <ms/MeasurementSets/MSSelection.h> +#include <ms/MeasurementSets/MSAntennaParse.h> +#include <ms/MeasurementSets/MSAntennaGram.h> #endif +#include <tables/Tables/Table.h> +#include <tables/Tables/SetupNewTab.h> +#include <tables/Tables/TableRecord.h> #include <tables/Tables/TableParse.h> #include <tables/Tables/ScalarColumn.h> +#include <tables/Tables/ScaColDesc.h> +#include <measures/Measures/MPosition.h> #include <casa/Arrays/Matrix.h> #include <casa/Arrays/Vector.h> + using namespace casa; namespace LOFAR { Matrix<bool> BaselineSelect::convert (const string& msName, - const string& baselineSelection) + const string& baselineSelection, + std::ostream& os) { // Find the unique baselines in the MS. // Do not use unique sort, because that is slow for a large MS. @@ -60,23 +75,91 @@ namespace LOFAR { } bltab = tab(Vector<uInt>(rows)); } - MeasurementSet ms(bltab); - MSSelection select; - // Set given selection strings. - select.setAntennaExpr (baselineSelection); - // Create a table expression over a MS representing the selection - TableExprNode node = select.toTableExprNode (&ms); - Table seltab = ms(node); - // Get the antenna numbers. - Vector<Int> a1 = ROScalarColumn<Int>(seltab, "ANTENNA1").getColumn(); - Vector<Int> a2 = ROScalarColumn<Int>(seltab, "ANTENNA2").getColumn(); - int nant = ms.antenna().nrow(); - Matrix<bool> bl(nant, nant, false); - for (uint i=0; i<a1.size(); ++i) { - bl(a1[i], a2[i]) = true; - bl(a2[i], a1[i]) = true; + TableExprNode a1 (bltab.col("ANTENNA1")); + TableExprNode a2 (bltab.col("ANTENNA2")); + Table anttab (bltab.keywordSet().asTable("ANTENNA")); + return convert (anttab, a1, a2, baselineSelection, os); + } + + casa::Matrix<bool> BaselineSelect::convert (const Vector<String>& names, + const vector<MPosition>& pos, + const Vector<Int>& ant1, + const Vector<Int>& ant2, + const string& baselineSelection, + std::ostream& os) + { + ASSERT (names.size() == pos.size()); + // Create a temporary MSAntenna table in memory for parsing purposes. + SetupNewTable antNew(String(), MSAntenna::requiredTableDesc(), + Table::New); + Table anttab(antNew, Table::Memory, names.size()); + MSAntenna msant(anttab); + MSAntennaColumns antcol(msant); + antcol.name().putColumn (names); + for (size_t i=0; i<pos.size(); ++i) { + antcol.positionMeas().put (i, pos[i]); } - return bl; + // Create a temporary table holding the antenna numbers of the baselines. + TableDesc td; + td.addColumn (ScalarColumnDesc<Int>("ANTENNA1")); + td.addColumn (ScalarColumnDesc<Int>("ANTENNA2")); + SetupNewTable tabNew(String(), td, Table::New); + Table tab(tabNew, Table::Memory, ant1.size()); + ScalarColumn<Int> ac1(tab, "ANTENNA1"); + ScalarColumn<Int> ac2(tab, "ANTENNA2"); + ac1.putColumn (ant1); + ac2.putColumn (ant2); + // Do the selection using the temporary tables. + TableExprNode a1 (tab.col("ANTENNA1")); + TableExprNode a2 (tab.col("ANTENNA2")); + return convert (anttab, a1, a2, baselineSelection, os); + } + + Matrix<bool> BaselineSelect::convert (Table& anttab, + TableExprNode& a1, + TableExprNode& a2, + const string& baselineSelection, + std::ostream& os) + { + // Overwrite the error handler to ignore errors for unknown antennas. + // First construct MSSelection, because it resets the error handler. + Vector<Int> selectedAnts1; + Vector<Int> selectedAnts2; + Matrix<Int> selectedBaselines; + MSSelectionErrorHandler* curHandler = MSAntennaParse::thisMSAErrorHandler; + BaselineSelectErrorHandler errorHandler (os); + MSAntennaParse::thisMSAErrorHandler = &errorHandler; + try { + // Create a table expression representing the selection. + TableExprNode node = msAntennaGramParseCommand + (anttab, a1, a2, baselineSelection, + selectedAnts1, selectedAnts2, selectedBaselines); + // Get the antenna numbers. + Table seltab = node.table()(node); + Vector<Int> a1 = ROScalarColumn<Int>(seltab, "ANTENNA1").getColumn(); + Vector<Int> a2 = ROScalarColumn<Int>(seltab, "ANTENNA2").getColumn(); + int nant = anttab.nrow(); + Matrix<bool> bl(nant, nant, false); + for (uint i=0; i<a1.size(); ++i) { + bl(a1[i], a2[i]) = true; + bl(a2[i], a1[i]) = true; + } + MSAntennaParse::thisMSAErrorHandler = curHandler; + return bl; + } catch (const std::exception&) { + MSAntennaParse::thisMSAErrorHandler = curHandler; + throw; + } + } + + + BaselineSelectErrorHandler::~BaselineSelectErrorHandler() + {} + + void BaselineSelectErrorHandler::reportError (const char* token, + const String message) + { + itsStream << message << token << endl; } } // end namespace diff --git a/CEP/MS/src/MSCreate.cc b/CEP/MS/src/MSCreate.cc index 95028ba30e66aabdd35f9bb4234a76398300fb42..b0aca89eb9f578a725f076e6101b6456a0e51a54 100644 --- a/CEP/MS/src/MSCreate.cc +++ b/CEP/MS/src/MSCreate.cc @@ -263,6 +263,28 @@ void MSCreate::createMS (const String& msName, } } +void MSCreate::closeSubTables() +{ + itsMS->antenna() = MSAntenna(); + itsMS->dataDescription() = MSDataDescription(); + itsMS->doppler() = MSDoppler(); + itsMS->feed() = MSFeed(); + itsMS->field() = MSField(); + itsMS->flagCmd() = MSFlagCmd(); + itsMS->freqOffset() = MSFreqOffset(); + itsMS->history() = MSHistory(); + itsMS->observation() = MSObservation(); + itsMS->pointing() = MSPointing(); + itsMS->polarization() = MSPolarization(); + itsMS->processor() = MSProcessor(); + itsMS->source() = MSSource(); + itsMS->spectralWindow() = MSSpectralWindow(); + itsMS->state() = MSState(); + itsMS->sysCal() = MSSysCal(); + itsMS->weather() = MSWeather(); + itsMS->keywordSet().closeTables(); +} + int MSCreate::addBand (int npolarizations, int nchannels, double refFreq, double chanWidth) { @@ -305,7 +327,7 @@ int MSCreate::addBand (int npolarizations, int nchannels, if (polnr < 0) { polnr = addPolarization (npolarizations); } - // Add a row to the DATADESCRIPTION subtable. + // Add a row to the DATA_DESCRIPTION subtable. MSDataDescription msdd = itsMS->dataDescription(); MSDataDescColumns msddCol(msdd); uInt rownr = msdd.nrow(); @@ -313,7 +335,7 @@ int MSCreate::addBand (int npolarizations, int nchannels, msddCol.spectralWindowId().put (rownr, rownr); msddCol.polarizationId().put (rownr, polnr); msddCol.flagRow().put (rownr, False); - // Add a row to the SPECTRALWINDOW subtable. + // Add a row to the SPECTRAL_WINDOW subtable. // Find the total bandwidth from the minimum and maximum. Vector<double> stFreqs = chanFreqs - chanWidths/2.; Vector<double> endFreqs = chanFreqs + chanWidths/2.; @@ -588,7 +610,7 @@ void MSCreate::updateTimes() Double midTime = (itsStartTime + endTime) / 2; // Update all rows in FEED subtable. { - MSFeed mssub = itsMS->feed(); + MSFeed mssub (itsMS->keywordSet().asTable("FEED")); MSFeedColumns mssubCol(mssub); Vector<Double> val(mssub.nrow()); val = midTime; @@ -598,7 +620,7 @@ void MSCreate::updateTimes() } // Update all rows in POINTING subtable. { - MSPointing mssub = itsMS->pointing(); + MSPointing mssub (itsMS->keywordSet().asTable("POINTING")); MSPointingColumns mssubCol(mssub); Vector<Double> val(mssub.nrow()); val = midTime; @@ -608,7 +630,7 @@ void MSCreate::updateTimes() } // Update all rows in OBSERVATION subtable. { - MSObservation msobs = itsMS->observation(); + MSObservation msobs (itsMS->keywordSet().asTable("OBSERVATION")); MSObservationColumns msobsCol(msobs); Vector<Double> timeRange(2); timeRange(0) = itsStartTime; @@ -843,6 +865,10 @@ void MSCreate::addImagerColumns (MeasurementSet& ms) ms.addColumn (td, stMan); // Set MODEL_DATA keyword for casa::VisSet. // Sort out the channel selection. + if (ms.spectralWindow().isNull()) { + ms.spectralWindow() = + MSSpectralWindow(ms.keywordSet().asTable("SPECTRAL_WINDOW")); + } MSSpWindowColumns msSpW(ms.spectralWindow()); Matrix<Int> selection(2, msSpW.nrow()); // Fill in default selection (all bands and channels). diff --git a/CEP/MS/src/makems.cc b/CEP/MS/src/makems.cc index 7f9b48afa3ca57e900483a36b53e98a53531dad3..2989ace240fe35eea19779a9839105edd49a99a8 100644 --- a/CEP/MS/src/makems.cc +++ b/CEP/MS/src/makems.cc @@ -37,6 +37,7 @@ #include <casa/Arrays/Matrix.h> #include <casa/Arrays/Vector.h> #include <casa/OS/Path.h> +#include <casa/OS/Timer.h> #include <tables/Tables/Table.h> #include <tables/Tables/ScalarColumn.h> #include <tables/Tables/ArrayColumn.h> @@ -173,6 +174,7 @@ void readParms (const string& parset) void createMS (int nband, int bandnr, const string& msName) { + Timer timer; int nfpb = itsNFreq/itsNBand; MSCreate msmaker(msName, itsStartTime, itsStepTime, nfpb, itsNCorr, itsAntPos, itsAntennaTableName, itsWriteAutoCorr, @@ -187,6 +189,9 @@ void createMS (int nband, int bandnr, const string& msName) for (uint i=0; i<itsRa.size(); ++i) { msmaker.addField (itsRa[i], itsDec[i]); } + // Close all subtables to reduce nr of open files. + msmaker.closeSubTables(); + timer.show ("Create MS"); for (int i=0; i<itsNTime; ++i) { msmaker.writeTimeStep(); } diff --git a/CEP/MS/src/makems.cfg b/CEP/MS/src/makems.cfg index 7330fa8aae890e24002ca7ea7539ff5b8608f926..8cd40b02cf9ad8850e8d9f5d416b78fb5c9cd11d 100644 --- a/CEP/MS/src/makems.cfg +++ b/CEP/MS/src/makems.cfg @@ -7,13 +7,13 @@ Declination=62.34.44.313606568 NBands=4 NFrequencies=64 NTimes=14 -NParts=2 +NParts=4 Dirs=["node1:abc", 'node1:def'] TileSizeFreq=8 TileSizeRest=10 WriteAutoCorr=T -AntennaTableName=~/WSRT_ANTENNA -MSName=tst/test.MS +AntennaTableName=~/data/WSRT_ANTENNA +MSName=test.MS VDSPath=. FlagColumn=LOFAR_FLAG NFlagBits=8 diff --git a/CEP/MS/src/msoverview.cc b/CEP/MS/src/msoverview.cc index d4d110fdf011724b8183bd6a360df4d34c3d8cd0..8d2e3df428d3e08352743a599318025330e5f269 100644 --- a/CEP/MS/src/msoverview.cc +++ b/CEP/MS/src/msoverview.cc @@ -112,7 +112,8 @@ int main (int argc, char* argv[]) logio << LogIO::NORMAL << endl << LogIO::POST; summ.listField (logio, False); logio << LogIO::NORMAL << endl << LogIO::POST; - summ.listSpectralAndPolInfo (logio, verbose); + ///summ.listSpectralAndPolInfo (logio, verbose, False); // new version + summ.listSpectralAndPolInfo (logio, verbose); // old version if (verbose) { logio << LogIO::NORMAL << endl << LogIO::POST; summ.listAntenna (logio, True); diff --git a/CEP/MS/test/tBaselineSelect.cc b/CEP/MS/test/tBaselineSelect.cc index b335d4696139c19e8c76cf23fe0e984edfe62389..29e6089085bf40b18cdf668304571d19b138382e 100644 --- a/CEP/MS/test/tBaselineSelect.cc +++ b/CEP/MS/test/tBaselineSelect.cc @@ -22,21 +22,60 @@ //# Includes #include <MS/BaselineSelect.h> +#include <Common/LofarLogger.h> + +#include <measures/Measures/MPosition.h> +#include <measures/Measures/MeasTable.h> #include <casa/Arrays/ArrayIO.h> +#include <casa/BasicSL/STLIO.h> #include <iostream> using namespace LOFAR; using namespace casa; using namespace std; +void testTemp() +{ + Vector<String> names(5); + names[0] = "CS013HBA0"; + names[1] = "CS013HBA1"; + names[2] = "CS014HBA0"; + names[3] = "RS015"; + names[4] = "DE013"; + Vector<Int> a1(15); + Vector<Int> a2(15); + int inx=0; + for (int i=0; i<5; ++i) { + for (int j=i; j<5; ++j) { + a1[inx] = i; + a2[inx] = j; + inx++; + } + } + vector<MPosition> pos(5); + ASSERT (MeasTable::Observatory (pos[0], "WSRT")); + ASSERT (MeasTable::Observatory (pos[1], "WSRT")); + ASSERT (MeasTable::Observatory (pos[2], "LOFAR")); + ASSERT (MeasTable::Observatory (pos[3], "WSRT")); + ASSERT (MeasTable::Observatory (pos[4], "LOFAR")); + // Try a few selection strings. + cout << BaselineSelect::convert (names, pos, a1, a2, "CS013HBA[01]", cout); + cout << BaselineSelect::convert (names, pos, a1, a2, "CS013HBA[23]", cout); + cout << BaselineSelect::convert (names, pos, a1, a2, "CS*&&[CR]S*", cout); + cout << BaselineSelect::convert (names, pos, a1, a2, "<10000", cout); +} + int main(int argc, char* argv[]) { try { - if (argc < 3) { + if (argc == 1) { + testTemp(); + } else if (argc < 3) { cout << "Run as: tBaselineSelect msname expr" << endl; return 0; + } else { + cout << BaselineSelect::convert (argv[1], argv[2], cout); } - cout << BaselineSelect::convert (argv[1], argv[2]); } catch (exception& x) { cout << "Unexpected expection: " << x.what() << endl; return 1;