PolygonTools.cpp

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <sstream>
#include <iomanip>
#include <vector>
#include <cmath>

#include <QDebug>
#include "geos/geom/BinaryOp.h"
#include "geos/geom/CoordinateArraySequence.h"
#include "geos/geom/CoordinateSequence.h"
#include "geos/geom/LinearRing.h"
#include "geos/geom/Point.h"
#include "geos/geom/Polygon.h"
#include "geos/operation/distance/DistanceOp.h"
#include "geos/opOverlay.h"
#include "geos/operation/overlay/snap/GeometrySnapper.h"

#include "SpecialPixel.h"
#include "PolygonTools.h"
#include "TProjection.h"
#include "ProjectionFactory.h"
#include "UniversalGroundMap.h"

using namespace std;
namespace Isis {
  geos::geom::MultiPolygon *PolygonTools::LatLonToXY(
      const geos::geom::MultiPolygon &lonLatPolygon, TProjection *projection) {
    if (projection == NULL) {
      string msg = "Unable to convert Lon/Lat polygon to X/Y. ";
      msg += "No projection has was supplied";
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }

    // Convert the Lat/Lon poly coordinates to X/Y coordinates
    if(lonLatPolygon.isEmpty()) {
      return globalFactory.createMultiPolygon();
    }
    else {
      vector<geos::geom::Geometry *> *xyPolys = new vector<geos::geom::Geometry *>;
      // Convert each polygon in this multi-polygon
      for(unsigned int g = 0; g < lonLatPolygon.getNumGeometries(); ++g) {
        const geos::geom::Polygon *poly =
            dynamic_cast<const geos::geom::Polygon *>(
              lonLatPolygon.getGeometryN(g));

        // Convert each hole inside this polygon
        vector<geos::geom::Geometry *> *holes = new vector<geos::geom::Geometry *>;
        for(unsigned int h = 0; h < poly->getNumInteriorRing(); ++h) {
          geos::geom::CoordinateSequence *xycoords = new geos::geom::CoordinateArraySequence();
          geos::geom::CoordinateSequence *llcoords =
            poly->getInteriorRingN(h)->getCoordinates();

          // Convert each coordinate in this hole
          for(unsigned int cord = 0; cord < llcoords->getSize(); ++cord) {
            projection->SetGround(llcoords->getAt(cord).y,
                                  llcoords->getAt(cord).x);
            xycoords->add(geos::geom::Coordinate(projection->XCoord(),
                                                 projection->YCoord()));
          } // end num coords in hole loop

          geos::geom::LinearRing *hole = globalFactory.createLinearRing(xycoords);

          if(hole->isValid() && !hole->isEmpty()) {
            holes->push_back(hole);
          }
          else {
            delete hole;
          }
        } // end num holes in polygon loop

        // Convert the exterior ring of this polygon
        geos::geom::CoordinateSequence *xycoords = new geos::geom::CoordinateArraySequence();
        geos::geom::CoordinateSequence *llcoords =
          poly->getExteriorRing()->getCoordinates();

        // Convert each coordinate in the exterior ring of this polygon
        for (unsigned int cord = 0; cord < llcoords->getSize(); ++cord) {
          projection->SetGround(llcoords->getAt(cord).y,
                                llcoords->getAt(cord).x);
          xycoords->add(geos::geom::Coordinate(projection->XCoord(),
                                               projection->YCoord()));
        } // end exterior ring coordinate loop

        geos::geom::Polygon *newPoly = globalFactory.createPolygon(
                                         globalFactory.createLinearRing(xycoords), holes);

        if(newPoly->isValid() && !newPoly->isEmpty() && newPoly->getArea() > 1.0e-14) {
          xyPolys->push_back(newPoly);
        }
        else {
          delete newPoly;
        }
      } // end num geometry in multi-poly

      // Create a new multipoly from all the new X/Y polygon(s)
      geos::geom::MultiPolygon *spikedPoly = globalFactory.createMultiPolygon(xyPolys);

      if(spikedPoly->isValid() && !spikedPoly->isEmpty()) {
        return spikedPoly;
      }
      else {
        try {
          geos::geom::MultiPolygon *despikedPoly = Despike(spikedPoly);

          delete spikedPoly;
          spikedPoly = NULL;

          return despikedPoly;
        }
        catch(IException &e) {
          IString msg = "Unable to convert polygon from Lat/Lon to X/Y";
          throw IException(IException::Programmer, msg, _FILEINFO_);
        }
      }

    } // end else
  }


  geos::geom::MultiPolygon *PolygonTools::XYToLatLon(
    const geos::geom::MultiPolygon &xYPolygon, TProjection *projection) {

    if(projection == NULL) {
      string msg = "Unable to convert X/Y polygon to Lon/Lat. ";
      msg += "No projection was supplied";
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }

    // Convert the X/Y poly coordinates to Lat/Lon coordinates
    if(xYPolygon.isEmpty()) {
      return globalFactory.createMultiPolygon();
    }
    else {
      vector<geos::geom::Geometry *> *llPolys = new vector<geos::geom::Geometry *>;
      // Convert each polygon in this multi-polygon
      for(unsigned int g = 0; g < xYPolygon.getNumGeometries(); ++g) {
        const geos::geom::Polygon *poly =
            dynamic_cast<const geos::geom::Polygon *>(
              xYPolygon.getGeometryN(g));

        // Convert each hole inside this polygon
        vector<geos::geom::Geometry *> *holes = new vector<geos::geom::Geometry *>;
        for(unsigned int h = 0; h < poly->getNumInteriorRing(); ++h) {
          geos::geom::CoordinateSequence *llcoords = new geos::geom::CoordinateArraySequence();
          geos::geom::CoordinateSequence *xycoords =
            poly->getInteriorRingN(h)->getCoordinates();

          // Convert each coordinate in this hole
          for(unsigned int cord = 0; cord < xycoords->getSize(); ++cord) {
            projection->SetWorld(xycoords->getAt(cord).x,
                                 xycoords->getAt(cord).y);
            llcoords->add(geos::geom::Coordinate(projection->Longitude(),
                                                 projection->Latitude()));
          } // end num coords in hole loop
          holes->push_back(globalFactory.createLinearRing(llcoords));
        } // end num holes in polygon loop

        // Convert the exterior ring of this polygon
        geos::geom::CoordinateSequence *llcoords = new geos::geom::DefaultCoordinateSequence();
        geos::geom::CoordinateSequence *xycoords =
          poly->getExteriorRing()->getCoordinates();

        // Convert each coordinate in the exterior ring of this polygon
        for(unsigned int cord = 0; cord < xycoords->getSize(); ++cord) {
          projection->SetWorld(xycoords->getAt(cord).x,
                               xycoords->getAt(cord).y);
          llcoords->add(geos::geom::Coordinate(projection->Longitude(),
                                               projection->Latitude()));
        } // end exterior ring coordinate loop

        llPolys->push_back(globalFactory.createPolygon(
                             globalFactory.createLinearRing(llcoords), holes));
      } // end num geometry in multi-poly


      // Create a new multipoly from all the new Lat/Lon polygon(s)
      geos::geom::MultiPolygon *spikedPoly = globalFactory.createMultiPolygon(llPolys);

      if(spikedPoly->isValid() && !spikedPoly->isEmpty()) {
        return spikedPoly;
      }
      else {
        try {
          geos::geom::MultiPolygon *despikedPoly = Despike(spikedPoly);

          delete spikedPoly;
          spikedPoly = NULL;

          return despikedPoly;
        }
        catch(IException &e) {
          IString msg = "Unable to convert polygon from X/Y to Lat/Lon";
          throw IException(IException::Programmer, msg, _FILEINFO_);
        }
      }
    } // end else
  }


  geos::geom::MultiPolygon *PolygonTools::LatLonToSampleLine(
    const geos::geom::MultiPolygon &lonLatPolygon, UniversalGroundMap *ugm) {

    if(ugm == NULL) {
      string msg = "Unable to convert Lon/Lat polygon to Sample/Line. ";
      msg += "No UniversalGroundMap was supplied";
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }

    // Convert the Lon/Lat poly coordinates to Sample/Line coordinates
    if(lonLatPolygon.isEmpty()) {
      return globalFactory.createMultiPolygon();
    }
    else {
      vector<geos::geom::Geometry *> *slPolys = new vector<geos::geom::Geometry *>;
      // Convert each polygon in this multi-polygon
      for(unsigned int g = 0; g < lonLatPolygon.getNumGeometries(); g++) {
        const geos::geom::Polygon *poly =
            dynamic_cast<const geos::geom::Polygon *>(
              lonLatPolygon.getGeometryN(g));

        // Convert each hole inside this polygon
        vector<geos::geom::Geometry *> *holes = new vector<geos::geom::Geometry *>;
        for(unsigned int h = 0; h < poly->getNumInteriorRing(); ++h) {
          geos::geom::CoordinateSequence *slcoords = new geos::geom::DefaultCoordinateSequence();
          geos::geom::CoordinateSequence *llcoords =
            poly->getInteriorRingN(h)->getCoordinates();

          // Convert each coordinate in this hole
          for(unsigned int cord = 0; cord < llcoords->getSize(); ++cord) {
            ugm->SetUniversalGround(llcoords->getAt(cord).y,
                                    llcoords->getAt(cord).x);
            slcoords->add(geos::geom::Coordinate(ugm->Sample(),
                                                 ugm->Line()));
          } // end num coords in hole loop
          holes->push_back(globalFactory.createLinearRing(slcoords));
          delete slcoords;
          delete llcoords;
        } // end num holes in polygon loop

        // Convert the exterior ring of this polygon
        geos::geom::CoordinateSequence *slcoords = new geos::geom::CoordinateArraySequence();
        geos::geom::CoordinateSequence *llcoords =
          poly->getExteriorRing()->getCoordinates();

        // Convert each coordinate in the exterior ring of this polygon
        for(unsigned int cord = 0; cord < llcoords->getSize(); ++cord) {
          if (ugm->SetUniversalGround(llcoords->getAt(cord).y,
                                      llcoords->getAt(cord).x)) {
            slcoords->add(geos::geom::Coordinate(ugm->Sample(),
                                                ugm->Line()));
          }
        } // end exterior ring coordinate loop
       
        // Make sure that the line string is closed.
        if (slcoords->getSize() > 0 && !slcoords->front().equals(slcoords->back())) {
          slcoords->add(slcoords->front());
        }

        try {
          slPolys->push_back(globalFactory.createPolygon(
                              globalFactory.createLinearRing(slcoords), holes));
        }
        catch (std::exception &e) {
          throw IException(IException::Unknown,
                           QObject::tr("Unable to convert polygon from Lat/Lon to Sample/Line. The "
                                       "error given was [%1].").arg(e.what()),
                           _FILEINFO_);
        }
        
        delete llcoords;
      } // end num geometry in multi-poly

      // Create a new multipoly from all the new Sample/Line polygon(s)
      geos::geom::MultiPolygon *spikedPoly = globalFactory.createMultiPolygon(slPolys);

      if(spikedPoly->isValid() && !spikedPoly->isEmpty()) {
        return spikedPoly;
      }
      else {
        try {
          geos::geom::MultiPolygon *despikedPoly = Despike(spikedPoly);

          delete spikedPoly;
          spikedPoly = NULL;

          return despikedPoly;
        }
        catch (IException &e) {
          IString msg = "Unable to convert polygon from Lat/Lon to Sample/Line";
          throw IException(IException::Programmer, msg, _FILEINFO_);
        }
      }
    } // end else
  }


  geos::geom::MultiPolygon *PolygonTools::CopyMultiPolygon(const geos::geom::MultiPolygon *mpolygon) {

    vector<geos::geom::Geometry *> *polys = new vector<geos::geom::Geometry *>;
    for(unsigned int i = 0; i < mpolygon->getNumGeometries(); ++i) {
      polys->push_back((mpolygon->getGeometryN(i))->clone());
    }
    return globalFactory.createMultiPolygon(polys);
  }


  geos::geom::MultiPolygon *PolygonTools::CopyMultiPolygon(const geos::geom::MultiPolygon &mpolygon) {

    vector<geos::geom::Geometry *> *polys = new vector<geos::geom::Geometry *>;
    for(unsigned int i = 0; i < mpolygon.getNumGeometries(); ++i) {
      polys->push_back((mpolygon.getGeometryN(i))->clone());
    }
    return globalFactory.createMultiPolygon(polys);
  }


  QString PolygonTools::ToGML(const geos::geom::MultiPolygon *mpolygon, QString idString,
                              QString schema) {

    ostringstream os;

    //Write out the GML header
    os << "<?xml version=\"1.0\" encoding=\"utf-8\" ?>" << endl;
    os << "<ogr:FeatureCollection" << endl;
    os << "    xmlns:xsi=\"http://www.w3.org/2001/XMLSchema-instance\"" << endl;
    os << "    xsi:schemaLocation=\"http://ogr.maptools.org/ " << schema << "\"" << endl;
    os << "    xmlns:ogr=\"http://ogr.maptools.org/\"" << endl;
    os << "    xmlns:gml=\"http://www.opengis.net/gml\">" << endl;
    os << "  <gml:boundedBy>" << endl;
    os << "    <gml:Box>" << endl;
    os << "      <gml:coord>";
    os <<          "<gml:X>" << 
                       setprecision(15) << mpolygon->getEnvelopeInternal()->getMinX() << 
                       "</gml:X>";
    os <<          "<gml:Y>" << 
                       setprecision(15) << mpolygon->getEnvelopeInternal()->getMinY() << 
                       "</gml:Y>";
    os <<      "</gml:coord>" << endl;
    os << "      <gml:coord>";
    os <<          "<gml:X>" << 
                       setprecision(15) << mpolygon->getEnvelopeInternal()->getMaxX() <<
                       "</gml:X>";
    os <<          "<gml:Y>" << 
                       setprecision(15) << mpolygon->getEnvelopeInternal()->getMaxY() << 
                       "</gml:Y>";
    os <<        "</gml:coord>" << endl;
    os << "    </gml:Box>" << endl;
    os << "  </gml:boundedBy>" << endl << endl;
    os << "  <gml:featureMember>" << endl;
    os << "   <ogr:multi_polygon fid=\"0\">" << endl;
    os << "      <ogr:ID>" << idString << "</ogr:ID>" << endl;
    os << "      <ogr:geometryProperty><gml:MultiPolygon><gml:polygonMember>" <<
                     "<gml:Polygon><gml:outerBoundaryIs><gml:LinearRing><gml:coordinates>";


    for(unsigned int polyN = 0; polyN < mpolygon->getNumGeometries(); polyN++) {
      geos::geom::CoordinateSequence *pts = mpolygon->getGeometryN(polyN)->getCoordinates();

      for(unsigned int i = 0; i < pts->getSize(); i++) {
        double lon = pts->getAt(i).x;
        double lat = pts->getAt(i).y;

        os << setprecision(15) << lon << "," << setprecision(15) << lat << " ";
      }
    }

    os << "</gml:coordinates></gml:LinearRing></gml:outerBoundaryIs>" <<
              "</gml:Polygon></gml:polygonMember></gml:MultiPolygon>" <<
              "</ogr:geometryProperty>" << endl;
    os << "    </ogr:multi_polygon>" << endl;
    os << "  </gml:featureMember>" << endl;
    os << "</ogr:FeatureCollection>";

    return os.str().c_str();
  }


  QString PolygonTools::GMLSchema() {

    ostringstream os;

    os << "<?xml version=\"1.0\" encoding=\"UTF-8\"?>" << endl;
    os << "<xs:schema targetNamespace=\"http://ogr.maptools.org/\" "
              "xmlns:ogr=\"http://ogr.maptools.org/\" "
              "xmlns:xs=\"http://www.w3.org/2001/XMLSchema\" "
              "xmlns:gml=\"http://www.opengis.net/gml\" "
              "elementFormDefault=\"qualified\" "
              "version=\"1.0\">" << endl;
    os << "  <xs:import namespace=\"http://www.opengis.net/gml\" "
              "schemaLocation=\"http://schemas.opengis.net/gml/2.1.2/feature.xsd\"/>" << endl;
    os << "  <xs:element name=\"FeatureCollection\" "
                 "type=\"ogr:FeatureCollectionType\" "
                 "substitutionGroup=\"gml:_FeatureCollection\"/>" << endl;
    os << "  <xs:complexType name=\"FeatureCollectionType\">" << endl;
    os << "    <xs:complexContent>" << endl;
    os << "      <xs:extension base=\"gml:AbstractFeatureCollectionType\">" << endl;
    os << "        <xs:attribute name=\"lockId\" type=\"xs:string\" use=\"optional\"/>" << endl;
    os << "        <xs:attribute name=\"scope\" type=\"xs:string\" use=\"optional\"/>" << endl;
    os << "      </xs:extension>" << endl;
    os << "    </xs:complexContent>" << endl;
    os << "  </xs:complexType>" << endl;
    os << "  <xs:element name=\"multi_polygon\" "
                 "type=\"ogr:multi_polygon_Type\" "
                 "substitutionGroup=\"gml:_Feature\"/>" << endl;
    os << "  <xs:complexType name=\"multi_polygon_Type\">" << endl;
    os << "    <xs:complexContent>" << endl;
    os << "      <xs:extension base=\"gml:AbstractFeatureType\">" << endl;
    os << "        <xs:sequence>" << endl;
    os << "          <xs:element name=\"geometryProperty\" "
                         "type=\"gml:MultiPolygonPropertyType\" "
                         "nillable=\"true\" minOccurs=\"0\" maxOccurs=\"1\"/>" << endl;
    os << "          <xs:element name=\"fid\" nillable=\"true\" minOccurs=\"0\" "
                         "maxOccurs=\"1\">" << endl;
    os << "            <xs:simpleType>" << endl;
    os << "              <xs:restriction base=\"xs:string\">" << endl;
    os << "              </xs:restriction>" << endl;
    os << "            </xs:simpleType>" << endl;
    os << "          </xs:element>" << endl;
    os << "          <xs:element name=\"ID\" nillable=\"true\" minOccurs=\"0\" "
                         "maxOccurs=\"1\">" << endl;
    os << "            <xs:simpleType>" << endl;
    os << "              <xs:restriction base=\"xs:integer\">" << endl;
    os << "                <xs:totalDigits value=\"16\"/>" << endl;
    os << "              </xs:restriction>" << endl;
    os << "            </xs:simpleType>" << endl;
    os << "          </xs:element>" << endl;
    os << "        </xs:sequence>" << endl;
    os << "      </xs:extension>" << endl;
    os << "    </xs:complexContent>" << endl;
    os << "  </xs:complexType>" << endl;
    os << "</xs:schema>";

    return os.str().c_str();
  }


  geos::geom::MultiPolygon *PolygonTools::To180(geos::geom::MultiPolygon *poly360) {
    try {
      // Let's take the 360 pieces that sit between 180 and 360 and move them
      //   to -180 to 0. To accomplish this, make a poly that fits 180 -> 360
      //   degrees longitude and intersect (this gives us the pieces that sit
      //   >180). Move this intersection to the left. Then make a poly that fits
      //   0 to 180 and intersect with the original. These two combined are the
      //   result.
      geos::geom::CoordinateArraySequence *leftOf180Pts =
          new geos::geom::CoordinateArraySequence();
      leftOf180Pts->add(geos::geom::Coordinate(0, -90));
      leftOf180Pts->add(geos::geom::Coordinate(0, 90));
      leftOf180Pts->add(geos::geom::Coordinate(180, 90));
      leftOf180Pts->add(geos::geom::Coordinate(180, -90));
      leftOf180Pts->add(geos::geom::Coordinate(0, -90));

      geos::geom::LinearRing *leftOf180Geom =
          globalFactory.createLinearRing(leftOf180Pts);

      geos::geom::Polygon *leftOf180Poly =
          globalFactory.createPolygon(leftOf180Geom, NULL);

      geos::geom::CoordinateArraySequence *rightOf180Pts =
          new geos::geom::CoordinateArraySequence();
      rightOf180Pts->add(geos::geom::Coordinate(180, -90));
      rightOf180Pts->add(geos::geom::Coordinate(180, 90));
      rightOf180Pts->add(geos::geom::Coordinate(360, 90));
      rightOf180Pts->add(geos::geom::Coordinate(360, -90));
      rightOf180Pts->add(geos::geom::Coordinate(180, -90));

      geos::geom::LinearRing *rightOf180Geom =
          globalFactory.createLinearRing(rightOf180Pts);

      geos::geom::Polygon *rightOf180Poly =
          globalFactory.createPolygon(rightOf180Geom, NULL);

      geos::geom::Geometry *preserved = Intersect(leftOf180Poly, poly360);
      geos::geom::Geometry *moving = Intersect(rightOf180Poly, poly360);

      geos::geom::CoordinateSequence *movingPts = moving->getCoordinates();
      geos::geom::CoordinateSequence *movedPts =
          new geos::geom::CoordinateArraySequence();

      for(unsigned int i = 0; i < movingPts->getSize(); i ++) {
        movedPts->add(geos::geom::Coordinate(movingPts->getAt(i).x - 360.0,
                                             movingPts->getAt(i).y));
      }

      if(movedPts->getSize()) {
        movedPts->add(geos::geom::Coordinate(movedPts->getAt(0).x,
                                            movedPts->getAt(0).y));
      }

      geos::geom::Geometry *moved = globalFactory.createPolygon(
          globalFactory.createLinearRing(movedPts), NULL);

      std::vector<geos::geom::Geometry *> *geomsForCollection = new
          std::vector<geos::geom::Geometry *>;
      geomsForCollection->push_back(preserved);
      geomsForCollection->push_back(moved);

      geos::geom::GeometryCollection *the180Polys =
          Isis::globalFactory.createGeometryCollection(geomsForCollection);

      geos::geom::MultiPolygon *result = MakeMultiPolygon(the180Polys);
      delete the180Polys;
      the180Polys = NULL;

      geos::geom::MultiPolygon *fixedResult = FixSeam(result);
      delete result;
      result = NULL;

      return fixedResult;
    }
    catch(geos::util::GEOSException *exc) {
      IString msg = "Conversion to 180 failed. The reason given was [" +
          IString(exc->what()) + "]";
      delete exc;
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }
    catch(...) {
      IString msg = "Conversion to 180 failed. Could not determine the reason";
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }
  }


  double PolygonTools::Thickness(const geos::geom::MultiPolygon *mpolygon) {
    const geos::geom::Envelope *envelope = mpolygon->getEnvelopeInternal();

    double x = fabs(envelope->getMaxX() - envelope->getMinX());
    double y = fabs(envelope->getMaxY() - envelope->getMinY());
    double extent = max(x, y);

    return mpolygon->getArea() / (extent * extent);
  }


  geos::geom::MultiPolygon *PolygonTools::Despike(const geos::geom::Geometry *geom) {
    geos::geom::MultiPolygon *spiked = MakeMultiPolygon(geom);
    geos::geom::MultiPolygon *despiked = Despike(spiked);

    delete spiked;
    spiked = NULL;

    return despiked;
  }


  geos::geom::MultiPolygon *PolygonTools::Despike(const geos::geom::MultiPolygon *multiPoly) {
    // Despike each polygon in the multipolygon
    vector<geos::geom::Geometry *> *newPolys = new vector<geos::geom::Geometry *>;
    for(unsigned int g = 0; g < multiPoly->getNumGeometries(); ++g) {
      const geos::geom::Polygon *poly =
          dynamic_cast<const geos::geom::Polygon *>(multiPoly->getGeometryN(g));

      // Despike each hole inside this polygon
      vector<geos::geom::Geometry *> *holes = new vector<geos::geom::Geometry *>;
      for(unsigned int h = 0; h < poly->getNumInteriorRing(); ++h) {
        const geos::geom::LineString *ls = poly->getInteriorRingN(h);
        geos::geom::LinearRing *lr;

        // If the hole is not valid fix it
        // If the hole is NOT valid despike it
        lr = Despike(ls);

        if(!lr->isValid()) {
          geos::geom::LinearRing *fixed = FixGeometry(lr);
          delete lr;
          lr = fixed;
        }

        // Save this hole if it is not too small
        if(!lr->isEmpty()) {
          holes->push_back(lr);
        }
        else {
          delete lr;
        }
      } // End holes loop

      // Work on the main polygon
      const geos::geom::LineString *ls = poly->getExteriorRing();
      geos::geom::LinearRing *lr;

      lr = Despike(ls);

      try {
        if(!lr->isValid()) {
          geos::geom::LinearRing *fixed = FixGeometry(lr);
          delete lr;
          lr = fixed;
        }
      }
      catch(IException &e) {
        // Sometimes despike and fix fail, but the input is really valid. We can just go
        // with the non-despiked polygon.
        if(ls->isValid() && ls->getGeometryTypeId() == geos::geom::GEOS_LINEARRING) {
          lr = dynamic_cast<geos::geom::LinearRing *>(ls->clone());
        }
        else {
          throw;
        }
      }

      // Create a new polygon with the holes and save it
      if(!lr->isEmpty()) {
        geos::geom::Polygon *tp = Isis::globalFactory.createPolygon(lr, holes);

        if(tp->isEmpty() || !tp->isValid()) {
          delete tp;
          newPolys->push_back(poly->clone());
        }
        else {
          newPolys->push_back(tp);
        }
      }
    } // End polygons loop

    // Create a new multipoly from the polygon(s)
    geos::geom::MultiPolygon *mp = Isis::globalFactory.createMultiPolygon(newPolys);

    if(!mp->isValid() || mp->isEmpty()) {
      delete mp;
      mp = NULL;
      IString msg = "Despike failed to correct the polygon";
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }

    // if multipoly area changes more than 25% we did something bad to the multipolygon
    if(fabs((mp->getArea() / multiPoly->getArea()) - 1.0) > 0.50) {
      IString msg = "Despike failed to correct the polygon " + mp->toString();
      delete mp;

      throw IException(IException::Programmer, msg, _FILEINFO_);
    }

    return mp;
  }


  geos::geom::LinearRing *PolygonTools::Despike(const geos::geom::LineString *lineString) {
    geos::geom::CoordinateSequence *vertices = lineString->getCoordinates();

    // We need a full polygon to despike = 3 points (+1 for end==first) = at least 4 points
    if(vertices->getSize() < 4) {
      delete vertices;
      return Isis::globalFactory.createLinearRing(geos::geom::CoordinateArraySequence());
    }

    // delete one of the duplicate first/end coordinates,
    //   spikes can occur here and the duplicate points throws off the test
    vertices->deleteAt(vertices->getSize() - 1);

    // Index will become negative if our first points are spiked, so we need
    //   to make it big enough to encapsulate the (unsigned int) from getSize()
    //   and allow negative numbers.
    for(long index = 0; index < (long)vertices->getSize(); index++) {
      // If we're under 3 vertices, we've despiked all the points out of the polygon :(
      if(vertices->getSize() < 3) {
        break;
      }

      // These are the array indices that we're going to test for a spike,
      //   the middle one is the one that's spiked if IsSpiked(...) is true.
      long testCoords[3] = {
        index - 1,
        index,
        index + 1
      };

      // Make sure the index is inside of our coordinates array (we'll have both too small/too 
      // big of an index)
      for(int j = 0; j < 3; j++) {
        while(testCoords[j] < 0) {
          testCoords[j] += vertices->getSize();
        }

        while(testCoords[j] >= (long)vertices->getSize()) {
          testCoords[j] -= vertices->getSize();
        }
      }

      // Test the middle point for a spike
      if(IsSpiked(vertices->getAt(testCoords[0]),
                  vertices->getAt(testCoords[1]),
                  vertices->getAt(testCoords[2]))) {
        // It's spiked, delete it
        vertices->deleteAt(testCoords[1]);

        // Back up to the first test that is affected by this change
        index -= 2;
      }
    }

    if(vertices->getSize() < 3) {
      delete vertices;
      vertices = NULL;

      return Isis::globalFactory.createLinearRing(geos::geom::CoordinateArraySequence());
    }
    else {
      // Duplicate the first vertex as the last to close the polygon
      vertices->add(vertices->getAt(0));
      return Isis::globalFactory.createLinearRing(vertices);
    }
  }

  bool PolygonTools::IsSpiked(geos::geom::Coordinate first, 
                              geos::geom::Coordinate middle, 
                              geos::geom::Coordinate last) {
    return TestSpiked(first, middle, last) || TestSpiked(last, middle, first);
  }

  bool PolygonTools::TestSpiked(geos::geom::Coordinate first, geos::geom::Coordinate middle, 
                                geos::geom::Coordinate last) {
    geos::geom::Point *firstPt = Isis::globalFactory.createPoint(first);
    geos::geom::Point *middlePt = Isis::globalFactory.createPoint(middle);
    geos::geom::Point *lastPt = Isis::globalFactory.createPoint(last);

    geos::geom::CoordinateSequence *coords = new geos::geom::CoordinateArraySequence();
    coords->add(first);
    coords->add(middle);
    geos::geom::LineString *line = Isis::globalFactory.createLineString(coords); // line takes ownership

    // The lower the tolerance, the less this algorithm removes and thus
    //   the better chance of success in findimageoverlaps. However, if you
    //   lower the tolerance then there is also a greater chance of programs
    //   such as autoseed failing. 1% is the current tolerance.
    double tolerance = line->getLength() / 100.0;

    bool spiked = true;

    double distanceLastMiddle = geos::operation::distance::DistanceOp::distance(lastPt, middlePt);
    double distanceLastLine = geos::operation::distance::DistanceOp::distance(lastPt, line);

    if(distanceLastMiddle == 0.0) return true; // get rid of same point

    // Checks the ratio of the distance between the last point and the line, and the last point
    // and the middle point if the ratio is very small then there is a spike
    if(distanceLastLine / distanceLastMiddle >= .05) {
      spiked = false;
    }

    // If the point is away from the line, keep it
    if(spiked && distanceLastLine > tolerance) {
      spiked = false;
    }

    if(!spiked) {
      geos::geom::CoordinateSequence *coords = new geos::geom::CoordinateArraySequence();
      coords->add(first);
      coords->add(middle);
      coords->add(last);
      coords->add(first);

      // shell takes ownership of coords
      geos::geom::LinearRing *shell = Isis::globalFactory.createLinearRing(coords);
      std::vector<geos::geom::Geometry *> *empty = new std::vector<geos::geom::Geometry *>;

      // poly takes ownership of shell and empty
      geos::geom::Polygon *poly = Isis::globalFactory.createPolygon(shell, empty);


      // if these 3 points define a straight line then the middle is worthless (defines nothing) 
      // or problematic
      if(poly->getArea() < 1.0e-10) {
        spiked = true;
      }

      delete poly;
    }


    delete firstPt;
    delete middlePt;
    delete lastPt;
    delete line;

    return spiked;
  }


  geos::geom::Geometry *PolygonTools::Intersect(const geos::geom::Geometry *geom1, 
                                                const geos::geom::Geometry *geom2) {
    try {
      return Operate(geom1, geom2, (unsigned int)geos::operation::overlay::OverlayOp::opINTERSECTION);
    }
    catch(geos::util::GEOSException *exc) {
      IString msg = "Intersect operation failed. The reason given was [" + IString(exc->what()) + "]";
      delete exc;
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }
    catch(IException &e) {
      IString msg = "Intersect operation failed";
      throw IException(e, IException::Programmer, msg, _FILEINFO_);
    }
    catch(...) {
      IString msg = "Intersect operation failed for an unknown reason";
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }
  }


  geos::geom::Geometry *PolygonTools::Operate(const geos::geom::Geometry *geom1,
                                              const geos::geom::Geometry *geom2, 
                                              unsigned int opcode) {

    geos::operation::overlay::OverlayOp::OpCode code =
      (geos::operation::overlay::OverlayOp::OpCode)opcode;

    geos::geom::Geometry *result = NULL;
    bool failed = true;
    geos::geom::Geometry *geomFirst  = MakeMultiPolygon(geom1);
    geos::geom::Geometry *geomSecond = MakeMultiPolygon(geom2);

    geos::operation::overlay::snap::GeometrySnapper snap(*geomFirst);
    geos::geom::Geometry *geomSnapped = snap.snapTo(*geomSecond, 1.0e-10)->clone();
    if(!geomSnapped->isValid()) {
      delete geomSnapped;
    }
    else {
      delete geomFirst;
      geomFirst = geomSnapped;
    }

    unsigned int precision = 15;
    unsigned int minPrecision = 13;
    while(failed) {
      try {
        // C++11: the geos BinaryOp returns an auto_ptr, we use release() to create a unique_ptr.
        std::unique_ptr< geos::geom::Geometry > resultAuto(
          BinaryOp(geomFirst, geomSecond, geos::operation::overlay::overlayOp(code)).release());
        failed = false;
        result = resultAuto->clone();
      }
      catch(geos::util::GEOSException *exc) {
        // Just in case the clone failed....
        if(!failed || precision == minPrecision) throw;

        delete exc;
      }
      catch(...) {
        if(precision == minPrecision) {
          IString msg = "An unknown geos error occurred";
          throw IException(IException::Programmer, msg, _FILEINFO_);
        }

        if(!failed) {
          IString msg = "An unknown geos error occurred when attempting to clone a geometry";
          throw IException(IException::Programmer, msg, _FILEINFO_);
        }
      }

      // try reducing the precision
      if(failed) {
        precision --;
        geos::geom::Geometry *tmp = ReducePrecision(geomFirst, precision);
        delete geomFirst;
        geomFirst = tmp;

        tmp = ReducePrecision(geomSecond, precision);
        delete geomSecond;
        geomSecond = tmp;
      }
    }

    delete geomFirst;
    geomFirst = NULL;

    delete geomSecond;
    geomSecond = NULL;

    if(result && !result->isValid()) {
      try {
        geos::geom::Geometry *newResult = FixGeometry(result);

        if(fabs(newResult->getArea() / result->getArea() - 1.0) > 0.50) {
          delete newResult;
          delete result;

          IString msg = "Operation [" + IString((int)opcode) + "] failed";
          throw IException(IException::Programmer, msg, _FILEINFO_);
        }

        delete result;
        result = newResult;
      }
      catch(IException &e) {
        IString msg = "Operation [" + IString((int)opcode) + "] failed";
        throw IException(IException::Programmer, msg, _FILEINFO_);
      }
    }

    if(result == NULL) {
      IString msg = "Operation [" + IString((int)opcode) + " failed";
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }

    return result;
  }

  geos::geom::Geometry *PolygonTools::FixGeometry(const geos::geom::Geometry *geom) {
    if(geom->getGeometryTypeId() == geos::geom::GEOS_MULTIPOLYGON) {
      return FixGeometry(dynamic_cast<const geos::geom::MultiPolygon *>(geom));
    }
    if(geom->getGeometryTypeId() == geos::geom::GEOS_LINEARRING) {
      return FixGeometry(dynamic_cast<const geos::geom::LinearRing *>(geom));
    }
    if(geom->getGeometryTypeId() == geos::geom::GEOS_POLYGON) {
      return FixGeometry(dynamic_cast<const geos::geom::Polygon *>(geom));
    }
    if(geom->getGeometryTypeId() == geos::geom::GEOS_GEOMETRYCOLLECTION) {
      return FixGeometry(MakeMultiPolygon(geom));
    }
    else {
      IString msg = "PolygonTools::FixGeometry does not support [" + GetGeometryName(geom) + "]";
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }
  }

  geos::geom::MultiPolygon *PolygonTools::FixGeometry(const geos::geom::MultiPolygon *poly) {
    // Maybe two points are on top of each other
    vector<geos::geom::Geometry *> *newPolys = new vector<geos::geom::Geometry *>;

    // Convert each polygon in this multi-polygon
    for(unsigned int geomIndex = 0; geomIndex < poly->getNumGeometries(); geomIndex ++) {
      geos::geom::Polygon *fixedpoly = FixGeometry(
          dynamic_cast<const geos::geom::Polygon *>(
            poly->getGeometryN(geomIndex)));
      if(fixedpoly->isValid()) {
        newPolys->push_back(fixedpoly);
      }
      else {
        delete fixedpoly;
      }
      fixedpoly = NULL;
    }

    geos::geom::MultiPolygon *mp = Isis::globalFactory.createMultiPolygon(newPolys);
    return mp;
  }


  geos::geom::Polygon *PolygonTools::FixGeometry(const geos::geom::Polygon *poly) {

    // Convert each hole inside this polygon
    vector<geos::geom::Geometry *> *holes = new vector<geos::geom::Geometry *>;
    for(unsigned int holeIndex = 0; holeIndex < poly->getNumInteriorRing(); holeIndex ++) {
      const geos::geom::LineString *thisHole = poly->getInteriorRingN(holeIndex);

      // We hope they are all linear rings (closed/simple), but if not just leave it be
      if(thisHole->getGeometryTypeId() != geos::geom::GEOS_LINEARRING) {
        holes->push_back(thisHole->clone());
        continue;
      }

      try {
        geos::geom::LinearRing *newHole = FixGeometry((geos::geom::LinearRing *)thisHole);
        holes->push_back(newHole);
      }
      catch (IException &e) {
        IString msg = "Failed when attempting to fix interior ring of multipolygon";
        throw IException(IException::Programmer, msg, _FILEINFO_);
      }
    } // end num holes in polygon loop


    const geos::geom::LineString *exterior = poly->getExteriorRing();

    try {
      geos::geom::LinearRing *newExterior = NULL;

      if(exterior->getGeometryTypeId() == geos::geom::GEOS_LINEARRING) {
        newExterior = FixGeometry((geos::geom::LinearRing *)exterior);
      }
      else {
        IString msg = "Failed when attempting to fix exterior ring of polygon. The exterior "
                      "ring is not simple and closed";
        throw IException(IException::Programmer, msg, _FILEINFO_);
      }

      return globalFactory.createPolygon(newExterior, holes);
    }
    catch (IException &e) {
      IString msg = "Failed when attempting to fix exterior ring of polygon";
      throw IException(e, IException::Programmer, msg, _FILEINFO_);
    }
  }


  geos::geom::LinearRing *PolygonTools::FixGeometry(const geos::geom::LinearRing *ring) {

    geos::geom::CoordinateSequence *coords = ring->getCoordinates();

    // Check this, just in case
    if(coords->getSize() < 4) {
      return globalFactory.createLinearRing(new geos::geom::DefaultCoordinateSequence());
    }

    geos::geom::CoordinateSequence *newCoords = new geos::geom::DefaultCoordinateSequence();
    const geos::geom::Coordinate *lastCoordinate = &coords->getAt(0);
    newCoords->add(*lastCoordinate);

    // Convert each coordinate in this hole
    for(unsigned int coordIndex = 1; coordIndex < coords->getSize() - 1; coordIndex ++) {
      const geos::geom::Coordinate *thisCoordinate = &coords->getAt(coordIndex);

      // we're going to compare the decimal place of the current point to the decimal place
      //   of the difference, if they are drastically different then geos might not be seeing them
      //   correctly.
      double difference[2] = {
        lastCoordinate->x - thisCoordinate->x,
        lastCoordinate->y - thisCoordinate->y,
      };

      // geos isnt differentiating between points this close
      double minDiff = fabs(DecimalPlace(thisCoordinate->x) - DecimalPlace(difference[0]));

      minDiff = min(minDiff, fabs(DecimalPlace(thisCoordinate->y) - DecimalPlace(difference[1])));

      // Cases where the difference in one direction is exactly zero, and the other direction is 
      // next to zero appear often enough (especially in despike).
      if(difference[0] == 0.0 && difference[1] != 0.0) {
        // subtracted the two points, got deltaX = 0.0, use the y difference decimal place
        minDiff = fabs(DecimalPlace(thisCoordinate->y) - DecimalPlace(difference[1]));
      }
      else if(difference[1] == 0.0 && difference[0] != 0.0) {
        // subtracted the two points, got deltaY = 0.0, use the x difference decimal place
        minDiff = fabs(DecimalPlace(thisCoordinate->x) - DecimalPlace(difference[0]));
      }
      else if(difference[0] == 0.0 && difference[1] == 0.0) {
        // subtracted the two points, got 0.0, so it's same point... make sure it gets ignored!
        minDiff = 1E99;
      }

      // geos has a hard time differentiating when points get too close...
      if(minDiff <= 15) {
        newCoords->add(*thisCoordinate);
        lastCoordinate = thisCoordinate;
      }
    } // end num coords in hole loop

    newCoords->add(geos::geom::Coordinate(newCoords->getAt(0).x, newCoords->getAt(0).y));
    geos::geom::LinearRing *newRing = NULL;

    // Now that we've weeded out any bad coordinates, let's rebuild the geometry
    try {
      if(newCoords->getSize() > 3) {
        newRing = globalFactory.createLinearRing(newCoords);
      }
      else {
        delete newCoords;
        newCoords = NULL;
      }
    }
    catch(geos::util::GEOSException *exc) {
      delete exc;
      exc = NULL;

      IString msg = "Error when attempting to fix linear ring";
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }

    if(newRing && !newRing->isValid() && ring->isValid()) {
      IString msg = "Failed when attempting to fix linear ring";
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }
    else if(!newRing || !newRing->isValid()) {
      if(newRing) {
        delete newRing;
      }

      newRing = dynamic_cast<geos::geom::LinearRing *>(ring->clone());
    }

    return newRing;
  }


  int PolygonTools::DecimalPlace(double num) {
    // 0.0 = decimal place 0
    if(num == 0.0) return 0;

    num = fabs(num);

    int decimalPlace = 1;
    while(num < 1.0) {
      num *= 10.0;
      decimalPlace --;
    }

    while(num > 10.0) {
      num /= 10.0;
      decimalPlace ++;
    }

    return decimalPlace;
  }

  geos::geom::Geometry *PolygonTools::Difference(const geos::geom::Geometry *geom1, 
                                                 const geos::geom::Geometry *geom2) {
    try {
      return Operate(geom1, geom2, (unsigned int)geos::operation::overlay::OverlayOp::opDIFFERENCE);
    }
    catch(geos::util::GEOSException *exc) {
      IString msg = "Difference operation failed. The reason given was [" + 
                    IString(exc->what()) + "]";
      delete exc;
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }
    catch(IException &e) {
      IString msg = "Difference operation failed";
      throw IException(e, IException::Programmer, msg, _FILEINFO_);
    }
    catch(...) {
      IString msg = "Difference operation failed for an unknown reason";
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }
  }


  geos::geom::MultiPolygon *PolygonTools::MakeMultiPolygon(const geos::geom::Geometry *geom) {
    // The area of the geometry is too small, so just ignore it.
    if(geom->isEmpty()) {
      return Isis::globalFactory.createMultiPolygon();
    }

    else if(geom->getArea() - DBL_EPSILON <= DBL_EPSILON) {
      return Isis::globalFactory.createMultiPolygon();
    }

    else if(geom->getGeometryTypeId() == geos::geom::GEOS_MULTIPOLYGON) {
      return dynamic_cast<geos::geom::MultiPolygon *>(geom->clone());
    }

    else if(geom->getGeometryTypeId() == geos::geom::GEOS_POLYGON) {
      vector<geos::geom::Geometry *> *polys = new vector<geos::geom::Geometry *>;
      polys->push_back(geom->clone());
      geos::geom::MultiPolygon *mp = Isis::globalFactory.createMultiPolygon(polys);
      return mp;
    }

    else if(geom->getGeometryTypeId() == geos::geom::GEOS_GEOMETRYCOLLECTION) {
      vector<geos::geom::Geometry *> * polys =
          new vector<geos::geom::Geometry *>;
      const geos::geom::GeometryCollection *gc =
          dynamic_cast<const geos::geom::GeometryCollection *>(geom);
      for(unsigned int i = 0; i < gc->getNumGeometries(); ++i) {
        geos::geom::MultiPolygon *subMultiPoly =
            MakeMultiPolygon(gc->getGeometryN(i));

        for(unsigned int subPoly = 0;
            subPoly < subMultiPoly->getNumGeometries();
            subPoly ++) {
          const geos::geom::Polygon *poly =
              dynamic_cast<const geos::geom::Polygon *>(
                subMultiPoly->getGeometryN(subPoly));
          polys->push_back(dynamic_cast<geos::geom::Polygon *>(poly->clone()));
        }
      }

      geos::geom::MultiPolygon *mp =
          Isis::globalFactory.createMultiPolygon(polys);
      if(mp->getArea() - DBL_EPSILON <= DBL_EPSILON) {
        delete mp;
        mp = Isis::globalFactory.createMultiPolygon();
      }

      return mp;
    }

    // All other geometry types are invalid so ignore them
    else {
      return Isis::globalFactory.createMultiPolygon();
    }
  }


  geos::geom::MultiPolygon *PolygonTools::FixSeam(
      const geos::geom::Polygon *polyA, const geos::geom::Polygon *polyB) {
    geos::geom::CoordinateSequence *polyAPoints = polyA->getCoordinates();
    geos::geom::CoordinateSequence *polyBPoints = polyB->getCoordinates();

    unsigned int aIntersectionBegin = 0;
    unsigned int aIntersectionEnd = 0;
    unsigned int bIntersectionBegin = 0;
    unsigned int bIntersectionEnd = 0;

    bool intersectionStarted = false;
    bool intersectionEnded = false;

    unsigned int lastBMatch  = 0;
    for (unsigned int i = 0;
        !intersectionEnded && i < polyAPoints->getSize();
        i++) {

      bool foundEquivalent = false;

      geos::geom::Coordinate coordA = polyAPoints->getAt(i);
      coordA = *ReducePrecision(&coordA, 13);

      for (unsigned int j = 0;
           !foundEquivalent && j < polyBPoints->getSize();
           j++) {
        geos::geom::Coordinate coordB = polyBPoints->getAt(j);
        coordB = *ReducePrecision(&coordB, 13);

        foundEquivalent = coordA.equals(coordB);

        if (foundEquivalent) lastBMatch = j;

        if (foundEquivalent && !intersectionStarted) {
          intersectionStarted = true;
          aIntersectionBegin = i;
          bIntersectionBegin = j;
        }
      }

      if (!foundEquivalent && intersectionStarted && !intersectionEnded) {
        intersectionEnded = true;
        aIntersectionEnd = i;
        bIntersectionEnd = lastBMatch;
      }
    }

    geos::geom::MultiPolygon * result = NULL;
    if (intersectionStarted && intersectionEnded) {
      geos::geom::CoordinateSequence *merged =
          new geos::geom::CoordinateArraySequence;

      unsigned int i = 0;
      for (i = 0; i < aIntersectionBegin; i ++) {
        merged->add(polyAPoints->getAt(i));
      }

      i = bIntersectionBegin;
      while (i != bIntersectionEnd) {
        merged->add(polyBPoints->getAt(i));
        i++;
        if (i >= polyBPoints->getSize()) i = 0;
      }

      for (i = aIntersectionEnd; i < polyAPoints->getSize() - 1; i++) {
        merged->add(polyAPoints->getAt(i));
      }

      merged->add(merged->getAt(0));
      result = MakeMultiPolygon(globalFactory.createPolygon(
          globalFactory.createLinearRing(merged), NULL));
    }

    return result;
  }


  geos::geom::MultiPolygon *PolygonTools::FixSeam(
      const geos::geom::MultiPolygon *poly) {

    std::vector<geos::geom::Geometry *> *polys =
        new std::vector<geos::geom::Geometry *>;


    for(unsigned int copyIndex = 0;
        copyIndex < poly->getNumGeometries();
        copyIndex ++) {
      polys->push_back(poly->getGeometryN(copyIndex)->clone());
    }

    unsigned int outerPolyIndex = 0;

    while(outerPolyIndex + 1 < polys->size()) {
      unsigned int innerPolyIndex = outerPolyIndex + 1;

      while(innerPolyIndex < polys->size()) {
        geos::geom::MultiPolygon *fixedPair = FixSeam(
            dynamic_cast<geos::geom::Polygon *>(polys->at(outerPolyIndex)),
            dynamic_cast<geos::geom::Polygon *>(polys->at(innerPolyIndex)));

        if(fixedPair != NULL) {
          geos::geom::Geometry *oldInnerPoly = polys->at(innerPolyIndex);
          geos::geom::Geometry *oldOuterPoly = polys->at(outerPolyIndex);

          polys->erase(polys->begin() + innerPolyIndex);
          (*polys)[outerPolyIndex] = fixedPair->getGeometryN(0)->clone();
          innerPolyIndex = outerPolyIndex + 1;

          delete fixedPair;
          fixedPair = NULL;

          delete oldInnerPoly;
          oldInnerPoly = NULL;

          delete oldOuterPoly;
          oldOuterPoly = NULL;
        }
        else {
          innerPolyIndex ++;
        }
      }

      outerPolyIndex ++;
    }

    return globalFactory.createMultiPolygon(polys);
  }


  geos::geom::Geometry *PolygonTools::ReducePrecision(const geos::geom::Geometry *geom, 
                                                      unsigned int precision) {
    if(geom->getGeometryTypeId() == geos::geom::GEOS_MULTIPOLYGON) {
      return ReducePrecision(
          dynamic_cast<const geos::geom::MultiPolygon *>(geom), precision);
    }
    if(geom->getGeometryTypeId() == geos::geom::GEOS_LINEARRING) {
      return ReducePrecision(
          dynamic_cast<const geos::geom::LinearRing *>(geom), precision);
    }
    if(geom->getGeometryTypeId() == geos::geom::GEOS_POLYGON) {
      return ReducePrecision(
          dynamic_cast<const geos::geom::Polygon *>(geom), precision);
    }
    if(geom->getGeometryTypeId() == geos::geom::GEOS_GEOMETRYCOLLECTION) {
      return ReducePrecision(MakeMultiPolygon(geom), precision);
    }
    else {
      IString msg = "PolygonTools::ReducePrecision does not support [" + 
                    GetGeometryName(geom) + "]";
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }
  }


  geos::geom::MultiPolygon *PolygonTools::ReducePrecision(const geos::geom::MultiPolygon *poly, 
                                                          unsigned int precision) {
    // Maybe two points are on top of each other
    vector<geos::geom::Geometry *> *newPolys = new vector<geos::geom::Geometry *>;

    // Convert each polygon in this multi-polygon
    for(unsigned int geomIndex = 0; geomIndex < poly->getNumGeometries(); geomIndex ++) {
      geos::geom::Geometry *lowerPrecision = ReducePrecision(
          dynamic_cast<const geos::geom::Polygon *>(
            poly->getGeometryN(geomIndex)),
          precision);

      if(!lowerPrecision->isEmpty()) {
        newPolys->push_back(lowerPrecision);
      }
      else {
        delete lowerPrecision;
      }
    }

    geos::geom::MultiPolygon *mp = Isis::globalFactory.createMultiPolygon(newPolys);
    return mp;
  }


  geos::geom::Polygon *PolygonTools::ReducePrecision(const geos::geom::Polygon *poly, 
                                                     unsigned int precision) {
    // Convert each hole inside this polygon
    vector<geos::geom::Geometry *> *holes = new vector<geos::geom::Geometry *>;
    for(unsigned int holeIndex = 0; holeIndex < poly->getNumInteriorRing(); holeIndex ++) {
      const geos::geom::LineString *thisHole = poly->getInteriorRingN(holeIndex);

      // We hope they are all linear rings (closed/simple), but if not just leave it be
      if(thisHole->getGeometryTypeId() != geos::geom::GEOS_LINEARRING) {
        holes->push_back(thisHole->clone());
        continue;
      }

      try {
        geos::geom::LinearRing *newHole = ReducePrecision((geos::geom::LinearRing *)thisHole, 
                                                          precision);

        if(!newHole->isEmpty()) {
          holes->push_back(newHole);
        }
        else {
          delete newHole;
        }

      }
      catch(IException &e) {
        IString msg = "Failed when attempting to fix interior ring of multipolygon";
        throw IException(e, IException::Programmer, msg, _FILEINFO_);
      }
    } // end num holes in polygon loop


    const geos::geom::LineString *exterior = poly->getExteriorRing();

    try {
      geos::geom::LinearRing *newExterior = NULL;

      if(exterior->getGeometryTypeId() == geos::geom::GEOS_LINEARRING) {
        newExterior = ReducePrecision((geos::geom::LinearRing *)exterior, precision);
      }
      else {
        IString msg = "Failed when attempting to fix exterior ring of polygon. The exterior "
                      "ring is not simple and closed";
        throw IException(IException::Programmer, msg, _FILEINFO_);
      }

      return globalFactory.createPolygon(newExterior, holes);
    }
    catch(IException &e) {
      IString msg = "Failed when attempting to fix exterior ring of polygon";
      throw IException(e, IException::Programmer, msg, _FILEINFO_);
    }
  }


  geos::geom::LinearRing *PolygonTools::ReducePrecision(const geos::geom::LinearRing *ring, 
                                                        unsigned int precision) {
    geos::geom::CoordinateSequence *coords = ring->getCoordinates();

    // Check this, just in case
    if(coords->getSize() <= 0) return dynamic_cast<geos::geom::LinearRing *>(
        ring->clone());

    geos::geom::CoordinateSequence *newCoords = new geos::geom::DefaultCoordinateSequence();
    geos::geom::Coordinate *coord = ReducePrecision(&coords->getAt(0), precision);
    newCoords->add(*coord);
    delete coord;
    coord = NULL;

    // Convert each coordinate in this ring
    for(unsigned int coordIndex = 1; coordIndex < coords->getSize() - 1; coordIndex ++) {
      const geos::geom::Coordinate *thisCoordinate = &coords->getAt(coordIndex);
      coord = ReducePrecision(thisCoordinate, precision);
      newCoords->add(*coord);
      delete coord;
      coord = NULL;
    }

    newCoords->add(geos::geom::Coordinate(newCoords->getAt(0).x, newCoords->getAt(0).y));
    geos::geom::LinearRing *newRing = NULL;

    // Now that we've weeded out any bad coordinates, let's rebuild the geometry
    try {
      newRing = globalFactory.createLinearRing(newCoords);
    }
    catch(geos::util::GEOSException *exc) {
      delete exc;
      exc = NULL;

      IString msg = "Error when attempting to reduce precision of linear ring";
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }

    // try to despike
    try {
      geos::geom::LinearRing *tmp = Despike(newRing);
      delete newRing;
      newRing = tmp;
    }
    catch(IException &e) {
    }

    if(!newRing->isValid()) {
      IString msg = "Failed when attempting to reduce precision of linear ring";
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }

    return newRing;
  }


  geos::geom::Coordinate *PolygonTools::ReducePrecision(const geos::geom::Coordinate *coord, 
                                                        unsigned int precision) {
    return new geos::geom::Coordinate(
             ReducePrecision(coord->x, precision),
             ReducePrecision(coord->y, precision),
             ReducePrecision(coord->z, precision)
           );
  }


  double PolygonTools::ReducePrecision(double num, unsigned int precision) {
    double result = num;

    // If num == nan then this test will fail
    if (num == num) {
      int decimalPlace = DecimalPlace(num);
      double factor = pow(10.0, (int)decimalPlace);

      // reduced num is in the form 0.nnnnnnnnnn...
      double reducedNum = num / factor;

      double cutoff = pow(10.0, (int)precision);
      double round_offset = (num < 0) ? -0.5 : 0.5;

      // cast off the digits past the precision's place
      reducedNum = ((long long)(reducedNum * cutoff + round_offset)) / cutoff;
      result = reducedNum * factor;
    }

    return result;
  }


  QString PolygonTools::GetGeometryName(const geos::geom::Geometry *geom) {
    switch(geom->getGeometryTypeId()) {
      case geos::geom::GEOS_POINT:
        return "Point";
      case geos::geom::GEOS_LINESTRING:
        return "Line String";
      case geos::geom::GEOS_LINEARRING:
        return "Linear Ring";
      case geos::geom::GEOS_POLYGON:
        return "Polygon";
      case geos::geom::GEOS_MULTIPOINT:
        return "Multi Point";
      case geos::geom::GEOS_MULTILINESTRING:
        return "Multi Line String";
      case geos::geom::GEOS_MULTIPOLYGON:
        return "Multi Polygon";
      case geos::geom::GEOS_GEOMETRYCOLLECTION:
        return "Geometry Collection";
      default:
        return "UNKNOWN";
    }
  }


  bool PolygonTools::Equal(const geos::geom::MultiPolygon *poly1, 
                           const geos::geom::MultiPolygon *poly2) {

    vector<const geos::geom::Polygon *> polys1;
    vector<const geos::geom::Polygon *> polys2;

    if(poly1->getNumGeometries() != poly2->getNumGeometries())  return false;

    // Convert each polygon in this multi-polygon
    for(unsigned int geomIndex = 0; geomIndex < poly1->getNumGeometries(); geomIndex ++) {
      polys1.push_back(dynamic_cast<const geos::geom::Polygon *>(
          poly1->getGeometryN(geomIndex)));
      polys2.push_back(dynamic_cast<const geos::geom::Polygon *>(
          poly2->getGeometryN(geomIndex)));
    }

    for(int p1 = polys1.size() - 1; (p1 >= 0) && polys1.size(); p1 --) {
      for(int p2 = polys2.size() - 1; (p2 >= 0) && polys2.size(); p2 --) {
        if(Equal(polys1[p1], polys2[p2])) {
          // Delete polys1[p1] by replacing it with the last Polygon in polys1
          polys1[p1] = polys1[polys1.size()-1];
          polys1.resize(polys1.size() - 1);
          // Delete polys2[p2] by replacing it with the last Polygon in polys2
          polys2[p2] = polys2[polys2.size()-1];
          polys2.resize(polys2.size() - 1);
        }
      }
    }

    return (polys1.size() == 0) && (polys2.size() == 0);
  }


  bool PolygonTools::Equal(const geos::geom::Polygon *poly1, const geos::geom::Polygon *poly2) {
    vector<const geos::geom::LineString *> holes1;
    vector<const geos::geom::LineString *> holes2;

    if(poly1->getNumInteriorRing() != poly2->getNumInteriorRing())  return false;

    if(!Equal(poly1->getExteriorRing(), poly2->getExteriorRing()))  return false;

    // Convert each hole inside this polygon
    for(unsigned int holeIndex = 0; holeIndex < poly1->getNumInteriorRing(); holeIndex ++) {

      // We hope they are all linear rings (closed/simple), but if not just leave it be
      if(poly1->getInteriorRingN(holeIndex)->getGeometryTypeId() == geos::geom::GEOS_LINESTRING) {
        holes1.push_back(poly1->getInteriorRingN(holeIndex));
      }

      if(poly2->getInteriorRingN(holeIndex)->getGeometryTypeId() == geos::geom::GEOS_LINESTRING) {
        holes2.push_back(poly2->getInteriorRingN(holeIndex));
      }

    }

    if(holes1.size() != holes2.size())  return false;

    for(int h1 = holes1.size() - 1; (h1 >= 0) && holes1.size(); h1 --) {
      for(int h2 = holes2.size() - 1; (h2 >= 0) && holes2.size(); h2 --) {
        if(Equal(holes1[h1], holes2[h2])) {
          // Delete holes1[h1] by replacing it with the last Polygon in holes1
          holes1[h1] = holes1[holes1.size()-1];
          holes1.resize(holes1.size() - 1);
          // Delete holes2[h2] by replacing it with the last Polygon in holes2
          holes2[h2] = holes2[holes2.size()-1];
          holes2.resize(holes2.size() - 1);
        }
      }
    }

    return (holes1.size() == 0) && (holes2.size() == 0);
  }


  bool PolygonTools::Equal(const geos::geom::LineString *lineString1, 
                           const geos::geom::LineString *lineString2) {

    geos::geom::CoordinateSequence *coords1 = lineString1->getCoordinates();
    geos::geom::CoordinateSequence *coords2 = lineString2->getCoordinates();
    bool isEqual = true;

    if(coords1->getSize() != coords2->getSize()) isEqual = false;

    unsigned int index1 = 0;
    unsigned int index2 = 0;

    // -1 extra for dupicate start/end coordinates
    for(; index2 < coords2->getSize() - 1 && isEqual; index2 ++) {
      if(Equal(coords1->getAt(index1), coords2->getAt(index2)))  break;
    }

    if(index2 == coords2->getSize() - 1) isEqual = false;

    for(; index1 < coords1->getSize() - 1 && isEqual; index1 ++, index2 ++) {
      if(!Equal(coords1->getAt(index1), coords2->getAt(index2 % (coords2->getSize() - 1)))) {
        isEqual = false;
      }
    }

    delete coords1;
    delete coords2;
    return isEqual;
  }


  bool PolygonTools::Equal(const geos::geom::Coordinate &coord1, 
                           const geos::geom::Coordinate &coord2) {

    if(!Equal(coord1.x, coord2.x))  return false;
    if(!Equal(coord1.y, coord2.y))  return false;
    if(!Equal(coord1.y, coord2.y))  return false;

    return true;
  }


  bool PolygonTools::Equal(const double d1, const double d2) {
    const double cutoff = 1e15;

    if(DecimalPlace(d1) != DecimalPlace(d2)) return false;

    int decimalPlace = DecimalPlace(d1);
    double factor = pow(10.0, (int)decimalPlace);

    // reduced num is in the form 0.nnnnnnnnnn...
    double reducedNum = d1 / factor;

    double round_offset = (d1 < 0) ? -0.5 : 0.5;

    // cast off the digits past the precision's place
    long long num1 = ((long long)(reducedNum * cutoff + round_offset));

    factor = pow(10.0, (int)decimalPlace);

    // reduced num is in the form 0.nnnnnnnnnn...
    reducedNum = d2 / factor;

    round_offset = (d2 < 0) ? -0.5 : 0.5;

    // cast off the digits past the precision's place
    long long num2 = ((long long)(reducedNum * cutoff + round_offset));


    return (num1 == num2);
  }


  geos::geom::MultiPolygon *PolygonTools::SplitPolygonOn360(const geos::geom::Polygon *inPoly) {
    bool convertLon = false;
    bool negAdjust = false;
    bool newCoords = false;  //  coordinates have been adjusted
    geos::geom::CoordinateSequence *newLonLatPts = new geos::geom::CoordinateArraySequence();
    double lon, lat;
    double lonOffset = 0;
    geos::geom::CoordinateSequence *inPolyCoords = inPoly->getCoordinates();
    double prevLon = inPolyCoords->getAt(0).x;
    double prevLat = inPolyCoords->getAt(0).y;

    newLonLatPts->add(geos::geom::Coordinate(prevLon, prevLat));
    double dist = 0.0;
    for (unsigned int i = 1; i < inPolyCoords->getSize(); i++) {
      lon = inPolyCoords->getAt(i).x;
      lat = inPolyCoords->getAt(i).y;
      // check to see if you just crossed the Meridian
      if (abs(lon - prevLon) > 180 && (prevLat != 90 && prevLat != -90)) {
        newCoords = true;
        // if you were already converting then stop (crossed Meridian even number of times)
        if (convertLon) {
          convertLon = false;
          lonOffset = 0;
        }
        else {   // Need to start converting again, deside how to adjust coordinates
          if ((lon - prevLon) > 0) {
            lonOffset = -360.;
            negAdjust = true;
          }
          else if ((lon - prevLon) < 0) {
            lonOffset = 360.;
            negAdjust = false;
          }
          convertLon = true;
        }
      }

      // Change to a minimum calculation
      if (newCoords  &&  dist == 0.0) {
        double longitude = (lon + lonOffset) - prevLon;
        double latitude = lat - prevLat;
        dist = std::sqrt((longitude * longitude) + (latitude * latitude));
      }

      // add coord
      newLonLatPts->add(geos::geom::Coordinate(lon + lonOffset, lat));

      // set current to old
      prevLon = lon;
      prevLat = lat;
    }

    delete inPolyCoords;

    // Nothing was done so return
    if (!newCoords) {
      geos::geom::Polygon *newPoly = globalFactory.createPolygon
                                     (globalFactory.createLinearRing(newLonLatPts), NULL);
      geos::geom::MultiPolygon *multi_polygon = PolygonTools::MakeMultiPolygon(newPoly);
      delete newLonLatPts;
      return multi_polygon;
    }

    // bisect into seperate polygons
    try {
      geos::geom::Polygon *newPoly = globalFactory.createPolygon
                                     (globalFactory.createLinearRing(newLonLatPts), NULL);

      geos::geom::CoordinateSequence *pts = new geos::geom::CoordinateArraySequence();
      geos::geom::CoordinateSequence *pts2 = new geos::geom::CoordinateArraySequence();

      // Depending on direction of compensation bound accordingly
      //***************************************************

      // please verify correct if you change these values
      //***************************************************
      if (negAdjust) {
        pts->add(geos::geom::Coordinate(0., 90.));
        pts->add(geos::geom::Coordinate(-360., 90.));
        pts->add(geos::geom::Coordinate(-360., -90.));
        pts->add(geos::geom::Coordinate(0., -90.));
        for (double lat = -90.0 + dist; lat < 90.0; lat += dist) {
          pts->add(geos::geom::Coordinate(0.0, lat));
        }
        pts->add(geos::geom::Coordinate(0., 90.));
        pts2->add(geos::geom::Coordinate(0., 90.));
        pts2->add(geos::geom::Coordinate(360., 90.));
        pts2->add(geos::geom::Coordinate(360., -90.));
        pts2->add(geos::geom::Coordinate(0., -90.));
        for (double lat = -90.0 + dist; lat < 90.0; lat += dist) {
          pts2->add(geos::geom::Coordinate(0.0, lat));
        }
        pts2->add(geos::geom::Coordinate(0., 90.));
      }
      else {
        pts->add(geos::geom::Coordinate(360., 90.));
        pts->add(geos::geom::Coordinate(720., 90.));
        pts->add(geos::geom::Coordinate(720., -90.));
        pts->add(geos::geom::Coordinate(360., -90.));
        for (double lat = -90.0 + dist; lat < 90.0; lat += dist) {
          pts->add(geos::geom::Coordinate(360.0, lat));
        }
        pts->add(geos::geom::Coordinate(360., 90.));
        pts2->add(geos::geom::Coordinate(360., 90.));
        pts2->add(geos::geom::Coordinate(0., 90.));
        pts2->add(geos::geom::Coordinate(0., -90.));
        pts2->add(geos::geom::Coordinate(360., -90.));
        for (double lat = -90.0 + dist; lat < 90.0; lat += dist) {
          pts2->add(geos::geom::Coordinate(360.0, lat));
        }
        pts2->add(geos::geom::Coordinate(360., 90.));
      }

      geos::geom::Polygon *boundaryPoly = globalFactory.createPolygon
                                          (globalFactory.createLinearRing(pts), NULL);
      geos::geom::Polygon *boundaryPoly2 = globalFactory.createPolygon
                                           (globalFactory.createLinearRing(pts2), NULL);
      /*------------------------------------------------------------------------
      /  Intersecting the original polygon (converted coordinates) with the
      /  boundary polygons will create the multi polygons with the converted coordinates.
      /  These will need to be converted back to the original coordinates.
      /-----------------------------------------------------------------------*/
      geos::geom::Geometry *intersection = PolygonTools::Intersect(newPoly, boundaryPoly);
      geos::geom::MultiPolygon *convertPoly = PolygonTools::MakeMultiPolygon(intersection);
      delete intersection;

      intersection = PolygonTools::Intersect(newPoly, boundaryPoly2);
      geos::geom::MultiPolygon *convertPoly2 = PolygonTools::MakeMultiPolygon(intersection);
      delete intersection;

      /*------------------------------------------------------------------------
      / Adjust points created in the negative space or >360 space to be back in
      / the 0-360 world.  This will always only need to be done on convertPoly.
      / Then add geometries to finalpolys.
      /-----------------------------------------------------------------------*/
      vector<geos::geom::Geometry *> *finalpolys = new vector<geos::geom::Geometry *>;
      geos::geom::Geometry *newGeom = NULL;

      for (unsigned int i = 0; i < convertPoly->getNumGeometries(); i++) {
        newGeom = (convertPoly->getGeometryN(i))->clone();
        pts = convertPoly->getGeometryN(i)->getCoordinates();
        geos::geom::CoordinateSequence *newLonLatPts = new geos::geom::CoordinateArraySequence();
        // fix the points

        for (unsigned int k = 0; k < pts->getSize() ; k++) {
          double lon = pts->getAt(k).x;
          double lat = pts->getAt(k).y;
          if (negAdjust) {
            lon = lon + 360;
          }
          else {
            lon = lon - 360;
          }
          newLonLatPts->add(geos::geom::Coordinate(lon, lat), k);

        }
        // Add the points to polys
        finalpolys->push_back(globalFactory.createPolygon
                              (globalFactory.createLinearRing(newLonLatPts), NULL));
      }

      // This loop is over polygons that will always be in 0-360 space no need to convert
      for (unsigned int i = 0; i < convertPoly2->getNumGeometries(); i++) {
        newGeom = (convertPoly2->getGeometryN(i))->clone();
        finalpolys->push_back(newGeom);
      }

      geos::geom::MultiPolygon *multi_polygon = globalFactory.createMultiPolygon(*finalpolys);

      delete finalpolys;
      delete newGeom;
      delete newLonLatPts;
      delete pts;
      delete pts2;
      return multi_polygon;
    }
    catch(geos::util::IllegalArgumentException *geosIll) {
      std::string msg = "Unable to split polygon on longitude boundry (SplitPolygonOn360) due to ";
      msg += "geos illegal argument [" + IString(geosIll->what()) + "]";
      delete geosIll;
      throw IException(IException::Unknown, msg, _FILEINFO_);
    }
    catch(geos::util::GEOSException *geosExc) {
      std::string msg = "Unable to split polygon on longitude boundry (SplitPolygonOn360) due to ";
      msg += "geos exception [" + IString(geosExc->what()) + "]";
      delete geosExc;
      throw IException(IException::Unknown, msg, _FILEINFO_);
    }
    catch(IException &e) {
      std::string msg = "Unable to split polygon on longitude boundry (SplitPolygonOn360) due to ";
      msg += "isis operation exception [" + IString(e.what()) + "]";
      throw IException(e, IException::Unknown, msg, _FILEINFO_);
    }
    catch(...) {
      std::string msg = "Unable to split polygon on longitude boundry (SplitPolygonOn360) due to ";
      msg += "unknown exception";
      throw IException(IException::Unknown, msg, _FILEINFO_);
    }
  }
} // end namespace isis