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

GeoRectangularProjection.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 GeoRectangularProjection.cc
    \brief Implementation of GeoRectangularProjection.
    \author Graphics Section, ECMWF

    Started: Thu Jan 10 17:24:36 2008

*/

#include <GeoRectangularProjection.h>

#include <Polyline.h>
#include <Text.h>

#include <GeoPoint.h>
#include <SceneVisitor.h>
#include <GridPlotting.h>
#include <LabelPlotting.h>


using namespace magics;

/*!
  \brief Constructor
*/
00045 GeoRectangularProjection::GeoRectangularProjection() : projection_(0)
{
      init();
}

/*!
  \brief Destructor

  \todo do we need here a "delete projection_;" as in MercatorProjection::~MercatorProjection() ?
*/
00055 GeoRectangularProjection::~GeoRectangularProjection() 
{
}

00059 void GeoRectangularProjection::print(ostream& out) const
{
      out << "GeoRectangularProjection[";
      GeoRectangularProjectionAttributes::print(out);
      out << "]"; 
} 

00066 PaperPoint GeoRectangularProjection::operator()(const UserPoint& point)  const
{
      Log::dev() << "GeoRectangularProjection::operator()(...) needs implementing." << endl;
      return Transformation::operator()(point);
}

00072 PaperPoint GeoRectangularProjection::operator()(const GeoPoint& point)  const
{
      if ( !projection_ ) 
            return PaperPoint(point.x(), point.y(), point.value(), point.missing(), point.border());
      
      TeCoord2D geo = TeCoord2D(point.longitude()*TeCDR, point.latitude()*TeCDR);
      TeCoord2D xy = projection_->LL2PC(geo);

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

00083 PaperPoint GeoRectangularProjection::operator()(const PaperPoint& point)  const
{
      Log::dev() << "GeoRectangularProjection::operator()(...) needs implementing." << endl;
      return Transformation::operator()(point);
}

00089 void GeoRectangularProjection::revert(const PaperPoint& xy, GeoPoint& point)  const
{
      if ( !projection_ ) {
            point = GeoPoint(xy.x(), xy.y());
            return;
      }
      TeCoord2D texy = TeCoord2D(xy.x(), xy.y());
      TeCoord2D geo = projection_->PC2LL(texy);  
      point = GeoPoint(geo.x()*TeCRD, geo.y()*TeCRD); 
}

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

00106 bool GeoRectangularProjection::needShiftedCoastlines()  const
{
      return true;
}

00111 void GeoRectangularProjection::boundingBox(double& xmin, double& ymin, 
                  double& xmax, double& ymax)  const
{
      xmin = std::min(min_longitude_, max_longitude_);
      xmax = std::max(min_longitude_, max_longitude_);
      ymin = std::min(min_latitude_, max_latitude_);
      ymax = std::max(min_latitude_, max_latitude_);
      
      xmin = xmin - 10;
      xmax = xmax +10;
      ymin = (ymin < - 80 ) ? -90. : ymin -10;
      ymax  = ( ymax >  80 ) ? 90. :  ymax + 10;
}

00125 void GeoRectangularProjection::aspectRatio(double& width, double& height) 
{
      init();
      Transformation::aspectRatio(width, height);
}


00132 double GeoRectangularProjection::getMinX()  const
{
       return min_longitude_;
}

00137 double GeoRectangularProjection::getMinY()  const
{
      return min_latitude_;
}

00142 double GeoRectangularProjection::getMaxX()  const
{
      return max_longitude_;
}

00147 double GeoRectangularProjection::getMaxY()  const
{
      return max_latitude_;
}

00152 void GeoRectangularProjection::setMinX(double x)  
{
      min_longitude_ = x;
}

00157 void GeoRectangularProjection::setMinY(double y)  
{
      min_latitude_ = y;
}

00162 void GeoRectangularProjection::setMaxX(double x)  
{
      max_longitude_ = x;
}

00167 void GeoRectangularProjection::setMaxY(double y)  
{
      max_latitude_ = y;
}

00172 double GeoRectangularProjection::getMinPCX()  const
{
      return xpcmin_;
}

00177 double GeoRectangularProjection::getMinPCY()  const
{
      return ypcmin_;
}

00182 double GeoRectangularProjection::getMaxPCX()  const
{
      return xpcmax_;
}

00187 double GeoRectangularProjection::getMaxPCY()  const
{
      return ypcmax_;
}

00192 void GeoRectangularProjection::gridLongitudes(const GridPlotting& grid)  const
{
      const vector<double>& longitudes = grid.longitudes();
      double min = std::max(min_latitude_, -90.);
      double max = std::min(max_latitude_, 90.);
      double step = (max - min)/20;
      for (vector<double>::const_iterator lon = longitudes.begin(); lon != longitudes.end(); ++lon)
      {
            Polyline poly;
            poly.setAntiAliasing(false);
      
            for (double lat = min; lat <= max+step; lat += step) {
                  ( lat >  max ) ?  
                              poly.push_back((*this)(GeoPoint(*lon, max))) :
                              poly.push_back((*this)(GeoPoint(*lon,lat)));
            }
            grid.add(poly);
      }
}

00212 void GeoRectangularProjection::gridLatitudes(const GridPlotting& grid)  const
{
      const vector<double>& latitudes = grid.latitudes();
      
      double step = (max_longitude_ - min_longitude_)/20;
      for(vector<double>::const_iterator lat = latitudes.begin(); lat != latitudes.end(); ++lat)
      {
            if ( *lat < -90 ) continue;
            if ( *lat > 90 ) continue;
            Polyline poly;        
            poly.setAntiAliasing(false);
            for (double lon = getMinX(); lon <= getMaxX() + step; lon += step) {
                  poly.push_back((*this)(GeoPoint(lon,*lat)));
            }
            grid.add(poly);
      }
}

/*!
 \brief generates text to mark longitudes at the top
 
 \sa Text
*/
00235 void GeoRectangularProjection::labels(const LabelPlotting& label, TopAxisVisitor&)  const
{
      Text *text;
      double y = min_latitude_ + (max_latitude_-min_latitude_)*.2;
      const vector<double>& longitudes = label.longitudes();
      for(vector<double>::const_iterator lon = longitudes.begin(); lon != longitudes.end(); ++lon)
      {     
            if ( *lon > min_longitude_ &&  *lon < max_longitude_ ) {
                  GeoPoint point(*lon, y);
                  text = new Text();      
                  text->setText(point.writeLongitude());    
                  text->setJustification(MCENTRE);
                  text->setVerticalAlign(MBOTTOM);
                  text->push_back((*this)(point));
                  label.add(text);
            }
      }           
}

/*!
 \brief generates text to mark longitudes at the bottom

 \sa Text
*/
00259 void GeoRectangularProjection::labels(const LabelPlotting& label, BottomAxisVisitor&)  const
{
      Text *text;
      double y = min_latitude_ + (max_latitude_-min_latitude_)*.8;
      const vector<double>& longitudes = label.longitudes();
      for(vector<double>::const_iterator lon = longitudes.begin(); lon != longitudes.end(); ++lon)
      {     
            if ( *lon > min_longitude_ &&  *lon < max_longitude_ ) {
                  GeoPoint point(*lon, y);
                  text = new Text();      
                  text->setText(point.writeLongitude());    
                  text->setJustification(MCENTRE);
                  text->setVerticalAlign(MTOP);
                  text->push_back((*this)(point));
                  label.add(text);
            }
      }           
}

/*!
 \brief generates text to mark latitudes at the left
  
 \sa Text
*/
00283 void GeoRectangularProjection::labels(const LabelPlotting& label, LeftAxisVisitor&)  const
{
      Text *text;
      const vector<double>& latitudes = label.latitudes();
      for(vector<double>::const_iterator lat = latitudes.begin(); lat != latitudes.end(); ++lat)
      {     
            if ( *lat > min_latitude_ &&  *lat < max_latitude_ ) {
                  const double lon = max_longitude_ - ((max_longitude_-min_longitude_)*.1);
                  GeoPoint point(lon, *lat);
                  text = new Text();      
                  text->setText(point.writeLatitude());     
                  text->setJustification(MRIGHT);
                  text->setVerticalAlign(MHALF);
                  text->push_back((*this)(point));
                  label.add(text);
            }
      }           
}
00301 void GeoRectangularProjection::labels(const LabelPlotting& label, DrawingVisitor&)  const
{
}

/*!
 \brief generates text to mark latitudes at the right
 
 \sa Text
*/
00310 void GeoRectangularProjection::labels(const LabelPlotting& label, RightAxisVisitor&)  const
{
      Text *text;
      const vector<double>& latitudes = label.latitudes();
      for(vector<double>::const_iterator lat = latitudes.begin(); lat != latitudes.end(); ++lat)
      {     
            if ( *lat > min_latitude_ &&  *lat < max_latitude_ ) {
                  const double lon = min_longitude_ + ((max_longitude_-min_longitude_)*.1);
                  GeoPoint point(lon, *lat);
                  text = new Text();      
                  text->setText(point.writeLatitude());     
                  text->setJustification(MLEFT);
                  text->setVerticalAlign(MHALF); 
                  text->push_back((*this)(point));
                  label.add(text);
            }
      }           
}


void GeoRectangularProjection::init()
{
      // make sure min < max! 
      if (min_longitude_ > max_longitude_) {
          min_longitude_-=360;
      }
      if (min_latitude_ > max_latitude_) {
            Log::warning() << "lower_left_latitude > upper_left_latitude---> swap" << endl;
            std::swap(min_latitude_, max_latitude_);
      }
      
      //if ( min_latitude_ < -90 ) min_latitude_ = -90;
      //if ( max_latitude_ > 90 ) max_latitude_ = 90;
      if ( min_longitude_ == max_longitude_) 
            max_longitude_ += 5;
      if ( min_latitude_ == max_latitude_) 
            max_latitude_ += 5;
      
      xpcmin_ = min_longitude_;
      ypcmin_ = min_latitude_;
      xpcmax_ = max_longitude_;
      ypcmax_ = max_latitude_;      
}

MercatorProjection::MercatorProjection()
{
}

MercatorProjection::~MercatorProjection()
{
      delete projection_;
}

void MercatorProjection::print(ostream& out) const
{
      out << " mercator[";
      GeoRectangularProjection::print(out);
      out << "]";
} 

void MercatorProjection::init()
{
      if (!projection_) 
            projection_ = new TeMercator(TeDatum(), 0);
      // make sure min < max! 
      if (min_longitude_ > max_longitude_) {
            Log::warning() << "lower_left_longitude > upper_left_longitude---> swap" << endl;
            std::swap(min_longitude_, max_longitude_);
      }
      if (min_latitude_ > max_latitude_) {
            Log::warning() << "lower_left_latitude > upper_left_latitude---> swap" << endl;
            std::swap(min_latitude_, max_latitude_);
      }
      
      min_latitude_ = std::max(min_latitude_, -85.);
      max_latitude_ = std::min(max_latitude_, 85.);

      
      GeoPoint   ll(min_longitude_, min_latitude_);
      GeoPoint   ur(max_longitude_, max_latitude_);

      PaperPoint xy;    

      xy = (*this)(ll); 
      xpcmin_ = xy.x();
      ypcmin_ = xy.y();

      xy = (*this)(ur);
      xpcmax_ = xy.x();
      ypcmax_ = xy.y();             
} 

double round( double x) {
        double prec = 100;
        return floor(x*prec+0.5)/prec;
    }

void GeoRectangularProjection::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();
      
         init();
}
/*
void GeoRectangularProjection::operator()(const vector<GeoPoint>& from, vector<PaperPoint>& to) const
{
      to.reserve(from.size());
      for (vector<GeoPoint>::const_iterator point = from.begin(); point != from.end(); ++point)
      {
            if ( in(*point) )
                  to.push_back(PaperPoint((*point).x(), (*point).y()));

            if ( min_longitude_ < -180 ) {
                  GeoPoint shift = point->left();           
                  if ( in(shift) )
                        to.push_back(PaperPoint(shift.x(), shift.y()));
            }
            if ( max_longitude_ > 360 ) {
                  GeoPoint shift = point->right();
                  if ( in(shift) )
                        to.push_back(PaperPoint(shift.x(), shift.y()));
            }
            if ( max_longitude_ > 720 ) {
                  GeoPoint shift = point->right().right();
                  if ( in(shift) )
                              to.push_back(PaperPoint(shift.x(), shift.y()));
            }
      }
}
*/

Generated by  Doxygen 1.6.0   Back to index