Changeset 2717 for trunk


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
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • trunk/python/asapplotter.py

    r2715 r2717  
    17671767    # plot spectra by pointing
    17681768    @asaplog_post_dec
    1769     def plotgrid(self, scan=None,center=None,spacing=None,rows=None,cols=None):
     1769    def plotgrid(self, scan=None,center="",spacing=[],rows=None,cols=None):
    17701770        """
    17711771        Plot spectra based on direction.
    … …  
    17731773        Parameters:
    17741774            scan:      a scantable to plot
    1775             center:    the grid center direction (a list) of plots in the
    1776                        unit of DIRECTION column.
     1775            center:    the grid center direction (a string)
    17771776                       (default) the center of map region
     1777                       (example) 'J2000 19h30m00s -25d00m00s'
    17781778            spacing:   a list of horizontal (R.A.) and vertical (Dec.)
    1779                        spacing in the unit of DIRECTION column.
     1779                       spacing.
    17801780                       (default) Calculated by the extent of map region and
     1781                       (example) ['1arcmin', '1arcmin']
    17811782                       the number of rows and cols to cover
    17821783            rows:      number of panels (grid points) in horizontal direction
    … …  
    18021803
    18031804        # Rows and cols
    1804         if rows:
    1805             self._rows = int(rows)
    1806         else:
    1807             msg = "Number of rows to plot are not specified. "
    1808             if self._rows:
    1809                 msg += "Using previous value = %d" % (self._rows)
    1810                 asaplog.push(msg)
    1811             else:
    1812                 self._rows = 1
    1813                 msg += "Setting rows = %d" % (self._rows)
    1814                 asaplog.post()
    1815                 asaplog.push(msg)
    1816                 asaplog.post("WARN")
    1817         if cols:
    1818             self._cols = int(cols)
    1819         else:
    1820             msg = "Number of cols to plot are not specified. "
    1821             if self._cols:
    1822                 msg += "Using previous value = %d" % (self._cols)
    1823                 asaplog.push(msg)
    1824             else:
    1825                 self._cols = 1
    1826                 msg += "Setting cols = %d" % (self._cols)
    1827                 asaplog.post()
    1828                 asaplog.push(msg)
    1829                 asaplog.post("WARN")
    1830 
    1831         # Center and spacing
    1832         dirarr = array(self._data.get_directionval()).transpose()
    1833         print "Pointing range: (x, y) = (%f - %f, %f - %f)" %\
    1834               (dirarr[0].min(),dirarr[0].max(),dirarr[1].min(),dirarr[1].max())
    1835         dircent = [0.5*(dirarr[0].max() + dirarr[0].min()),
    1836                    0.5*(dirarr[1].max() + dirarr[1].min())]
    1837         del dirarr
    1838         if center is None:
    1839             #asaplog.post()
    1840             asaplog.push("Grid center is not specified. Automatically calculated from pointing center.")
    1841             #asaplog.post("WARN")
    1842             #center = [dirarr[0].mean(), dirarr[1].mean()]
    1843             center = dircent
    1844         elif (type(center) in (list, tuple)) and len(center) > 1:
    1845             from numpy import pi
    1846             # make sure center_x is in +-pi of pointing center
    1847             # (assumes dirs are in rad)
    1848             rotnum = round(abs(center[0] - dircent[0])/(2*pi))
    1849             if center[0] < dircent[0]: rotnum *= -1
    1850             cenx = center[0] - rotnum*2*pi
    1851             center = [cenx, center[1]]
    1852         else:
    1853             msg = "Direction of grid center should be a list of float (R.A., Dec.)"
    1854             raise ValueError, msg
    1855         asaplog.push("Grid center: (%f, %f) " % (center[0],center[1]))
    1856 
    1857         if spacing is None:
    1858             #asaplog.post()
    1859             asaplog.push("Grid spacing not specified. Automatically calculated from map coverage")
    1860             #asaplog.post("WARN")
    1861             # automatically get spacing
    1862             dirarr = array(self._data.get_directionval()).transpose()
    1863             wx = 2. * max(abs(dirarr[0].max()-center[0]),
    1864                           abs(dirarr[0].min()-center[0]))
    1865             wy = 2. * max(abs(dirarr[1].max()-center[1]),
    1866                           abs(dirarr[1].min()-center[1]))
    1867             ## slightly expand area to plot the edges
    1868             #wx *= 1.1
    1869             #wy *= 1.1
    1870             xgrid = wx/max(self._cols-1.,1.)
    1871             #xgrid = wx/float(max(self._cols,1.))
    1872             xgrid *= cos(center[1])
    1873             ygrid = wy/max(self._rows-1.,1.)
    1874             #ygrid = wy/float(max(self._rows,1.))
    1875             # single pointing (identical R.A. and/or Dec. for all spectra.)
    1876             if xgrid == 0:
    1877                 xgrid = 1.
    1878             if ygrid == 0:
    1879                 ygrid = 1.
    1880             # spacing should be negative to transpose plot
    1881             spacing = [- xgrid, - ygrid]
    1882             del dirarr, xgrid, ygrid
    1883         #elif isinstance(spacing, str):
    1884         #    # spacing is a quantity
    1885         elif (type(spacing) in (list, tuple)) and len(spacing) > 1:
    1886             for i in xrange(2):
    1887                 val = spacing[i]
    1888                 if not isinstance(val, float):
    1889                     raise TypeError("spacing should be a list of float")
    1890                 if val > 0.:
    1891                     spacing[i] = -val
    1892             spacing = spacing[0:2]
    1893         else:
    1894             msg = "Invalid spacing."
    1895             raise TypeError(msg)
    1896         asaplog.push("Spacing: (%f, %f) (projected)" % (spacing[0],spacing[1]))
    1897 
     1805        if (self._rows is None):
     1806            rows = max(1, rows)
     1807        if (self._cols is None):
     1808            cols = max(1, cols)
     1809        self.set_layout(rows,cols,False)
     1810
     1811        # Select the first IF, POL, and BEAM for plotting
    18981812        ntotpl = self._rows * self._cols
    18991813        ifs = self._data.getifnos()
    … …  
    19221836            asaplog.push(msg)
    19231837            asaplog.post("WARN")
    1924        
     1838
     1839#         # Center and spacing
     1840#         dirarr = array(self._data.get_directionval()).transpose()
     1841#         print "Pointing range: (x, y) = (%f - %f, %f - %f)" %\
     1842#               (dirarr[0].min(),dirarr[0].max(),dirarr[1].min(),dirarr[1].max())
     1843#         dircent = [0.5*(dirarr[0].max() + dirarr[0].min()),
     1844#                    0.5*(dirarr[1].max() + dirarr[1].min())]
     1845#         del dirarr
     1846#         if center is None:
     1847#             #asaplog.post()
     1848#             asaplog.push("Grid center is not specified. Automatically calculated from pointing center.")
     1849#             #asaplog.post("WARN")
     1850#             #center = [dirarr[0].mean(), dirarr[1].mean()]
     1851#             center = dircent
     1852#         elif (type(center) in (list, tuple)) and len(center) > 1:
     1853#             from numpy import pi
     1854#             # make sure center_x is in +-pi of pointing center
     1855#             # (assumes dirs are in rad)
     1856#             rotnum = round(abs(center[0] - dircent[0])/(2*pi))
     1857#             if center[0] < dircent[0]: rotnum *= -1
     1858#             cenx = center[0] - rotnum*2*pi
     1859#             center = [cenx, center[1]]
     1860#         else:
     1861#             msg = "Direction of grid center should be a list of float (R.A., Dec.)"
     1862#             raise ValueError, msg
     1863#         asaplog.push("Grid center: (%f, %f) " % (center[0],center[1]))
     1864
     1865#         if spacing is None:
     1866#             #asaplog.post()
     1867#             asaplog.push("Grid spacing not specified. Automatically calculated from map coverage")
     1868#             #asaplog.post("WARN")
     1869#             # automatically get spacing
     1870#             dirarr = array(self._data.get_directionval()).transpose()
     1871#             wx = 2. * max(abs(dirarr[0].max()-center[0]),
     1872#                           abs(dirarr[0].min()-center[0]))
     1873#             wy = 2. * max(abs(dirarr[1].max()-center[1]),
     1874#                           abs(dirarr[1].min()-center[1]))
     1875#             ## slightly expand area to plot the edges
     1876#             #wx *= 1.1
     1877#             #wy *= 1.1
     1878#             xgrid = wx/max(self._cols-1.,1.)
     1879#             #xgrid = wx/float(max(self._cols,1.))
     1880#             xgrid *= cos(center[1])
     1881#             ygrid = wy/max(self._rows-1.,1.)
     1882#             #ygrid = wy/float(max(self._rows,1.))
     1883#             # single pointing (identical R.A. and/or Dec. for all spectra.)
     1884#             if xgrid == 0:
     1885#                 xgrid = 1.
     1886#             if ygrid == 0:
     1887#                 ygrid = 1.
     1888#             # spacing should be negative to transpose plot
     1889#             spacing = [- xgrid, - ygrid]
     1890#             del dirarr, xgrid, ygrid
     1891#         #elif isinstance(spacing, str):
     1892#         #    # spacing is a quantity
     1893#         elif (type(spacing) in (list, tuple)) and len(spacing) > 1:
     1894#             for i in xrange(2):
     1895#                 val = spacing[i]
     1896#                 if not isinstance(val, float):
     1897#                     raise TypeError("spacing should be a list of float")
     1898#                 if val > 0.:
     1899#                     spacing[i] = -val
     1900#             spacing = spacing[0:2]
     1901#         else:
     1902#             msg = "Invalid spacing."
     1903#             raise TypeError(msg)
     1904#         asaplog.push("Spacing: (%f, %f) (projected)" % (spacing[0],spacing[1]))
     1905
     1906        # Prepare plotter
    19251907        self._assert_plotter(action="reload")
    19261908        self._plotter.hold()
    … …  
    19381920        from asap._asap import plothelper as plhelper
    19391921        ph = plhelper(self._data)
    1940         ph.set_gridval(self._cols, self._rows, spacing[0], spacing[1],
    1941                           center[0], center[1], epoch="J2000", projname="SIN")
     1922        #ph.set_gridval(self._cols, self._rows, spacing[0], spacing[1],
     1923        #                  center[0], center[1], epoch="J2000", projname="SIN")
     1924        if type(spacing) in (list, tuple, array):
     1925            if len(spacing) == 0:
     1926                spacing = ["", ""]
     1927            elif len(spacing) == 1:
     1928                spacing = [spacing[0], spacing[0]]
     1929        else:
     1930            spacing = [spacing, spacing]
     1931        ph.set_grid(self._cols, self._rows, spacing[0], spacing[1], \
     1932                    center, projname="SIN")
     1933
    19421934        # Actual plot
    19431935        npl = 0

  • 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.