source: trunk/src/STFrequencies.cpp@ 3003

Last change on this file since 3003 was 2900, checked in by Takeshi Nakazato, 11 years ago

New Development: No

JIRA Issue: Yes CAS-5875

Ready for Test: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: test_sdcoadd

Put in Release Notes: No

Module(s): sd

Description: Describe your changes here...

Added freq_tol parameter to asapmath.merge.
The freq_tol allows to specify frequency tolerance as
numeric value (1.0e6) or string ('1MHz'). The value will
be used to merge FREQUENCIES rows, FREQ_ID, and IFNO.


File size: 13.0 KB
Line 
1//
2// C++ Implementation: STFrequencies
3//
4// Description:
5//
6//
7// Author: Malte Marquarding <asap@atnf.csiro.au>, (C) 2006
8//
9// Copyright: See COPYING file that comes with this distribution
10//
11//
12#include <casa/iostream.h>
13#include <casa/iomanip.h>
14#include <casa/Exceptions/Error.h>
15#include <casa/Containers/RecordField.h>
16#include <casa/Arrays/IPosition.h>
17#include <casa/Logging/LogIO.h>
18#include <casa/BasicMath/Math.h>
19
20#include <tables/Tables/TableDesc.h>
21#include <tables/Tables/SetupNewTab.h>
22#include <tables/Tables/ScaColDesc.h>
23#include <tables/Tables/TableRecord.h>
24#include <tables/Tables/TableParse.h>
25#include <tables/Tables/TableRow.h>
26
27#include <coordinates/Coordinates/CoordinateSystem.h>
28#include <coordinates/Coordinates/CoordinateUtil.h>
29
30#include "STFrequencies.h"
31
32
33using namespace casa;
34
35namespace asap {
36
37const String STFrequencies::name_ = "FREQUENCIES";
38
39STFrequencies::STFrequencies(const Scantable& parent) :
40 STSubTable(parent, name_)
41{
42 setup();
43}
44
45STFrequencies::STFrequencies( casa::Table tab ) :
46 STSubTable(tab, name_)
47{
48 refpixCol_.attach(table_,"REFPIX");
49 refvalCol_.attach(table_,"REFVAL");
50 incrCol_.attach(table_,"INCREMENT");
51
52}
53
54STFrequencies::~STFrequencies()
55{
56}
57
58STFrequencies & STFrequencies::operator=( const STFrequencies & other )
59{
60 if ( this != &other ) {
61 static_cast<STSubTable&>(*this) = other;
62 refpixCol_.attach(table_,"REFPIX");
63 refvalCol_.attach(table_,"REFVAL");
64 incrCol_.attach(table_,"INCREMENT");
65 }
66 return *this;
67}
68
69void STFrequencies::setup( )
70{
71 // add to base class table
72 table_.addColumn(ScalarColumnDesc<Double>("REFPIX"));
73 table_.addColumn(ScalarColumnDesc<Double>("REFVAL"));
74 table_.addColumn(ScalarColumnDesc<Double>("INCREMENT"));
75
76 table_.rwKeywordSet().define("FRAME", String("TOPO"));
77 table_.rwKeywordSet().define("BASEFRAME", String("TOPO"));
78 table_.rwKeywordSet().define("EQUINOX",String( "J2000"));
79 table_.rwKeywordSet().define("UNIT", String(""));
80 table_.rwKeywordSet().define("DOPPLER", String("RADIO"));
81
82 // new cached columns
83 refpixCol_.attach(table_,"REFPIX");
84 refvalCol_.attach(table_,"REFVAL");
85 incrCol_.attach(table_,"INCREMENT");
86}
87
88uInt STFrequencies::addEntry( Double refpix, Double refval, Double inc )
89{
90 // test if this already exists
91 Table result = table_( near(table_.col("REFVAL"), refval)
92 && near(table_.col("REFPIX"), refpix)
93 && near(table_.col("INCREMENT"), inc), 1 );
94 uInt resultid = 0;
95 if ( result.nrow() > 0) {
96 ROScalarColumn<uInt> c(result, "ID");
97 c.get(0, resultid);
98
99 } else {
100 uInt rno = table_.nrow();
101 table_.addRow();
102 // get last assigned freq_id and increment
103 if ( rno > 0 ) {
104 idCol_.get(rno-1, resultid);
105 resultid++;
106 }
107 refpixCol_.put(rno, refpix);
108 refvalCol_.put(rno, refval);
109 incrCol_.put(rno, inc);
110 idCol_.put(rno, resultid);
111 }
112 return resultid;
113}
114
115
116
117void STFrequencies::getEntry( Double& refpix, Double& refval, Double& inc,
118 uInt id )
119{
120 Table t = table_(table_.col("ID") == Int(id), 1 );
121 if (t.nrow() == 0 ) {
122 throw(AipsError("STFrequencies::getEntry - freqID out of range"));
123 }
124 ROTableRow row(t);
125 // get first row - there should only be one matching id
126 const TableRecord& rec = row.get(0);
127 refpix = rec.asDouble("REFPIX");
128 refval = rec.asDouble("REFVAL");
129 inc = rec.asDouble("INCREMENT");
130}
131
132void STFrequencies::setEntry( Double refpix, Double refval, Double inc, uInt id )
133{
134 Table t = table_(table_.col("ID") == Int(id), 1 );
135 if (t.nrow() == 0 ) {
136 throw(AipsError("STFrequencies::getEntry - freqID out of range"));
137 }
138 for ( uInt i = 0 ; i < table_.nrow() ; i++ ) {
139 uInt fid ;
140 idCol_.get( i, fid ) ;
141 if ( fid == id ) {
142 refpixCol_.put( i, refpix ) ;
143 refvalCol_.put( i, refval ) ;
144 incrCol_.put( i, inc ) ;
145 }
146 }
147}
148
149SpectralCoordinate STFrequencies::getSpectralCoordinate( uInt id ) const
150{
151 Table t = table_(table_.col("ID") == Int(id), 1 );
152
153 if (t.nrow() == 0 ) {
154 throw(AipsError("STFrequencies::getSpectralCoordinate - ID out of range"));
155 }
156
157 // get the data
158 ROTableRow row(t);
159 // get first row - there should only be one matching id
160 const TableRecord& rec = row.get(0);
161 return SpectralCoordinate( getFrame(true), rec.asDouble("REFVAL"),
162 rec.asDouble("INCREMENT"),
163 rec.asDouble("REFPIX"));
164}
165
166/**
167SpectralCoordinate
168 STFrequencies::getSpectralCoordinate( const MDirection& md,
169 const MPosition& mp,
170 const MEpoch& me,
171 Double restfreq, uInt id ) const
172**/
173SpectralCoordinate
174 STFrequencies::getSpectralCoordinate( const MDirection& md,
175 const MPosition& mp,
176 const MEpoch& me,
177 Vector<Double> restfreq, uInt id ) const
178{
179 SpectralCoordinate spc = getSpectralCoordinate(id);
180 //spc.setRestFrequency(restfreq, True);
181 // for now just use the first rest frequency
182 if (restfreq.nelements()==0 ) {
183 restfreq.resize(1);
184 restfreq[0] = 0;
185 }
186 spc.setRestFrequency(restfreq[0], True);
187 if ( !spc.setReferenceConversion(getFrame(), me, mp, md) ) {
188 throw(AipsError("Couldn't convert frequency frame."));
189 }
190 String unitstr = getUnitString();
191 if ( !unitstr.empty() ) {
192 Unit unitu(unitstr);
193 if ( unitu == Unit("Hz") ) {
194 Vector<String> wau(1); wau = unitu.getName();
195 spc.setWorldAxisUnits(wau);
196 } else {
197 spc.setVelocity(unitstr, getDoppler());
198 }
199 }
200 return spc;
201}
202
203
204void STFrequencies::rescale( Float factor, const std::string& mode )
205{
206 TableRow row(table_);
207 TableRecord& outrec = row.record();
208 RecordFieldPtr<Double> rv(outrec, "REFVAL");
209 RecordFieldPtr<Double> rp(outrec, "REFPIX");
210 RecordFieldPtr<Double> inc(outrec, "INCREMENT");
211 for (uInt i=0; i<table_.nrow(); ++i) {
212
213 const TableRecord& rec = row.get(i);
214
215 SpectralCoordinate sc ( getFrame(true), rec.asDouble("REFVAL"),
216 rec.asDouble("INCREMENT"), rec.asDouble("REFPIX") );
217
218 SpectralCoordinate scout;
219 if (mode == "BIN") {
220 scout = binCsys(sc, Int(factor));
221 } else if (mode == "RESAMPLE") {
222 scout = resampleCsys(sc, factor);
223 }
224 *rv = scout.referenceValue()[0];
225 *rp = scout.referencePixel()[0];
226 *inc = scout.increment()[0];
227 row.put(i);
228 }
229}
230
231SpectralCoordinate STFrequencies::binCsys(const SpectralCoordinate& sc,
232 Int factor)
233{
234 CoordinateSystem csys;
235 csys.addCoordinate(sc);
236 IPosition factors(1, factor);
237 CoordinateSystem binnedcs =
238 CoordinateUtil::makeBinnedCoordinateSystem(factors, csys, False);
239 return binnedcs.spectralCoordinate(0);
240}
241
242SpectralCoordinate STFrequencies::resampleCsys(const SpectralCoordinate& sc,
243 Float width)
244{
245 Vector<Float> offset(1,0.0);
246 Vector<Float> factors(1,width);
247 Vector<Int> newshape;
248 CoordinateSystem csys;
249 csys.addCoordinate(sc);
250 CoordinateSystem csys2 = csys.subImage(offset, factors, newshape);
251 return csys2.spectralCoordinate(0);
252}
253
254
255MFrequency::Types STFrequencies::getFrame(bool base) const
256{
257 // get the ref frame
258 String rf;
259 if ( base )
260 rf = table_.keywordSet().asString("BASEFRAME");
261 else
262 rf = table_.keywordSet().asString("FRAME");
263
264 // Create SpectralCoordinate (units Hz)
265 MFrequency::Types mft;
266 if (!MFrequency::getType(mft, rf)) {
267 ostringstream oss;
268 LogIO os( casa::LogOrigin( "STFrequencies", "getFrame") );
269 os << LogIO::WARN << "WARNING: Frequency type unknown assuming TOPO"
270 << LogIO::POST;
271 mft = MFrequency::TOPO;
272 }
273
274 return mft;
275}
276
277std::string STFrequencies::getFrameString( bool base ) const
278{
279 if ( base ) return table_.keywordSet().asString("BASEFRAME");
280 else return table_.keywordSet().asString("FRAME");
281}
282
283std::string STFrequencies::getUnitString( ) const
284{
285 return table_.keywordSet().asString("UNIT");
286}
287
288Unit STFrequencies::getUnit( ) const
289{
290 return Unit(table_.keywordSet().asString("UNIT"));
291}
292
293std::string STFrequencies::getDopplerString( ) const
294{
295 return table_.keywordSet().asString("DOPPLER");
296}
297
298MDoppler::Types STFrequencies::getDoppler( ) const
299{
300 String dpl = table_.keywordSet().asString("DOPPLER");
301
302 // Create SpectralCoordinate (units Hz)
303 MDoppler::Types mdt;
304 if (!MDoppler::getType(mdt, dpl)) {
305 throw(AipsError("Doppler type unknown"));
306 }
307 return mdt;
308}
309
310std::string STFrequencies::print( int id, Bool strip ) const
311{
312 Table t;
313 ostringstream oss;
314 if ( id < 0 ) t = table_;
315 else t = table_(table_.col("ID") == Int(id), 1 );
316 ROTableRow row(t);
317 for (uInt i=0; i<t.nrow(); ++i) {
318 const TableRecord& rec = row.get(i);
319 oss << setw(8)
320 << t.keywordSet().asString("BASEFRAME") << setw(16) << setprecision(8)
321 << rec.asDouble("REFVAL") << setw(7)
322 << rec.asDouble("REFPIX")
323 << setw(15)
324 << rec.asDouble("INCREMENT");
325 }
326 String outstr(oss);
327 if (strip) {
328 int f = outstr.find_first_not_of(' ');
329 int l = outstr.find_last_not_of(' ', outstr.size());
330 if (f < 0) {
331 f = 0;
332 }
333 if ( l < f || l < f ) {
334 l = outstr.size();
335 }
336 return outstr.substr(f,l);
337 }
338 return outstr;
339}
340
341float STFrequencies::getRefFreq( uInt id, uInt channel )
342{
343 Table t = table_(table_.col("ID") == Int(id), 1 );
344 if ( t.nrow() == 0 ) throw(AipsError("Selected Illegal frequency id"));
345 ROTableRow row(t);
346 const TableRecord& rec = row.get(0);
347 return (Double(channel/2) - rec.asDouble("REFPIX"))
348 * rec.asDouble("INCREMENT") + rec.asDouble("REFVAL");
349}
350
351bool STFrequencies::conformant( const STFrequencies& other ) const
352{
353 const Record& r = table_.keywordSet();
354 const Record& ro = other.table_.keywordSet();
355 return ( r.asString("FRAME") == ro.asString("FRAME") &&
356 r.asString("EQUINOX") == ro.asString("EQUINOX") &&
357 r.asString("UNIT") == ro.asString("UNIT") &&
358 r.asString("DOPPLER") == ro.asString("DOPPLER")
359 );
360}
361
362std::vector< std::string > STFrequencies::getInfo( ) const
363{
364 const Record& r = table_.keywordSet();
365 std::vector<std::string> out;
366 out.push_back(r.asString("UNIT"));
367 out.push_back(r.asString("FRAME"));
368 out.push_back(r.asString("DOPPLER"));
369 return out;
370}
371
372void STFrequencies::setInfo( const std::vector< std::string >& theinfo )
373{
374 if ( theinfo.size() != 3 ) throw(AipsError("setInfo needs three parameters"));
375 try {
376 setUnit(theinfo[0]);
377 setFrame(theinfo[1]);
378 setDoppler(theinfo[2]);
379 } catch (AipsError& e) {
380 throw(e);
381 }
382}
383
384void STFrequencies::setUnit( const std::string & unit )
385{
386 if (unit == "" || unit == "pixel" || unit == "channel" ) {
387 table_.rwKeywordSet().define("UNIT", "");
388 } else {
389 Unit u(unit);
390 if ( u == Unit("km/s") || u == Unit("Hz") )
391 table_.rwKeywordSet().define("UNIT", unit);
392 else {
393 throw(AipsError("Illegal spectral unit."));
394 }
395 }
396}
397
398void STFrequencies::setFrame(MFrequency::Types frame, bool base )
399{
400 String f = MFrequency::showType(frame);
401 if (base)
402 table_.rwKeywordSet().define("BASEFRAME", f);
403 else
404 table_.rwKeywordSet().define("FRAME", f);
405
406}
407
408void STFrequencies::setFrame( const std::string & frame, bool base )
409{
410 MFrequency::Types mdr;
411 if (!MFrequency::getType(mdr, frame)) {
412 Int a,b;const uInt* c;
413 const String* valid = MFrequency::allMyTypes(a, b, c);
414 Vector<String> ftypes(IPosition(1,a), valid);
415 ostringstream oss;
416 oss << String("Please specify a legal frequency type. Types are\n");
417 oss << ftypes;
418 String msg(oss);
419 throw(AipsError(msg));
420 } else {
421 if (base)
422 table_.rwKeywordSet().define("BASEFRAME", frame);
423 else
424 table_.rwKeywordSet().define("FRAME", frame);
425 }
426}
427
428void STFrequencies::setDoppler( const std::string & doppler )
429{
430 MDoppler::Types mdt;
431 if (!MDoppler::getType(mdt, doppler)) {
432 Int a,b;const uInt* c;
433 const String* valid = MDoppler::allMyTypes(a, b, c);
434 Vector<String> ftypes(IPosition(1,a), valid);
435 ostringstream oss;
436 oss << String("Please specify a legal doppler type. Types are\n");
437 oss << ftypes;
438 String msg(oss);
439 throw(AipsError(msg));
440 } else {
441 table_.rwKeywordSet().define("DOPPLER", doppler);
442 }
443}
444
445void STFrequencies::shiftRefPix(int npix, uInt id)
446{
447 Table t = table_(table_.col("ID") == Int(id), 1 );
448 if ( t.nrow() == 0 ) throw(AipsError("Selected Illegal frequency id"));
449 ScalarColumn<Double> tcol(t, "REFPIX");
450 tcol.put(0, tcol(0)+Double(npix));
451}
452
453bool STFrequencies::match( Double refpix, Double refval,
454 Double inc, Double freqTolInHz,
455 uInt &id)
456{
457 ROScalarColumn<uInt> idCol(table_, "ID");
458 ROScalarColumn<Double> refPixCol(table_, "REFPIX");
459 ROScalarColumn<Double> refValCol(table_, "REFVAL");
460 ROScalarColumn<Double> incCol(table_, "INCREMENT");
461 for (uInt irow = 0; irow < table_.nrow(); ++irow) {
462 Double refInc = incCol(irow);
463 Double refFreq = refValCol(irow) - refPixCol(irow) * refInc;
464 Double freq0 = refval - refpix * inc;
465 if (nearAbs(inc, refInc, freqTolInHz) &&
466 nearAbs(refFreq, freq0, freqTolInHz)) {
467 id = irow;
468 return true;
469 }
470 }
471 return false;
472}
473
474} // namespace
Note: See TracBrowser for help on using the repository browser.