source: trunk/src/PlotHelper.cpp@ 2695

Last change on this file since 2695 was 2690, checked in by Kana Sugimoto, 12 years ago

New Development: No

JIRA Issue: No

Ready for Test: Yes

Interface Changes: No

What Interface Changed:

Test Programs: sdplot unit tests

Put in Release Notes: Yes

Module(s): asapplotter and sdplot

Description:

Proper WCS projection in plottype = "grid" in sdplot.
Removed debug flag.


File size: 6.3 KB
Line 
1//
2// C++ Interface: PlotHelper
3//
4// Description:
5// A small helper class to handle direction coordinate in asapplotter
6//
7// Author: Kana Sugimoto <kana.sugi@nao.ac.jp>, (C) 2012
8//
9// Copyright: See COPYING file that comes with this distribution
10//
11//
12
13// casacore
14#include <casa/Arrays/Vector.h>
15#include <casa/Quanta/Quantum.h>
16#include <casa/Quanta/QuantumHolder.h>
17#include <measures/Measures/MDirection.h>
18#include <casa/Logging/LogIO.h>
19
20#include "PlotHelper.h"
21
22//#define KS_DEBUG
23
24using namespace std ;
25using namespace casa ;
26using namespace asap ;
27
28namespace asap {
29
30PlotHelper::PlotHelper() : dircoord_(0)
31{
32#ifdef KS_DEBUG
33 cout << "Default constructor nrow = " << data_.nrow() << endl;
34#endif
35};
36
37PlotHelper::PlotHelper( const ScantableWrapper &s) : dircoord_(0)
38{
39#ifdef KS_DEBUG
40 cout << "Constructing PlotHelper with scantable wrapper" << endl;
41#endif
42 setScantable(s);
43};
44
45PlotHelper::~PlotHelper(){
46#ifdef KS_DEBUG
47 cout << "Called PlotHelper destructor" << endl;
48#endif
49 if (dircoord_){
50#ifdef KS_DEBUG
51 cout << "Destructing dircoord_" << endl;
52#endif
53 delete dircoord_;
54 dircoord_ = 0;
55 }
56};
57
58void PlotHelper::setScantable( const ScantableWrapper &s )
59{
60#ifdef KS_DEBUG
61 cout << "Setting scantable" << endl;
62#endif
63 data_ = s ;
64};
65
66// void PlotHelper::setGridParam(const int nx, const int ny, const string cellx, const string celly, string center, const string projection){
67// // Value check of nx and ny
68// if (nx < 1)
69// throw(AipsError("nx should be > 0"));
70// if (ny < 1)
71// throw(AipsError("ny should be > 0"));
72// // Destroy old coord
73// if (dircoord_){
74// cout << "Destructing dircoord_" << endl;
75// delete dircoord_;
76// dircoord_ = 0;
77// }
78
79// //
80// Double incX, incY;
81// Double centX, centY;
82// MDirection::Types mdt;
83// // projection
84// Projection::Type projType(Projection::type(String(projname)));
85// ///
86// incX = cellx;
87// incY = celly;
88// centX = centx;
89// centY = centy;
90// MDirection::getType(mdt, String(epoch));
91// ///
92// // Get map extent
93// Double xmax, xmin, ymax, ymin;
94// ROArrayColumn<Double> dirCol;
95// MDirection::Types stdt;
96// const bool stset = (data_ ? true : false);
97
98// if (stset) {
99// dirCol_.attach( data_, "DIRECTION" );
100// Matrix<Double> direction = dirCol.getColumn();
101// minMax(xmin, xmax, direction.row( 0 ));
102// minMax(ymin, ymax, direction.row( 1 ));
103// if (!MDirection::getType(stdt, data_.getDirectionRefString()))
104// throw AipsError("Failed to get direction reference from scantable.");
105// }
106// // center
107// if (center.size() == 0){
108// if (!stset)
109// throw AipsError("Scantable is not set. Could not resolve map center.");
110// cenX = 0.5 * (xmin + xmax);
111// cenY = 0.5 * (ymin + ymax);
112// mdt = stdt;
113// } else {
114// if (!MDirection::getType(mdt,center))
115// throw AipsError("Invalid direction reference in center");
116// if (stset && mdt != stdt)
117// throw AipsError("Direction reference of center should be the same as input scantable");
118// MDirection centdir;
119// centdir.setRefString(String(center));
120// ///
121// Vector<Double> centrad;
122// centrad = centdir.getAngle("rad").getValue();
123// cenX = centrad[0];
124// cenY = centrad[1];
125// // rotaion
126// if (stset) {
127// Double rotnum = round(abs(cenX - )/(C::))
128// }
129// }
130// // cell
131// Quantum<Double> qcellx, qcelly;
132// if (cellx.size() == 0 && celly.size() == 0){
133// // Need resolution
134// if (!stset)
135// throw AipsError("Scantable is not set. Could not resolve cell size.");
136// } else {
137// if (cellx.size() != 0){
138// readQuantity(qcellx, String(cellx));
139// qcellx.convert("rad");
140// } else if (celly.size() == 0) {
141// celly = cellx;
142//
143// }
144
145void PlotHelper::setGridParamVal(const int nx, const int ny, const double cellx, const double celly, const double centx, const double centy, const string epoch, const string projname){
146 // Value check of nx and ny
147 if (nx < 1)
148 throw(AipsError("nx should be > 0"));
149 if (ny < 1)
150 throw(AipsError("ny should be > 0"));
151 // Destroy old coord
152 if (dircoord_){
153#ifdef KS_DEBUG
154 cout << "Destructing dircoord_" << endl;
155#endif
156 delete dircoord_;
157 dircoord_ = 0;
158 }
159
160 // center (in rad)
161 Double centX(centx), centY(centy);
162 // cell size (in rad)
163 Double incX(cellx), incY(celly);
164 // epoch
165 MDirection::Types mdt;
166 MDirection::getType(mdt, String(epoch));
167 // projection
168 Projection::Type projType(Projection::type(String(projname)));
169
170 Matrix<Double> xform(2,2) ;
171 xform = 0.0 ;
172 xform.diagonal() = 1.0 ;
173 dircoord_ = new DirectionCoordinate(mdt, projType,
174 centX, centY, incX, incY,
175 xform,
176 0.5*Double(nx),
177 0.5*Double(ny)) ; // pixel at center
178// 0.5*Double(nx-1),
179// 0.5*Double(ny-1)) ; // pixel at grid
180#ifdef KS_DEBUG
181 {//Debug outputs
182 cout << "Direction coordinate is set: " << endl;
183 Vector<String> units = dircoord_->worldAxisUnits();
184 Vector<Double> refv = dircoord_->referenceValue();
185 cout <<"Reference: " << MDirection::showType(dircoord_->directionType()) << " " << refv[0] << units[0] << " " << refv[1] << units[1] << endl;
186 Vector<Double> refpix = dircoord_->referencePixel();
187 cout <<"Reference Pixel: [" << refpix[0] << ", " << refpix[1] << "]" << endl;
188 Vector<Double> inc = dircoord_->increment();
189 cout <<"Increments: [" << inc[0] << ", " << inc[1] << "]" << endl;
190 cout <<"Projection: " << dircoord_->projection().name() << endl;
191 }
192#endif
193};
194
195vector<double> PlotHelper::getGridPixel(const int whichrow){
196 if (data_.nrow() < 1)
197 throw AipsError("Scantable is not set. Could not get direction.");
198 else if (whichrow > int(data_.nrow()) - 1)
199 throw AipsError("Row index out of range.");
200 if (!dircoord_)
201 throw AipsError("Direction coordinate is not defined.");
202
203 Vector<Double> pixel;
204 MDirection world;
205 vector<double> outvec;
206 world = data_.getCP()->getDirection(whichrow);
207#ifdef KS_DEBUG
208 cout << "searching pixel position (world = " << data_.getCP()->getDirectionString(whichrow) << " = [" << world.getAngle("rad").getValue()[0] << ", " << world.getAngle("rad").getValue()[1] << "])" << endl;
209#endif
210 dircoord_->toPixel(pixel, world);
211#ifdef KS_DEBUG
212 cout << "got pixel = [" << pixel[0] << ", " << pixel[1] << "]" << endl;
213#endif
214 // convert pixel to std vector
215 pixel.tovector(outvec);
216 return outvec;
217};
218
219} //namespace asap
Note: See TracBrowser for help on using the repository browser.