source: trunk/src/PlotHelper.cpp@ 2784

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

New Development: No

JIRA Issue: No

Ready for Test: Yes

Interface Changes: Yes

What Interface Changed: added two new methods, get_gref() and get_gcellval() in plothelper class.

Test Programs: test_sdplot[sdplot_gridTest]

Put in Release Notes: No

Module(s): asapplotter and sdplot

Description:

Added two new methods in plothelper class to get grid center
and spacings of plot grids.

  • get_gref(): returns the reference direction (grid center) of grid (a string)
  • get_gcellval(): returns the absolute increment (in radian) of grid (a double vector)


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