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

CoastPlotting.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 CoastPlotting.cc
    \brief Implementation of the Template class CoastPlotting.
    
    Magics Team - ECMWF 2004
    
    Started: Mon 2-Feb-2004
    
    Changes:
    
*/

#include "Netcdf.h"
#include "CoastPlotting.h"
#include "Timer.h"
#include "SceneVisitor.h"
#include "PolyCoast.h"
#include "ViewNode.h"

using namespace magics;

VectorOfPointers<vector<PolyCoast* >  > CoastPlotting::shiftedGlobal_;
VectorOfPointers<vector<PolyCoast* >  > CoastPlotting::global_;
string CoastPlotting::resolution_;

CoastPlotting::CoastPlotting() : layer_(0)
{}


void CoastPlotting::operator()(PreviewVisitor& parent)
{
      
      const Transformation& transformation = parent.transformation();
      
      CoastPlotting& preview = parent.coastlines();

      LandShading* land = new LandShading();
      SeaShading* sea = new SeaShading();
      land->setColour(new Colour("cream"));
      sea->setColour(new Colour("white"));
      land->setSea(sea);
      preview.setLand(land);
      preview.setColour(new Colour("tan"));
      preview.decode(parent);

      for (vector<PolyCoast* >::iterator poly = preview.coastlines_.begin(); poly != preview.coastlines_.end(); ++poly)
      {
            PolyCoast* coast = static_cast<PolyCoast*>(*poly);
            
            coast->setThickness(thickness_);
            coast->setColour(*colour_);
            coast->setLineStyle(style_);  

            bool in = false;

            vector<GeoPoint>& points = coast->coastlines();
            vector<GeoPoint>::const_iterator point = points.begin();

            while ( in == false &&  point != points.end() )
            {
                  in = transformation.in(*point); 
                  ++point;                                  
            }

            if ( in )
            {
                  (*preview.land_)(coast);                        
                  transformation(*coast, parent);                             
            }
            
      }
      //Now we had the frame! 
      Polyline* frame = new Polyline();
      frame->setAntiAliasing(false);
      frame->setThickness(thickness_);
      frame->setColour(Colour("tan"));
      frame->setLineStyle(style_);  
      frame->push_back(PaperPoint(transformation.getMinX(), transformation.getMinY()));
      frame->push_back(PaperPoint(transformation.getMaxX(), transformation.getMinY()));
      frame->push_back(PaperPoint(transformation.getMaxX(), transformation.getMaxY()));
      frame->push_back(PaperPoint(transformation.getMinX(), transformation.getMaxY()));
      frame->push_back(PaperPoint(transformation.getMinX(), transformation.getMinY()));

      parent.push_back(frame);
}


struct SortHelper
{
      bool operator() (PolyCoast*first, PolyCoast* second) {
            return first->coastlines().size() > second->coastlines().size();
      }
};

void CoastPlotting::operator()(DrawingVisitor& parent)
{
      Timer timer("coastlines", "prepare the coastlines");
      if ( !layer_) layer_ = &parent;
      const Transformation& transformation = parent.transformation();
      decode(parent);
      land_->setBoundingBox(transformation);
      
      Polyline* sea = new Polyline();
      sea->setThickness(0);
      sea->setColour(Colour("none"));
      sea->setLineStyle(style_);    
      
      land_->sea(sea, *layer_);

      std::sort(coastlines_.begin(), coastlines_.end(), SortHelper());

      for (vector<PolyCoast* >::iterator poly = coastlines_.begin(); poly != coastlines_.end(); ++poly)
      {
            PolyCoast* coast = static_cast<PolyCoast*>(*poly);
      
            coast->setThickness(thickness_);
            coast->setColour(*colour_);
            coast->setLineStyle(style_);  
            
            bool in = true;
/*
                  vector<GeoPoint>& points = coast->coastlines();
                              vector<GeoPoint>::const_iterator point = points.begin();

                              while ( in == false &&  point != points.end() )
                              {
                                    in = transformation.in(*point); 
                                    ++point;                                  
                              }

                  if ( in )
                  {
                        (*land_)(coast);                    
                        transformation(*coast, parent);                             
                  }
                  else {
                        delete  coast;
                  }
*/
            
            
            //Calculate Bounding of the caostline...
            vector<GeoPoint>& points = coast->coastlines();
            vector<GeoPoint>::const_iterator point = points.begin();
            double xmin = point->x();
            double ymin = point->y();
            double xmax = point->x();
            double ymax = point->y();
            
            while ( point != points.end() ) {
                        xmin = std::min(xmin, point->x());
                        xmax = std::max(xmax, point->x());
                        ymin = std::min(ymin, point->y());
                        ymax = std::max(ymax, point->y());
                        ++point;
            }
      
      
            if ( transformation.getMaxX() < xmin ) in = false;
            if ( transformation.getMaxY() < ymin ) in = false;
            if ( transformation.getMinX() > xmax ) in = false;
            if ( transformation.getMinY() >  ymax ) in = false;
      
            if ( in )
            {
                  (*land_)(coast);                    
                  transformation(*coast, *layer_);          
            }
      }
      // Steps to do ... recuperer tous les polycoast et les trier suivant leur taille ( et leurs levels???)
      
      (*boundaries_)(*layer_);
      /*
            Polyline* frame = new Polyline();
            frame->setAntiAliasing(false);
            frame->setThickness(thickness_);
            frame->setColour(*colour_);
            frame->setLineStyle(style_);  
            frame->push_back(PaperPoint(transformation.getMinX(), transformation.getMinY()));
            frame->push_back(PaperPoint(transformation.getMaxX(), transformation.getMinY()));
            frame->push_back(PaperPoint(transformation.getMaxX(), transformation.getMaxY()));
            frame->push_back(PaperPoint(transformation.getMinX(), transformation.getMaxY()));
            frame->push_back(PaperPoint(transformation.getMinX(), transformation.getMinY()));
            
            parent.push_back(frame);
      */
}

CoastPlotting::~CoastPlotting() 
{}

/*!
 Class information are given to the output-stream.
*/
void CoastPlotting::print(ostream& out)  const
{
      out << "CoastPlotting[";
      CoastPlottingAttributes::print(out);
      out << "] ";
}

void CoastPlotting::decode(const Layout& parent )
{
      prepareGlobal(parent);
      
      const Transformation& transformation = parent.transformation();
      
      coastlines_.clear();
      coastlines_.reserve(global_.size());

      int i = 0;

      for (vector<PolyCoast* >::iterator coast = global_.begin(); coast != global_.end(); ++coast)
      {
            //if ( transformation.out(**coast) ) continue;
            //if (i < 1 )
            coastlines_.push_back((*coast)->clone());
            i++;
      }

      if ( !transformation.needShiftedCoastlines() ) return;

      for (vector<PolyCoast* >::iterator coast = shiftedGlobal_.begin(); coast != shiftedGlobal_.end(); ++coast)
      {
            //if ( transformation.out(**coast) ) continue;
            
             coastlines_.push_back((*coast)->clone());
            i++;
      }
}


/*!
      ShiftPolyCoastCoordinates was written as a relpacement for using std::transform to
      shift the coastline polygon coordinates. This makes the class ShiftHelper redundant.
      The problem with using std::transform was that it did not call the puch_back()
      function and therefore the polygon's bounding box was not updated (minX, maxX, minY, maxY).
      Without the bounding box, clipping is much harder. Note also that since push_back() is
      now being used, the new polygon does not have to be resized before the shifting operation.
*/
00258 void ShiftPolyCoastCoordinates (PolyCoast* src_poly, PolyCoast* dest_poly, double shift)
{
    double previous = 0;
    vector<GeoPoint>& from = src_poly->coastlines();
    vector<GeoPoint>& to = dest_poly->coastlines();

    for (vector<GeoPoint>::iterator point = from.begin(); point != from.end(); ++point)
    {
        if ( shift == 0 && previous < 175 )
        {
            previous = point->x() - 360;
            if ( point->x() > 180 )
            {
                  to.push_back (GeoPoint(point->x()-360, point->y()));
                  continue;
            }
      }
      previous = point->x();
      to.push_back (GeoPoint(point->x()+shift, point->y()));
    }
}

/*!
  general logic in determining which coastline file to use:

  - first, map_coastline_path can override anything else
  - but, if the path is invalid then we pretend the user did not set it
  - then, map_coastline_resolution is checked - if set to one of the predefined
    strings such as 'low' and 'high', then use the corresponding coast file.
  - but if it is set to something silly, then pretend it was set to 'automatic'.
  - last, if there is no user-specified file, and the resolution is set to
    'automatic', then the coastline file is automatically deduced.
*/
void CoastPlotting::prepareGlobal(const Layout& owner)
{
    
    const Transformation& transformation = owner.transformation();
    string netcdf;
    string strSpecifiedResolution = lowerCase (getResolution());
    string strSpecifiedCoastFile  =  getFile();
    bool bPathSupplied = true;

    // check - has the user tried to specify their own file to read?
    if ( magCompare(strSpecifiedCoastFile, "off") ) 
      strSpecifiedCoastFile = ""; // Metview default is strange!
    if (strSpecifiedCoastFile != ""  )
    {
        // yes - but can we read it?
        if (fileReadable (strSpecifiedCoastFile))
        {
            netcdf = strSpecifiedCoastFile;  // yes - let's use it
        }
        else
        {
            // no - if it's just a filename with no path, then try to find the
            // file in the MAGICS shared directory
            bool bReadable = false;
            if (strSpecifiedCoastFile.find_first_of("/") == string::npos)
            {
                netcdf = getEnvVariable("MAGPLUS_HOME") + MAGPLUS_PATH_TO_SHARE_ + strSpecifiedCoastFile;
                bReadable = fileReadable (netcdf);
            }
            if (bReadable)
            {
                Log::debug() << "User-specified coastline file found in " << netcdf << "\n";
            }
            else
            {
                // nope, the user supplied a path to a file which we cannot read
                // - issue an error message and revert to the default coastline
                Log::error() << "Coastline file '" << strSpecifiedCoastFile
                             << "' cannot be found/read. Using default coastline file." << "\n";
                strSpecifiedCoastFile = "";
            }
        }
    }

    // check - are we NOT using a a user-defined coastline file?
    if (strSpecifiedCoastFile == "")
    {
        // use the resolution defined by the user?
        if (strSpecifiedResolution != "automatic")
        {
            if      (strSpecifiedResolution == "tiny")        netcdf = "coast_tiny.nc";
            else if (strSpecifiedResolution == "very_low")    netcdf = "coast_very_low.nc";
            else if (strSpecifiedResolution == "low")         netcdf = "coast_low.nc";
            else if (strSpecifiedResolution == "medium")      netcdf = "coast_medium.nc";
            else if (strSpecifiedResolution == "medium_high") netcdf = "coast_medium_high.nc";
            else if (strSpecifiedResolution == "high")        netcdf = "coast_high.nc";
            else 
            {
                Log::error() << "Coastline resolution '" << strSpecifiedResolution
                             << "' is not valid. Use 'tiny', 'very_low', 'low', 'medium', 'medium_high' or 'high'. "
                             << "Using default coastline file." << "\n";

                strSpecifiedResolution = "automatic";
            }
            bPathSupplied = false;
        }

        // if resolution set to automatic, then find the right file ourselves

        if (strSpecifiedResolution == "automatic")
        {
            // get the bounds of the plot in geographical coordinates
            double fMinX = transformation.getMinX();
            double fMaxX = transformation.getMaxX();
            double fMinY = transformation.getMinY();
            double fMaxY = transformation.getMaxY();

            // get the geographical and paper dimensions
            double fGeoAreaWidth  = fMaxX - fMinX;
            double fGeoAreaHeight = fMaxY - fMinY;

            double fPaperAreaWidth  = owner.absoluteWidth();
            double fPaperAreaHeight = owner.absoluteHeight();

            // work out the ratios of geographical to paper lengths
            double fGeoOverPaperWidth  = fGeoAreaWidth  / fPaperAreaWidth;
            double fGeoOverPaperHeight = fGeoAreaHeight / fPaperAreaHeight;

            // choose the smallest (smaller ratio means more detail required)
            double fRatio = min (fGeoOverPaperWidth, fGeoOverPaperHeight);


            // calculate the lowest resolution that will adequately plot the coastline
            //   - this logic was derived from trial and error, and subjective judgement
                 if (fRatio > 25)  netcdf = "coast_very_low.nc";      // lowest-res
            else if (fRatio >  6)  netcdf = "coast_low.nc";           // low-res
            else if (fRatio >  3)  netcdf = "coast_medium.nc";        // medium-res
            else if (fRatio >  1)  netcdf = "coast_medium_high.nc";   // medium-high-res
            else                   netcdf = "coast_high.nc";          // highest-res

            bPathSupplied = false;

            Log::debug() << "Automatic coastlines selection: " << netcdf << "\n";
        }
    }


    // we do not need to read the coastline file again if it is already in memory
    if ( netcdf == resolution_) return;
    resolution_ = netcdf;

    string filename = (bPathSupplied)
                         ? resolution_
                         : getEnvVariable("MAGPLUS_HOME") + MAGPLUS_PATH_TO_SHARE_ + resolution_;

    //filename = getEnvVariable("MAGPLUS_HOME") + MAGPLUS_PATH_TO_SHARE_ + file_;

    try {
            global_.clear();
            shiftedGlobal_.clear();
            vector<PolyCoast* > tmp;
            {
                  Netcdf netcdf(filename);
                  vector<double> first;
                  vector<double> nb;
                  vector<double> lat;
                  vector<double> lon;
                  vector<double> level;
                  vector<double> area;
                  vector<double> greenwich;
                  netcdf.get("position", first);
                  netcdf.get("size", nb);
                  netcdf.get("longitude", lon);
                  netcdf.get("latitude", lat);
                  netcdf.get("level", level);
                  netcdf.get("area", area);
                  netcdf.get("greenwich", greenwich);

                  global_.reserve(first.size()*2);
                  
                  for (unsigned int p1 = 0; p1 < first.size(); p1++ )
                  {
                        int i = (int)first[p1];
                  
                        PolyCoast* poly = new PolyCoast();
                        poly->setThickness(thickness_);
                        poly->setColour(*colour_);
                        poly->setLineStyle(style_);
                        poly->setFilled(false);
                        poly->level(int(level[p1]));
                        poly->greenwich(int(greenwich[p1]));

                        i = (int)first[p1];
                        vector<GeoPoint>& coast = poly->coastlines(); 
                        vector<double> lons;
                        vector<double> lats;
                        
                                                
                        for (int p2 = 0; p2 < (int)nb[p1]; p2++ )
                        {
                              coast.push_back(GeoPoint(lon[i], lat[i]));
                              lons.push_back(lon[i]);
                              lats.push_back(lat[i]);
                              
                              i++;
                        }

                        global_.push_back(poly);
                        // Now we study the polyline! What is the bouding box?

                        double xmin = *std::min_element(lons.begin(), lons.end());
                        double xmax =*std::max_element(lons.begin(), lons.end());
                        double ymin = *std::min_element(lats.begin(), lats.end());
                        double ymax =*std::max_element(lats.begin(), lats.end());

                        
                        PolyCoast* shift;
                        PolyCoast* shift2;

                        if (xmin < 0 && xmax > 180) {
                              //poly->push_back(poly->front());
                              Log::dev() << "case 1--> on triple le poly+360/-360" << endl;
                              shift = new PolyCoast();
                              shift->level(int(level[p1]));
                              ShiftPolyCoastCoordinates (poly, shift, 360);
                              shiftedGlobal_.push_back(shift);
                              
                              shift = new PolyCoast();
                              shift->level(int(level[p1]));
                              ShiftPolyCoastCoordinates (poly, shift, 720);
                              shiftedGlobal_.push_back(shift);

                              shift = new PolyCoast();
                              ShiftPolyCoastCoordinates (poly, shift, -360);
                              shift->level(int(poly->level()));
                              shiftedGlobal_.push_back(shift);
                              
                        }
                        else  if (xmax-xmin > 350)  {
                              Log::dev() << "Bounding box = [" << xmin << "/" << ymin << "/" << xmax << "/" << ymax << "]\n";       
                              Log::dev() << "greenwich box = [" << greenwich[p1] << "]\n";      

                              Log::dev() << "WRAP AROUND to be done for " << p1 << "[" << area[p1] << "]" << endl;
                              if ( greenwich[p1] ) {
                                    shift2 = new PolyCoast();
                                    shift2->level(int(level[p1]));
                                    ShiftPolyCoastCoordinates (poly, shift2, 0);
                                    global_.pop_back();
                                    global_.push_back(shift2);

                                    // on shift et on triple!
                                    shift = new PolyCoast();
                                    shift->level(int(level[p1]));
                                    ShiftPolyCoastCoordinates (shift2, shift, 360);
                                    shiftedGlobal_.push_back(shift);
                                    // on shift et on triple!
                                    
                                    shift = new PolyCoast();
                                    shift->level(int(level[p1]));
                                    ShiftPolyCoastCoordinates (shift2, shift, 720);
                                    shiftedGlobal_.push_back(shift);

                                    shift = new PolyCoast();
                                    shift->level(int(level[p1]));
                                    ShiftPolyCoastCoordinates (shift2, shift, -360);
                                    shiftedGlobal_.push_back(shift);

                                    // delete poly;  // poly was defined, but not added to the global list
                              }
                              else {
                                    // On fait une grande ligne de -360 a 360
                                    // First find the gap ... then add 

                                    vector<GeoPoint>::iterator previous = coast.begin();
                                    
                                    // try to reorganise thenit goest from one side to another!
                                    for (vector<GeoPoint>::iterator point = coast.begin(); point != coast.end(); ++point) {                       
                                          if (abs(previous->x() - point->x()) > 100) {                            
                                                std::rotate(coast.begin(), point, coast.end());
                                          }
                                          
                                          previous = point;
                                    }
                                     // in our coastlines file we are always going from east to west!
                                    // Then we reverse to be sure that we rae going from west to east!
                                    //.. but this has perhaps to be checked! 
                                    std::reverse(coast.begin(), coast.end()); 
                                 
                                    PolyCoast *keep = new PolyCoast();
                                    vector<GeoPoint>& kc = keep->coastlines();
                                    for ( vector<GeoPoint>::iterator point = coast.begin(); point != coast.end(); ++point) {                      
                                          kc.push_back(*point);
                                    }
                                    
                                    
                                    shift = new PolyCoast();
                                    shift->level(int(level[p1]));

                                    ShiftPolyCoastCoordinates (keep, shift, -360);
                                    
                                    coast.insert(coast.begin(), shift->coastlines().begin(), shift->coastlines().end());
                                    
                                    if ( shift->getMinX() < poly->getMinX() ) poly->setMinX(shift->getMinX());
                                    if ( shift->getMaxX() > poly->getMaxX() ) poly->setMaxX(shift->getMaxX());
                                    if ( shift->getMinY() < poly->getMinY() ) poly->setMinY(shift->getMinY());
                                    if ( shift->getMaxY() > poly->getMaxY() ) poly->setMaxY(shift->getMaxY());

                                    delete shift;
                                    
                                    
                                    
                                    shift = new PolyCoast();
                                    shift->level(int(level[p1]));
                              
                                    ShiftPolyCoastCoordinates (keep, shift, 360);
                                    coast.insert(coast.end(), shift->coastlines().begin(), shift->coastlines().end());
                              
                                    poly->level(0);
                                    
                                    // we've used insert(), but that does not update the min/max values
                                    
                                    if (shift->getMinX() < poly->getMinX()) poly->setMinX(shift->getMinX());
                                    if (shift->getMaxX() > poly->getMaxX()) poly->setMaxX(shift->getMaxX());
                                    if (shift->getMinY() < poly->getMinY()) poly->setMinY(shift->getMinY());
                                    if (shift->getMaxY() > poly->getMaxY()) poly->setMaxY(shift->getMaxY());

                                    delete shift;
                                    delete keep;
                              }      
                        }
                        else if ( xmin < 0 && xmax < 180) {
                                    //poly->push_back(poly->front());
                              Log::dev() << "case 4/2 --> on double le poly+360" << endl;
                              shift = new PolyCoast();
                              shift->level(int(level[p1]));
                              ShiftPolyCoastCoordinates (poly, shift, 360);
                              shiftedGlobal_.push_back(shift);
                              shift = new PolyCoast();
                              shift->level(int(level[p1]));
                              ShiftPolyCoastCoordinates (poly, shift, 720);
                              shiftedGlobal_.push_back(shift);
                        }
                        else if ( xmin > 0  && xmax < 180 ) {
                              //poly->push_back(poly->front());
                               //  Log::dev() << "Bounding box = [" << xmin << "/" << ymin << "/" << xmax << "/" << ymax << "]\n";  
                             //Log::dev() << "greenwich box = [" << greenwich[p1] << "]\n";     
                              //Log::dev() << "case 3 --> on fait rien le poly" << endl;
                              shift = new PolyCoast();
                                                            shift->level(int(level[p1]));
                                                            ShiftPolyCoastCoordinates (poly, shift, 360);
                                                            shiftedGlobal_.push_back(shift);
                                                            shift = new PolyCoast();
                                                                                                                        shift->level(int(level[p1]));
                                                                                                                        ShiftPolyCoastCoordinates (poly, shift, 720);
                                                                                                                        shiftedGlobal_.push_back(shift);
                        }
                        else if ( xmin > 0 && xmax > 180) {
                              //poly->push_back(poly->front());
                              //Log::dev() << "case 5/6 --> on double le poly-360" << endl;
                              shift = new PolyCoast();;
                              shift->level(int(level[p1]));
                              ShiftPolyCoastCoordinates (poly, shift, -360);
                              global_.push_back(shift);
                              shift = new PolyCoast();
                                                            shift->level(int(level[p1]));
                                                            ShiftPolyCoastCoordinates (poly, shift, 360);
                                                            shiftedGlobal_.push_back(shift);
                        }
//                      else { 
//                            Log::dev() << "Bounding box = [" << xmin << "/" << ymin << "/" << xmax << "/" << ymax << "]\n";       
//                            Log::dev() << "greenwich box = [" << greenwich[p1] << "]\n";      
//
//                            Log::dev() << "bof![" << xmin << ", " << xmax << "]" << endl; 
//                      }
                  }
            }
      }

      catch (exception& e) 
      {
            Log::error() << "Problem reading Coastlines file : " << filename << "\n";
            throw e;
      }
}

Generated by  Doxygen 1.6.0   Back to index