Line data Source code
1 : //# Copyright (C) 2009
2 : //# Associated Universities, Inc. Washington DC, USA.
3 : //#
4 : //# This library is free software; you can redistribute it and/or modify it
5 : //# under the terms of the GNU Library General Public License as published by
6 : //# the Free Software Foundation; either version 2 of the License, or (at your
7 : //# option) any later version.
8 : //#
9 : //# This library is distributed in the hope that it will be useful, but WITHOUT
10 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
11 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
12 : //# License for more details.
13 : //#
14 : //# You should have received a copy of the GNU Library General Public License
15 : //# along with this library; if not, write to the Free Software Foundation,
16 : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
17 : //#
18 : //# Correspondence concerning AIPS++ should be addressed as follows:
19 : //# Internet email: casa-feedback@nrao.edu.
20 : //# Postal address: AIPS++ Project Office
21 : //# National Radio Astronomy Observatory
22 : //# 520 Edgemont Road
23 : //# Charlottesville, VA 22903-2475 USA
24 : //#
25 :
26 : #ifndef IMAGEANALYSIS_IMAGEMETADATABASE_TCC
27 : #define IMAGEANALYSIS_IMAGEMETADATABASE_TCC
28 :
29 : #include <imageanalysis/ImageAnalysis/ImageMetaDataBase.h>
30 :
31 : #include <casacore/casa/aips.h>
32 :
33 : #include <casacore/casa/Arrays/ArrayLogical.h>
34 : #include <casacore/casa/Containers/ValueHolder.h>
35 : #include <casacore/casa/Quanta/QuantumHolder.h>
36 : #include <casacore/casa/Utilities/DataType.h>
37 : #include <casacore/images/Images/ImageSummary.h>
38 : #include <casacore/images/Images/ImageStatistics.h>
39 : #include <casacore/measures/Measures/MeasureHolder.h>
40 : #include <casacore/casa/Utilities/Regex.h>
41 :
42 : #include <iostream>
43 : #include <iomanip>
44 :
45 : #define _ORIGINB LogOrigin("ImageMetaDataBase", __func__, WHERE)
46 :
47 : using namespace std;
48 :
49 : using namespace casacore;
50 :
51 : namespace casa {
52 :
53 0 : template <class T> ImageMetaDataBase<T>::ImageMetaDataBase(SPCIIT image)
54 0 : : _image(image), _log(), _shape() {
55 0 : ThrowIf(! _image, "image cannot be NULL");
56 0 : _shape = _image->shape();
57 0 : }
58 :
59 0 : template <class T> Record ImageMetaDataBase<T>::_makeHeader() const {
60 0 : Record header;
61 0 : header.define(ImageMetaDataConstants::_IMTYPE, _getImType());
62 0 : header.define(ImageMetaDataConstants::_OBJECT, _getObject());
63 0 : const auto& csys = _getCoords();
64 :
65 0 : if (csys.hasDirectionCoordinate()) {
66 0 : const DirectionCoordinate& dc = csys.directionCoordinate();
67 0 : String equinox = MDirection::showType(
68 : dc.directionType()
69 : );
70 0 : header.define(ImageMetaDataConstants::_EQUINOX, _getEquinox());
71 0 : header.define(ImageMetaDataConstants::_PROJECTION, _getProjection());
72 0 : }
73 0 : header.define(ImageMetaDataConstants::_OBSDATE, _getEpochString());
74 0 : header.define(ImageMetaDataConstants::MASKS, _getMasks());
75 0 : header.define(ImageMetaDataConstants::_OBSERVER, _getObserver());
76 0 : header.define(ImageMetaDataConstants::_SHAPE, _getShape().asVector());
77 0 : header.define(ImageMetaDataConstants::_TELESCOPE, _getTelescope());
78 0 : header.define(ImageMetaDataConstants::_BUNIT, _getBrightnessUnit());
79 0 : if (csys.hasSpectralAxis()) {
80 0 : const SpectralCoordinate& sc = csys.spectralCoordinate();
81 0 : header.define(ImageMetaDataConstants::_RESTFREQ, sc.restFrequencies());
82 0 : header.define(
83 0 : ImageMetaDataConstants::_REFFREQTYPE , _getRefFreqType()
84 : );
85 : }
86 0 : const auto& info = _getInfo();
87 0 : if (info.hasSingleBeam()) {
88 0 : GaussianBeam beam = _getBeam();
89 0 : header.defineRecord(
90 : ImageMetaDataConstants::_BEAMMAJOR,
91 0 : QuantumHolder(beam.getMajor()).toRecord()
92 : );
93 0 : header.defineRecord(
94 : ImageMetaDataConstants::_BEAMMINOR,
95 0 : QuantumHolder(beam.getMinor()).toRecord()
96 : );
97 0 : header.defineRecord(
98 : ImageMetaDataConstants::_BEAMPA,
99 0 : QuantumHolder(beam.getPA(true)).toRecord()
100 : );
101 0 : }
102 0 : else if (info.hasMultipleBeams()) {
103 0 : String error;
104 0 : Record rec;
105 0 : info.toRecord(error, rec);
106 0 : static const String recName = "perplanebeams";
107 0 : Record beamRec = rec.asRecord(recName);
108 0 : beamRec.defineRecord(
109 0 : "median area beam", info.getBeamSet().getMedianAreaBeam().toRecord()
110 : );
111 0 : header.defineRecord(recName, beamRec);
112 0 : }
113 0 : auto cdelt = _getIncrements();
114 0 : auto units = _getAxisUnits();
115 0 : auto crpix = _getRefPixel();
116 0 : auto crval = _getRefValue();
117 0 : auto types = _getAxisNames();
118 0 : header.merge(_getStatistics());
119 0 : for (uInt i=0; i<cdelt.size(); ++i) {
120 0 : auto iString = String::toString(i + 1);
121 0 : auto delt = ImageMetaDataConstants::_CDELT + iString;
122 0 : header.define(delt, cdelt[i].getValue());
123 0 : auto unit = ImageMetaDataConstants::_CUNIT + iString;
124 0 : header.define(unit, units[i]);
125 0 : auto pix = ImageMetaDataConstants::_CRPIX + iString;
126 0 : header.define(pix, crpix[i]);
127 0 : auto val = ImageMetaDataConstants::_CRVAL + iString;
128 0 : header.define(val, crval[i].getValue());
129 0 : auto type = ImageMetaDataConstants::_CTYPE + iString;
130 0 : header.define(type, types[i]);
131 : }
132 0 : return header;
133 0 : }
134 :
135 : template <class T> const TableRecord ImageMetaDataBase<T>::_miscInfo() const {
136 : return _image->miscInfo();
137 : }
138 :
139 0 : template <class T> CoordinateSystem ImageMetaDataBase<T>::coordsys(
140 : const std::vector<Int>& pixelAxes
141 : ) const {
142 : // Recover CoordinateSytem into a Record
143 0 : auto cSys = _getCoords();
144 0 : if (pixelAxes.empty()) {
145 0 : return cSys;
146 : }
147 0 : Record rec;
148 0 : CoordinateSystem cSys2;
149 : // Fish out the coordinate of the desired axes
150 0 : uInt j = 0;
151 0 : const Int nPixelAxes = cSys.nPixelAxes();
152 0 : Vector<uInt> coordinates(cSys.nCoordinates(), uInt(0));
153 : Int coord, axisInCoord;
154 0 : for (const auto& axis: pixelAxes) {
155 0 : ThrowIf (
156 : axis < 0 || axis >= nPixelAxes,
157 : "Specified zero-based pixel axis " + String::toString(axis)
158 : + " is not a valid pixel axis"
159 : );
160 0 : cSys.findPixelAxis(coord, axisInCoord, uInt(axis));
161 0 : ThrowIf(
162 : coord < 0,
163 : "Zero-based pixel axis " + String::toString(axis)
164 : + " has been removed"
165 : );
166 0 : coordinates(coord)++;
167 : // Copy desired coordinate (once)
168 0 : if (coordinates(coord) == 1) {
169 0 : cSys2.addCoordinate(cSys.coordinate(coord));
170 : }
171 : }
172 : // Find mapping. Says where world axis i in cSys is in cSys2
173 0 : Vector<Int> worldAxisMap, worldAxisTranspose;
174 0 : Vector<Bool> refChange;
175 0 : ThrowIf(
176 : ! cSys2.worldMap(worldAxisMap, worldAxisTranspose, refChange, cSys),
177 : "Error finding world map because " + cSys2.errorMessage()
178 : );
179 : // Generate list of world axes to keep
180 0 : Vector<Int> keepList(cSys.nWorldAxes());
181 0 : Vector<Double> worldReplace;
182 0 : j = 0;
183 0 : for (const auto& axis: pixelAxes) {
184 0 : if (axis >= 0 && axis < nPixelAxes) {
185 0 : Int worldAxis = cSys.pixelAxisToWorldAxis(uInt(axis));
186 0 : ThrowIf(
187 : worldAxis < 0,
188 : "World axis corresponding to zero-based pixel axis "
189 : + String::toString(axis) + " has been removed"
190 : );
191 0 : keepList(j++) = worldAxisMap(worldAxis);
192 : }
193 : }
194 : // Remove unwanted world (and pixel) axes. Better would be to just
195 : // remove the pixel axes and leave the world axes there...
196 0 : if (j > 0) {
197 0 : keepList.resize(j, true);
198 0 : CoordinateUtil::removeAxes(cSys2, worldReplace, keepList, false);
199 : }
200 : // Copy the ObsInfo
201 0 : cSys2.setObsInfo(cSys.obsInfo());
202 0 : return cSys2;
203 0 : }
204 :
205 : template <class T> DataType ImageMetaDataBase<T>::dataType() const {
206 : return _image->dataType();
207 : }
208 :
209 : template <class T> Record* ImageMetaDataBase<T>::getBoundingBox(
210 : const Record& region
211 : ) const {
212 : const auto& csys = _getCoords();
213 : const auto shape = _getShape();
214 : const unique_ptr<ImageRegion> pRegion(
215 : ImageRegion::fromRecord(
216 : nullptr, csys, shape, region
217 : )
218 : );
219 : LatticeRegion latRegion = pRegion->toLatticeRegion(
220 : csys, shape
221 : );
222 : Slicer sl = latRegion.slicer();
223 : IPosition blc(sl.start());
224 : IPosition trc(sl.end());
225 : IPosition inc(sl.stride());
226 : IPosition length(sl.length());
227 : std::unique_ptr<Record> outRec(new Record());
228 : outRec->define("blc", blc.asVector());
229 : outRec->define("trc", trc.asVector());
230 : outRec->define("inc", inc.asVector());
231 : outRec->define("bbShape", (trc - blc + 1).asVector());
232 : outRec->define("regionShape", length.asVector());
233 : outRec->define("imageShape", shape.asVector());
234 : outRec->define("blcf", CoordinateUtil::formatCoordinate(blc, csys)); // 0-rel for use in C++
235 : outRec->define("trcf", CoordinateUtil::formatCoordinate(trc, csys));
236 : return outRec.release();
237 : }
238 :
239 : template <class T> ValueHolder ImageMetaDataBase<T>::getFITSValue(const String& key) const {
240 : String c = key;
241 : c.downcase();
242 : const TableRecord info = _miscInfo();
243 : if (c == ImageMetaDataConstants::_BUNIT) {
244 : return ValueHolder(_getBrightnessUnit());
245 : }
246 : else if (
247 : c.startsWith(ImageMetaDataConstants::_CDELT) || c.startsWith(ImageMetaDataConstants::_CRPIX)
248 : || c.startsWith(ImageMetaDataConstants::_CRVAL) || c.startsWith(ImageMetaDataConstants::_CTYPE)
249 : || c.startsWith(ImageMetaDataConstants::_CUNIT)
250 : ) {
251 : String prefix = c.substr(0, 5);
252 : uInt n = _getAxisNumber(c);
253 : if (prefix == ImageMetaDataConstants::_CDELT) {
254 : return ValueHolder(
255 : QuantumHolder(
256 : _getIncrements()[n-1]
257 : ).toRecord()
258 : );
259 : }
260 : else if (prefix == ImageMetaDataConstants::_CRPIX) {
261 : return ValueHolder(_getRefPixel()[n-1]);
262 : }
263 : else if (prefix == ImageMetaDataConstants::_CRVAL) {
264 : if (_getCoords().polarizationAxisNumber(false) == (Int)(n-1)) {
265 : return ValueHolder(
266 : _getStokes()
267 : );
268 : }
269 : else {
270 : return ValueHolder(
271 : QuantumHolder(_getRefValue()[n-1]).toRecord()
272 : );
273 : }
274 : }
275 : else if (prefix == ImageMetaDataConstants::_CTYPE) {
276 : return ValueHolder(_getAxisNames()[n-1]);
277 : }
278 : else if (prefix == ImageMetaDataConstants::_CUNIT) {
279 : return ValueHolder(_getAxisUnits()[n-1]);
280 : }
281 : }
282 : else if (c == ImageMetaDataConstants::_EQUINOX) {
283 : return ValueHolder(_getEquinox());
284 : }
285 : else if (c == ImageMetaDataConstants::_IMTYPE) {
286 : return ValueHolder(_getImType());
287 : }
288 : else if (c == ImageMetaDataConstants::MASKS) {
289 : return ValueHolder(_getMasks());
290 : }
291 : else if (c == ImageMetaDataConstants::_OBJECT) {
292 : return ValueHolder(_getObject());
293 : }
294 : else if (c == ImageMetaDataConstants::_OBSDATE || c == ImageMetaDataConstants::_EPOCH) {
295 : return ValueHolder(_getEpochString());
296 : }
297 : else if (c == ImageMetaDataConstants::_OBSERVER) {
298 : return ValueHolder(_getObserver());
299 : }
300 : else if (c == ImageMetaDataConstants::_PROJECTION) {
301 : return ValueHolder(_getProjection());
302 : }
303 : else if (c == ImageMetaDataConstants::_REFFREQTYPE) {
304 : return ValueHolder(_getRefFreqType());
305 : }
306 : else if (c == ImageMetaDataConstants::_RESTFREQ) {
307 : return ValueHolder(
308 : QuantumHolder(_getRestFrequency()).toRecord()
309 : );
310 : }
311 : else if (c == ImageMetaDataConstants::_SHAPE) {
312 : return ValueHolder(_getShape().asVector());
313 : }
314 : else if (c == ImageMetaDataConstants::_TELESCOPE) {
315 : return ValueHolder(_getTelescope());
316 : }
317 : else if (
318 : c == ImageMetaDataConstants::_BEAMMAJOR || c == ImageMetaDataConstants::_BEAMMINOR || c == ImageMetaDataConstants::_BEAMPA
319 : || c == ImageMetaDataConstants::_BMAJ || c == ImageMetaDataConstants::_BMIN || c == ImageMetaDataConstants::_BPA
320 : ) {
321 : GaussianBeam beam = _getBeam();
322 : if (c == ImageMetaDataConstants::_BEAMMAJOR || c == ImageMetaDataConstants::_BMAJ) {
323 : return ValueHolder(QuantumHolder(beam.getMajor()).toRecord());
324 : }
325 : else if (c == ImageMetaDataConstants::_BEAMMINOR || c == ImageMetaDataConstants::_BMIN) {
326 : return ValueHolder(QuantumHolder(beam.getMinor()).toRecord());
327 : }
328 : else {
329 : return ValueHolder(QuantumHolder(beam.getPA()).toRecord());
330 : }
331 : }
332 : else if (
333 : c == ImageMetaDataConstants::_DATAMIN || c == ImageMetaDataConstants::_DATAMAX || c == ImageMetaDataConstants::_MINPIXPOS
334 : || c == ImageMetaDataConstants::_MINPOS || c == ImageMetaDataConstants::_MAXPIXPOS || c == ImageMetaDataConstants::_MAXPOS
335 : ) {
336 : Record x = _getStatistics();
337 : if (c == ImageMetaDataConstants::_DATAMIN || c == ImageMetaDataConstants::_DATAMAX) {
338 : auto dt = dataType();
339 : if (dt == TpFloat) {
340 : Float val;
341 : x.get(c, val);
342 : return ValueHolder(val);
343 : }
344 : else if (dt == TpComplex) {
345 : Complex val;
346 : x.get(c, val);
347 : return ValueHolder(val);
348 : }
349 : else if (dt == TpDouble) {
350 : Double val;
351 : x.get(c, val);
352 : return ValueHolder(val);
353 : }
354 : else if (dt == TpDComplex) {
355 : DComplex val;
356 : x.get(c, val);
357 : return ValueHolder(val);
358 : }
359 : else {
360 : ThrowCc("Logic error");
361 : }
362 : }
363 : else if (c == ImageMetaDataConstants::_MINPOS || c == ImageMetaDataConstants::_MAXPOS) {
364 : return ValueHolder(x.asString(c));
365 : }
366 : else if (c == ImageMetaDataConstants::_MINPIXPOS || c == ImageMetaDataConstants::_MAXPIXPOS) {
367 : return ValueHolder(x.asArrayInt(c));
368 : }
369 : }
370 : else if (
371 : info.isDefined(key) || info.isDefined(c)
372 : ) {
373 : String x = info.isDefined(key) ? key : c;
374 : switch (info.type(info.fieldNumber(x))) {
375 : case TpString:
376 : return ValueHolder(info.asString(x));
377 : break;
378 : case TpInt:
379 : return ValueHolder(info.asInt(x));
380 : break;
381 : case TpDouble:
382 : return ValueHolder(info.asDouble(x));
383 : break;
384 : case TpRecord:
385 : // allow fall through
386 : case TpQuantity: {
387 : return ValueHolder(info.asRecord(x));
388 : break;
389 : }
390 : default:
391 : ostringstream os;
392 : os << info.type(info.fieldNumber(x));
393 : ThrowCc(
394 : "Unhandled data type "
395 : + os.str()
396 : + " for user defined type. Send us a bug report"
397 : );
398 : }
399 : }
400 : ThrowCc(
401 : "Unknown keyword " + key + ". If you are trying to use a key name from "
402 : "the imhead summary dictionary, note that some keys in "
403 : "mode='put'/'get' are different from mode='summary'. Please see imhead "
404 : "description in the CASA online documentation for complete details."
405 : );
406 : return ValueHolder();
407 : }
408 :
409 : template <class T> uInt ImageMetaDataBase<T>::_ndim() const {
410 : return _image->ndim();
411 : }
412 :
413 : template <class T> uInt ImageMetaDataBase<T>::_getAxisNumber(
414 : const String& key
415 : ) const {
416 : uInt n = 0;
417 : string sre = key.substr(0, 5) + "[0-9]+";
418 : Regex myre(Regex::makeCaseInsensitive(sre));
419 : if (key.find(myre) != String::npos) {
420 : n = String::toInt(key.substr(5));
421 : uInt ndim = _ndim();
422 : ThrowIf(
423 : n == 0,
424 : "The FITS convention is that axes "
425 : "are 1-based. Therefore, " + key + " is not a valid "
426 : "FITS keyword specification"
427 : );
428 : ThrowIf(
429 : n > ndim,
430 : "This image only has " + String::toString(ndim)
431 : + " axes."
432 : );
433 : }
434 : else {
435 : ThrowCc("Unsupported key " + key);
436 : }
437 : return n;
438 : }
439 :
440 0 : template <class T> String ImageMetaDataBase<T>::_getEpochString() const {
441 0 : return MVTime(_getObsDate().getValue()).string(MVTime::YMD, 12);
442 : }
443 :
444 0 : template <class T> IPosition ImageMetaDataBase<T>::_getShape() const {
445 0 : if (_shape.empty()) {
446 0 : _shape = _image->shape();
447 : }
448 0 : return _shape;
449 : }
450 :
451 0 : template <class T> void ImageMetaDataBase<T>::_fieldToLog(
452 : const Record& header,const String& field, Int precision
453 : ) const {
454 0 : _log << " -- " << field << ": ";
455 0 : if (header.isDefined(field)) {
456 0 : DataType type = header.type(header.idToNumber(field));
457 0 : if (precision >= 0) {
458 0 : _log.output() << setprecision(precision);
459 : }
460 0 : switch (type) {
461 0 : case TpArrayDouble: {
462 0 : _log << header.asArrayDouble(field);
463 0 : break;
464 : }
465 0 : case TpArrayInt: {
466 0 : _log << header.asArrayInt(field);
467 0 : break;
468 : }
469 0 : case TpArrayString: {
470 0 : _log << header.asArrayString(field);
471 0 : break;
472 : }
473 0 : case TpDouble: {
474 0 : _log << header.asDouble(field);
475 0 : break;
476 : }
477 0 : case TpRecord: {
478 0 : Record r = header.asRecord(field);
479 0 : QuantumHolder qh;
480 0 : String error;
481 0 : if (qh.fromRecord(error, r) && qh.isQuantity()) {
482 0 : Quantity q = qh.asQuantity();
483 0 : _log << q.getValue() << q.getUnit();
484 0 : }
485 : else {
486 0 : _log << "Logic Error: Don't know how to deal with records of this type "
487 0 : << LogIO::EXCEPTION;
488 : }
489 0 : break;
490 0 : }
491 0 : case TpString: {
492 0 : _log << header.asString(field);
493 0 : break;
494 : }
495 0 : default: {
496 0 : _log << "Logic Error: Unsupported type "
497 0 : << type << LogIO::EXCEPTION;
498 0 : break;
499 : }
500 : }
501 : }
502 : else {
503 0 : _log << "Not found";
504 : }
505 0 : _log << LogIO::POST;
506 0 : }
507 :
508 0 : template <class T> void ImageMetaDataBase<T>::_toLog(const Record& header) const {
509 0 : _log << _ORIGINB << "General --" << LogIO::POST;
510 0 : _fieldToLog(header, ImageMetaDataConstants::_IMTYPE);
511 0 : _fieldToLog(header, ImageMetaDataConstants::_OBJECT);
512 0 : _fieldToLog(header, ImageMetaDataConstants::_EQUINOX);
513 0 : _fieldToLog(header, ImageMetaDataConstants::_OBSDATE);
514 0 : _fieldToLog(header, ImageMetaDataConstants::_OBSERVER);
515 0 : _fieldToLog(header, ImageMetaDataConstants::_PROJECTION);
516 0 : if (header.isDefined(ImageMetaDataConstants::_RESTFREQ)) {
517 0 : _log << " -- " << ImageMetaDataConstants::_RESTFREQ << ": ";
518 0 : _log.output() << std::fixed << std::setprecision(1);
519 0 : _log << header.asArrayDouble(ImageMetaDataConstants::_RESTFREQ) << LogIO::POST;
520 : }
521 0 : _fieldToLog(header, ImageMetaDataConstants::_REFFREQTYPE);
522 0 : _fieldToLog(header, ImageMetaDataConstants::_TELESCOPE);
523 0 : _fieldToLog(header, ImageMetaDataConstants::_BEAMMAJOR, 12);
524 0 : _fieldToLog(header, ImageMetaDataConstants::_BEAMMINOR, 12);
525 0 : _fieldToLog(header, ImageMetaDataConstants::_BEAMPA, 12);
526 0 : _fieldToLog(header, ImageMetaDataConstants::_BUNIT);
527 0 : _fieldToLog(header, ImageMetaDataConstants::MASKS);
528 0 : _fieldToLog(header, ImageMetaDataConstants::_SHAPE);
529 0 : _fieldToLog(header, ImageMetaDataConstants::_DATAMIN);
530 0 : _fieldToLog(header, ImageMetaDataConstants::_DATAMAX);
531 0 : _fieldToLog(header, ImageMetaDataConstants::_MINPOS);
532 0 : _fieldToLog(header, ImageMetaDataConstants::_MINPIXPOS);
533 0 : _fieldToLog(header, ImageMetaDataConstants::_MAXPOS);
534 0 : _fieldToLog(header, ImageMetaDataConstants::_MAXPIXPOS);
535 :
536 0 : uInt i = 1;
537 0 : _log << LogIO::NORMAL << "axes --" << LogIO::POST;
538 0 : while (true) {
539 0 : String iString = String::toString(i);
540 0 : String key = ImageMetaDataConstants::_CTYPE + iString;
541 0 : if (! header.isDefined(key)) {
542 0 : break;
543 : }
544 0 : _log << " -- " << key << ": "
545 0 : << header.asString(key) << LogIO::POST;
546 0 : String unit = ImageMetaDataConstants::_CUNIT + iString;
547 0 : i++;
548 : }
549 0 : i = 1;
550 0 : _log << LogIO::NORMAL << ImageMetaDataConstants::_CRPIX << " --" << LogIO::POST;
551 0 : while (true) {
552 0 : String iString = String::toString(i);
553 0 : String key = ImageMetaDataConstants::_CRPIX + iString;
554 0 : if (! header.isDefined(key)) {
555 0 : break;
556 : }
557 0 : _log.output() << std::fixed << std::setprecision(1);
558 0 : _log << " -- " << key << ": " << header.asDouble(key)
559 0 : << LogIO::POST;
560 0 : i++;
561 : }
562 0 : i = 1;
563 0 : _log << LogIO::NORMAL << ImageMetaDataConstants::_CRVAL << " --" << LogIO::POST;
564 0 : while (true) {
565 0 : String iString = String::toString(i);
566 0 : String key = ImageMetaDataConstants::_CRVAL + iString;
567 0 : if (! header.isDefined(key)) {
568 0 : break;
569 : }
570 0 : _log << " -- " << key << ": ";
571 0 : ostringstream x;
572 0 : Double val = header.asDouble(key);
573 0 : x << val;
574 0 : String unit = ImageMetaDataConstants::_CUNIT + iString;
575 0 : if (header.isDefined(unit)) {
576 0 : x << header.asString(unit);
577 : }
578 0 : String valunit = x.str();
579 0 : if (header.isDefined(unit)) {
580 0 : String myunit = header.asString(unit);
581 0 : if (header.asString(unit).empty()) {
582 0 : String ctype = ImageMetaDataConstants::_CTYPE + iString;
583 0 : if (
584 0 : header.isDefined(ctype)
585 0 : && header.asString(ctype) == "Stokes"
586 : ) {
587 0 : valunit = "['" + Stokes::name((Stokes::StokesTypes)((Int)val)) + "']";
588 : }
589 0 : }
590 : else {
591 0 : String tmp = _doStandardFormat(val, myunit);
592 0 : if (! tmp.empty()) {
593 0 : valunit = tmp;
594 : }
595 0 : }
596 0 : }
597 0 : _log << valunit << LogIO::POST;
598 0 : i++;
599 : }
600 0 : i = 1;
601 0 : _log << LogIO::NORMAL << ImageMetaDataConstants::_CDELT << " --" << LogIO::POST;
602 0 : while (true) {
603 0 : String iString = String::toString(i);
604 0 : String key = ImageMetaDataConstants::_CDELT + iString;
605 0 : if (! header.isDefined(key)) {
606 0 : break;
607 : }
608 0 : _log << " -- " << key << ": ";
609 0 : Double val = header.asDouble(key);
610 0 : String unit = ImageMetaDataConstants::_CUNIT + iString;
611 0 : String myunit;
612 0 : if (header.isDefined(unit)) {
613 0 : myunit = header.asString(unit);
614 : }
615 0 : ostringstream x;
616 0 : x << val << myunit;
617 0 : String valunit = x.str();
618 0 : if (header.isDefined(unit)) {
619 0 : String myunit = header.asString(unit);
620 0 : if (! header.asString(unit).empty()) {
621 0 : String tmp = _doStandardFormat(val, myunit);
622 0 : if (! tmp.empty()) {
623 0 : valunit = tmp;
624 : }
625 0 : }
626 0 : }
627 0 : _log << valunit << LogIO::POST;
628 0 : i++;
629 : }
630 0 : i = 1;
631 0 : _log << LogIO::NORMAL << "units --" << LogIO::POST;
632 0 : while (true) {
633 0 : String iString = String::toString(i);
634 0 : String key = ImageMetaDataConstants::_CUNIT + iString;
635 0 : if (! header.isDefined(key)) {
636 0 : break;
637 : }
638 0 : _log << " -- " << key << ": "
639 0 : << header.asString(key) << LogIO::POST;
640 0 : String unit = ImageMetaDataConstants::_CUNIT + iString;
641 0 : i++;
642 : }
643 0 : }
644 :
645 0 : template <class T> String ImageMetaDataBase<T>::_doStandardFormat(
646 : Double value, const String& unit
647 : ) const {
648 0 : String valunit;
649 : try {
650 0 : Quantity q(1, unit);
651 0 : if (q.isConform(Quantity(1, "rad"))) {
652 : // to dms
653 0 : valunit = MVAngle(Quantity(value, unit)).string(MVAngle::CLEAN, 9) + "deg.min.sec";
654 : }
655 0 : else if (unit == "Hz") {
656 0 : ostringstream x;
657 0 : x << std::fixed << std::setprecision(1);
658 0 : x << value << "Hz";
659 0 : valunit = x.str();
660 0 : }
661 0 : }
662 0 : catch (const AipsError& x) {}
663 0 : return valunit;
664 0 : }
665 :
666 0 : template <class T> uInt ImageMetaDataBase<T>::nChannels() const {
667 0 : const CoordinateSystem csys = _getCoords();
668 0 : if (! csys.hasSpectralAxis()) {
669 0 : return 0;
670 : }
671 0 : return _getShape()[csys.spectralAxisNumber()];
672 0 : }
673 :
674 : template <class T> Bool ImageMetaDataBase<T>::isChannelNumberValid(
675 : const uInt chan
676 : ) const {
677 : if (! _getCoords().hasSpectralAxis()) {
678 : return false;
679 : }
680 : return (chan < nChannels());
681 : }
682 :
683 0 : template <class T> uInt ImageMetaDataBase<T>::nStokes() const {
684 0 : const CoordinateSystem& csys = _getCoords();
685 0 : if (! csys.hasPolarizationCoordinate()) {
686 0 : return 0;
687 : }
688 0 : return _getShape()[csys.polarizationAxisNumber()];
689 : }
690 :
691 : template <class T> Int ImageMetaDataBase<T>::stokesPixelNumber(
692 : const String& stokesString) const {
693 : Int pixNum = _getCoords().stokesPixelNumber(stokesString);
694 : if (pixNum >= (Int)nStokes()) {
695 : pixNum = -1;
696 : }
697 : return pixNum;
698 : }
699 :
700 0 : template <class T> String ImageMetaDataBase<T>::_getProjection() const {
701 0 : const CoordinateSystem csys = _getCoords();
702 0 : if (! csys.hasDirectionCoordinate()) {
703 0 : return "";
704 : }
705 0 : const DirectionCoordinate dc = csys.directionCoordinate();
706 0 : Projection proj = dc.projection();
707 0 : if (proj.type() == Projection::SIN) {
708 0 : Vector<Double> pars = proj.parameters();
709 0 : if (dc.isNCP()) {
710 0 : ostringstream os;
711 0 : os << "SIN (" << pars << "): NCP";
712 0 : return os.str();
713 0 : }
714 0 : else if(pars.size() == 2 && (anyNE(pars, 0.0))) {
715 : // modified SIN
716 0 : ostringstream os;
717 0 : os << "SIN (" << pars << ")";
718 0 : return os.str();
719 0 : }
720 0 : }
721 0 : return proj.name();
722 0 : }
723 :
724 : template <class T> String ImageMetaDataBase<T>::stokesAtPixel(
725 : const uInt pixel
726 : ) const {
727 : const CoordinateSystem& csys = _getCoords();
728 : if (! csys.hasPolarizationCoordinate() || pixel >= nStokes()) {
729 : return "";
730 : }
731 : return csys.stokesAtPixel(pixel);
732 : }
733 :
734 : template <class T> Bool ImageMetaDataBase<T>::isStokesValid(
735 : const String& stokesString
736 : ) const {
737 : if (! _getCoords().hasPolarizationCoordinate()) {
738 : return false;
739 : }
740 : Int stokesPixNum = stokesPixelNumber(stokesString);
741 : return stokesPixNum >= 0 && stokesPixNum < (Int)nStokes();
742 : }
743 :
744 0 : template <class T> Vector<Int> ImageMetaDataBase<T>::directionShape() const {
745 0 : Vector<Int> dirAxesNums = _getCoords().directionAxesNumbers();
746 0 : if (dirAxesNums.nelements() == 0) {
747 0 : return Vector<Int>();
748 : }
749 0 : Vector<Int> dirShape(2);
750 0 : IPosition shape = _getShape();
751 0 : dirShape[0] = shape[dirAxesNums[0]];
752 0 : dirShape[1] = shape[dirAxesNums[1]];
753 0 : return dirShape;
754 0 : }
755 :
756 : template <class T> Bool ImageMetaDataBase<T>::areChannelAndStokesValid(
757 : String& message, const uInt chan, const String& stokesString
758 : ) const {
759 : ostringstream os;
760 : Bool areValid = true;
761 : if (! isChannelNumberValid(chan)) {
762 : os << "Zero-based channel number " << chan << " is too large. There are only "
763 : << nChannels() << " spectral channels in this image.";
764 : areValid = false;
765 : }
766 : if (! isStokesValid(stokesString)) {
767 : if (! areValid) {
768 : os << " and ";
769 : }
770 : os << "Stokes parameter " << stokesString << " is not in image";
771 : areValid = false;
772 : }
773 : if (! areValid) {
774 : message = os.str();
775 : }
776 : return areValid;
777 : }
778 :
779 0 : template <class T> Record ImageMetaDataBase<T>::_calcStats() const {
780 0 : return _calcStatsT(_image);
781 : }
782 :
783 0 : template <class T> Record ImageMetaDataBase<T>::_calcStatsT(
784 : std::shared_ptr<const ImageInterface<T> > image
785 : ) const {
786 0 : Record x;
787 0 : if (! isReal(image->dataType())) {
788 : // the min and max and associated positions
789 : // cannot be calculated for complex images
790 0 : return x;
791 : }
792 0 : ImageStatistics<T> stats(*image);
793 0 : Array<typename NumericTraits<T>::PrecisionType> min;
794 0 : stats.getStatistic(min, LatticeStatsBase::MIN);
795 0 : if (min.size() == 0) {
796 : // image is completely masked
797 0 : return x;
798 : }
799 0 : x.define(ImageMetaDataConstants::_DATAMIN, min(IPosition(min.ndim(), 0)));
800 0 : Array<typename NumericTraits<T>::PrecisionType> max;
801 0 : stats.getStatistic(max, LatticeStatsBase::MAX);
802 0 : x.define(ImageMetaDataConstants::_DATAMAX, max(IPosition(max.ndim(), 0)));
803 0 : IPosition minPixPos, maxPixPos;
804 0 : stats.getMinMaxPos(minPixPos, maxPixPos);
805 0 : x.define(ImageMetaDataConstants::_MINPIXPOS, minPixPos.asVector());
806 0 : x.define(ImageMetaDataConstants::_MAXPIXPOS, maxPixPos.asVector());
807 0 : const auto& csys = _getCoords();
808 0 : Vector<Double> minPos = csys.toWorld(minPixPos);
809 0 : Vector<Double> maxPos = csys.toWorld(maxPixPos);
810 :
811 0 : String minFormat, maxFormat;
812 0 : uInt ndim = csys.nPixelAxes();
813 0 : Int spAxis = csys.spectralAxisNumber();
814 :
815 0 : for (uInt i=0; i<ndim; i++) {
816 0 : Int worldAxis = csys.pixelAxisToWorldAxis(i);
817 0 : String foundUnit;
818 0 : minFormat += csys.format(
819 : foundUnit, Coordinate::DEFAULT,
820 : minPos[i], worldAxis
821 : );
822 0 : maxFormat += csys.format(
823 : foundUnit, Coordinate::DEFAULT,
824 : maxPos[i], worldAxis
825 : );
826 0 : if ((Int)i == spAxis) {
827 0 : minFormat += foundUnit;
828 0 : maxFormat += foundUnit;
829 : }
830 0 : if (i != ndim-1) {
831 0 : minFormat += " ";
832 0 : maxFormat += " ";
833 : }
834 : }
835 0 : x.define(ImageMetaDataConstants::_MINPOS, minFormat);
836 0 : x.define(ImageMetaDataConstants::_MAXPOS, maxFormat);
837 0 : return x;
838 0 : }
839 :
840 : template <class T> Record ImageMetaDataBase<T>::toWorld(
841 : const Vector<Double>& pixel, const String& format, Bool doVelocity,
842 : const String& dirFrame, const String& freqFrame
843 : ) const {
844 : Vector<Double> pixel2 = pixel.copy();
845 : auto csys = _getCoords();
846 : {
847 : Vector<Double> replace = csys.referencePixel();
848 : const Int nIn = pixel2.size();
849 : const Int nOut = replace.size();
850 : Vector<Double> out(nOut);
851 : for (Int i = 0; i < nOut; ++i) {
852 : if (i > nIn - 1) {
853 : out(i) = replace(i);
854 : }
855 : else {
856 : out(i) = pixel2(i);
857 : }
858 : }
859 : pixel2.assign(out);
860 : }
861 :
862 : // Convert to world
863 :
864 : Vector<Double> world;
865 : Record rec;
866 : String dFrame = dirFrame;
867 : dFrame.upcase();
868 : String fFrame = freqFrame;
869 : fFrame.upcase();
870 : MDirection::Types dirType = csys.hasDirectionCoordinate()
871 : ? csys.directionCoordinate().directionType(dFrame == "CL")
872 : : MDirection::J2000;
873 : MFrequency::Types freqType = csys.hasSpectralAxis()
874 : ? csys.spectralCoordinate().frequencySystem(fFrame == "CL")
875 : : MFrequency::LSRK;
876 : if (
877 : (! csys.hasDirectionCoordinate() || dFrame == "CL")
878 : && (! csys.hasSpectralAxis() || fFrame == "CL")
879 : ) {
880 : ThrowIf(
881 : ! csys.toWorld(world, pixel2, true),
882 : "Error converting to world coordinates: " + csys.errorMessage()
883 : );
884 : }
885 : else if (
886 : (! csys.hasDirectionCoordinate() || dFrame == "NATIVE")
887 : && (! csys.hasSpectralAxis() || fFrame == "NATIVE")
888 : ) {
889 : ThrowIf(
890 : ! csys.toWorld(world, pixel2, false),
891 : "Error converting to world coordinates: " + csys.errorMessage()
892 : );
893 : }
894 : else {
895 : if (csys.hasDirectionCoordinate() && dFrame != "CL") {
896 : if (dFrame == "NATIVE") {
897 : dirType = csys.directionCoordinate().directionType(false);
898 : }
899 : else {
900 : ThrowIf(
901 : ! MDirection::getType(dirType, dFrame),
902 : "Unknown direction reference frame " + dirFrame
903 : );
904 : }
905 : auto dirCoord = csys.directionCoordinate();
906 : dirCoord.setReferenceConversion(dirType);
907 : csys.replaceCoordinate(dirCoord, csys.directionCoordinateNumber());
908 : }
909 : if (csys.hasSpectralAxis() && fFrame != "CL") {
910 : if (fFrame == "NATIVE") {
911 : freqType = csys.spectralCoordinate().frequencySystem(false);
912 : }
913 : else {
914 : ThrowIf(
915 : ! MFrequency::getType(freqType, fFrame),
916 : "Unknown frequency reference frame " + freqFrame
917 : );
918 : }
919 : auto specCoord = csys.spectralCoordinate();
920 : MFrequency::Types clFrame;
921 : MEpoch epoch;
922 : MPosition pos;
923 : MDirection dir;
924 : specCoord.getReferenceConversion(clFrame, epoch, pos, dir);
925 : specCoord.setReferenceConversion(freqType, epoch, pos, dir);
926 : csys.replaceCoordinate(specCoord, csys.spectralCoordinateNumber());
927 : }
928 : ThrowIf(
929 : ! csys.toWorld(world, pixel2, true),
930 : "Error converting to world coordinates: " + csys.errorMessage()
931 : );
932 :
933 : }
934 : return _worldVectorToRecord(
935 : csys, world, -1, format, true, true,
936 : doVelocity, dirType, freqType
937 : );
938 : }
939 :
940 : template <class T> Record ImageMetaDataBase<T>::_worldVectorToRecord(
941 : const CoordinateSystem& csys, const Vector<Double>& world, Int c,
942 : const String& format, Bool isAbsolute, Bool showAsAbsolute, Bool doVelocity,
943 : MDirection::Types dirFrame, MFrequency::Types freqFrame
944 : ) {
945 : // World vector must be in the native units of cSys
946 : // c = -1 means world must be length cSys.nWorldAxes
947 : // c > 0 means world must be length cSys.coordinate(c).nWorldAxes()
948 : // format from 'n,q,s,m'
949 :
950 : auto ct = upcase(format);
951 : Vector<String> units;
952 : if (c < 0) {
953 : units = csys.worldAxisUnits();
954 : }
955 : else {
956 : units = csys.coordinate(c).worldAxisUnits();
957 : }
958 : AlwaysAssert(world.size() == units.size(),AipsError);
959 : Record rec;
960 : if (ct.contains("N")) {
961 : rec.define("numeric", world);
962 : }
963 : if (ct.contains("Q")) {
964 : String error;
965 : Record recQ1, recQ2;
966 : for (uInt i = 0; i < world.size(); ++i) {
967 : Quantum<Double> worldQ(world(i), Unit(units(i)));
968 : QuantumHolder h(worldQ);
969 : ThrowIf(! h.toRecord(error, recQ1), error);
970 : recQ2.defineRecord(i, recQ1);
971 : }
972 : rec.defineRecord("quantity", recQ2);
973 : }
974 : if (ct.contains("S")) {
975 : Vector<Int> worldAxes;
976 : if (c < 0) {
977 : worldAxes.resize(world.size());
978 : indgen(worldAxes);
979 : }
980 : else {
981 : worldAxes = csys.worldAxes(c);
982 : }
983 : Coordinate::formatType fType = Coordinate::SCIENTIFIC;
984 : Int prec = 8;
985 : String u;
986 : Int coord, axisInCoord;
987 : Vector<String> fs(world.nelements());
988 : for (uInt i = 0; i < world.size(); ++i) {
989 : csys.findWorldAxis(coord, axisInCoord, i);
990 : if (
991 : csys.type(coord) == Coordinate::DIRECTION
992 : || csys.type(coord) == Coordinate::STOKES
993 : ) {
994 : fType = Coordinate::DEFAULT;
995 : }
996 : else {
997 : fType = Coordinate::SCIENTIFIC;
998 : }
999 : u = "";
1000 : fs(i) = csys.format(
1001 : u, fType, world(i), worldAxes(i),
1002 : isAbsolute, showAsAbsolute, prec
1003 : );
1004 : if (! u.empty() && (u != " ")) {
1005 : fs(i) += " " + u;
1006 : }
1007 : }
1008 : rec.define("string", fs);
1009 : }
1010 : if (ct.contains(String("M"))) {
1011 : Record recM = _worldVectorToMeasures(
1012 : csys, world, c, isAbsolute, doVelocity,
1013 : dirFrame, freqFrame
1014 : );
1015 : rec.defineRecord("measure", recM);
1016 : }
1017 : return rec;
1018 : }
1019 :
1020 : template <class T> Record ImageMetaDataBase<T>::_worldVectorToMeasures(
1021 : const CoordinateSystem& csys, const Vector<Double>& world,
1022 : Int c, Bool abs, Bool doVelocity, MDirection::Types dirFrame,
1023 : MFrequency::Types freqFrame
1024 : ) {
1025 : LogIO log;
1026 : log << LogOrigin("ImageMetaDataBase", __func__);
1027 : uInt directionCount, spectralCount, linearCount, stokesCount, tabularCount;
1028 : directionCount = spectralCount = linearCount = stokesCount = tabularCount
1029 : = 0;
1030 : // Loop over desired Coordinates
1031 : Record rec;
1032 : String error;
1033 : uInt s, e;
1034 : if (c < 0) {
1035 : AlwaysAssert(world.nelements()==csys.nWorldAxes(), AipsError);
1036 : s = 0;
1037 : e = csys.nCoordinates();
1038 : }
1039 : else {
1040 : AlwaysAssert(world.nelements()==csys.coordinate(c).nWorldAxes(), AipsError);
1041 : s = c;
1042 : e = c + 1;
1043 : }
1044 : for (uInt i = s; i < e; ++i) {
1045 : // Find the world axes in the CoordinateSystem that this coordinate belongs to
1046 :
1047 : const auto& worldAxes = csys.worldAxes(i);
1048 : const auto nWorldAxes = worldAxes.size();
1049 : Vector<Double> world2(nWorldAxes);
1050 : const auto& coord = csys.coordinate(i);
1051 : auto units = coord.worldAxisUnits();
1052 : Bool none = true;
1053 :
1054 : // Fill in missing world axes if all coordinates specified
1055 :
1056 : if (c < 0) {
1057 : for (uInt j = 0; j < nWorldAxes; ++j) {
1058 : if (worldAxes(j) < 0) {
1059 : world2[j] = coord.referenceValue()[j];
1060 : }
1061 : else {
1062 : world2(j) = world(worldAxes[j]);
1063 : none = false;
1064 : }
1065 : }
1066 : }
1067 : else {
1068 : world2 = world;
1069 : none = false;
1070 : }
1071 : if (
1072 : csys.type(i) == Coordinate::LINEAR
1073 : || csys.type(i) == Coordinate::TABULAR
1074 : ) {
1075 : if (!none) {
1076 : Record linRec1, linRec2;
1077 : for (uInt k = 0; k < world2.size(); ++k) {
1078 : Quantum<Double> value(world2(k), units(k));
1079 : QuantumHolder h(value);
1080 : ThrowIf(
1081 : ! h.toRecord(error, linRec1), error
1082 : );
1083 : linRec2.defineRecord(k, linRec1);
1084 : }
1085 : if (csys.type(i) == Coordinate::LINEAR) {
1086 : rec.defineRecord("linear", linRec2);
1087 : }
1088 : else if (csys.type(i) == Coordinate::TABULAR) {
1089 : rec.defineRecord("tabular", linRec2);
1090 : }
1091 : }
1092 : if (csys.type(i) == Coordinate::LINEAR) {
1093 : ++linearCount;
1094 : }
1095 : if (csys.type(i) == Coordinate::TABULAR) {
1096 : ++tabularCount;
1097 : }
1098 : }
1099 : else if (csys.type(i) == Coordinate::DIRECTION) {
1100 : ThrowIf(
1101 : ! abs,
1102 : "It is not possible to have a relative MDirection measure"
1103 : );
1104 : AlwaysAssert(worldAxes.nelements() == 2,AipsError);
1105 : if (!none) {
1106 : // Make an MDirection and stick in record
1107 :
1108 : Quantum<Double> t1(world2(0), units(0));
1109 : Quantum<Double> t2(world2(1), units(1));
1110 : MDirection direction(
1111 : t1, t2, dirFrame
1112 : );
1113 : MeasureHolder h(direction);
1114 : Record dirRec;
1115 : ThrowIf(
1116 : ! h.toRecord(error, dirRec), error
1117 : );
1118 : rec.defineRecord("direction", dirRec);
1119 : }
1120 : directionCount++;
1121 : }
1122 : else if (csys.type(i) == Coordinate::SPECTRAL) {
1123 : ThrowIf(
1124 : ! abs,
1125 : "It is not possible to have a relative MFrequency measure"
1126 : );
1127 : AlwaysAssert(worldAxes.nelements()==1,AipsError);
1128 : if (!none) {
1129 : // Make an MFrequency and stick in record
1130 :
1131 : Record specRec, specRec1;
1132 : Quantum<Double> t1(world2(0), units(0));
1133 : const auto& sc0 = csys.spectralCoordinate(i);
1134 : MFrequency frequency(t1, freqFrame);
1135 : MeasureHolder h(frequency);
1136 : ThrowIf(
1137 : ! h.toRecord(error, specRec1), error
1138 : );
1139 : specRec.defineRecord("frequency", specRec1);
1140 : if (doVelocity) {
1141 : SpectralCoordinate sc(sc0);
1142 :
1143 : // Do velocity conversions and stick in MDOppler
1144 : // Radio
1145 :
1146 : sc.setVelocity(String("km/s"), MDoppler::RADIO);
1147 : Quantum<Double> velocity;
1148 : ThrowIf(
1149 : ! sc.frequencyToVelocity(velocity, frequency),
1150 : sc.errorMessage()
1151 : );
1152 : MDoppler v(velocity, MDoppler::RADIO);
1153 : MeasureHolder h(v);
1154 : ThrowIf(
1155 : ! h.toRecord(error, specRec1), error
1156 : );
1157 : specRec.defineRecord("radiovelocity", specRec1);
1158 : // Optical
1159 :
1160 : sc.setVelocity(String("km/s"), MDoppler::OPTICAL);
1161 : ThrowIf(
1162 : ! sc.frequencyToVelocity(velocity, frequency),
1163 : sc.errorMessage()
1164 : );
1165 : v = MDoppler(velocity, MDoppler::OPTICAL);
1166 : h = MeasureHolder(v);
1167 : ThrowIf(
1168 : ! h.toRecord(error, specRec1), error
1169 : );
1170 : specRec.defineRecord("opticalvelocity", specRec1);
1171 :
1172 : // beta (relativistic/true)
1173 :
1174 : sc.setVelocity(String("km/s"), MDoppler::BETA);
1175 : ThrowIf(
1176 : ! sc.frequencyToVelocity(velocity, frequency),
1177 : sc.errorMessage()
1178 : );
1179 : v = MDoppler(velocity, MDoppler::BETA);
1180 : h = MeasureHolder(v);
1181 : ThrowIf(
1182 : ! h.toRecord(error, specRec1), error
1183 : );
1184 : specRec.defineRecord("betavelocity", specRec1);
1185 : }
1186 : rec.defineRecord("spectral", specRec);
1187 : }
1188 : ++spectralCount;
1189 : }
1190 : else if (csys.type(i) == Coordinate::STOKES) {
1191 : ThrowIf (
1192 : ! abs,
1193 : "It makes no sense to have a relative Stokes measure"
1194 : );
1195 : AlwaysAssert(worldAxes.size() == 1, AipsError);
1196 : if (!none) {
1197 : const auto& coord0 = csys.stokesCoordinate(i);
1198 : StokesCoordinate coord(coord0); // non-const
1199 : String u;
1200 : auto s = coord.format(
1201 : u, Coordinate::DEFAULT, world2(0),
1202 : 0, true, true, -1
1203 : );
1204 : rec.define("stokes", s);
1205 : }
1206 : ++stokesCount;
1207 : }
1208 : else {
1209 : ThrowCc("Cannot handle Coordinates of type " + csys.showType(i));
1210 : }
1211 : }
1212 : if (directionCount > 1) {
1213 : log << LogIO::WARN
1214 : << "There was more than one DirectionCoordinate in the "
1215 : << LogIO::POST;
1216 : log << LogIO::WARN << "CoordinateSystem. Only the last one is returned"
1217 : << LogIO::POST;
1218 : }
1219 : if (spectralCount > 1) {
1220 : log << LogIO::WARN
1221 : << "There was more than one SpectralCoordinate in the "
1222 : << LogIO::POST;
1223 : log << LogIO::WARN << "CoordinateSystem. Only the last one is returned"
1224 : << LogIO::POST;
1225 : }
1226 : if (stokesCount > 1) {
1227 : log << LogIO::WARN << "There was more than one StokesCoordinate in the "
1228 : << LogIO::POST;
1229 : log << LogIO::WARN << "CoordinateSystem. Only the last one is returned"
1230 : << LogIO::POST;
1231 : }
1232 : if (linearCount > 1) {
1233 : log << LogIO::WARN << "There was more than one LinearCoordinate in the "
1234 : << LogIO::POST;
1235 : log << LogIO::WARN << "CoordinateSystem. Only the last one is returned"
1236 : << LogIO::POST;
1237 : }
1238 : if (tabularCount > 1) {
1239 : log << LogIO::WARN
1240 : << "There was more than one TabularCoordinate in the "
1241 : << LogIO::POST;
1242 : log << LogIO::WARN << "CoordinateSystem. Only the last one is returned"
1243 : << LogIO::POST;
1244 : }
1245 : return rec;
1246 : }
1247 :
1248 : }
1249 :
1250 : #endif
|
