source: trunk/src/PlotHelper.cpp@ 2718

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

New Development: No (a bug fix)

JIRA Issue: No

Ready for Test: Yes

Interface Changes: No

What Interface Changed:

Test Programs:

Put in Release Notes: No

Module(s): asapplotter and sdplot

Description:

fixed a bug in increment.
better log messages.


File size: 13.2 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 <casa/Logging/LogIO.h>
18
19#include <measures/Measures/MDirection.h>
20#include <tables/Tables/TableRecord.h>
21
22
23#include "PlotHelper.h"
24
25//#ifndef KS_DEBUG
26//#define KS_DEBUG
27//#endif
28
29using namespace std ;
30using namespace casa ;
31using namespace asap ;
32
33namespace asap {
34
35PlotHelper::PlotHelper() : dircoord_(0)
36{
37#ifdef KS_DEBUG
38 cout << "Default constructor nrow = " << data_p->nrow() << endl;
39#endif
40};
41
42PlotHelper::PlotHelper( const ScantableWrapper &s) : dircoord_(0)
43{
44#ifdef KS_DEBUG
45 cout << "Constructing PlotHelper with scantable wrapper" << endl;
46#endif
47 setScantable(s);
48};
49
50PlotHelper::~PlotHelper(){
51#ifdef KS_DEBUG
52 cout << "Called PlotHelper destructor" << endl;
53#endif
54 if (dircoord_){
55#ifdef KS_DEBUG
56 cout << "Destructing dircoord_" << endl;
57#endif
58 delete dircoord_;
59 dircoord_ = 0;
60 }
61};
62
63void PlotHelper::setScantable( const ScantableWrapper &s )
64{
65#ifdef KS_DEBUG
66 cout << "Setting scantable" << endl;
67#endif
68 data_p = s.getCP();
69};
70
71
72DirectionCoordinate PlotHelper::getSTCoord(const int nx, const int ny,
73 const Projection::Type ptype)
74{
75 LogIO os(LogOrigin("PlotHelper","getSTCoord()", WHERE));
76 os << "Getting pointing information of the scantable." << LogIO::POST;
77 if (data_p->nrow() < 1)
78 throw AipsError("Scantable is not set. Please set a scantable first.");
79
80 // First, generate rough direction coordinate.
81 DirectionCoordinate coord;
82 Double incx, incy;
83 MDirection::Types mdt;
84 ROArrayColumn<Double> dircol;
85 Double xmax, xmin, ymax, ymin;
86 Double centx, centy;
87 Matrix<Double> xform(2,2);
88 xform = 0.0;
89 xform.diagonal() = 1.0;
90 // Rough estimates of center and cell from scantable DIRECTIONs.
91 dircol.attach( data_p->table(), "DIRECTION" );
92 const Vector<String> udir = dircol.keywordSet().asArrayString("QuantumUnits");
93 const Matrix<Double> direction = dircol.getColumn();
94 minMax(xmin, xmax, direction.row(0));
95 minMax(ymin, ymax, direction.row(1));
96 if (!MDirection::getType(mdt, data_p->getDirectionRefString()))
97 throw AipsError("Failed to get direction reference from scantable.");
98 centx = 0.5 * (xmin + xmax);
99 centy = 0.5 * (ymin + ymax);
100 incx = abs(xmax - xmin) / (double) nx * cos(centy);
101 incy = abs(ymax - ymin) / (double) ny;
102 // Generate a temporal direction coordinte
103 coord = DirectionCoordinate(mdt, ptype,
104 centx, centy, incx, incy, xform,
105 0.5*Double(nx), 0.5*Double(ny));
106 coord.setWorldAxisUnits(udir);
107#ifdef KS_DEBUG
108 {//Debug outputs
109 cout << "Generated a temporal direction coordinate of a scantable: " << endl;
110 Vector<String> units = coord.worldAxisUnits();
111 Vector<Double> refv = coord.referenceValue();
112 cout <<"- Reference: " << MDirection::showType(coord.directionType()) << " " << refv[0] << units[0] << " " << refv[1] << units[1] << endl;
113 Vector<Double> refpix = coord.referencePixel();
114 cout <<"- Reference Pixel: [" << refpix[0] << ", " << refpix[1] << "]" << endl;
115 Vector<Double> inc = coord.increment();
116 cout <<"- Increments: [" << inc[0] << ", " << inc[1] << "]" << endl;
117 cout <<"- Projection: " << coord.projection().name() << endl;
118 }
119#endif
120 return coord;
121}
122
123void PlotHelper::setGridParam(const int nx, const int ny,
124 const string cellx, const string celly,
125 string center, const string projname)
126{
127 LogIO os(LogOrigin("PlotHelper","setGridParam()", WHERE));
128 // Value check of nx and ny
129 if (nx < 1)
130 throw(AipsError("nx should be > 0"));
131 if (ny < 1)
132 throw(AipsError("ny should be > 0"));
133 // Destroy old coord
134 if (dircoord_){
135#ifdef KS_DEBUG
136 cout << "Destructing dircoord_" << endl;
137#endif
138 delete dircoord_;
139 dircoord_ = 0;
140 }
141
142 // Check for availability of data_p
143 const bool stset = ((data_p->nrow() > 0) ? true : false);
144 const bool needst = center.empty() || (cellx.empty() && celly.empty());
145 if (needst && !stset)
146 throw AipsError("Could not resolve grid parameters. Please set a scantable first.");
147
148 // projection
149 Projection::Type projtype(Projection::type(String(projname)));
150
151 // Calculate projected map center (in world coordinate)
152 // and extent (in pixel coordinate) from scantable DIRECIONs (if necessary).
153 MDirection stcent;
154 Double stxmax, stxmin, stymax, stymin;
155 MDirection::Types stdt;
156 DirectionCoordinate stcoord;
157 Quantum<Double> stincx, stincy;
158 if (needst && stset) {
159 stcoord = getSTCoord(nx, ny, projtype);
160 Vector<Double> inc = stcoord.increment();
161 Vector<String> units = stcoord.worldAxisUnits();
162 stincx = Quantum<Double>(inc[0], units[0]);
163 stincy = Quantum<Double>(inc[1], units[0]);
164 stdt = stcoord.directionType();
165 // Get the extent of directions in Pixel coordinate
166 ROArrayColumn<Double> dircol;
167 dircol.attach( data_p->table(), "DIRECTION" );
168 const Matrix<Double> direction = dircol.getColumn();
169 Matrix<Double> dirpix;
170 Vector<Bool> failures;
171 if (!stcoord.toPixelMany(dirpix, direction, failures))
172 throw AipsError("Failed to get directions in pixel coordinate.");
173 minMax(stxmin, stxmax, dirpix.row(0));
174 minMax(stymin, stymax, dirpix.row(1));
175 // Get the direction center in World coordinate.
176 Vector<Double> centpix(2);
177 centpix[0] = 0.5 * (stxmin + stxmax);
178 centpix[1] = 0.5 * (stymin + stymax);
179 stcoord.toWorld(stcent, centpix);
180#ifdef KS_DEBUG
181 {//Debug output
182 Quantum< Vector<Double> > qcent;
183 cout << "Got the center and map extent of scantable." << endl;
184 qcent = stcent.getAngle();
185 cout << "- Center of DIRECTIONs: " << stcent.getRefString() << " "
186 << qcent.getValue()[0] << qcent.getUnit() << " "
187 << qcent.getValue()[1] << qcent.getUnit() << endl;
188 cout << "- Map extent: [" << (stxmax-stxmin)*stincx.getValue() << ", "
189 << (stymax-stymin)*stincy.getValue() << "] ( " << stincx.getUnit() << ")" << endl;
190 }
191#endif
192 }
193
194 os << "Setting grid parameters" << LogIO::POST;
195 // Now, define direction coordinate from input parameters (in radian)
196 Double incx, incy;
197 Double centx, centy;
198 MDirection::Types mdt;
199 // center
200 if (center.empty()){
201 os << "center is not specified. Using pointing center of the scantable." << LogIO::POST;
202 if (!stset)
203 throw AipsError("Scantable is not set. Could not resolve map center.");
204
205 centx = stcent.getAngle("rad").getValue()[0];
206 centy = stcent.getAngle("rad").getValue()[1];
207 mdt = stdt;
208 } else {
209 os << "Using user defined center" << LogIO::POST;
210 // Parse center string
211 string::size_type pos0 = center.find(" ");
212
213 if (pos0 == string::npos)
214 throw AipsError("bad string format in center direction");
215
216 string::size_type pos1 = center.find(" ", pos0+1);
217 String sepoch, sra, sdec;
218 if (pos1 != string::npos) {
219 sepoch = center.substr(0, pos0);
220 sra = center.substr(pos0+1, pos1-pos0);
221 sdec = center.substr(pos1+1);
222 } else {
223 sepoch = "J2000";
224 sra = center.substr(0, pos0);
225 sdec = center.substr(pos0+1);
226 }
227 if (!MDirection::getType(mdt,sepoch))
228 throw AipsError("Invalid direction reference in center");
229 if (stset && (mdt != stdt))
230 throw AipsError("Direction reference of center should be the same as input scantable");
231 QuantumHolder qh ;
232 String err ;
233 qh.fromString(err, sra);
234 Quantum<Double> ra = qh.asQuantumDouble();
235 qh.fromString(err, sdec) ;
236 Quantum<Double> dec = qh.asQuantumDouble();
237 centx = ra.getValue("rad");
238 centy = dec.getValue("rad");
239 // rotaion
240 if (stset) {
241 Double stcentx = stcent.getAngle("rad").getValue()[0];
242 Double rotnum = round( abs(centx - stcentx) / (C::_2pi) );
243 if (centx < stcentx) rotnum *= -1;
244 centx -= (rotnum * C::_2pi);
245 }
246 }
247 os << "Grid center: ( " << centx << " rad , " << centy << " rad ) " << LogIO::POST;
248
249 // cell
250 if (cellx.empty() && celly.empty()){
251 os << "cell size is not specified. Using cell size to cover all pointings in the scantable." << LogIO::POST;
252 if (!stset)
253 throw AipsError("Scantable is not set. Could not resolve cell size.");
254
255 Vector<Double> centpix;
256 MDirection centmd = MDirection(Quantum<Double>(centx, "rad"),
257 Quantum<Double>(centy, "rad"), mdt);
258 stcoord.toPixel(centpix, centmd);
259#ifdef KS_DEBUG
260 cout << "- got centpix [" << centpix[0] << ", " << centpix[1] << "]" <<endl;
261#endif
262 Double wx = max( abs(stxmax-centpix[0]), abs(stxmin-centpix[0]) )
263 * 2 * stincx.getValue("rad");
264 Double wy = max( abs(stymax-centpix[1]), abs(stymin-centpix[1]) )
265 * 2 * stincy.getValue("rad");
266 incx = wx / max(nx - 1., 1.);
267 incy = wy / max(ny - 1., 1.);
268 } else {
269 os << "Using user defined cell size" << LogIO::POST;
270 Quantum<Double> qcellx, qcelly;
271 if (!cellx.empty() && !celly.empty()){
272 readQuantity(qcellx, String(cellx));
273 readQuantity(qcelly, String(celly));
274 } else if (celly.empty()) {
275 readQuantity(qcellx, String(cellx));
276 qcelly = qcellx;
277 } else { //only celly is defined
278 readQuantity(qcelly, String(celly));
279 qcellx = qcelly;
280 }
281 incx = qcellx.getValue("rad");
282 incy = qcelly.getValue("rad");
283 }
284 // inc should be negative to transpose plot
285 incx = -abs(incx);
286 incy = -abs(incy);
287
288 os << "Spacing: ( " << abs(incx) << " rad , " << abs(incy) << " rad )" <<endl;
289
290 Matrix<Double> xform(2,2) ;
291 xform = 0.0 ;
292 xform.diagonal() = 1.0 ;
293 dircoord_ = new DirectionCoordinate(mdt, projtype,
294 centx, centy, incx, incy,
295 xform,
296 0.5*Double(nx),
297 0.5*Double(ny)) ; // pixel at center
298 {//Summary
299 os << "Successfully generated grid coordinate:" << LogIO::POST;
300 Vector<String> units = dircoord_->worldAxisUnits();
301 Vector<Double> refv = dircoord_->referenceValue();
302 os <<"- Reference Direction : " << MDirection::showType(dircoord_->directionType())
303 << " " << refv[0] << units[0] << " " << refv[1] << units[1] << LogIO::POST;
304 Vector<Double> refpix = dircoord_->referencePixel();
305 os <<"- Reference Pixel : [" << refpix[0] << ", " << refpix[1] << "]" << LogIO::POST;
306 Vector<Double> inc = dircoord_->increment();
307 os <<"- Increments : [" << inc[0] << ", " << inc[1] << "]" << LogIO::POST;
308 os <<"- Projection Type : " << dircoord_->projection().name() << LogIO::POST;
309 }
310};
311
312void PlotHelper::setGridParamVal(const int nx, const int ny,
313 const double cellx, const double celly,
314 const double centx, const double centy,
315 const string epoch, const string projname){
316 LogIO os(LogOrigin("PlotHelper","setGridParamVal()", WHERE));
317 // Value check of nx and ny
318 if (nx < 1)
319 throw(AipsError("nx should be > 0"));
320 if (ny < 1)
321 throw(AipsError("ny should be > 0"));
322 // Destroy old coord
323 if (dircoord_){
324#ifdef KS_DEBUG
325 cout << "Destructing dircoord_" << endl;
326#endif
327 delete dircoord_;
328 dircoord_ = 0;
329 }
330
331 // center (in rad)
332 Double centX(centx), centY(centy);
333 // cell size (in rad)
334 Double incX(cellx), incY(celly);
335 // epoch
336 MDirection::Types mdt;
337 MDirection::getType(mdt, String(epoch));
338 // projection
339 Projection::Type projType(Projection::type(String(projname)));
340
341 Matrix<Double> xform(2,2) ;
342 xform = 0.0 ;
343 xform.diagonal() = 1.0 ;
344 dircoord_ = new DirectionCoordinate(mdt, projType,
345 centX, centY, incX, incY,
346 xform,
347 0.5*Double(nx),
348 0.5*Double(ny)) ; // pixel at center
349// 0.5*Double(nx-1),
350// 0.5*Double(ny-1)) ; // pixel at grid
351 {//Summary
352 os << "Successfully generated grid coordinate:" << LogIO::POST;
353 Vector<String> units = dircoord_->worldAxisUnits();
354 Vector<Double> refv = dircoord_->referenceValue();
355 os <<"- Reference Direction : " << MDirection::showType(dircoord_->directionType())
356 << " " << refv[0] << units[0] << " " << refv[1] << units[1] << LogIO::POST;
357 Vector<Double> refpix = dircoord_->referencePixel();
358 os <<"- Reference Pixel : [" << refpix[0] << ", " << refpix[1] << "]" << LogIO::POST;
359 Vector<Double> inc = dircoord_->increment();
360 os <<"- Increments : [" << inc[0] << ", " << inc[1] << "]" << LogIO::POST;
361 os <<"- Projection Type : " << dircoord_->projection().name() << LogIO::POST;
362 }
363};
364
365vector<double> PlotHelper::getGridPixel(const int whichrow){
366 if (data_p->nrow() < 1)
367 throw AipsError("Scantable is not set. Could not get direction.");
368 else if (whichrow > int(data_p->nrow()) - 1)
369 throw AipsError("Row index out of range.");
370 if (!dircoord_)
371 throw AipsError("Direction coordinate is not defined.");
372
373 Vector<Double> pixel;
374 MDirection world;
375 vector<double> outvec;
376 world = data_p->getDirection(whichrow);
377#ifdef KS_DEBUG
378 cout << "searching pixel position (world = " << data_p->getDirectionString(whichrow) << " = [" << world.getAngle("rad").getValue()[0] << ", " << world.getAngle("rad").getValue()[1] << "])" << endl;
379#endif
380 dircoord_->toPixel(pixel, world);
381#ifdef KS_DEBUG
382 cout << "got pixel = [" << pixel[0] << ", " << pixel[1] << "]" << endl;
383#endif
384 // convert pixel to std vector
385 pixel.tovector(outvec);
386 return outvec;
387};
388
389} //namespace asap
Note: See TracBrowser for help on using the repository browser.