Merge branch 'pan-195-dp3-pol-info' into 'master'

PAN-195: Add polarization information to DPInfo

See merge request !1250

CMakeLists.txt

@@ -773,6 +773,7 @@ if(BUILD_TESTING)
       steps/test/unit/tMSBDAReader.cc
       steps/test/unit/tMSBDAWriter.cc
       steps/test/unit/tMsColumnReader.cc
+      steps/test/unit/tMSReader.cc
       steps/test/unit/tMSUpdater.cc
       steps/test/unit/tMSWriter.cc
       steps/test/unit/tNullStokes.cc

base/DPInfo.cc

@@ -54,7 +54,8 @@ DPInfo::DPInfo(unsigned int n_correlations, unsigned int original_n_channels,
       resolutions_(1),          // can retrieve a first list.
       effective_bandwidth_(1),
       total_bandwidth_(0.0),
-      spectral_window_(0) {}
+      spectral_window_(0),
+      polarizations_() {}
 
 void DPInfo::setTimes(double first_time, double last_time,
                       double time_interval) {

include/dp3/base/DPInfo.h

@@ -9,6 +9,10 @@
 #ifndef DP3_BASE_DPINFO_H_
 #define DP3_BASE_DPINFO_H_
 
+#include <set>
+
+#include <aocommon/polarization.h>
+
 #include <casacore/measures/Measures/MDirection.h>
 #include <casacore/measures/Measures/MPosition.h>
 #include <casacore/measures/Measures/MeasureHolder.h>
@@ -129,6 +133,11 @@ class DPInfo {
     phase_center_ = phase_center;
   }
 
+  void setPolarizations(
+      const std::set<aocommon::PolarizationEnum>& polarizations) {
+    polarizations_ = polarizations;
+  }
+
   /// Get the info.
   ///@{
   const std::string& msName() const { return ms_name_; }
@@ -258,6 +267,10 @@ class DPInfo {
   /// Determine if the channels have a regular layout.
   bool channelsAreRegular() const;
 
+  const std::set<aocommon::PolarizationEnum>& polarizations() const {
+    return polarizations_;
+  }
+
  private:
   /// Set which antennae are actually used.
   void setAntUsed();
@@ -317,6 +330,7 @@ class DPInfo {
   mutable std::vector<double> baseline_lengths_;
   /// For each antenna, the auto correlation index.
   mutable std::vector<int> auto_correlation_indices_;
+  std::set<aocommon::PolarizationEnum> polarizations_;
 };
 
 }  // namespace base

pythondp3/PyDpInfo.cc

@@ -173,7 +173,10 @@ The start time equals the first time minus half a time interval)")
 
       /* Other properties */
       .def_property("ms_name", &DPInfo::msName, &DPInfo::setMsName,
-                    "The name of the measurement set");
+                    "The name of the measurement set")
+      .def_property("polarizations", &DPInfo::polarizations,
+                    &DPInfo::setPolarizations,
+                    "Polarizations of the measurement set");
 }
 
 }  // namespace pythondp3

steps/MSReader.cc

@@ -686,6 +686,29 @@ void MSReader::prepare2(int spectralWindow) {
                        std::vector<double>(effbwBegin, effbwBegin + itsNrChan),
                        refFreq, spectralWindow);
   }
+
+  std::set<aocommon::PolarizationEnum> polarizations;
+  casacore::MSDataDescription data_description_table = itsMS.dataDescription();
+  casacore::ScalarColumn<int> polarization_index_column(
+      data_description_table,
+      casacore::MSDataDescription::columnName(
+          casacore::MSDataDescription::POLARIZATION_ID));
+  const size_t polarization_index = polarization_index_column(spectralWindow);
+  // Populate the polarization
+  casacore::MSPolarization pol_table = itsMS.polarization();
+  casacore::ArrayColumn<int> corr_type_column(
+      pol_table, casacore::MSPolarization::columnName(
+                     casacore::MSPolarizationEnums::CORR_TYPE));
+  casacore::Array<int> corr_type_vec(corr_type_column(polarization_index));
+  for (casacore::Array<int>::const_contiter p = corr_type_vec.cbegin();
+       p != corr_type_vec.cend(); ++p) {
+    polarizations.emplace(aocommon::Polarization::AipsIndexToEnum(*p));
+  }
+  if (polarizations.size() != 4) {
+    throw std::runtime_error(
+        "DP3 expects a measurement set with 4 polarizations");
+  }
+  info().setPolarizations(polarizations);
 }
 
 void MSReader::skipFirstTimes() {

steps/test/unit/tMSReader.cc

+// Copyright (C) 2024 ASTRON (Netherlands Institute for Radio Astronomy)
+// SPDX-License-Identifier: GPL-3.0-or-later
+
+#include "../../MSReader.h"
+
+#include <memory>
+
+#include <boost/test/unit_test.hpp>
+
+#include <casacore/tables/TaQL/TableParse.h>
+#include <casacore/tables/Tables/TableCopy.h>
+
+#include "../../../common/ParameterSet.h"
+#include "../../../common/test/unit/fixtures/fDirectory.h"
+
+using casacore::Table;
+
+using dp3::common::ParameterSet;
+using dp3::common::test::FixtureDirectory;
+using dp3::steps::MSReader;
+
+const std::string kInputMs = "../tNDPPP_tmp.MS";
+const std::string kCopyMs = "tNDPPP_tmp.copy.MS";
+const std::string kCopyMsPol = "tNDPPP_tmp.copy.MS/POLARIZATION";
+
+BOOST_AUTO_TEST_SUITE(msreader)
+
+class FixtureCopyAndUpdatePol : FixtureDirectory {
+ public:
+  FixtureCopyAndUpdatePol() : FixtureDirectory() {
+    Table(kInputMs).deepCopy(kCopyMs, Table::New);
+    casacore::tableCommand("update " + kCopyMsPol +
+                           " set CORR_TYPE=[9, 12], "
+                           "CORR_PRODUCT=[[0,0],[1,1]], NUM_CORR=2");
+  }
+};
+
+BOOST_FIXTURE_TEST_CASE(polarization_initialization, FixtureDirectory) {
+  const casacore::MeasurementSet ms(kInputMs,
+                                    casacore::TableLock::AutoNoReadLocking);
+  ParameterSet parset;
+  parset.add("msin", kInputMs);
+  const MSReader reader(ms, parset, "msin.");
+
+  const std::set<aocommon::PolarizationEnum> expected_polarizations{
+      aocommon::PolarizationEnum::XX, aocommon::PolarizationEnum::XY,
+      aocommon::PolarizationEnum::YX, aocommon::PolarizationEnum::YY};
+
+  const std::set<aocommon::PolarizationEnum> actual_polarizations =
+      reader.getInfo().polarizations();
+
+  BOOST_CHECK_EQUAL_COLLECTIONS(
+      actual_polarizations.begin(), actual_polarizations.end(),
+      expected_polarizations.begin(), expected_polarizations.end());
+}
+
+BOOST_FIXTURE_TEST_CASE(expect_four_polarizations, FixtureCopyAndUpdatePol) {
+  casacore::MeasurementSet ms(kCopyMs, casacore::TableLock::AutoNoReadLocking);
+  ParameterSet parset;
+  parset.add("msin", kCopyMs);
+
+  BOOST_CHECK_THROW(std::make_unique<MSReader>(ms, parset, "msin."),
+                    std::runtime_error);
+}
+
+BOOST_AUTO_TEST_SUITE_END()