Logo Search packages:      
Sourcecode: magics++ version File versions  Download package

PolarStereographicProjection.cc

Go to the documentation of this file.
/******************************** LICENSE ********************************


 Copyright 2007 European Centre for Medium-Range Weather Forecasts (ECMWF)
 
 Licensed under the Apache License, Version 2.0 (the "License"); 
 you may not use this file except in compliance with the License. 
 You may obtain a copy of the License at 
 
      http://www.apache.org/licenses/LICENSE-2.0
 
 Unless required by applicable law or agreed to in writing, software 
 distributed under the License is distributed on an "AS IS" BASIS, 
 WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 
 See the License for the specific language governing permissions and 
 limitations under the License.


 ******************************** LICENSE ********************************/

/*! \file PolarStereographicProjection.cc
    \brief Implementation of PolarStereographicProjection.
    \author Graphics Section, ECMWF

    Started: Fri Jan 11 15:01:45 2008

*/

#include <PolarStereographicProjection.h>
#include <BasicSceneObject.h>
#include <GridPlotting.h>
#include <LabelPlotting.h>
#include <Text.h>
#include <MatrixHandler.h>
#include <MetaData.h>
using namespace magics;

/*!
  \brief Constructor
*/
00041 PolarStereographicProjection::PolarStereographicProjection()  : projection_(0)
{
}

/*!
  \brief Destructor
*/
00048 PolarStereographicProjection::~PolarStereographicProjection()
{
      delete projection_;
}

00053 void PolarStereographicProjection::print(ostream& out) const
{
      out << "PolarStereographicProjection[";
      PolarStereographicProjectionAttributes::print(out);
      out << "]"; 
} 

double PolarStereographicProjection::unitToCm(double width, double height) const
{
      // Reproject 2 point and find the unit
      int hemis = (hemisphere_ ==  NORTH) ? 1 : -1;

      GeoPoint ll1(0,  hemis*50);
      GeoPoint ll2(0,  hemis*51);
      PaperPoint xy1 = (*this)(ll1);
      PaperPoint xy2 = (*this)(ll2);
      
      double unit = abs(xy1.y() - xy2.y());
      
      return unit*(height/(getMaxY() -getMinY()) ); 
}


00076 PaperPoint PolarStereographicProjection::operator()(const UserPoint& point)  const
{
      Log::warning() << "PolarStereographicProjection::operator()(const UserPoint& point) should not be called!" << endl;
      return Transformation::operator()(point);
}

00082 PaperPoint PolarStereographicProjection::operator()(const GeoPoint& point)  const
{
      assert(projection_);
      TeCoord2D geo = TeCoord2D(point.longitude()*TeCDR, point.latitude()*TeCDR);
      TeCoord2D xy = projection_->LL2PC(geo);

      return PaperPoint(xy.x(), xy.y(), point.value(), point.missing(), point.border());
}

00091 PaperPoint PolarStereographicProjection::operator()(const PaperPoint& point)  const
{
      return Transformation::operator()(point);
}

00096 void PolarStereographicProjection::revert(const PaperPoint& xy, GeoPoint& point)  const
{
      assert(projection_);
      TeCoord2D texy = TeCoord2D(xy.x(), xy.y());
      TeCoord2D geo = projection_->PC2LL(texy);
  
      point = GeoPoint(geo.x()*TeCRD, geo.y()*TeCRD);
}

00105 void PolarStereographicProjection::revert(const PaperPoint& xy, UserPoint& point)  const
{
      Log::dev() << "PolarStereographicProjection::revert(...) needs implementing." << endl;
      Transformation::revert(xy, point);
}

00111 bool PolarStereographicProjection::needShiftedCoastlines()  const
{
      return false;
}

00116 void PolarStereographicProjection::init(double width, double height)
{
      if ( !projection_ ) 
            projection_ = new TePolarStereographic(TeDatum(), vertical_longitude_*TeCDR, 0., 0., "Meters", (hemisphere_ == NORTH) ? TeNORTH_HEM : TeSOUTH_HEM);
      
      magCompare(area_, "corners" ) ?  corners() : centre(width, height);
      double llx, urx , lly , ury;
            
      if ( magCompare(system_, "projection") == false )
      {
            TeCoord2D ll = TeCoord2D(xmin_*TeCDR, ymin_*TeCDR);
            TeCoord2D ur = TeCoord2D(xmax_*TeCDR, ymax_*TeCDR);

            TeCoord2D llxy = projection_->LL2PC(ll);
            TeCoord2D urxy = projection_->LL2PC(ur);
    
            TeCoord2D ell = TeCoord2D(-20*TeCDR, 40*TeCDR);
            TeCoord2D exy = projection_->LL2PC(ell);
    
            ell = projection_->PC2LL(exy);
    
            // Create a grid 100/100 to calculate the corners ...
     
            llx = ::min(urxy.x(), llxy.x());
            urx = ::max(urxy.x(), llxy.x());

            if (llx == urx)
            urx = llx+1000;

            lly = ::min(urxy.y(), llxy.y());
            ury = ::max(urxy.y(), llxy.y());

            if (lly == ury)
            ury = lly+1000;
            xmin_ = int_MAX;
            ymin_ = int_MAX;
            xmax_ = int_MIN;
            ymax_ = int_MIN;

            double stepx= (urx - llx)/100;
            double stepy= (ury - lly)/100;
            for (double x = llx; x <= urx; x += stepx )
            {
               for (double y = lly; y <= ury; y += stepy )
               {
                   TeCoord2D xy(x, y);
                   TeCoord2D ll = projection_->PC2LL(xy);
                   double xx = ll.x()*TeCRD;
                   double yy = ll.y()*TeCRD;

                   if (xx < xmin_) xmin_ = xx;
                   if (xx > xmax_) xmax_ = xx;
                   if (yy < ymin_) ymin_ = yy;
                   if (yy > ymax_) ymax_ = yy;
               }
            }

            if (xmax_ - xmin_ > 358)
            {
                  // Wrap around!
                  xmax_=180;
                  xmin_=-180; 
            }
      
            xpcmin_ = llx;
            ypcmin_ = lly;
            xpcmax_ = urx;
            ypcmax_ = ury; 
      }
      else {
            llx = ::min(xmin_, xmax_);
            urx = ::max(xmin_, xmax_);
            lly = ::min(ymin_, ymax_);
            ury = ::max(ymin_, ymax_);
            TeCoord2D llxy(llx, lly);
            TeCoord2D ll = projection_->PC2LL(llxy); 
            TeCoord2D urxy(urx, ury);
            TeCoord2D  ur = projection_->PC2LL(urxy); 
                     
            xmin_ =  ::min(ll.x()*TeCRD, ur.x()*TeCRD);
            xmax_ = ::max(ll.x()*TeCRD, ur.x()*TeCRD); 
            ymin_ =  ::min(ll.y()*TeCRD, ur.y()*TeCRD);
            ymax_ =  ::max(ll.y()*TeCRD, ur.y()*TeCRD);
            double stepx= (urx - llx)/100;
            double stepy= (ury - lly)/100;
            for (double x = llx; x <= urx; x += stepx )
            {
               for (double y = lly; y <= ury; y += stepy )
               {
                    TeCoord2D xy(x, y);
                    TeCoord2D ll = projection_->PC2LL(xy);
                    double xx = ll.x()*TeCRD;
                    double yy = ll.y()*TeCRD;

                    if (xx < xmin_) xmin_ = xx;
                    if (xx > xmax_) xmax_ = xx;
                    if (yy < ymin_) ymin_ = yy;
                    if (yy > ymax_) ymax_ = yy;
               }
            }

            if (xmax_ - xmin_ > 358)
            {
                  // Wrap around!
                  xmax_=180;
                  xmin_=-180; 
            }

            xpcmin_ = llx;
            ypcmin_ = lly;
            xpcmax_ = urx;
            ypcmax_ = ury; 
            Log::dev() << " Projection definition-->[" << ymin_ << ", " << xmin_ << ", " << xmax_ << ", " << ymax_ << "]" << endl;
      }
}
void PolarStereographicProjection::fill(double& width, double& height)  
{  
      init(width, height);
    Transformation::fill(width, height);
    
    

}

00240 void PolarStereographicProjection::aspectRatio(double& width, double& height)  
{  
      init(width, height);
      Transformation::aspectRatio(width, height);
}

00246 void PolarStereographicProjection::boundingBox(double& xmin, double& ymin, double& xmax, double& ymax)  const
{
      xmin=-180;
      ymin=-90;
      xmax=180;
      ymax=90;
}

00254 double PolarStereographicProjection::getMinX()  const
{
      return xmin_;
}

00259 double PolarStereographicProjection::getMinY()  const
{
      return ymin_;
}

00264 double PolarStereographicProjection::getMaxX()  const
{
      return xmax_;
}

00269 double PolarStereographicProjection::getMaxY()  const
{
      return ymax_;
}

00274 void PolarStereographicProjection::setMinX(double x)  
{
      Log::dev() << "PolarStereographicProjection::setMinX(...) needs implementing." << endl;
      Transformation::setMinX(x);
}

00280 void PolarStereographicProjection::setMinY(double y)  
{
      Log::dev() << "PolarStereographicProjection::setMinY(...) needs implementing." << endl;
      Transformation::setMinY(y);
}

00286 void PolarStereographicProjection::setMaxX(double x)  
{
      Log::dev() << "PolarStereographicProjection::setMaxX(...) needs implementing." << endl;
      Transformation::setMaxX(x);
}

00292 void PolarStereographicProjection::setMaxY(double y)  
{
      Log::dev() << "PolarStereographicProjection::setMaxY(...) needs implementing." << endl;
      Transformation::setMaxY(y);
}

00298 double PolarStereographicProjection::getMinPCX()  const
{
      return xpcmin_;
}

00303 double PolarStereographicProjection::getMinPCY()  const
{
      return ypcmin_;
}

00308 double PolarStereographicProjection::getMaxPCX()  const
{
      return xpcmax_;
}

00313 double PolarStereographicProjection::getMaxPCY()  const
{
      return ypcmax_;
}

00318 void PolarStereographicProjection::gridLongitudes(const GridPlotting& grid)  const
{
      const vector<double>& longitudes = grid.longitudes();
      const vector<double>& latitudes = grid.latitudes();
      for (vector<double>::const_iterator lon = longitudes.begin(); lon != longitudes.end(); ++lon)
      {
            Polyline poly;

            double min = ::min(latitudes.front(), latitudes.back());
            double max = 80;
            
            for (double lat = min; lat <= max; lat += 1) { 
                  poly.push_back((*this)(GeoPoint(*lon,lat)));
            }  
            grid.add(poly);     
      }
}

00336 void PolarStereographicProjection::gridLatitudes(const GridPlotting& grid)  const
{
      const vector<double>& latitudes = grid.latitudes();
      for (vector<double>::const_iterator lat = latitudes.begin(); lat != latitudes.end(); ++lat)
      {
          Polyline poly;        
            for (int lon = -180; lon <= 180; lon += 1) { 
                poly.push_back((*this)(GeoPoint(lon,*lat)));
                
            }     
          grid.add(poly);        
      }
}

inline double CA(TeCoord2D& p1, TeCoord2D& p2)
{
    return (p2.y() - p1.y())/(p2.x() - p1.x());
}

inline double CB(double a, TeCoord2D& p)
{
    return p.y() - a * p.x();
}

inline double CX(double a, double b, double y)
{
    return (y - b)/a;
}

inline double CY(double a, double b, double x)
{
    return (a * x) + b;
}

void PolarStereographicProjection::horizontalLabels(const LabelPlotting& label, double y, double yy, bool top) const
{
    const vector<double>& longitudes = label.longitudes();
    for (vector<double>::const_iterator lon = longitudes.begin(); lon != longitudes.end(); ++lon)
    {
        // find the equation of the line using 2 points : lon/-20 -->lon/ +20
        TeCoord2D p1 = TeCoord2D((*lon)*TeCDR, 20*TeCDR);
        TeCoord2D p2 = TeCoord2D((*lon)*TeCDR, -20*TeCDR);

        TeCoord2D p1xy = projection_->LL2PC(p1);
        TeCoord2D p2xy = projection_->LL2PC(p2);
    
        double a = CA(p1xy, p2xy);
        double b = CB(a, p1xy);
            
        TeCoord2D xy = TeCoord2D(CX(a, b, y), y);   

        PaperPoint point(CX(a, b, y), y);
        if ( !in(point) ) continue;
            GeoPoint geo;
            revert(point, geo);
            point.y(yy);                  
        Text* text = new Text();
        text->setJustification(MCENTRE);
        text->setVerticalAlign(MBOTTOM);
      if(top) text->setVerticalAlign(MBOTTOM);
      else    text->setVerticalAlign(MTOP);     
        text->setText(geo.writeLongitude());
        text->push_back(point);
                  
        label.add(text);         
     }
}

/*!
 Method to draw vertical labels.
 
 \sa LeftAxisVisitor RightAxisVisitor
 
*/
00410 void PolarStereographicProjection::verticalLabels(const LabelPlotting& label, double x, double xx, bool left)  const
{
    const vector<double>& longitudes = label.longitudes();
    for (vector<double>::const_iterator lon = longitudes.begin(); lon != longitudes.end(); ++lon)
    {
        // find the equation of the line using 2 points : lon/-20 -->lon/ +20
        TeCoord2D p1 = TeCoord2D((*lon)*TeCDR, 20*TeCDR);
        TeCoord2D p2 = TeCoord2D((*lon)*TeCDR, -20*TeCDR);

        TeCoord2D p1xy = projection_->LL2PC(p1);
        TeCoord2D p2xy = projection_->LL2PC(p2);
    
        double a = CA(p1xy, p2xy);
        double b = CB(a, p1xy);
        
        TeCoord2D xy = TeCoord2D(x, CY(a, b, x) );       

        PaperPoint point(x, CY(a, b, x));
            if ( !in(point) ) continue;
        GeoPoint geo;
        revert(point, geo);
        point.x(xx);
            
        Text* text = new Text();
      if(left) text->setJustification(MRIGHT);
      else     text->setJustification(MLEFT);
        text->setVerticalAlign(MHALF);    
        text->setText(geo.writeLongitude());
        text->push_back(point);
      
        label.add(text);
     }
}


00445 void PolarStereographicProjection::labels(const LabelPlotting& label, DrawingVisitor& )  const
{
      Text *text;
      const vector<double>& latitudes = label.latitudes();
      const vector<double>& longitudes = label.longitudes();
      
      unsigned int flat = (unsigned int) std::max(1, (int) maground(latitudes.size()/4));
      unsigned int flon = (unsigned int) std::max(1, (int) maground(longitudes.size()/4));
      
      
      
      for (unsigned int lat = 0; lat < latitudes.size(); lat += flat)
      {  
          for (unsigned int lon = 0 ; lon < longitudes.size(); lon += flon)
          {          
               GeoPoint point(longitudes[lon],latitudes[lat]);
               PaperPoint xy = (*this)(point);
               
               if ( !in(xy) ) continue;      
         
               text = new Text();
               text->setText(point.writeLatitude());
                 text->push_back(xy);
                 text->setBlanking(false);
               label.add(text);
          }
      }
}
void PolarStereographicProjection::labels(const LabelPlotting& label, TopAxisVisitor&)  const
{
      // Find intersection of the latitude line  with the top!
      // Y = top

      const double yy = ypcmin_ + ( ( ypcmax_ - ypcmin_  )*0.1);
      horizontalLabels(label, ypcmax_, yy, true);
}

void PolarStereographicProjection::labels(const LabelPlotting& label, BottomAxisVisitor&)  const
{
      // Find intersection of the latitude line  with the bottom!
      // Y = bottom
     
      const double yy = ypcmax_ - ( ( ypcmax_ - ypcmin_  )*0.1);
      horizontalLabels(label, ypcmin_, yy, false);
}

void PolarStereographicProjection::labels(const LabelPlotting& label, LeftAxisVisitor&)  const
{
      // Find intersection of the latitude line  with the left!
      // X = left
    
      const double xx = xpcmax_ - ( ( xpcmax_ - xpcmin_  )*0.1);
      verticalLabels(label, xpcmin_, xx, true);
}

void PolarStereographicProjection::labels(const LabelPlotting& label, RightAxisVisitor& )  const
{
      // Find intersection of the latitude line  with the right!
      // X = right
      
      const double xx = xpcmin_ + ( ( xpcmax_ - xpcmin_  )*0.1);
      verticalLabels(label, xpcmax_, xx, false);
}

/*!
 This method generates meta output for web interaction through JavaScript...
*/
00512 void  PolarStereographicProjection::visit(MetaDataVisitor& visitor, double left, double top, double width, double height) 
{
      ostringstream java;

      double w = getMaxPCX() - getMinPCX();
      double h = getMaxPCY() - getMinPCY();

      java << "{";
      projection_->LL2PC(java);
            
      java << "\"top\" : " << top <<  ",";            
      java << "\"left\" : " << left <<  ",";          
      java << "\"width\" : " << width <<  ",";  
      java << "\"height\" : " << height <<  ",";      
      
      java << "\"pcxmin\" : " << getMinPCX() <<  ",";       
      java << "\"pcymin\" : " << getMinPCY() <<  ",";       
      java << "\"pcwidth\" : " << w <<  ",";    
      java << "\"pcheight\" : " << h <<  "";    
      java << "}";      
      visitor.add("projection", java.str());
}

void PolarStereographicProjection::corners()
{
      xmin_ = min_longitude_;
      xmax_ = max_longitude_;
      ymin_ = min_latitude_;
      ymax_ = max_latitude_;
            
      if ( magCompare(system_, "projection") ) 
            return; 
      
      double xmin = std::min(min_longitude_, max_longitude_);
      double xmax = std::max(min_longitude_, max_longitude_);
      
      double ymin = std::min(min_latitude_, max_latitude_);
      double ymax = std::max(min_latitude_, max_latitude_);
   
      // change the default ...
       if ( ymin == -90 ) 
             ymin_ = ( hemisphere_ == NORTH ) ? -20. : 20.;
      if ( ymax == 90 ) 
             ymax_ = ( hemisphere_ == NORTH ) ? -20. : 20.;
      if ( xmin == -180 ) 
             xmin_ = ( hemisphere_ == NORTH ) ? -45. : 45.;
      if ( xmax == 180 ) 
             xmax_ = ( hemisphere_ == NORTH ) ? 135. : -135.;   
}



void PolarStereographicProjection::centre(double width, double height)
{
    PaperPoint centre = (*this)(GeoPoint(centre_longitude_, centre_latitude_));
      
    double x = (width* map_scale_)/200;
    double y = (height * map_scale_)/200;
    
           
    PaperPoint llxy(centre.x() - x, centre.y() - y);
    PaperPoint urxy(centre.x() + x, centre.y() + y);
    
    GeoPoint ll;
    revert(PaperPoint(centre.x() - x, centre.y() - y), ll);
    GeoPoint ur;
    revert(PaperPoint(centre.x() + x, centre.y() + y), ur);    

    xmin_ = ll.x();
    ymin_ = ll.y();
    xmax_ = ur.x();
    ymax_ = ur.y();
}

#define ff(p1, p2) sqrt( ((p1.x()-p2.x())*(p1.x()-p2.x())) + ((p1.y()-p2.y())*(p1.y()-p2.y())) )



/*!
 Read in the documentation:   
 For Polar Stereographic projections, the thinning factor is the distance, 
 in both X and Y directions, corresponding to the projected INPUT_FIELD_LONGITUDE_STEP ,
 along 60 degrees latitude, multiplied by the value of WIND_THINNING_FACTOR . 
 After plotting at a grid point, all subsequent grid points, 
 whose distance from the current grid point is less than the thinning factor, will be ignored. 
 The default value is 2.0, e.g. the statement
*/

00600 void PolarStereographicProjection::thin(MatrixHandler<GeoPoint>& matrix, double x, double y, vector<GeoPoint>& out) const
{
      int xfactor = (int) ceil((float) x);
      int yfactor = (int) ceil((float) y);
    
      ThinningMatrixHandler<GeoPoint> original(matrix, 1, 1);
      double lat = hemisphere_ == NORTH ? 90 : -60;
   
      GeoPoint ll1(0, lat);
      GeoPoint ll2(original.XResolution(),  lat);
      PaperPoint xy1 = (*this)(ll1);
      PaperPoint xy2 = (*this)(ll2);
      
      double factor = ff(xy1, xy2)*xfactor;
      
      
      int columns = original.columns();
      int rows = original.rows();

      PaperPoint last, first;
      double savelat = -1;
       
      for ( int row = 0; row < rows; row+=yfactor)
      {
            double savelon = -1;

            for ( int column = 0; column < columns; column++)
            {
                  double lat = original.row(row, column);   
                  double lon = original.column(row, column);
                  if (hemisphere_ == NORTH && lat < -20) continue;
                  if (hemisphere_ == SOUTH && lat > 20) continue;
                  if (savelat != -1)
                  {
                        PaperPoint previous = (*this)(GeoPoint(savelon, savelat));
                        PaperPoint current = (*this)(GeoPoint(lon, lat));
                        if ( ff(previous, current) < factor ) continue;
                        savelat = lat;
                  }                 
                  GeoPoint point(lon, lat, original(row, column));

                  PaperPoint xy = (*this)(point); 
                  if ( column == 0 ) first = xy;

                  if ( savelon == -1)
                  {
                        savelon = lon;
                        savelat = lat;
                        last = xy;
                        if ( in(xy) )
                              out.push_back(point);
                        continue;
                  }
                  if ( ff(last, xy) > factor && ff(first, xy) > factor)
                  {                                   
                        last = xy;
                        if ( in(xy) )
                              out.push_back(point); 
                  }     
            }
      }
}

bool PolarStereographicProjection::concatenate(vector<Polyline* >& lines, Polyline* poly) const 
{
      if (poly->empty()) return true;
      for ( vector<Polyline* >::iterator  line = lines.begin(); line != lines.end(); ++line)
      {
            Polyline* building = *line;
            // concatenate front???
            if ( building->empty() ) continue;
            
            if ( (*line)->front() == poly->back() )
            {
//                Log::dev()<< "concatenate to the front" << (*line)->front() << " == " <<  poly->back() << endl;
                  (*line)->insert((*line)->begin(), poly->begin(), poly->end());
                  lines.erase(line);
                  concatenate(lines, building);
                  return true;                  
            }
            // concatenate back???
            if ( (*line)->back() == poly->front()  )
            {
//                Log::dev()<< "concatenate to the back" << (*line)->back() << " == " <<  poly->front() << endl;
                  (*line)->insert((*line)->end(), poly->begin(), poly->end());
                  lines.erase(line);
                  concatenate(lines, building);             
                  return true;
            }
      }
      
      lines.push_back(poly);
      return true;
}


void PolarStereographicProjection::setNewPCBox(double minx, double miny, double maxx, double maxy)
{
      PaperPoint p1(minx, miny);
      PaperPoint p2(maxx, maxy);
      GeoPoint   ll, ur;

      revert(p1, ll);
      revert(p2, ur);

      cout << ll << endl;
      cout << ur << endl;

      min_longitude_ = ll.x();
      max_longitude_ = ur.x();
      min_latitude_ = ll.y();
      max_latitude_ = ur.y(); 

      corners();
}

void PolarStereographicProjection::operator()(const GeoPoint& geo, vector<PaperPoint>& out) const
{
      PaperPoint pp = (*this)(geo);
            if ( in(pp) ) 
                  out.push_back(pp);
}

void PolarStereographicProjection::reprojectComponents(const GeoPoint& point, pair<double, double>& components) const
{
      double s1 = sqrt((components.first * components.first) + (components.second * components.second));
      GeoPoint u(point.x() + components.first, point.y());
      GeoPoint v(point.x(), point.y() + components.second);

      PaperPoint up = (*this)(u);
      PaperPoint vp = (*this)(v);
      PaperPoint p = (*this)(point);

      components.first = (up.x() - p.x()) + (vp.x() - p.x());
      components.second = (up.y() - p.y()) + (vp.y() - p.y());

      double s2 = sqrt((components.first * components.first) + (components.second * components.second));
      double f = s1/s2;
      components.first *= f;
      components.second *= f; 
}

void PolarStereographicProjection::reprojectSpeedDirection(const PaperPoint& point, pair<double, double>& wind) const
{
      double a = 90 - (wind.second); 
      double speed =wind.first; 
      double x = 3.14/180.;
      a *= x;
    GeoPoint geo;
    wind.first = speed-1 * cos(a);
    wind.second = speed-1 * sin(a);
    reprojectComponents(geo, wind);
}



Generated by  Doxygen 1.6.0   Back to index