Changeset 2717 for trunk/src


Ignore:
Timestamp:
01/09/13 12:43:18 (12 years ago)
Author:
Kana Sugimoto
Message:

New Development: Yes

JIRA Issue: No

Ready for Test: Yes

Interface Changes: Yes

What Interface Changed: added a new method, set_grid(), in plothelper class.

Test Programs: test_sdplot[sdplot_gridTest]

Put in Release Notes: Yes

Module(s): asapplotter, sdplot(plottype='grid')

Description:

The method, asapplotter.plotgrid(), considers tangential projection
effects when getting grid center and extents from a scantable.
This affects behavior of plottype='grid' in sdplot task of CASA.

Location:
trunk/src
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/PlotHelper.cpp

    r2690 r2717  
    1515#include <casa/Quanta/Quantum.h>
    1616#include <casa/Quanta/QuantumHolder.h>
     17#include <casa/Logging/LogIO.h>
     18
    1719#include <measures/Measures/MDirection.h>
    18 #include <casa/Logging/LogIO.h>
     20#include <tables/Tables/TableRecord.h>
     21
    1922
    2023#include "PlotHelper.h"
    2124
     25//#ifndef KS_DEBUG
    2226//#define KS_DEBUG
     27//#endif
    2328
    2429using namespace std ;
    … …  
    3136{
    3237#ifdef KS_DEBUG
    33   cout << "Default constructor nrow = " << data_.nrow() << endl;
     38  cout << "Default constructor nrow = " << data_p->nrow() << endl;
    3439#endif
    3540};
    … …  
    6166  cout << "Setting scantable" << endl;
    6267#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 
    145 void 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){
     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 << "Start getSTCoord()" << 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 << "Generating a coordinate from DIRECTIONs in 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  os << "Start setGridParam()" << LogIO::POST;
     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_){
     136#ifdef KS_DEBUG
     137    cout << "Destructing dircoord_" << endl;
     138#endif
     139    delete dircoord_;
     140    dircoord_ = 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  // 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    if (!stset)
     202      throw AipsError("Scantable is not set. Could not resolve map center.");
     203#ifdef KS_DEBUG
     204    cout << "Using pointing center from DIRECTION column" << endl;
     205#endif
     206    centx = stcent.getAngle("rad").getValue()[0];
     207    centy = stcent.getAngle("rad").getValue()[1];
     208    mdt = stdt;
     209  } else {
     210#ifdef KS_DEBUG
     211    cout << "Using user defined grid center" << endl;
     212#endif
     213    // Parse center string
     214    string::size_type pos0 = center.find(" ");
     215   
     216    if (pos0 == string::npos)
     217      throw AipsError("bad string format in center direction");
     218
     219    string::size_type pos1 = center.find(" ", pos0+1);
     220    String sepoch, sra, sdec;
     221    if (pos1 != string::npos) {
     222      sepoch = center.substr(0, pos0);
     223      sra = center.substr(pos0+1, pos1-pos0);
     224      sdec = center.substr(pos1+1);
     225    } else {
     226      sepoch = "J2000";
     227      sra = center.substr(0, pos0);
     228      sdec = center.substr(pos0+1);
     229    }
     230    if (!MDirection::getType(mdt,sepoch))
     231      throw AipsError("Invalid direction reference in center");
     232    if (stset && mdt != stdt)
     233      throw AipsError("Direction reference of center should be the same as input scantable");
     234    QuantumHolder qh ;
     235    String err ;
     236    qh.fromString(err, sra);
     237    Quantum<Double> ra  = qh.asQuantumDouble();
     238    qh.fromString(err, sdec) ;
     239    Quantum<Double> dec = qh.asQuantumDouble();
     240    centx = ra.getValue("rad");
     241    centy = dec.getValue("rad");
     242    // rotaion
     243    if (stset) {
     244      Double stcentx = stcent.getAngle("rad").getValue()[0];
     245      Double rotnum = round( abs(centx - stcentx) / (C::_2pi) );
     246      if (centx < stcentx) rotnum *= -1;
     247      centx -= (rotnum * C::_2pi);
     248    }
     249  }
     250#ifdef KS_DEBUG
     251  cout << "The center direction of plotting grids: [" << centx  << ", " << centy << "] (rad)" <<endl;
     252#endif
     253
     254  // cell
     255  if (cellx.empty() && celly.empty()){
     256    if (!stset)
     257      throw AipsError("Scantable is not set. Could not resolve cell size.");
     258#ifdef KS_DEBUG
     259    cout << "Using cell size defined from DIRECTION column" << endl;
     260#endif
     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#ifdef KS_DEBUG
     276    cout << "Using user defined cell size" << endl;
     277#endif
     278    Quantum<Double> qcellx, qcelly;
     279    if (!cellx.empty() && !celly.empty()){
     280      readQuantity(qcellx, String(cellx));
     281      readQuantity(qcelly, String(celly));
     282    } else if (celly.empty()) {
     283      readQuantity(qcellx, String(cellx));
     284      qcelly = qcellx;
     285    } else { //only celly is defined
     286      readQuantity(qcelly, String(celly));
     287      qcellx = qcelly;
     288    }
     289    incx = qcellx.getValue("rad");
     290    incy = qcelly.getValue("rad");
     291  }
     292#ifdef KS_DEBUG
     293  cout << "The cell size of plotting grids: [" << incx << ", " << incy << "] (rad)" <<endl;
     294#endif
     295
     296  Matrix<Double> xform(2,2) ;
     297  xform = 0.0 ;
     298  xform.diagonal() = 1.0 ;
     299  dircoord_ = new DirectionCoordinate(mdt, projtype,
     300                                      centx, centy, incx, incy,
     301                                      xform,
     302                                      0.5*Double(nx),
     303                                      0.5*Double(ny)) ; // pixel at center
     304  os << "Successfully finished generation of Direction Coordinate" << LogIO::POST;
     305};
     306
     307void PlotHelper::setGridParamVal(const int nx, const int ny,
     308                                 const double cellx, const double celly,
     309                                 const double centx, const double centy,
     310                                 const string epoch, const string projname){
    146311  // Value check of nx and ny
    147312  if (nx < 1)
    … …  
    194359
    195360vector<double>  PlotHelper::getGridPixel(const int whichrow){
    196   if (data_.nrow() < 1)
     361  if (data_p->nrow() < 1)
    197362    throw AipsError("Scantable is not set. Could not get direction.");
    198   else if (whichrow > int(data_.nrow()) - 1)
     363  else if (whichrow > int(data_p->nrow()) - 1)
    199364    throw AipsError("Row index out of range.");
    200365  if (!dircoord_)
    … …  
    204369  MDirection world;
    205370  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;
     371  world = data_p->getDirection(whichrow);
     372#ifdef KS_DEBUG
     373  cout << "searching pixel position (world = " << data_p->getDirectionString(whichrow) << " = [" << world.getAngle("rad").getValue()[0] << ", " << world.getAngle("rad").getValue()[1] << "])" << endl;
    209374#endif
    210375  dircoord_->toPixel(pixel, world);

  • trunk/src/PlotHelper.h

    r2689 r2717  
    2323
    2424#include "ScantableWrapper.h"
     25#include "Scantable.h"
    2526
    2627namespace asap {
    … …  
    5253   * Set grid parameters for plot panel (by quantity)
    5354   **/
    54 /*   void setGridParam(const int nx, const int ny,  */
    55 /*                  const string cellx="", const string celly="", */
    56 /*                  string center="", const string projname="SIN") ; */
     55  void setGridParam(const int nx, const int ny,
     56                    const string cellx="", const string celly="",
     57                    string center="", const string projname="SIN") ;
    5758
    5859  /**
    … …  
    6162  vector<double> getGridPixel(const int whichrow=0);
    6263
    63  private:
     64private:
     65  casa::DirectionCoordinate getSTCoord(const int nx, const int ny,
     66                                       const casa::Projection::Type ptype);
     67
    6468  casa::DirectionCoordinate *dircoord_;
    65   ScantableWrapper data_;
     69  casa::CountedPtr<Scantable> data_p;
    6670
    6771};

  • trunk/src/python_PlotHelper.cpp

    r2711 r2717  
    2626         (boost::python::arg("epoch")="J2000",
    2727          boost::python::arg("projname")="SIN"))
    28 //     .def("set_gird", &PlotHelper::setGridParam,
    29 //          (boost::python::arg("scellx")="",
    30 //           boost::python::arg("scelly")="",
    31 //           boost::python::arg("scenter")=""))
     28    .def("set_grid", &PlotHelper::setGridParam,
     29         (boost::python::arg("cellx")="",
     30          boost::python::arg("celly")="",
     31          boost::python::arg("center")="",
     32          boost::python::arg("projname")="SIN") )
    3233    .def("get_gpos", &PlotHelper::getGridPixel,
    3334         (boost::python::arg("whichrow")=0) )
Note: See TracChangeset for help on using the changeset viewer.