source: trunk/src/PlotHelper.cpp@ 3137

Last change on this file since 3137 was 3106, checked in by Takeshi Nakazato, 8 years ago

New Development: No

JIRA Issue: No

Ready for Test: Yes/No

Interface Changes: Yes/No

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes/No

Module(s): Module Names change impacts.

Description: Describe your changes here...


Check-in asap modifications from Jim regarding casacore namespace conversion.

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