ProcessGroundPolygons.cpp

#include "ProcessGroundPolygons.h"

#include <vector>

#include "geos/geom/CoordinateSequence.h"
#include "geos/geom/CoordinateArraySequence.h"
#include "geos/geom/LineString.h"
#include "geos/geosAlgorithm.h"

#include "Application.h"
#include "BoxcarCachingAlgorithm.h"
#include "IException.h"
#include "PolygonTools.h"
#include "Projection.h"

using namespace std;
namespace Isis {

  ProcessGroundPolygons::ProcessGroundPolygons() {
    p_groundMap = NULL;
  }


  void ProcessGroundPolygons::Rasterize(std::vector<double> &lat,
                                        std::vector<double> &lon,
                                        std::vector<double> &values) {

    // Decide if we need to split the poly on the 360 boundry
    // Yes, we can do this better. The lat and lon vectors should be passed in as polygons, but
    // for now this is reusing older code.
    bool crosses = false;
    for (unsigned int i = 0; i < lon.size() - 1; i++) {
      if (fabs(lon[i] - lon[i+1]) > 180.0) {
        crosses = true;
        break;
      }
    }

    if (crosses) {
      // Make a polygon from the lat/lon vectors and split it on 360
      geos::geom::CoordinateSequence *pts = new geos::geom::CoordinateArraySequence();
      for (unsigned int i = 0; i < lat.size(); i++) {
        pts->add(geos::geom::Coordinate(lon[i], lat[i]));
      }
      pts->add(geos::geom::Coordinate(lon[0], lat[0]));

      geos::geom::Polygon *crossingPoly = Isis::globalFactory.createPolygon(
          globalFactory.createLinearRing(pts), NULL);

      geos::geom::MultiPolygon *splitPoly = NULL;
      try {
        splitPoly = PolygonTools::SplitPolygonOn360(crossingPoly);
      }
      // Ignore any pixel footprints that could not be split. This should only be pixels
      // that contain the pole.
      catch (IException &) {
        // See leading comment
      }

      delete crossingPoly;

      if (splitPoly != NULL) {
        // Process the polygons in the split multipolygon as if we were still using the lat/lon vectors
        for (unsigned int g = 0; g < splitPoly->getNumGeometries(); ++g) {
          const geos::geom::Polygon *poly = 
              dynamic_cast<const geos::geom::Polygon *>(splitPoly->getGeometryN(g));

          geos::geom::CoordinateSequence *llcoords = poly->getExteriorRing()->getCoordinates();

          // Move each coordinate in the exterior ring of this lat/lon polygon to vectors
          // Ignore any holes in the polygon
          std::vector<double> tlat;
          std::vector<double> tlon;
          for (unsigned int cord = 0; cord < llcoords->getSize() - 1; ++cord) {
            tlon.push_back(llcoords->getAt(cord).x);
            tlat.push_back(llcoords->getAt(cord).y);
          }

          Convert(tlat, tlon);
          ProcessPolygons::Rasterize(p_samples, p_lines, values);
        }
        delete splitPoly;
      }
    }
    else {
      Convert(lat, lon);
      ProcessPolygons::Rasterize(p_samples, p_lines, values);
    }
  }


  void ProcessGroundPolygons::Rasterize(std::vector<double> &lat,
                                        std::vector<double> &lon,
                                        int &band, double &value) {

    Convert(lat, lon);
    ProcessPolygons::Rasterize(p_samples, p_lines, band, value);
  }


  void ProcessGroundPolygons::Convert(std::vector<double> &lat,
                                      std::vector<double> &lon) {

    p_samples.clear();
    p_lines.clear();

    for (unsigned int i = 0; i < lat.size(); i++) {
      if (p_groundMap->SetUniversalGround(lat[i], lon[i])) {
        p_samples.push_back(p_groundMap->Sample());
        p_lines.push_back(p_groundMap->Line());
      }
    }
  }


  void ProcessGroundPolygons::EndProcess() {

    if(p_groundMap != NULL) {
      delete p_groundMap;
    }

    ProcessPolygons::EndProcess();
  }


  void ProcessGroundPolygons::Finalize() {

    if(p_groundMap != NULL) {
      delete p_groundMap;
    }

    ProcessPolygons::Finalize();
  }


  void ProcessGroundPolygons::AppendOutputCube(QString &cubeStr,
      const QString &avgFileName,
      const QString &countFileName) {
    /*We need a ground map for converting lat/long to line/sample  see Convert()*/
    Cube cube(cubeStr, "r");
    p_groundMap = new UniversalGroundMap(cube);
    ProcessPolygons::AppendOutputCube(avgFileName, countFileName);

  }


  void ProcessGroundPolygons::SetStatCubes(const QString &avgFileName,
      const QString &countFileName,
      Isis::CubeAttributeOutput &outAtts,
      QString &cubeStr) {
    /*We need a ground map for converting lat/long to line/sample  see Convert()*/
    Cube cube(cubeStr);
    p_groundMap = new UniversalGroundMap(cube);

    /*setup input cube to transfer projection or camera labels*/
    CubeAttributeInput inAtts;
    Isis::Process::SetInputCube(cubeStr, inAtts, 0);
    int nBands = this->InputCubes[0]->bandCount();
    int nLines = this->InputCubes[0]->lineCount();
    int nSamples = this->InputCubes[0]->sampleCount();

    this->Process::SetOutputCube(avgFileName, outAtts, nSamples, nLines, nBands);
    this->Process::SetOutputCube(countFileName, outAtts, nSamples, nLines, nBands);

    OutputCubes[0]->addCachingAlgorithm(new BoxcarCachingAlgorithm());
    OutputCubes[1]->addCachingAlgorithm(new BoxcarCachingAlgorithm());

    ClearInputCubes();

  }


  void ProcessGroundPolygons::SetStatCubes(const QString &parameter,
      QString &cube) {

    QString avgString =
      Application::GetUserInterface().GetFileName(parameter);
    CubeAttributeOutput atts =
      Application::GetUserInterface().GetOutputAttribute(parameter);

    FileName file(avgString);
    QString path = file.path();
    QString filename = file.baseName();
    QString countString = path + "/" + filename + "-count-";

    SetStatCubes(avgString, countString, atts, cube);

  }


  void ProcessGroundPolygons::SetStatCubes(const QString &parameter,
      Isis::Pvl &map, int bands) {

    QString avgString =
      Application::GetUserInterface().GetFileName(parameter);
    CubeAttributeOutput atts =
      Application::GetUserInterface().GetOutputAttribute(parameter);

    FileName file(avgString);
    QString path = file.path();
    QString filename = file.baseName();
    QString countString = path + "/" + filename + "-count-";

    SetStatCubes(avgString, countString, atts, map, bands);

  }


  void ProcessGroundPolygons::SetStatCubes(const QString &avgFileName,
      const QString &countFileName,
      Isis::CubeAttributeOutput &atts,
      Isis::Pvl &map, int bands) {
    int samples, lines;

    Projection *proj = ProjectionFactory::CreateForCube(map, samples, lines,
                       false);

    this->ProcessPolygons::SetStatCubes(avgFileName, countFileName, atts,
                                         samples, lines, bands);

    OutputCubes[0]->addCachingAlgorithm(new BoxcarCachingAlgorithm());
    OutputCubes[1]->addCachingAlgorithm(new BoxcarCachingAlgorithm());

    /*Write the pvl group to the cube files.*/

    PvlGroup group = map.findGroup("Mapping", Pvl::Traverse);

    OutputCubes[0]->putGroup(group);
    OutputCubes[1]->putGroup(group);

    // If there is an alpha cube in the label passed, attach to output cubes
    if (map.hasGroup("AlphaCube")) {
      PvlGroup alpha = map.findGroup("AlphaCube", Pvl::Traverse);
      OutputCubes[0]->putGroup(alpha);
      OutputCubes[1]->putGroup(alpha);
    }

    /*We need a ground map for converting lat/long to line/sample  see Convert()*/
    p_groundMap = new UniversalGroundMap(*OutputCubes[0]);

    delete proj;
  }

} /* end namespace isis*/