Skip to content

Commit

Permalink
Add UVFITS support for large antenna numbers
Browse files Browse the repository at this point in the history 
Follow the scheme documented in AIPS memo 117 (revised, 19 June 2017)
and use separate SUBARRAY, ANTENNA1 and ANTENNA2 random parameters
when antenna numbers larger than 255 are encountered when writing
UVFITS files.

Add support for reading such UVFITS files back in as well.
  • Loading branch information
@kettenis
kettenis committed Nov 30, 2021
1 parent e793b3d commit f72e8f0
Show file tree
Hide file tree
Showing 3 changed files with 133 additions and 51 deletions.
77 changes: 64 additions & 13 deletions msfits/MSFits/MSFitsInput.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1138,6 +1138,10 @@ void MSFitsInput::fillMSMainTableColWise(Int& nField, Int& nSpW) {
}
// get index for baseline
Int iBsln = getIndex(pType, "BASELINE");
// get indices for subarray, antenna1 and antenna2
Int iSubarr = getIndex(pType, "SUBARRAY");
Int iAnt1 = getIndex(pType, "ANTENNA1");
Int iAnt2 = getIndex(pType, "ANTENNA2");
// get indices for time
Int iTime0 = getIndex(pType, "DATE", 0);
Int iTime1 = getIndex(pType, "DATE", 1);
Expand Down Expand Up @@ -1231,8 +1235,18 @@ void MSFitsInput::fillMSMainTableColWise(Int& nField, Int& nSpW) {
// make 0-based
fieldId = (Int) _priGroup.parm(iSource) - 1;
}
Float baseline = _priGroup.parm(iBsln);
Int arrayId = Int(100.0 * (baseline - Int(baseline) + 0.001));
Int arrayId = 0;
std::pair<Int, Int> ants;
if (iBsln >= 0) {
Float baseline = _priGroup.parm(iBsln);
ants = _extractAntennas(baseline);
arrayId = Int(100.0 * (baseline - Int(baseline) + 0.001));
} else {
Int antenna1 = _priGroup.parm(iAnt1);
Int antenna2 = _priGroup.parm(iAnt2);
ants = _extractAntennas(antenna1, antenna2);
arrayId = _priGroup.parm(iSubarr);
}
_nArray = max(_nArray, arrayId + 1);
for (Int k = 0; k < nif; ++k) {
Int index = group * nif + k;
Expand All @@ -1241,7 +1255,6 @@ void MSFitsInput::fillMSMainTableColWise(Int& nField, Int& nSpW) {
uvw(1, index) = _priGroup.parm(iV) * C::c;
uvw(2, index) = _priGroup.parm(iW) * C::c;
// Convert from units of seconds to meters
std::pair<Int, Int> ants = _extractAntennas(baseline);
ant1[index] = ants.first;
ant2[index] = ants.second;
}
Expand Down Expand Up @@ -1461,6 +1474,10 @@ void MSFitsInput::fillMSMainTable(Int& nField, Int& nSpW) {
}
// get index for baseline
Int iBsln = getIndex(pType, "BASELINE");
// get indices for subarray, antenna1 and antenna2
Int iSubarr = getIndex(pType, "SUBARRAY");
Int iAnt1 = getIndex(pType, "ANTENNA1");
Int iAnt2 = getIndex(pType, "ANTENNA2");
// get indices for time
Int iTime0 = getIndex(pType, "DATE", 0);
Int iTime1 = getIndex(pType, "DATE", 1);
Expand Down Expand Up @@ -1556,10 +1573,19 @@ void MSFitsInput::fillMSMainTable(Int& nField, Int& nSpW) {
uvw *= C::c;

// Extract array/baseline/antenna info
Float baseline = _priGroup.parm(iBsln);
Int arrayId = Int(100.0 * (baseline - Int(baseline) + 0.001));
Int arrayId = 0;
std::pair<Int, Int> ants;
if (iBsln >= 0) {
Float baseline = _priGroup.parm(iBsln);
ants = _extractAntennas(baseline);
arrayId = Int(100.0 * (baseline - Int(baseline) + 0.001));
} else {
Int antenna1 = _priGroup.parm(iAnt1);
Int antenna2 = _priGroup.parm(iAnt2);
ants = _extractAntennas(antenna1, antenna2);
arrayId = _priGroup.parm(iSubarr);
}
_nArray = max(_nArray, arrayId + 1);
std::pair<Int, Int> ants = _extractAntennas(baseline);
Int ant1 = ants.first;
Int ant2 = ants.second;
// Ensure arrayId-specific params are of correct length:
Expand Down Expand Up @@ -3196,6 +3222,9 @@ void MSFitsInput::fillMSMainTable(BinaryTable& bt) {
}

Int iBsln = getIndex(TType, "BASELINE");
Int iSubarr = getIndex(TType, "SUBARRAY");
Int iAnt1 = getIndex(TType, "ANTENNA1");
Int iAnt2 = getIndex(TType, "ANTENNA2");
Int iTime0 = getIndex(TType, "DATE");
Int iSource = getIndex(TType, "SOURCE");
Int iFreq = getIndex(TType, "FREQSEL");
Expand Down Expand Up @@ -3239,9 +3268,18 @@ void MSFitsInput::fillMSMainTable(BinaryTable& bt) {
ScalarColumn<Float> colUU(tb, TType(iU));
ScalarColumn<Float> colVV(tb, TType(iV));
ScalarColumn<Float> colWW(tb, TType(iW));
ScalarColumn<Float> colBL(tb, TType(iBsln));
ScalarColumn<Float> colBL;
ScalarColumn<Float> colSubarr, colAnt1, colAnt2;
ArrayColumn<Float> colVis(tb, TType(iVis));

if (iBsln >= 0) {
colBL = ScalarColumn<Float>(tb, TType(iBsln));
} else {
colSubarr = ScalarColumn<Float>(tb, TType(iSubarr));
colAnt1 = ScalarColumn<Float>(tb, TType(iAnt1));
colAnt2 = ScalarColumn<Float>(tb, TType(iAnt2));
}

Int visL = 1;
for (uInt i = 0; i < _nPixel.nelements(); i++)
visL *= _nPixel(i);
Expand Down Expand Up @@ -3275,10 +3313,19 @@ void MSFitsInput::fillMSMainTable(BinaryTable& bt) {
uvw *= C::c;

// Extract array/baseline/antenna info
Float baseline = colBL.asfloat(0);
Int arrayId = Int(100.0 * (baseline - Int(baseline) + 0.001));
Int arrayId;
std::pair<Int, Int> ants;
if (iBsln >= 0) {
Float baseline = colBL.asfloat(0);
ants = _extractAntennas(baseline);
arrayId = Int(100.0 * (baseline - Int(baseline) + 0.001));
} else {
Int antenna1 = _priGroup.parm(iAnt1);
Int antenna2 = _priGroup.parm(iAnt2);
ants = _extractAntennas(antenna1, antenna2);
arrayId = _priGroup.parm(iSubarr) - 1;
}
_nArray = max(_nArray, arrayId + 1);
std::pair<Int, Int> ants = _extractAntennas(baseline);
Int ant1 = ants.first;
Int ant2 = ants.second;

Expand Down Expand Up @@ -3459,9 +3506,7 @@ void MSFitsInput::fillMSMainTable(BinaryTable& bt) {
}
}

std::pair<Int, Int> MSFitsInput::_extractAntennas(Float baseline) {
Int ant1 = Int(baseline) / 256;
Int ant2 = Int(baseline) - ant1 * 256;
std::pair<Int, Int> MSFitsInput::_extractAntennas(Int ant1, Int ant2) {
_nAntRow = max(_nAntRow, ant1);
_nAntRow = max(_nAntRow, ant2);
_uniqueAnts.insert(ant1);
Expand All @@ -3472,6 +3517,12 @@ std::pair<Int, Int> MSFitsInput::_extractAntennas(Float baseline) {
return make_pair(ant1, ant2);
}

std::pair<Int, Int> MSFitsInput::_extractAntennas(Float baseline) {
Int ant1 = Int(baseline) / 256;
Int ant2 = Int(baseline) - ant1 * 256;
return _extractAntennas(ant1, ant2);
}

void MSFitsInput::fillObservationTable(ConstFitsKeywordList& kwl) {
const FitsKeyword* kw;
const Regex trailing(" *$"); // trailing blanks
Expand Down
3 changes: 2 additions & 1 deletion msfits/MSFits/MSFitsInput.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -421,8 +421,9 @@ class MSFitsInput
void readRandomGroupUVFits(Int obsType);
void readPrimaryTableUVFits(Int obsType);

std::pair<Int, Int> _extractAntennas(Int antenna1, Int antenna2);
std::pair<Int, Int> _extractAntennas(Float baseline);

void _fillSysPowerTable(BinaryTable& bt);

void _doFillSysPowerSingleIF(
Expand Down
104 changes: 67 additions & 37 deletions msfits/MSFits/MSFitsOutput.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -366,6 +366,23 @@ uInt MSFitsOutput::get_tbf_end(const uInt rownr, const uInt nrow,
return tbfend;
}

// Define a FITS random group parameter with default scaling and offset
static void defineRandomParam(Record& ek, Int n, const String& name) {
String ptype = "ptype" + String::toString(n);
String pscal = "pscal" + String::toString(n);
String pzero = "pzero" + String::toString(n);
ek.define(ptype, name);
ek.define(pscal, 1.0);
ek.define(pzero, 0.0);
}

static void defineRandomParam(Record& ek, Int n, const String& name,
const String& comment) {
String ptype = "ptype" + String::toString(n);
defineRandomParam(ek, n, name);
ek.setComment(ptype, comment);
}

std::shared_ptr<FitsOutput> MSFitsOutput::_writeMain(Int& refPixelFreq, Double& refFreq,
Double& chanbw, const String &outFITSFile,
const Block<Int>& spwidMap, Int nrspw,
Expand Down Expand Up @@ -477,6 +494,10 @@ std::shared_ptr<FitsOutput> MSFitsOutput::_writeMain(Int& refPixelFreq, Double&
}
}

Vector<Int> antnumbers;
_handleAntNumbers(_ms, antnumbers);
Int maxant = max(antnumbers);

// Also find out what the Stokes are and make sure that they are the same
// throughout the MS. In principle we could handle the same stokes in
// different order by transposing, but this may well never happen.
Expand Down Expand Up @@ -678,8 +699,16 @@ std::shared_ptr<FitsOutput> MSFitsOutput::_writeMain(Int& refPixelFreq, Double&
// DATE
desc.addField("date1", TpFloat);
desc.addField("date2", TpFloat);
// BASELINE
desc.addField("baseline", TpFloat);
// BASELINE if maximum antenna number < 256
// SUBARRAY, ANTENNA1 and ANTENNA2 otherwise
// (see AIPS memo 117, section 3.1.2)
if (maxant < 256) {
desc.addField("baseline", TpFloat);
} else {
desc.addField("subarray", TpFloat);
desc.addField("antenna1", TpFloat);
desc.addField("antenna2", TpFloat);
}
// FREQSEL
ScalarColumn<Int> inddid(_ms, MS::columnName(MS::DATA_DESC_ID));
desc.addField("freqsel", TpFloat);
Expand Down Expand Up @@ -748,36 +777,23 @@ std::shared_ptr<FitsOutput> MSFitsOutput::_writeMain(Int& refPixelFreq, Double&
ek.define("crota7", 0.0);

// PTYPE PSCALE PZERO
ek.define("ptype1", "UU");
ek.define("pscal1", 1.0);
ek.define("pzero1", 0.0);
ek.define("ptype2", "VV");
ek.define("pscal2", 1.0);
ek.define("pzero2", 0.0);
ek.define("ptype3", "WW");
ek.define("pscal3", 1.0);
ek.define("pzero3", 0.0);
ek.define("ptype4", "DATE");
ek.define("pscal4", 1.0);
ek.define("pzero4", 0.0);
ek.setComment("ptype4", "Day number");
ek.define("ptype5", "DATE");
ek.define("pscal5", 1.0);
ek.define("pzero5", 0.0);
ek.setComment("ptype5", "Day fraction");
ek.define("ptype6", "BASELINE");
ek.define("pscal6", 1.0);
ek.define("pzero6", 0.0);
ek.define("ptype7", "FREQSEL");
ek.define("pscal7", 1.0);
ek.define("pzero7", 0.0);
Int idx = 1;
defineRandomParam(ek, idx++, "UU");
defineRandomParam(ek, idx++, "VV");
defineRandomParam(ek, idx++, "WW");
defineRandomParam(ek, idx++, "DATE", "Day number");
defineRandomParam(ek, idx++, "DATE", "Day fraction");
if (maxant < 256) {
defineRandomParam(ek, idx++, "BASELINE");
} else {
defineRandomParam(ek, idx++, "SUBARRAY");
defineRandomParam(ek, idx++, "ANTENNA1");
defineRandomParam(ek, idx++, "ANTENNA2");
}
defineRandomParam(ek, idx++, "FREQSEL");
if (asMultiSource) {
ek.define("ptype8", "SOURCE");
ek.define("pscal8", 1.0);
ek.define("pzero8", 0.0);
ek.define("ptype9", "INTTIM");
ek.define("pscal9", 1.0);
ek.define("pzero9", 0.0);
defineRandomParam(ek, idx++, "SOURCE");
defineRandomParam(ek, idx++, "INTTIM");
}

// EXTEND - already written by FITSGroupWriter
Expand Down Expand Up @@ -1046,7 +1062,17 @@ std::shared_ptr<FitsOutput> MSFitsOutput::_writeMain(Int& refPixelFreq, Double&
RecordFieldPtr<Float> oww(writer.row(), "w");
RecordFieldPtr<Float> odate1(writer.row(), "date1");
RecordFieldPtr<Float> odate2(writer.row(), "date2");
RecordFieldPtr<Float> obaseline(writer.row(), "baseline");
RecordFieldPtr<Float> obaseline;
RecordFieldPtr<Float> osubarray;
RecordFieldPtr<Float> oantenna1;
RecordFieldPtr<Float> oantenna2;
if (maxant < 256) {
obaseline = RecordFieldPtr<Float> (writer.row(), "baseline");
} else {
osubarray = RecordFieldPtr<Float> (writer.row(), "subarray");
oantenna1 = RecordFieldPtr<Float> (writer.row(), "antenna1");
oantenna2 = RecordFieldPtr<Float> (writer.row(), "antenna2");
}
RecordFieldPtr<Float> ofreqsel(writer.row(), "freqsel");
RecordFieldPtr<Float> osource;
RecordFieldPtr<Float> ointtim;
Expand All @@ -1065,9 +1091,6 @@ std::shared_ptr<FitsOutput> MSFitsOutput::_writeMain(Int& refPixelFreq, Double&
}
}

Vector<Int> antnumbers;
_handleAntNumbers(_ms, antnumbers);

// Loop through all rows.
ProgressMeter meter(0.0, nOutRow * 1.0, "UVFITS Writer", "Rows copied", "",
"", True, nOutRow / 100);
Expand Down Expand Up @@ -1343,8 +1366,15 @@ std::shared_ptr<FitsOutput> MSFitsOutput::_writeMain(Int& refPixelFreq, Double&
*odate2 = dayFraction;

// BASELINE
*obaseline = antnumbers(inant1(tbfrownr)) * 256 + antnumbers(inant2(
tbfrownr)) + inarray(tbfrownr) * 0.01;
if (maxant < 256) {
*obaseline = antnumbers(inant1(tbfrownr)) * 256 +
antnumbers(inant2(tbfrownr)) +
inarray(tbfrownr) * 0.01;
} else {
*osubarray = inarray(tbfrownr) + 1;
*oantenna1 = antnumbers(inant1(tbfrownr));
*oantenna2 = antnumbers(inant2(tbfrownr));
}

// FREQSEL (in the future it might be FREQ_GRP+1)
// *ofreqsel = inddid(i) + 1;
Expand Down

0 comments on commit f72e8f0

Please sign in to comment.