USGS

Isis 3.0 Object Programmers' Reference

 Home
PolygonTools.cpp
Go to the documentation of this file.
1 
25 #include <string>
26 #include <iostream>
27 #include <sstream>
28 #include <iomanip>
29 #include <vector>
30 #include <cmath>
31 
32 #include <QDebug>
33 #include "geos/geom/BinaryOp.h"
34 #include "geos/geom/CoordinateArraySequence.h"
35 #include "geos/geom/CoordinateSequence.h"
36 #include "geos/geom/LinearRing.h"
37 #include "geos/geom/Point.h"
38 #include "geos/geom/Polygon.h"
39 #include "geos/operation/distance/DistanceOp.h"
40 #include "geos/opOverlay.h"
41 #include "geos/operation/overlay/snap/GeometrySnapper.h"
42 
43 #include "SpecialPixel.h"
44 #include "PolygonTools.h"
45 #include "TProjection.h"
46 #include "ProjectionFactory.h"
47 #include "UniversalGroundMap.h"
48 
49 using namespace std;
50 namespace Isis {
63  geos::geom::MultiPolygon *PolygonTools::LatLonToXY(
64  const geos::geom::MultiPolygon &lonLatPolygon, TProjection *projection) {
65  if (projection == NULL) {
66  string msg = "Unable to convert Lon/Lat polygon to X/Y. ";
67  msg += "No projection has was supplied";
68  throw IException(IException::Programmer, msg, _FILEINFO_);
69  }
70 
71  // Convert the Lat/Lon poly coordinates to X/Y coordinates
72  if(lonLatPolygon.isEmpty()) {
73  return globalFactory.createMultiPolygon();
74  }
75  else {
76  vector<geos::geom::Geometry *> *xyPolys = new vector<geos::geom::Geometry *>;
77  // Convert each polygon in this multi-polygon
78  for(unsigned int g = 0; g < lonLatPolygon.getNumGeometries(); ++g) {
79  const geos::geom::Polygon *poly =
80  dynamic_cast<const geos::geom::Polygon *>(
81  lonLatPolygon.getGeometryN(g));
82 
83  // Convert each hole inside this polygon
84  vector<geos::geom::Geometry *> *holes = new vector<geos::geom::Geometry *>;
85  for(unsigned int h = 0; h < poly->getNumInteriorRing(); ++h) {
86  geos::geom::CoordinateSequence *xycoords = new geos::geom::CoordinateArraySequence();
87  geos::geom::CoordinateSequence *llcoords =
88  poly->getInteriorRingN(h)->getCoordinates();
89 
90  // Convert each coordinate in this hole
91  for(unsigned int cord = 0; cord < llcoords->getSize(); ++cord) {
92  projection->SetGround(llcoords->getAt(cord).y,
93  llcoords->getAt(cord).x);
94  xycoords->add(geos::geom::Coordinate(projection->XCoord(),
95  projection->YCoord()));
96  } // end num coords in hole loop
97 
98  geos::geom::LinearRing *hole = globalFactory.createLinearRing(xycoords);
99 
100  if(hole->isValid() && !hole->isEmpty()) {
101  holes->push_back(hole);
102  }
103  else {
104  delete hole;
105  }
106  } // end num holes in polygon loop
107 
108  // Convert the exterior ring of this polygon
109  geos::geom::CoordinateSequence *xycoords = new geos::geom::CoordinateArraySequence();
110  geos::geom::CoordinateSequence *llcoords =
111  poly->getExteriorRing()->getCoordinates();
112 
113  // Convert each coordinate in the exterior ring of this polygon
114  for (unsigned int cord = 0; cord < llcoords->getSize(); ++cord) {
115  projection->SetGround(llcoords->getAt(cord).y,
116  llcoords->getAt(cord).x);
117  xycoords->add(geos::geom::Coordinate(projection->XCoord(),
118  projection->YCoord()));
119  } // end exterior ring coordinate loop
120 
121  geos::geom::Polygon *newPoly = globalFactory.createPolygon(
122  globalFactory.createLinearRing(xycoords), holes);
123 
124  if(newPoly->isValid() && !newPoly->isEmpty() && newPoly->getArea() > 1.0e-14) {
125  xyPolys->push_back(newPoly);
126  }
127  else {
128  delete newPoly;
129  }
130  } // end num geometry in multi-poly
131 
132  // Create a new multipoly from all the new X/Y polygon(s)
133  geos::geom::MultiPolygon *spikedPoly = globalFactory.createMultiPolygon(xyPolys);
134 
135  if(spikedPoly->isValid() && !spikedPoly->isEmpty()) {
136  return spikedPoly;
137  }
138  else {
139  try {
140  geos::geom::MultiPolygon *despikedPoly = Despike(spikedPoly);
141 
142  delete spikedPoly;
143  spikedPoly = NULL;
144 
145  return despikedPoly;
146  }
147  catch(IException &e) {
148  IString msg = "Unable to convert polygon from Lat/Lon to X/Y";
149  throw IException(IException::Programmer, msg, _FILEINFO_);
150  }
151  }
152 
153  } // end else
154  }
155 
156 
169  geos::geom::MultiPolygon *PolygonTools::XYToLatLon(
170  const geos::geom::MultiPolygon &xYPolygon, TProjection *projection) {
171 
172  if(projection == NULL) {
173  string msg = "Unable to convert X/Y polygon to Lon/Lat. ";
174  msg += "No projection was supplied";
175  throw IException(IException::Programmer, msg, _FILEINFO_);
176  }
177 
178  // Convert the X/Y poly coordinates to Lat/Lon coordinates
179  if(xYPolygon.isEmpty()) {
180  return globalFactory.createMultiPolygon();
181  }
182  else {
183  vector<geos::geom::Geometry *> *llPolys = new vector<geos::geom::Geometry *>;
184  // Convert each polygon in this multi-polygon
185  for(unsigned int g = 0; g < xYPolygon.getNumGeometries(); ++g) {
186  const geos::geom::Polygon *poly =
187  dynamic_cast<const geos::geom::Polygon *>(
188  xYPolygon.getGeometryN(g));
189 
190  // Convert each hole inside this polygon
191  vector<geos::geom::Geometry *> *holes = new vector<geos::geom::Geometry *>;
192  for(unsigned int h = 0; h < poly->getNumInteriorRing(); ++h) {
193  geos::geom::CoordinateSequence *llcoords = new geos::geom::CoordinateArraySequence();
194  geos::geom::CoordinateSequence *xycoords =
195  poly->getInteriorRingN(h)->getCoordinates();
196 
197  // Convert each coordinate in this hole
198  for(unsigned int cord = 0; cord < xycoords->getSize(); ++cord) {
199  projection->SetWorld(xycoords->getAt(cord).x,
200  xycoords->getAt(cord).y);
201  llcoords->add(geos::geom::Coordinate(projection->Longitude(),
202  projection->Latitude()));
203  } // end num coords in hole loop
204  holes->push_back(globalFactory.createLinearRing(llcoords));
205  } // end num holes in polygon loop
206 
207  // Convert the exterior ring of this polygon
208  geos::geom::CoordinateSequence *llcoords = new geos::geom::DefaultCoordinateSequence();
209  geos::geom::CoordinateSequence *xycoords =
210  poly->getExteriorRing()->getCoordinates();
211 
212  // Convert each coordinate in the exterior ring of this polygon
213  for(unsigned int cord = 0; cord < xycoords->getSize(); ++cord) {
214  projection->SetWorld(xycoords->getAt(cord).x,
215  xycoords->getAt(cord).y);
216  llcoords->add(geos::geom::Coordinate(projection->Longitude(),
217  projection->Latitude()));
218  } // end exterior ring coordinate loop
219 
220  llPolys->push_back(globalFactory.createPolygon(
221  globalFactory.createLinearRing(llcoords), holes));
222  } // end num geometry in multi-poly
223 
224 
225  // Create a new multipoly from all the new Lat/Lon polygon(s)
226  geos::geom::MultiPolygon *spikedPoly = globalFactory.createMultiPolygon(llPolys);
227 
228  if(spikedPoly->isValid() && !spikedPoly->isEmpty()) {
229  return spikedPoly;
230  }
231  else {
232  try {
233  geos::geom::MultiPolygon *despikedPoly = Despike(spikedPoly);
234 
235  delete spikedPoly;
236  spikedPoly = NULL;
237 
238  return despikedPoly;
239  }
240  catch(IException &e) {
241  IString msg = "Unable to convert polygon from X/Y to Lat/Lon";
242  throw IException(IException::Programmer, msg, _FILEINFO_);
243  }
244  }
245  } // end else
246  }
247 
248 
261  geos::geom::MultiPolygon *PolygonTools::LatLonToSampleLine(
262  const geos::geom::MultiPolygon &lonLatPolygon, UniversalGroundMap *ugm) {
263 
264  if(ugm == NULL) {
265  string msg = "Unable to convert Lon/Lat polygon to Sample/Line. ";
266  msg += "No UniversalGroundMap was supplied";
267  throw IException(IException::Programmer, msg, _FILEINFO_);
268  }
269 
270  // Convert the Lon/Lat poly coordinates to Sample/Line coordinates
271  if(lonLatPolygon.isEmpty()) {
272  return globalFactory.createMultiPolygon();
273  }
274  else {
275  vector<geos::geom::Geometry *> *slPolys = new vector<geos::geom::Geometry *>;
276  // Convert each polygon in this multi-polygon
277  for(unsigned int g = 0; g < lonLatPolygon.getNumGeometries(); g++) {
278  const geos::geom::Polygon *poly =
279  dynamic_cast<const geos::geom::Polygon *>(
280  lonLatPolygon.getGeometryN(g));
281 
282  // Convert each hole inside this polygon
283  vector<geos::geom::Geometry *> *holes = new vector<geos::geom::Geometry *>;
284  for(unsigned int h = 0; h < poly->getNumInteriorRing(); ++h) {
285  geos::geom::CoordinateSequence *slcoords = new geos::geom::DefaultCoordinateSequence();
286  geos::geom::CoordinateSequence *llcoords =
287  poly->getInteriorRingN(h)->getCoordinates();
288 
289  // Convert each coordinate in this hole
290  for(unsigned int cord = 0; cord < llcoords->getSize(); ++cord) {
291  ugm->SetUniversalGround(llcoords->getAt(cord).y,
292  llcoords->getAt(cord).x);
293  slcoords->add(geos::geom::Coordinate(ugm->Sample(),
294  ugm->Line()));
295  } // end num coords in hole loop
296  holes->push_back(globalFactory.createLinearRing(slcoords));
297  delete slcoords;
298  delete llcoords;
299  } // end num holes in polygon loop
300 
301  // Convert the exterior ring of this polygon
302  geos::geom::CoordinateSequence *slcoords = new geos::geom::CoordinateArraySequence();
303  geos::geom::CoordinateSequence *llcoords =
304  poly->getExteriorRing()->getCoordinates();
305 
306  // Convert each coordinate in the exterior ring of this polygon
307  for(unsigned int cord = 0; cord < llcoords->getSize(); ++cord) {
308  if (ugm->SetUniversalGround(llcoords->getAt(cord).y,
309  llcoords->getAt(cord).x)) {
310  slcoords->add(geos::geom::Coordinate(ugm->Sample(),
311  ugm->Line()));
312  }
313  } // end exterior ring coordinate loop
314 
315  // Make sure that the line string is closed.
316  if (slcoords->getSize() > 0 && !slcoords->front().equals(slcoords->back())) {
317  slcoords->add(slcoords->front());
318  }
319 
320  try {
321  slPolys->push_back(globalFactory.createPolygon(
322  globalFactory.createLinearRing(slcoords), holes));
323  }
324  catch (std::exception &e) {
325  throw IException(IException::Unknown,
326  QObject::tr("Unable to convert polygon from Lat/Lon to Sample/Line. The "
327  "error given was [%1].").arg(e.what()),
328  _FILEINFO_);
329  }
330 
331  delete llcoords;
332  } // end num geometry in multi-poly
333 
334  // Create a new multipoly from all the new Sample/Line polygon(s)
335  geos::geom::MultiPolygon *spikedPoly = globalFactory.createMultiPolygon(slPolys);
336 
337  if(spikedPoly->isValid() && !spikedPoly->isEmpty()) {
338  return spikedPoly;
339  }
340  else {
341  try {
342  geos::geom::MultiPolygon *despikedPoly = Despike(spikedPoly);
343 
344  delete spikedPoly;
345  spikedPoly = NULL;
346 
347  return despikedPoly;
348  }
349  catch (IException &e) {
350  IString msg = "Unable to convert polygon from Lat/Lon to Sample/Line";
351  throw IException(IException::Programmer, msg, _FILEINFO_);
352  }
353  }
354  } // end else
355  }
356 
357 
370  geos::geom::MultiPolygon *PolygonTools::CopyMultiPolygon(const geos::geom::MultiPolygon *mpolygon) {
371 
372  vector<geos::geom::Geometry *> *polys = new vector<geos::geom::Geometry *>;
373  for(unsigned int i = 0; i < mpolygon->getNumGeometries(); ++i) {
374  polys->push_back((mpolygon->getGeometryN(i))->clone());
375  }
376  return globalFactory.createMultiPolygon(polys);
377  }
378 
379 
392  geos::geom::MultiPolygon *PolygonTools::CopyMultiPolygon(const geos::geom::MultiPolygon &mpolygon) {
393 
394  vector<geos::geom::Geometry *> *polys = new vector<geos::geom::Geometry *>;
395  for(unsigned int i = 0; i < mpolygon.getNumGeometries(); ++i) {
396  polys->push_back((mpolygon.getGeometryN(i))->clone());
397  }
398  return globalFactory.createMultiPolygon(polys);
399  }
400 
401 
410  QString PolygonTools::ToGML(const geos::geom::MultiPolygon *mpolygon, QString idString,
411  QString schema) {
412 
413  ostringstream os;
414 
415  //Write out the GML header
416  os << "<?xml version=\"1.0\" encoding=\"utf-8\" ?>" << endl;
417  os << "<ogr:FeatureCollection" << endl;
418  os << " xmlns:xsi=\"http://www.w3.org/2001/XMLSchema-instance\"" << endl;
419  os << " xsi:schemaLocation=\"http://ogr.maptools.org/ " << schema << "\"" << endl;
420  os << " xmlns:ogr=\"http://ogr.maptools.org/\"" << endl;
421  os << " xmlns:gml=\"http://www.opengis.net/gml\">" << endl;
422  os << " <gml:boundedBy>" << endl;
423  os << " <gml:Box>" << endl;
424  os << " <gml:coord>";
425  os << "<gml:X>" <<
426  setprecision(15) << mpolygon->getEnvelopeInternal()->getMinX() <<
427  "</gml:X>";
428  os << "<gml:Y>" <<
429  setprecision(15) << mpolygon->getEnvelopeInternal()->getMinY() <<
430  "</gml:Y>";
431  os << "</gml:coord>" << endl;
432  os << " <gml:coord>";
433  os << "<gml:X>" <<
434  setprecision(15) << mpolygon->getEnvelopeInternal()->getMaxX() <<
435  "</gml:X>";
436  os << "<gml:Y>" <<
437  setprecision(15) << mpolygon->getEnvelopeInternal()->getMaxY() <<
438  "</gml:Y>";
439  os << "</gml:coord>" << endl;
440  os << " </gml:Box>" << endl;
441  os << " </gml:boundedBy>" << endl << endl;
442  os << " <gml:featureMember>" << endl;
443  os << " <ogr:multi_polygon fid=\"0\">" << endl;
444  os << " <ogr:ID>" << idString << "</ogr:ID>" << endl;
445  os << " <ogr:geometryProperty><gml:MultiPolygon><gml:polygonMember>" <<
446  "<gml:Polygon><gml:outerBoundaryIs><gml:LinearRing><gml:coordinates>";
447 
448 
449  for(unsigned int polyN = 0; polyN < mpolygon->getNumGeometries(); polyN++) {
450  geos::geom::CoordinateSequence *pts = mpolygon->getGeometryN(polyN)->getCoordinates();
451 
452  for(unsigned int i = 0; i < pts->getSize(); i++) {
453  double lon = pts->getAt(i).x;
454  double lat = pts->getAt(i).y;
455 
456  os << setprecision(15) << lon << "," << setprecision(15) << lat << " ";
457  }
458  }
459 
460  os << "</gml:coordinates></gml:LinearRing></gml:outerBoundaryIs>" <<
461  "</gml:Polygon></gml:polygonMember></gml:MultiPolygon>" <<
462  "</ogr:geometryProperty>" << endl;
463  os << " </ogr:multi_polygon>" << endl;
464  os << " </gml:featureMember>" << endl;
465  os << "</ogr:FeatureCollection>";
466 
467  return os.str().c_str();
468  }
469 
470 
477  QString PolygonTools::GMLSchema() {
478 
479  ostringstream os;
480 
481  os << "<?xml version=\"1.0\" encoding=\"UTF-8\"?>" << endl;
482  os << "<xs:schema targetNamespace=\"http://ogr.maptools.org/\" "
483  "xmlns:ogr=\"http://ogr.maptools.org/\" "
484  "xmlns:xs=\"http://www.w3.org/2001/XMLSchema\" "
485  "xmlns:gml=\"http://www.opengis.net/gml\" "
486  "elementFormDefault=\"qualified\" "
487  "version=\"1.0\">" << endl;
488  os << " <xs:import namespace=\"http://www.opengis.net/gml\" "
489  "schemaLocation=\"http://schemas.opengis.net/gml/2.1.2/feature.xsd\"/>" << endl;
490  os << " <xs:element name=\"FeatureCollection\" "
491  "type=\"ogr:FeatureCollectionType\" "
492  "substitutionGroup=\"gml:_FeatureCollection\"/>" << endl;
493  os << " <xs:complexType name=\"FeatureCollectionType\">" << endl;
494  os << " <xs:complexContent>" << endl;
495  os << " <xs:extension base=\"gml:AbstractFeatureCollectionType\">" << endl;
496  os << " <xs:attribute name=\"lockId\" type=\"xs:string\" use=\"optional\"/>" << endl;
497  os << " <xs:attribute name=\"scope\" type=\"xs:string\" use=\"optional\"/>" << endl;
498  os << " </xs:extension>" << endl;
499  os << " </xs:complexContent>" << endl;
500  os << " </xs:complexType>" << endl;
501  os << " <xs:element name=\"multi_polygon\" "
502  "type=\"ogr:multi_polygon_Type\" "
503  "substitutionGroup=\"gml:_Feature\"/>" << endl;
504  os << " <xs:complexType name=\"multi_polygon_Type\">" << endl;
505  os << " <xs:complexContent>" << endl;
506  os << " <xs:extension base=\"gml:AbstractFeatureType\">" << endl;
507  os << " <xs:sequence>" << endl;
508  os << " <xs:element name=\"geometryProperty\" "
509  "type=\"gml:MultiPolygonPropertyType\" "
510  "nillable=\"true\" minOccurs=\"0\" maxOccurs=\"1\"/>" << endl;
511  os << " <xs:element name=\"fid\" nillable=\"true\" minOccurs=\"0\" "
512  "maxOccurs=\"1\">" << endl;
513  os << " <xs:simpleType>" << endl;
514  os << " <xs:restriction base=\"xs:string\">" << endl;
515  os << " </xs:restriction>" << endl;
516  os << " </xs:simpleType>" << endl;
517  os << " </xs:element>" << endl;
518  os << " <xs:element name=\"ID\" nillable=\"true\" minOccurs=\"0\" "
519  "maxOccurs=\"1\">" << endl;
520  os << " <xs:simpleType>" << endl;
521  os << " <xs:restriction base=\"xs:integer\">" << endl;
522  os << " <xs:totalDigits value=\"16\"/>" << endl;
523  os << " </xs:restriction>" << endl;
524  os << " </xs:simpleType>" << endl;
525  os << " </xs:element>" << endl;
526  os << " </xs:sequence>" << endl;
527  os << " </xs:extension>" << endl;
528  os << " </xs:complexContent>" << endl;
529  os << " </xs:complexType>" << endl;
530  os << "</xs:schema>";
531 
532  return os.str().c_str();
533  }
534 
535 
543  geos::geom::MultiPolygon *PolygonTools::To180(geos::geom::MultiPolygon *poly360) {
544  try {
545  // Let's take the 360 pieces that sit between 180 and 360 and move them
546  // to -180 to 0. To accomplish this, make a poly that fits 180 -> 360
547  // degrees longitude and intersect (this gives us the pieces that sit
548  // >180). Move this intersection to the left. Then make a poly that fits
549  // 0 to 180 and intersect with the original. These two combined are the
550  // result.
551  geos::geom::CoordinateArraySequence *leftOf180Pts =
552  new geos::geom::CoordinateArraySequence();
553  leftOf180Pts->add(geos::geom::Coordinate(0, -90));
554  leftOf180Pts->add(geos::geom::Coordinate(0, 90));
555  leftOf180Pts->add(geos::geom::Coordinate(180, 90));
556  leftOf180Pts->add(geos::geom::Coordinate(180, -90));
557  leftOf180Pts->add(geos::geom::Coordinate(0, -90));
558 
559  geos::geom::LinearRing *leftOf180Geom =
560  globalFactory.createLinearRing(leftOf180Pts);
561 
562  geos::geom::Polygon *leftOf180Poly =
563  globalFactory.createPolygon(leftOf180Geom, NULL);
564 
565  geos::geom::CoordinateArraySequence *rightOf180Pts =
566  new geos::geom::CoordinateArraySequence();
567  rightOf180Pts->add(geos::geom::Coordinate(180, -90));
568  rightOf180Pts->add(geos::geom::Coordinate(180, 90));
569  rightOf180Pts->add(geos::geom::Coordinate(360, 90));
570  rightOf180Pts->add(geos::geom::Coordinate(360, -90));
571  rightOf180Pts->add(geos::geom::Coordinate(180, -90));
572 
573  geos::geom::LinearRing *rightOf180Geom =
574  globalFactory.createLinearRing(rightOf180Pts);
575 
576  geos::geom::Polygon *rightOf180Poly =
577  globalFactory.createPolygon(rightOf180Geom, NULL);
578 
579  geos::geom::Geometry *preserved = Intersect(leftOf180Poly, poly360);
580  geos::geom::Geometry *moving = Intersect(rightOf180Poly, poly360);
581 
582  geos::geom::CoordinateSequence *movingPts = moving->getCoordinates();
583  geos::geom::CoordinateSequence *movedPts =
584  new geos::geom::CoordinateArraySequence();
585 
586  for(unsigned int i = 0; i < movingPts->getSize(); i ++) {
587  movedPts->add(geos::geom::Coordinate(movingPts->getAt(i).x - 360.0,
588  movingPts->getAt(i).y));
589  }
590 
591  if(movedPts->getSize()) {
592  movedPts->add(geos::geom::Coordinate(movedPts->getAt(0).x,
593  movedPts->getAt(0).y));
594  }
595 
596  geos::geom::Geometry *moved = globalFactory.createPolygon(
597  globalFactory.createLinearRing(movedPts), NULL);
598 
599  std::vector<geos::geom::Geometry *> *geomsForCollection = new
600  std::vector<geos::geom::Geometry *>;
601  geomsForCollection->push_back(preserved);
602  geomsForCollection->push_back(moved);
603 
604  geos::geom::GeometryCollection *the180Polys =
605  Isis::globalFactory.createGeometryCollection(geomsForCollection);
606 
607  geos::geom::MultiPolygon *result = MakeMultiPolygon(the180Polys);
608  delete the180Polys;
609  the180Polys = NULL;
610 
611  geos::geom::MultiPolygon *fixedResult = FixSeam(result);
612  delete result;
613  result = NULL;
614 
615  return fixedResult;
616  }
617  catch(geos::util::GEOSException *exc) {
618  IString msg = "Conversion to 180 failed. The reason given was [" +
619  IString(exc->what()) + "]";
620  delete exc;
621  throw IException(IException::Programmer, msg, _FILEINFO_);
622  }
623  catch(...) {
624  IString msg = "Conversion to 180 failed. Could not determine the reason";
625  throw IException(IException::Programmer, msg, _FILEINFO_);
626  }
627  }
628 
629 
640  double PolygonTools::Thickness(const geos::geom::MultiPolygon *mpolygon) {
641  const geos::geom::Envelope *envelope = mpolygon->getEnvelopeInternal();
642 
643  double x = fabs(envelope->getMaxX() - envelope->getMinX());
644  double y = fabs(envelope->getMaxY() - envelope->getMinY());
645  double extent = max(x, y);
646 
647  return mpolygon->getArea() / (extent * extent);
648  }
649 
650 
662  geos::geom::MultiPolygon *PolygonTools::Despike(const geos::geom::Geometry *geom) {
663  geos::geom::MultiPolygon *spiked = MakeMultiPolygon(geom);
664  geos::geom::MultiPolygon *despiked = Despike(spiked);
665 
666  delete spiked;
667  spiked = NULL;
668 
669  return despiked;
670  }
671 
672 
683  geos::geom::MultiPolygon *PolygonTools::Despike(const geos::geom::MultiPolygon *multiPoly) {
684  // Despike each polygon in the multipolygon
685  vector<geos::geom::Geometry *> *newPolys = new vector<geos::geom::Geometry *>;
686  for(unsigned int g = 0; g < multiPoly->getNumGeometries(); ++g) {
687  const geos::geom::Polygon *poly =
688  dynamic_cast<const geos::geom::Polygon *>(multiPoly->getGeometryN(g));
689 
690  // Despike each hole inside this polygon
691  vector<geos::geom::Geometry *> *holes = new vector<geos::geom::Geometry *>;
692  for(unsigned int h = 0; h < poly->getNumInteriorRing(); ++h) {
693  const geos::geom::LineString *ls = poly->getInteriorRingN(h);
694  geos::geom::LinearRing *lr;
695 
696  // If the hole is not valid fix it
697  // If the hole is NOT valid despike it
698  lr = Despike(ls);
699 
700  if(!lr->isValid()) {
701  geos::geom::LinearRing *fixed = FixGeometry(lr);
702  delete lr;
703  lr = fixed;
704  }
705 
706  // Save this hole if it is not too small
707  if(!lr->isEmpty()) {
708  holes->push_back(lr);
709  }
710  else {
711  delete lr;
712  }
713  } // End holes loop
714 
715  // Work on the main polygon
716  const geos::geom::LineString *ls = poly->getExteriorRing();
717  geos::geom::LinearRing *lr;
718 
719  lr = Despike(ls);
720 
721  try {
722  if(!lr->isValid()) {
723  geos::geom::LinearRing *fixed = FixGeometry(lr);
724  delete lr;
725  lr = fixed;
726  }
727  }
728  catch(IException &e) {
729  // Sometimes despike and fix fail, but the input is really valid. We can just go
730  // with the non-despiked polygon.
731  if(ls->isValid() && ls->getGeometryTypeId() == geos::geom::GEOS_LINEARRING) {
732  lr = dynamic_cast<geos::geom::LinearRing *>(ls->clone());
733  }
734  else {
735  throw;
736  }
737  }
738 
739  // Create a new polygon with the holes and save it
740  if(!lr->isEmpty()) {
741  geos::geom::Polygon *tp = Isis::globalFactory.createPolygon(lr, holes);
742 
743  if(tp->isEmpty() || !tp->isValid()) {
744  delete tp;
745  newPolys->push_back(poly->clone());
746  }
747  else {
748  newPolys->push_back(tp);
749  }
750  }
751  } // End polygons loop
752 
753  // Create a new multipoly from the polygon(s)
754  geos::geom::MultiPolygon *mp = Isis::globalFactory.createMultiPolygon(newPolys);
755 
756  if(!mp->isValid() || mp->isEmpty()) {
757  delete mp;
758  mp = NULL;
759  IString msg = "Despike failed to correct the polygon";
760  throw IException(IException::Programmer, msg, _FILEINFO_);
761  }
762 
763  // if multipoly area changes more than 25% we did something bad to the multipolygon
764  if(fabs((mp->getArea() / multiPoly->getArea()) - 1.0) > 0.50) {
765  IString msg = "Despike failed to correct the polygon " + mp->toString();
766  delete mp;
767 
768  throw IException(IException::Programmer, msg, _FILEINFO_);
769  }
770 
771  return mp;
772  }
773 
774 
788  geos::geom::LinearRing *PolygonTools::Despike(const geos::geom::LineString *lineString) {
789  geos::geom::CoordinateSequence *vertices = lineString->getCoordinates();
790 
791  // We need a full polygon to despike = 3 points (+1 for end==first) = at least 4 points
792  if(vertices->getSize() < 4) {
793  delete vertices;
794  return Isis::globalFactory.createLinearRing(geos::geom::CoordinateArraySequence());
795  }
796 
797  // delete one of the duplicate first/end coordinates,
798  // spikes can occur here and the duplicate points throws off the test
799  vertices->deleteAt(vertices->getSize() - 1);
800 
801  // Index will become negative if our first points are spiked, so we need
802  // to make it big enough to encapsulate the (unsigned int) from getSize()
803  // and allow negative numbers.
804  for(long index = 0; index < (long)vertices->getSize(); index++) {
805  // If we're under 3 vertices, we've despiked all the points out of the polygon :(
806  if(vertices->getSize() < 3) {
807  break;
808  }
809 
810  // These are the array indices that we're going to test for a spike,
811  // the middle one is the one that's spiked if IsSpiked(...) is true.
812  long testCoords[3] = {
813  index - 1,
814  index,
815  index + 1
816  };
817 
818  // Make sure the index is inside of our coordinates array (we'll have both too small/too
819  // big of an index)
820  for(int j = 0; j < 3; j++) {
821  while(testCoords[j] < 0) {
822  testCoords[j] += vertices->getSize();
823  }
824 
825  while(testCoords[j] >= (long)vertices->getSize()) {
826  testCoords[j] -= vertices->getSize();
827  }
828  }
829 
830  // Test the middle point for a spike
831  if(IsSpiked(vertices->getAt(testCoords[0]),
832  vertices->getAt(testCoords[1]),
833  vertices->getAt(testCoords[2]))) {
834  // It's spiked, delete it
835  vertices->deleteAt(testCoords[1]);
836 
837  // Back up to the first test that is affected by this change
838  index -= 2;
839  }
840  }
841 
842  if(vertices->getSize() < 3) {
843  delete vertices;
844  vertices = NULL;
845 
846  return Isis::globalFactory.createLinearRing(geos::geom::CoordinateArraySequence());
847  }
848  else {
849  // Duplicate the first vertex as the last to close the polygon
850  vertices->add(vertices->getAt(0));
851  return Isis::globalFactory.createLinearRing(vertices);
852  }
853  }
854 
864  bool PolygonTools::IsSpiked(geos::geom::Coordinate first,
865  geos::geom::Coordinate middle,
866  geos::geom::Coordinate last) {
867  return TestSpiked(first, middle, last) || TestSpiked(last, middle, first);
868  }
869 
898  bool PolygonTools::TestSpiked(geos::geom::Coordinate first, geos::geom::Coordinate middle,
899  geos::geom::Coordinate last) {
900  geos::geom::Point *firstPt = Isis::globalFactory.createPoint(first);
901  geos::geom::Point *middlePt = Isis::globalFactory.createPoint(middle);
902  geos::geom::Point *lastPt = Isis::globalFactory.createPoint(last);
903 
904  geos::geom::CoordinateSequence *coords = new geos::geom::CoordinateArraySequence();
905  coords->add(first);
906  coords->add(middle);
907  geos::geom::LineString *line = Isis::globalFactory.createLineString(coords); // line takes ownership
908 
909  // The lower the tolerance, the less this algorithm removes and thus
910  // the better chance of success in findimageoverlaps. However, if you
911  // lower the tolerance then there is also a greater chance of programs
912  // such as autoseed failing. 1% is the current tolerance.
913  double tolerance = line->getLength() / 100.0;
914 
915  bool spiked = true;
916 
917  double distanceLastMiddle = geos::operation::distance::DistanceOp::distance(lastPt, middlePt);
918  double distanceLastLine = geos::operation::distance::DistanceOp::distance(lastPt, line);
919 
920  if(distanceLastMiddle == 0.0) return true; // get rid of same point
921 
922  // Checks the ratio of the distance between the last point and the line, and the last point
923  // and the middle point if the ratio is very small then there is a spike
924  if(distanceLastLine / distanceLastMiddle >= .05) {
925  spiked = false;
926  }
927 
928  // If the point is away from the line, keep it
929  if(spiked && distanceLastLine > tolerance) {
930  spiked = false;
931  }
932 
933  if(!spiked) {
934  geos::geom::CoordinateSequence *coords = new geos::geom::CoordinateArraySequence();
935  coords->add(first);
936  coords->add(middle);
937  coords->add(last);
938  coords->add(first);
939 
940  // shell takes ownership of coords
941  geos::geom::LinearRing *shell = Isis::globalFactory.createLinearRing(coords);
942  std::vector<geos::geom::Geometry *> *empty = new std::vector<geos::geom::Geometry *>;
943 
944  // poly takes ownership of shell and empty
945  geos::geom::Polygon *poly = Isis::globalFactory.createPolygon(shell, empty);
946 
947 
948  // if these 3 points define a straight line then the middle is worthless (defines nothing)
949  // or problematic
950  if(poly->getArea() < 1.0e-10) {
951  spiked = true;
952  }
953 
954  delete poly;
955  }
956 
957 
958  delete firstPt;
959  delete middlePt;
960  delete lastPt;
961  delete line;
962 
963  return spiked;
964  }
965 
966 
980  geos::geom::Geometry *PolygonTools::Intersect(const geos::geom::Geometry *geom1,
981  const geos::geom::Geometry *geom2) {
982  try {
983  return Operate(geom1, geom2, (unsigned int)geos::operation::overlay::OverlayOp::opINTERSECTION);
984  }
985  catch(geos::util::GEOSException *exc) {
986  IString msg = "Intersect operation failed. The reason given was [" + IString(exc->what()) + "]";
987  delete exc;
988  throw IException(IException::Programmer, msg, _FILEINFO_);
989  }
990  catch(IException &e) {
991  IString msg = "Intersect operation failed";
992  throw IException(e, IException::Programmer, msg, _FILEINFO_);
993  }
994  catch(...) {
995  IString msg = "Intersect operation failed for an unknown reason";
996  throw IException(IException::Programmer, msg, _FILEINFO_);
997  }
998  }
999 
1000 
1001  geos::geom::Geometry *PolygonTools::Operate(const geos::geom::Geometry *geom1,
1002  const geos::geom::Geometry *geom2,
1003  unsigned int opcode) {
1004 
1005  geos::operation::overlay::OverlayOp::OpCode code =
1006  (geos::operation::overlay::OverlayOp::OpCode)opcode;
1007 
1008  geos::geom::Geometry *result = NULL;
1009  bool failed = true;
1010  geos::geom::Geometry *geomFirst = MakeMultiPolygon(geom1);
1011  geos::geom::Geometry *geomSecond = MakeMultiPolygon(geom2);
1012 
1013  geos::operation::overlay::snap::GeometrySnapper snap(*geomFirst);
1014  geos::geom::Geometry *geomSnapped = snap.snapTo(*geomSecond, 1.0e-10)->clone();
1015  if(!geomSnapped->isValid()) {
1016  delete geomSnapped;
1017  }
1018  else {
1019  delete geomFirst;
1020  geomFirst = geomSnapped;
1021  }
1022 
1023  unsigned int precision = 15;
1024  unsigned int minPrecision = 13;
1025  while(failed) {
1026  try {
1027  std::auto_ptr< geos::geom::Geometry > resultAuto =
1028  BinaryOp(geomFirst, geomSecond, geos::operation::overlay::overlayOp(code));
1029  failed = false;
1030  result = resultAuto->clone();
1031  }
1032  catch(geos::util::GEOSException *exc) {
1033  // Just in case the clone failed....
1034  if(!failed || precision == minPrecision) throw;
1035 
1036  delete exc;
1037  }
1038  catch(...) {
1039  if(precision == minPrecision) {
1040  IString msg = "An unknown geos error occurred";
1041  throw IException(IException::Programmer, msg, _FILEINFO_);
1042  }
1043 
1044  if(!failed) {
1045  IString msg = "An unknown geos error occurred when attempting to clone a geometry";
1046  throw IException(IException::Programmer, msg, _FILEINFO_);
1047  }
1048  }
1049 
1050  // try reducing the precision
1051  if(failed) {
1052  precision --;
1053  geos::geom::Geometry *tmp = ReducePrecision(geomFirst, precision);
1054  delete geomFirst;
1055  geomFirst = tmp;
1056 
1057  tmp = ReducePrecision(geomSecond, precision);
1058  delete geomSecond;
1059  geomSecond = tmp;
1060  }
1061  }
1062 
1063  delete geomFirst;
1064  geomFirst = NULL;
1065 
1066  delete geomSecond;
1067  geomSecond = NULL;
1068 
1069  if(result && !result->isValid()) {
1070  try {
1071  geos::geom::Geometry *newResult = FixGeometry(result);
1072 
1073  if(fabs(newResult->getArea() / result->getArea() - 1.0) > 0.50) {
1074  delete newResult;
1075  delete result;
1076 
1077  IString msg = "Operation [" + IString((int)opcode) + "] failed";
1078  throw IException(IException::Programmer, msg, _FILEINFO_);
1079  }
1080 
1081  delete result;
1082  result = newResult;
1083  }
1084  catch(IException &e) {
1085  IString msg = "Operation [" + IString((int)opcode) + "] failed";
1086  throw IException(IException::Programmer, msg, _FILEINFO_);
1087  }
1088  }
1089 
1090  if(result == NULL) {
1091  IString msg = "Operation [" + IString((int)opcode) + " failed";
1092  throw IException(IException::Programmer, msg, _FILEINFO_);
1093  }
1094 
1095  return result;
1096  }
1097 
1108  geos::geom::Geometry *PolygonTools::FixGeometry(const geos::geom::Geometry *geom) {
1109  if(geom->getGeometryTypeId() == geos::geom::GEOS_MULTIPOLYGON) {
1110  return FixGeometry(dynamic_cast<const geos::geom::MultiPolygon *>(geom));
1111  }
1112  if(geom->getGeometryTypeId() == geos::geom::GEOS_LINEARRING) {
1113  return FixGeometry(dynamic_cast<const geos::geom::LinearRing *>(geom));
1114  }
1115  if(geom->getGeometryTypeId() == geos::geom::GEOS_POLYGON) {
1116  return FixGeometry(dynamic_cast<const geos::geom::Polygon *>(geom));
1117  }
1118  if(geom->getGeometryTypeId() == geos::geom::GEOS_GEOMETRYCOLLECTION) {
1119  return FixGeometry(MakeMultiPolygon(geom));
1120  }
1121  else {
1122  IString msg = "PolygonTools::FixGeometry does not support [" + GetGeometryName(geom) + "]";
1123  throw IException(IException::Programmer, msg, _FILEINFO_);
1124  }
1125  }
1126 
1136  geos::geom::MultiPolygon *PolygonTools::FixGeometry(const geos::geom::MultiPolygon *poly) {
1137  // Maybe two points are on top of each other
1138  vector<geos::geom::Geometry *> *newPolys = new vector<geos::geom::Geometry *>;
1139 
1140  // Convert each polygon in this multi-polygon
1141  for(unsigned int geomIndex = 0; geomIndex < poly->getNumGeometries(); geomIndex ++) {
1142  geos::geom::Polygon *fixedpoly = FixGeometry(
1143  dynamic_cast<const geos::geom::Polygon *>(
1144  poly->getGeometryN(geomIndex)));
1145  if(fixedpoly->isValid()) {
1146  newPolys->push_back(fixedpoly);
1147  }
1148  else {
1149  delete fixedpoly;
1150  }
1151  fixedpoly = NULL;
1152  }
1153 
1154  geos::geom::MultiPolygon *mp = Isis::globalFactory.createMultiPolygon(newPolys);
1155  return mp;
1156  }
1157 
1158 
1159  geos::geom::Polygon *PolygonTools::FixGeometry(const geos::geom::Polygon *poly) {
1160 
1161  // Convert each hole inside this polygon
1162  vector<geos::geom::Geometry *> *holes = new vector<geos::geom::Geometry *>;
1163  for(unsigned int holeIndex = 0; holeIndex < poly->getNumInteriorRing(); holeIndex ++) {
1164  const geos::geom::LineString *thisHole = poly->getInteriorRingN(holeIndex);
1165 
1166  // We hope they are all linear rings (closed/simple), but if not just leave it be
1167  if(thisHole->getGeometryTypeId() != geos::geom::GEOS_LINEARRING) {
1168  holes->push_back(thisHole->clone());
1169  continue;
1170  }
1171 
1172  try {
1173  geos::geom::LinearRing *newHole = FixGeometry((geos::geom::LinearRing *)thisHole);
1174  holes->push_back(newHole);
1175  }
1176  catch (IException &e) {
1177  IString msg = "Failed when attempting to fix interior ring of multipolygon";
1178  throw IException(IException::Programmer, msg, _FILEINFO_);
1179  }
1180  } // end num holes in polygon loop
1181 
1182 
1183  const geos::geom::LineString *exterior = poly->getExteriorRing();
1184 
1185  try {
1186  geos::geom::LinearRing *newExterior = NULL;
1187 
1188  if(exterior->getGeometryTypeId() == geos::geom::GEOS_LINEARRING) {
1189  newExterior = FixGeometry((geos::geom::LinearRing *)exterior);
1190  }
1191  else {
1192  IString msg = "Failed when attempting to fix exterior ring of polygon. The exterior "
1193  "ring is not simple and closed";
1194  throw IException(IException::Programmer, msg, _FILEINFO_);
1195  }
1196 
1197  return globalFactory.createPolygon(newExterior, holes);
1198  }
1199  catch (IException &e) {
1200  IString msg = "Failed when attempting to fix exterior ring of polygon";
1201  throw IException(e, IException::Programmer, msg, _FILEINFO_);
1202  }
1203  }
1204 
1205 
1221  geos::geom::LinearRing *PolygonTools::FixGeometry(const geos::geom::LinearRing *ring) {
1222 
1223  geos::geom::CoordinateSequence *coords = ring->getCoordinates();
1224 
1225  // Check this, just in case
1226  if(coords->getSize() < 4) {
1227  return globalFactory.createLinearRing(new geos::geom::DefaultCoordinateSequence());
1228  }
1229 
1230  geos::geom::CoordinateSequence *newCoords = new geos::geom::DefaultCoordinateSequence();
1231  const geos::geom::Coordinate *lastCoordinate = &coords->getAt(0);
1232  newCoords->add(*lastCoordinate);
1233 
1234  // Convert each coordinate in this hole
1235  for(unsigned int coordIndex = 1; coordIndex < coords->getSize() - 1; coordIndex ++) {
1236  const geos::geom::Coordinate *thisCoordinate = &coords->getAt(coordIndex);
1237 
1238  // we're going to compare the decimal place of the current point to the decimal place
1239  // of the difference, if they are drastically different then geos might not be seeing them
1240  // correctly.
1241  double difference[2] = {
1242  lastCoordinate->x - thisCoordinate->x,
1243  lastCoordinate->y - thisCoordinate->y,
1244  };
1245 
1246  // geos isnt differentiating between points this close
1247  double minDiff = fabs(DecimalPlace(thisCoordinate->x) - DecimalPlace(difference[0]));
1248 
1249  minDiff = min(minDiff, fabs(DecimalPlace(thisCoordinate->y) - DecimalPlace(difference[1])));
1250 
1251  // Cases where the difference in one direction is exactly zero, and the other direction is
1252  // next to zero appear often enough (especially in despike).
1253  if(difference[0] == 0.0 && difference[1] != 0.0) {
1254  // subtracted the two points, got deltaX = 0.0, use the y difference decimal place
1255  minDiff = fabs(DecimalPlace(thisCoordinate->y) - DecimalPlace(difference[1]));
1256  }
1257  else if(difference[1] == 0.0 && difference[0] != 0.0) {
1258  // subtracted the two points, got deltaY = 0.0, use the x difference decimal place
1259  minDiff = fabs(DecimalPlace(thisCoordinate->x) - DecimalPlace(difference[0]));
1260  }
1261  else if(difference[0] == 0.0 && difference[1] == 0.0) {
1262  // subtracted the two points, got 0.0, so it's same point... make sure it gets ignored!
1263  minDiff = 1E99;
1264  }
1265 
1266  // geos has a hard time differentiating when points get too close...
1267  if(minDiff <= 15) {
1268  newCoords->add(*thisCoordinate);
1269  lastCoordinate = thisCoordinate;
1270  }
1271  } // end num coords in hole loop
1272 
1273  newCoords->add(geos::geom::Coordinate(newCoords->getAt(0).x, newCoords->getAt(0).y));
1274  geos::geom::LinearRing *newRing = NULL;
1275 
1276  // Now that we've weeded out any bad coordinates, let's rebuild the geometry
1277  try {
1278  if(newCoords->getSize() > 3) {
1279  newRing = globalFactory.createLinearRing(newCoords);
1280  }
1281  else {
1282  delete newCoords;
1283  newCoords = NULL;
1284  }
1285  }
1286  catch(geos::util::GEOSException *exc) {
1287  delete exc;
1288  exc = NULL;
1289 
1290  IString msg = "Error when attempting to fix linear ring";
1291  throw IException(IException::Programmer, msg, _FILEINFO_);
1292  }
1293 
1294  if(newRing && !newRing->isValid() && ring->isValid()) {
1295  IString msg = "Failed when attempting to fix linear ring";
1296  throw IException(IException::Programmer, msg, _FILEINFO_);
1297  }
1298  else if(!newRing || !newRing->isValid()) {
1299  if(newRing) {
1300  delete newRing;
1301  }
1302 
1303  newRing = dynamic_cast<geos::geom::LinearRing *>(ring->clone());
1304  }
1305 
1306  return newRing;
1307  }
1308 
1309 
1321  int PolygonTools::DecimalPlace(double num) {
1322  // 0.0 = decimal place 0
1323  if(num == 0.0) return 0;
1324 
1325  num = fabs(num);
1326 
1327  int decimalPlace = 1;
1328  while(num < 1.0) {
1329  num *= 10.0;
1330  decimalPlace --;
1331  }
1332 
1333  while(num > 10.0) {
1334  num /= 10.0;
1335  decimalPlace ++;
1336  }
1337 
1338  return decimalPlace;
1339  }
1340 
1349  geos::geom::Geometry *PolygonTools::Difference(const geos::geom::Geometry *geom1,
1350  const geos::geom::Geometry *geom2) {
1351  try {
1352  return Operate(geom1, geom2, (unsigned int)geos::operation::overlay::OverlayOp::opDIFFERENCE);
1353  }
1354  catch(geos::util::GEOSException *exc) {
1355  IString msg = "Difference operation failed. The reason given was [" +
1356  IString(exc->what()) + "]";
1357  delete exc;
1358  throw IException(IException::Programmer, msg, _FILEINFO_);
1359  }
1360  catch(IException &e) {
1361  IString msg = "Difference operation failed";
1362  throw IException(e, IException::Programmer, msg, _FILEINFO_);
1363  }
1364  catch(...) {
1365  IString msg = "Difference operation failed for an unknown reason";
1366  throw IException(IException::Programmer, msg, _FILEINFO_);
1367  }
1368  }
1369 
1370 
1384  geos::geom::MultiPolygon *PolygonTools::MakeMultiPolygon(const geos::geom::Geometry *geom) {
1385  // The area of the geometry is too small, so just ignore it.
1386  if(geom->isEmpty()) {
1387  return Isis::globalFactory.createMultiPolygon();
1388  }
1389 
1390  else if(geom->getArea() - DBL_EPSILON <= DBL_EPSILON) {
1391  return Isis::globalFactory.createMultiPolygon();
1392  }
1393 
1394  else if(geom->getGeometryTypeId() == geos::geom::GEOS_MULTIPOLYGON) {
1395  return dynamic_cast<geos::geom::MultiPolygon *>(geom->clone());
1396  }
1397 
1398  else if(geom->getGeometryTypeId() == geos::geom::GEOS_POLYGON) {
1399  vector<geos::geom::Geometry *> *polys = new vector<geos::geom::Geometry *>;
1400  polys->push_back(geom->clone());
1401  geos::geom::MultiPolygon *mp = Isis::globalFactory.createMultiPolygon(polys);
1402  return mp;
1403  }
1404 
1405  else if(geom->getGeometryTypeId() == geos::geom::GEOS_GEOMETRYCOLLECTION) {
1406  vector<geos::geom::Geometry *> * polys =
1407  new vector<geos::geom::Geometry *>;
1408  const geos::geom::GeometryCollection *gc =
1409  dynamic_cast<const geos::geom::GeometryCollection *>(geom);
1410  for(unsigned int i = 0; i < gc->getNumGeometries(); ++i) {
1411  geos::geom::MultiPolygon *subMultiPoly =
1412  MakeMultiPolygon(gc->getGeometryN(i));
1413 
1414  for(unsigned int subPoly = 0;
1415  subPoly < subMultiPoly->getNumGeometries();
1416  subPoly ++) {
1417  const geos::geom::Polygon *poly =
1418  dynamic_cast<const geos::geom::Polygon *>(
1419  subMultiPoly->getGeometryN(subPoly));
1420  polys->push_back(dynamic_cast<geos::geom::Polygon *>(poly->clone()));
1421  }
1422  }
1423 
1424  geos::geom::MultiPolygon *mp =
1425  Isis::globalFactory.createMultiPolygon(polys);
1426  if(mp->getArea() - DBL_EPSILON <= DBL_EPSILON) {
1427  delete mp;
1428  mp = Isis::globalFactory.createMultiPolygon();
1429  }
1430 
1431  return mp;
1432  }
1433 
1434  // All other geometry types are invalid so ignore them
1435  else {
1436  return Isis::globalFactory.createMultiPolygon();
1437  }
1438  }
1439 
1440 
1441  geos::geom::MultiPolygon *PolygonTools::FixSeam(
1442  const geos::geom::Polygon *polyA, const geos::geom::Polygon *polyB) {
1443  geos::geom::CoordinateSequence *polyAPoints = polyA->getCoordinates();
1444  geos::geom::CoordinateSequence *polyBPoints = polyB->getCoordinates();
1445 
1446  unsigned int aIntersectionBegin = 0;
1447  unsigned int aIntersectionEnd = 0;
1448  unsigned int bIntersectionBegin = 0;
1449  unsigned int bIntersectionEnd = 0;
1450 
1451  bool intersectionStarted = false;
1452  bool intersectionEnded = false;
1453 
1454  unsigned int lastBMatch = 0;
1455  for (unsigned int i = 0;
1456  !intersectionEnded && i < polyAPoints->getSize();
1457  i++) {
1458 
1459  bool foundEquivalent = false;
1460 
1461  geos::geom::Coordinate coordA = polyAPoints->getAt(i);
1462  coordA = *ReducePrecision(&coordA, 13);
1463 
1464  for (unsigned int j = 0;
1465  !foundEquivalent && j < polyBPoints->getSize();
1466  j++) {
1467  geos::geom::Coordinate coordB = polyBPoints->getAt(j);
1468  coordB = *ReducePrecision(&coordB, 13);
1469 
1470  foundEquivalent = coordA.equals(coordB);
1471 
1472  if (foundEquivalent) lastBMatch = j;
1473 
1474  if (foundEquivalent && !intersectionStarted) {
1475  intersectionStarted = true;
1476  aIntersectionBegin = i;
1477  bIntersectionBegin = j;
1478  }
1479  }
1480 
1481  if (!foundEquivalent && intersectionStarted && !intersectionEnded) {
1482  intersectionEnded = true;
1483  aIntersectionEnd = i;
1484  bIntersectionEnd = lastBMatch;
1485  }
1486  }
1487 
1488  geos::geom::MultiPolygon * result = NULL;
1489  if (intersectionStarted && intersectionEnded) {
1490  geos::geom::CoordinateSequence *merged =
1491  new geos::geom::CoordinateArraySequence;
1492 
1493  unsigned int i = 0;
1494  for (i = 0; i < aIntersectionBegin; i ++) {
1495  merged->add(polyAPoints->getAt(i));
1496  }
1497 
1498  i = bIntersectionBegin;
1499  while (i != bIntersectionEnd) {
1500  merged->add(polyBPoints->getAt(i));
1501  i++;
1502  if (i >= polyBPoints->getSize()) i = 0;
1503  }
1504 
1505  for (i = aIntersectionEnd; i < polyAPoints->getSize() - 1; i++) {
1506  merged->add(polyAPoints->getAt(i));
1507  }
1508 
1509  merged->add(merged->getAt(0));
1510  result = MakeMultiPolygon(globalFactory.createPolygon(
1511  globalFactory.createLinearRing(merged), NULL));
1512  }
1513 
1514  return result;
1515  }
1516 
1517 
1518  geos::geom::MultiPolygon *PolygonTools::FixSeam(
1519  const geos::geom::MultiPolygon *poly) {
1520 
1521  std::vector<geos::geom::Geometry *> *polys =
1522  new std::vector<geos::geom::Geometry *>;
1523 
1524 
1525  for(unsigned int copyIndex = 0;
1526  copyIndex < poly->getNumGeometries();
1527  copyIndex ++) {
1528  polys->push_back(poly->getGeometryN(copyIndex)->clone());
1529  }
1530 
1531  unsigned int outerPolyIndex = 0;
1532 
1533  while(outerPolyIndex + 1 < polys->size()) {
1534  unsigned int innerPolyIndex = outerPolyIndex + 1;
1535 
1536  while(innerPolyIndex < polys->size()) {
1537  geos::geom::MultiPolygon *fixedPair = FixSeam(
1538  dynamic_cast<geos::geom::Polygon *>(polys->at(outerPolyIndex)),
1539  dynamic_cast<geos::geom::Polygon *>(polys->at(innerPolyIndex)));
1540 
1541  if(fixedPair != NULL) {
1542  geos::geom::Geometry *oldInnerPoly = polys->at(innerPolyIndex);
1543  geos::geom::Geometry *oldOuterPoly = polys->at(outerPolyIndex);
1544 
1545  polys->erase(polys->begin() + innerPolyIndex);
1546  (*polys)[outerPolyIndex] = fixedPair->getGeometryN(0)->clone();
1547  innerPolyIndex = outerPolyIndex + 1;
1548 
1549  delete fixedPair;
1550  fixedPair = NULL;
1551 
1552  delete oldInnerPoly;
1553  oldInnerPoly = NULL;
1554 
1555  delete oldOuterPoly;
1556  oldOuterPoly = NULL;
1557  }
1558  else {
1559  innerPolyIndex ++;
1560  }
1561  }
1562 
1563  outerPolyIndex ++;
1564  }
1565 
1566  return globalFactory.createMultiPolygon(polys);
1567  }
1568 
1569 
1579  geos::geom::Geometry *PolygonTools::ReducePrecision(const geos::geom::Geometry *geom,
1580  unsigned int precision) {
1581  if(geom->getGeometryTypeId() == geos::geom::GEOS_MULTIPOLYGON) {
1582  return ReducePrecision(
1583  dynamic_cast<const geos::geom::MultiPolygon *>(geom), precision);
1584  }
1585  if(geom->getGeometryTypeId() == geos::geom::GEOS_LINEARRING) {
1586  return ReducePrecision(
1587  dynamic_cast<const geos::geom::LinearRing *>(geom), precision);
1588  }
1589  if(geom->getGeometryTypeId() == geos::geom::GEOS_POLYGON) {
1590  return ReducePrecision(
1591  dynamic_cast<const geos::geom::Polygon *>(geom), precision);
1592  }
1593  if(geom->getGeometryTypeId() == geos::geom::GEOS_GEOMETRYCOLLECTION) {
1594  return ReducePrecision(MakeMultiPolygon(geom), precision);
1595  }
1596  else {
1597  IString msg = "PolygonTools::ReducePrecision does not support [" +
1598  GetGeometryName(geom) + "]";
1599  throw IException(IException::Programmer, msg, _FILEINFO_);
1600  }
1601  }
1602 
1603 
1613  geos::geom::MultiPolygon *PolygonTools::ReducePrecision(const geos::geom::MultiPolygon *poly,
1614  unsigned int precision) {
1615  // Maybe two points are on top of each other
1616  vector<geos::geom::Geometry *> *newPolys = new vector<geos::geom::Geometry *>;
1617 
1618  // Convert each polygon in this multi-polygon
1619  for(unsigned int geomIndex = 0; geomIndex < poly->getNumGeometries(); geomIndex ++) {
1620  geos::geom::Geometry *lowerPrecision = ReducePrecision(
1621  dynamic_cast<const geos::geom::Polygon *>(
1622  poly->getGeometryN(geomIndex)),
1623  precision);
1624 
1625  if(!lowerPrecision->isEmpty()) {
1626  newPolys->push_back(lowerPrecision);
1627  }
1628  else {
1629  delete lowerPrecision;
1630  }
1631  }
1632 
1633  geos::geom::MultiPolygon *mp = Isis::globalFactory.createMultiPolygon(newPolys);
1634  return mp;
1635  }
1636 
1637 
1647  geos::geom::Polygon *PolygonTools::ReducePrecision(const geos::geom::Polygon *poly,
1648  unsigned int precision) {
1649  // Convert each hole inside this polygon
1650  vector<geos::geom::Geometry *> *holes = new vector<geos::geom::Geometry *>;
1651  for(unsigned int holeIndex = 0; holeIndex < poly->getNumInteriorRing(); holeIndex ++) {
1652  const geos::geom::LineString *thisHole = poly->getInteriorRingN(holeIndex);
1653 
1654  // We hope they are all linear rings (closed/simple), but if not just leave it be
1655  if(thisHole->getGeometryTypeId() != geos::geom::GEOS_LINEARRING) {
1656  holes->push_back(thisHole->clone());
1657  continue;
1658  }
1659 
1660  try {
1661  geos::geom::LinearRing *newHole = ReducePrecision((geos::geom::LinearRing *)thisHole,
1662  precision);
1663 
1664  if(!newHole->isEmpty()) {
1665  holes->push_back(newHole);
1666  }
1667  else {
1668  delete newHole;
1669  }
1670 
1671  }
1672  catch(IException &e) {
1673  IString msg = "Failed when attempting to fix interior ring of multipolygon";
1674  throw IException(e, IException::Programmer, msg, _FILEINFO_);
1675  }
1676  } // end num holes in polygon loop
1677 
1678 
1679  const geos::geom::LineString *exterior = poly->getExteriorRing();
1680 
1681  try {
1682  geos::geom::LinearRing *newExterior = NULL;
1683 
1684  if(exterior->getGeometryTypeId() == geos::geom::GEOS_LINEARRING) {
1685  newExterior = ReducePrecision((geos::geom::LinearRing *)exterior, precision);
1686  }
1687  else {
1688  IString msg = "Failed when attempting to fix exterior ring of polygon. The exterior "
1689  "ring is not simple and closed";
1690  throw IException(IException::Programmer, msg, _FILEINFO_);
1691  }
1692 
1693  return globalFactory.createPolygon(newExterior, holes);
1694  }
1695  catch(IException &e) {
1696  IString msg = "Failed when attempting to fix exterior ring of polygon";
1697  throw IException(e, IException::Programmer, msg, _FILEINFO_);
1698  }
1699  }
1700 
1701 
1711  geos::geom::LinearRing *PolygonTools::ReducePrecision(const geos::geom::LinearRing *ring,
1712  unsigned int precision) {
1713  geos::geom::CoordinateSequence *coords = ring->getCoordinates();
1714 
1715  // Check this, just in case
1716  if(coords->getSize() <= 0) return dynamic_cast<geos::geom::LinearRing *>(
1717  ring->clone());
1718 
1719  geos::geom::CoordinateSequence *newCoords = new geos::geom::DefaultCoordinateSequence();
1720  geos::geom::Coordinate *coord = ReducePrecision(&coords->getAt(0), precision);
1721  newCoords->add(*coord);
1722  delete coord;
1723  coord = NULL;
1724 
1725  // Convert each coordinate in this ring
1726  for(unsigned int coordIndex = 1; coordIndex < coords->getSize() - 1; coordIndex ++) {
1727  const geos::geom::Coordinate *thisCoordinate = &coords->getAt(coordIndex);
1728  coord = ReducePrecision(thisCoordinate, precision);
1729  newCoords->add(*coord);
1730  delete coord;
1731  coord = NULL;
1732  }
1733 
1734  newCoords->add(geos::geom::Coordinate(newCoords->getAt(0).x, newCoords->getAt(0).y));
1735  geos::geom::LinearRing *newRing = NULL;
1736 
1737  // Now that we've weeded out any bad coordinates, let's rebuild the geometry
1738  try {
1739  newRing = globalFactory.createLinearRing(newCoords);
1740  }
1741  catch(geos::util::GEOSException *exc) {
1742  delete exc;
1743  exc = NULL;
1744 
1745  IString msg = "Error when attempting to reduce precision of linear ring";
1746  throw IException(IException::Programmer, msg, _FILEINFO_);
1747  }
1748 
1749  // try to despike
1750  try {
1751  geos::geom::LinearRing *tmp = Despike(newRing);
1752  delete newRing;
1753  newRing = tmp;
1754  }
1755  catch(IException &e) {
1756  }
1757 
1758  if(!newRing->isValid()) {
1759  IString msg = "Failed when attempting to reduce precision of linear ring";
1760  throw IException(IException::Programmer, msg, _FILEINFO_);
1761  }
1762 
1763  return newRing;
1764  }
1765 
1766 
1776  geos::geom::Coordinate *PolygonTools::ReducePrecision(const geos::geom::Coordinate *coord,
1777  unsigned int precision) {
1778  return new geos::geom::Coordinate(
1779  ReducePrecision(coord->x, precision),
1780  ReducePrecision(coord->y, precision),
1781  ReducePrecision(coord->z, precision)
1782  );
1783  }
1784 
1785 
1795  double PolygonTools::ReducePrecision(double num, unsigned int precision) {
1796  double result = num;
1797 
1798  // If num == nan then this test will fail
1799  if (num == num) {
1800  int decimalPlace = DecimalPlace(num);
1801  double factor = pow(10.0, (int)decimalPlace);
1802 
1803  // reduced num is in the form 0.nnnnnnnnnn...
1804  double reducedNum = num / factor;
1805 
1806  double cutoff = pow(10.0, (int)precision);
1807  double round_offset = (num < 0) ? -0.5 : 0.5;
1808 
1809  // cast off the digits past the precision's place
1810  reducedNum = ((long long)(reducedNum * cutoff + round_offset)) / cutoff;
1811  result = reducedNum * factor;
1812  }
1813 
1814  return result;
1815  }
1816 
1817 
1826  QString PolygonTools::GetGeometryName(const geos::geom::Geometry *geom) {
1827  switch(geom->getGeometryTypeId()) {
1828  case geos::geom::GEOS_POINT:
1829  return "Point";
1830  case geos::geom::GEOS_LINESTRING:
1831  return "Line String";
1832  case geos::geom::GEOS_LINEARRING:
1833  return "Linear Ring";
1834  case geos::geom::GEOS_POLYGON:
1835  return "Polygon";
1836  case geos::geom::GEOS_MULTIPOINT:
1837  return "Multi Point";
1838  case geos::geom::GEOS_MULTILINESTRING:
1839  return "Multi Line String";
1840  case geos::geom::GEOS_MULTIPOLYGON:
1841  return "Multi Polygon";
1842  case geos::geom::GEOS_GEOMETRYCOLLECTION:
1843  return "Geometry Collection";
1844  default:
1845  return "UNKNOWN";
1846  }
1847  }
1848 
1849 
1850  bool PolygonTools::Equal(const geos::geom::MultiPolygon *poly1,
1851  const geos::geom::MultiPolygon *poly2) {
1852 
1853  vector<const geos::geom::Polygon *> polys1;
1854  vector<const geos::geom::Polygon *> polys2;
1855 
1856  if(poly1->getNumGeometries() != poly2->getNumGeometries()) return false;
1857 
1858  // Convert each polygon in this multi-polygon
1859  for(unsigned int geomIndex = 0; geomIndex < poly1->getNumGeometries(); geomIndex ++) {
1860  polys1.push_back(dynamic_cast<const geos::geom::Polygon *>(
1861  poly1->getGeometryN(geomIndex)));
1862  polys2.push_back(dynamic_cast<const geos::geom::Polygon *>(
1863  poly2->getGeometryN(geomIndex)));
1864  }
1865 
1866  for(int p1 = polys1.size() - 1; (p1 >= 0) && polys1.size(); p1 --) {
1867  for(int p2 = polys2.size() - 1; (p2 >= 0) && polys2.size(); p2 --) {
1868  if(Equal(polys1[p1], polys2[p2])) {
1869  // Delete polys1[p1] by replacing it with the last Polygon in polys1
1870  polys1[p1] = polys1[polys1.size()-1];
1871  polys1.resize(polys1.size() - 1);
1872  // Delete polys2[p2] by replacing it with the last Polygon in polys2
1873  polys2[p2] = polys2[polys2.size()-1];
1874  polys2.resize(polys2.size() - 1);
1875  }
1876  }
1877  }
1878 
1879  return (polys1.size() == 0) && (polys2.size() == 0);
1880  }
1881 
1882 
1883  bool PolygonTools::Equal(const geos::geom::Polygon *poly1, const geos::geom::Polygon *poly2) {
1884  vector<const geos::geom::LineString *> holes1;
1885  vector<const geos::geom::LineString *> holes2;
1886 
1887  if(poly1->getNumInteriorRing() != poly2->getNumInteriorRing()) return false;
1888 
1889  if(!Equal(poly1->getExteriorRing(), poly2->getExteriorRing())) return false;
1890 
1891  // Convert each hole inside this polygon
1892  for(unsigned int holeIndex = 0; holeIndex < poly1->getNumInteriorRing(); holeIndex ++) {
1893 
1894  // We hope they are all linear rings (closed/simple), but if not just leave it be
1895  if(poly1->getInteriorRingN(holeIndex)->getGeometryTypeId() == geos::geom::GEOS_LINESTRING) {
1896  holes1.push_back(poly1->getInteriorRingN(holeIndex));
1897  }
1898 
1899  if(poly2->getInteriorRingN(holeIndex)->getGeometryTypeId() == geos::geom::GEOS_LINESTRING) {
1900  holes2.push_back(poly2->getInteriorRingN(holeIndex));
1901  }
1902 
1903  }
1904 
1905  if(holes1.size() != holes2.size()) return false;
1906 
1907  for(int h1 = holes1.size() - 1; (h1 >= 0) && holes1.size(); h1 --) {
1908  for(int h2 = holes2.size() - 1; (h2 >= 0) && holes2.size(); h2 --) {
1909  if(Equal(holes1[h1], holes2[h2])) {
1910  // Delete holes1[h1] by replacing it with the last Polygon in holes1
1911  holes1[h1] = holes1[holes1.size()-1];
1912  holes1.resize(holes1.size() - 1);
1913  // Delete holes2[h2] by replacing it with the last Polygon in holes2
1914  holes2[h2] = holes2[holes2.size()-1];
1915  holes2.resize(holes2.size() - 1);
1916  }
1917  }
1918  }
1919 
1920  return (holes1.size() == 0) && (holes2.size() == 0);
1921  }
1922 
1923 
1924  bool PolygonTools::Equal(const geos::geom::LineString *lineString1,
1925  const geos::geom::LineString *lineString2) {
1926 
1927  geos::geom::CoordinateSequence *coords1 = lineString1->getCoordinates();
1928  geos::geom::CoordinateSequence *coords2 = lineString2->getCoordinates();
1929  bool isEqual = true;
1930 
1931  if(coords1->getSize() != coords2->getSize()) isEqual = false;
1932 
1933  unsigned int index1 = 0;
1934  unsigned int index2 = 0;
1935 
1936  // -1 extra for dupicate start/end coordinates
1937  for(; index2 < coords2->getSize() - 1 && isEqual; index2 ++) {
1938  if(Equal(coords1->getAt(index1), coords2->getAt(index2))) break;
1939  }
1940 
1941  if(index2 == coords2->getSize() - 1) isEqual = false;
1942 
1943  for(; index1 < coords1->getSize() - 1 && isEqual; index1 ++, index2 ++) {
1944  if(!Equal(coords1->getAt(index1), coords2->getAt(index2 % (coords2->getSize() - 1)))) {
1945  isEqual = false;
1946  }
1947  }
1948 
1949  delete coords1;
1950  delete coords2;
1951  return isEqual;
1952  }
1953 
1954 
1955  bool PolygonTools::Equal(const geos::geom::Coordinate &coord1,
1956  const geos::geom::Coordinate &coord2) {
1957 
1958  if(!Equal(coord1.x, coord2.x)) return false;
1959  if(!Equal(coord1.y, coord2.y)) return false;
1960  if(!Equal(coord1.y, coord2.y)) return false;
1961 
1962  return true;
1963  }
1964 
1965 
1966  bool PolygonTools::Equal(const double d1, const double d2) {
1967  const double cutoff = 1e15;
1968 
1969  if(DecimalPlace(d1) != DecimalPlace(d2)) return false;
1970 
1971  int decimalPlace = DecimalPlace(d1);
1972  double factor = pow(10.0, (int)decimalPlace);
1973 
1974  // reduced num is in the form 0.nnnnnnnnnn...
1975  double reducedNum = d1 / factor;
1976 
1977  double round_offset = (d1 < 0) ? -0.5 : 0.5;
1978 
1979  // cast off the digits past the precision's place
1980  long long num1 = ((long long)(reducedNum * cutoff + round_offset));
1981 
1982  factor = pow(10.0, (int)decimalPlace);
1983 
1984  // reduced num is in the form 0.nnnnnnnnnn...
1985  reducedNum = d2 / factor;
1986 
1987  round_offset = (d2 < 0) ? -0.5 : 0.5;
1988 
1989  // cast off the digits past the precision's place
1990  long long num2 = ((long long)(reducedNum * cutoff + round_offset));
1991 
1992 
1993  return (num1 == num2);
1994  }
1995 
1996 
2011  geos::geom::MultiPolygon *PolygonTools::SplitPolygonOn360(const geos::geom::Polygon *inPoly) {
2012  bool convertLon = false;
2013  bool negAdjust = false;
2014  bool newCoords = false; // coordinates have been adjusted
2015  geos::geom::CoordinateSequence *newLonLatPts = new geos::geom::CoordinateArraySequence();
2016  double lon, lat;
2017  double lonOffset = 0;
2018  geos::geom::CoordinateSequence *inPolyCoords = inPoly->getCoordinates();
2019  double prevLon = inPolyCoords->getAt(0).x;
2020  double prevLat = inPolyCoords->getAt(0).y;
2021 
2022  newLonLatPts->add(geos::geom::Coordinate(prevLon, prevLat));
2023  double dist = 0.0;
2024  for (unsigned int i = 1; i < inPolyCoords->getSize(); i++) {
2025  lon = inPolyCoords->getAt(i).x;
2026  lat = inPolyCoords->getAt(i).y;
2027  // check to see if you just crossed the Meridian
2028  if (abs(lon - prevLon) > 180 && (prevLat != 90 && prevLat != -90)) {
2029  newCoords = true;
2030  // if you were already converting then stop (crossed Meridian even number of times)
2031  if (convertLon) {
2032  convertLon = false;
2033  lonOffset = 0;
2034  }
2035  else { // Need to start converting again, deside how to adjust coordinates
2036  if ((lon - prevLon) > 0) {
2037  lonOffset = -360.;
2038  negAdjust = true;
2039  }
2040  else if ((lon - prevLon) < 0) {
2041  lonOffset = 360.;
2042  negAdjust = false;
2043  }
2044  convertLon = true;
2045  }
2046  }
2047 
2048  // Change to a minimum calculation
2049  if (newCoords && dist == 0.0) {
2050  double longitude = (lon + lonOffset) - prevLon;
2051  double latitude = lat - prevLat;
2052  dist = std::sqrt((longitude * longitude) + (latitude * latitude));
2053  }
2054 
2055  // add coord
2056  newLonLatPts->add(geos::geom::Coordinate(lon + lonOffset, lat));
2057 
2058  // set current to old
2059  prevLon = lon;
2060  prevLat = lat;
2061  }
2062 
2063  delete inPolyCoords;
2064 
2065  // Nothing was done so return
2066  if (!newCoords) {
2067  geos::geom::Polygon *newPoly = globalFactory.createPolygon
2068  (globalFactory.createLinearRing(newLonLatPts), NULL);
2069  geos::geom::MultiPolygon *multi_polygon = PolygonTools::MakeMultiPolygon(newPoly);
2070  delete newLonLatPts;
2071  return multi_polygon;
2072  }
2073 
2074  // bisect into seperate polygons
2075  try {
2076  geos::geom::Polygon *newPoly = globalFactory.createPolygon
2077  (globalFactory.createLinearRing(newLonLatPts), NULL);
2078 
2079  geos::geom::CoordinateSequence *pts = new geos::geom::CoordinateArraySequence();
2080  geos::geom::CoordinateSequence *pts2 = new geos::geom::CoordinateArraySequence();
2081 
2082  // Depending on direction of compensation bound accordingly
2083  //***************************************************
2084 
2085  // please verify correct if you change these values
2086  //***************************************************
2087  if (negAdjust) {
2088  pts->add(geos::geom::Coordinate(0., 90.));
2089  pts->add(geos::geom::Coordinate(-360., 90.));
2090  pts->add(geos::geom::Coordinate(-360., -90.));
2091  pts->add(geos::geom::Coordinate(0., -90.));
2092  for (double lat = -90.0 + dist; lat < 90.0; lat += dist) {
2093  pts->add(geos::geom::Coordinate(0.0, lat));
2094  }
2095  pts->add(geos::geom::Coordinate(0., 90.));
2096  pts2->add(geos::geom::Coordinate(0., 90.));
2097  pts2->add(geos::geom::Coordinate(360., 90.));
2098  pts2->add(geos::geom::Coordinate(360., -90.));
2099  pts2->add(geos::geom::Coordinate(0., -90.));
2100  for (double lat = -90.0 + dist; lat < 90.0; lat += dist) {
2101  pts2->add(geos::geom::Coordinate(0.0, lat));
2102  }
2103  pts2->add(geos::geom::Coordinate(0., 90.));
2104  }
2105  else {
2106  pts->add(geos::geom::Coordinate(360., 90.));
2107  pts->add(geos::geom::Coordinate(720., 90.));
2108  pts->add(geos::geom::Coordinate(720., -90.));
2109  pts->add(geos::geom::Coordinate(360., -90.));
2110  for (double lat = -90.0 + dist; lat < 90.0; lat += dist) {
2111  pts->add(geos::geom::Coordinate(360.0, lat));
2112  }
2113  pts->add(geos::geom::Coordinate(360., 90.));
2114  pts2->add(geos::geom::Coordinate(360., 90.));
2115  pts2->add(geos::geom::Coordinate(0., 90.));
2116  pts2->add(geos::geom::Coordinate(0., -90.));
2117  pts2->add(geos::geom::Coordinate(360., -90.));
2118  for (double lat = -90.0 + dist; lat < 90.0; lat += dist) {
2119  pts2->add(geos::geom::Coordinate(360.0, lat));
2120  }
2121  pts2->add(geos::geom::Coordinate(360., 90.));
2122  }
2123 
2124  geos::geom::Polygon *boundaryPoly = globalFactory.createPolygon
2125  (globalFactory.createLinearRing(pts), NULL);
2126  geos::geom::Polygon *boundaryPoly2 = globalFactory.createPolygon
2127  (globalFactory.createLinearRing(pts2), NULL);
2128  /*------------------------------------------------------------------------
2129  / Intersecting the original polygon (converted coordinates) with the
2130  / boundary polygons will create the multi polygons with the converted coordinates.
2131  / These will need to be converted back to the original coordinates.
2132  /-----------------------------------------------------------------------*/
2133  geos::geom::Geometry *intersection = PolygonTools::Intersect(newPoly, boundaryPoly);
2134  geos::geom::MultiPolygon *convertPoly = PolygonTools::MakeMultiPolygon(intersection);
2135  delete intersection;
2136 
2137  intersection = PolygonTools::Intersect(newPoly, boundaryPoly2);
2138  geos::geom::MultiPolygon *convertPoly2 = PolygonTools::MakeMultiPolygon(intersection);
2139  delete intersection;
2140 
2141  /*------------------------------------------------------------------------
2142  / Adjust points created in the negative space or >360 space to be back in
2143  / the 0-360 world. This will always only need to be done on convertPoly.
2144  / Then add geometries to finalpolys.
2145  /-----------------------------------------------------------------------*/
2146  vector<geos::geom::Geometry *> *finalpolys = new vector<geos::geom::Geometry *>;
2147  geos::geom::Geometry *newGeom = NULL;
2148 
2149  for (unsigned int i = 0; i < convertPoly->getNumGeometries(); i++) {
2150  newGeom = (convertPoly->getGeometryN(i))->clone();
2151  pts = convertPoly->getGeometryN(i)->getCoordinates();
2152  geos::geom::CoordinateSequence *newLonLatPts = new geos::geom::CoordinateArraySequence();
2153  // fix the points
2154 
2155  for (unsigned int k = 0; k < pts->getSize() ; k++) {
2156  double lon = pts->getAt(k).x;
2157  double lat = pts->getAt(k).y;
2158  if (negAdjust) {
2159  lon = lon + 360;
2160  }
2161  else {
2162  lon = lon - 360;
2163  }
2164  newLonLatPts->add(geos::geom::Coordinate(lon, lat), k);
2165 
2166  }
2167  // Add the points to polys
2168  finalpolys->push_back(globalFactory.createPolygon
2169  (globalFactory.createLinearRing(newLonLatPts), NULL));
2170  }
2171 
2172  // This loop is over polygons that will always be in 0-360 space no need to convert
2173  for (unsigned int i = 0; i < convertPoly2->getNumGeometries(); i++) {
2174  newGeom = (convertPoly2->getGeometryN(i))->clone();
2175  finalpolys->push_back(newGeom);
2176  }
2177 
2178  geos::geom::MultiPolygon *multi_polygon = globalFactory.createMultiPolygon(*finalpolys);
2179 
2180  delete finalpolys;
2181  delete newGeom;
2182  delete newLonLatPts;
2183  delete pts;
2184  delete pts2;
2185  return multi_polygon;
2186  }
2187  catch(geos::util::IllegalArgumentException *geosIll) {
2188  std::string msg = "Unable to split polygon on longitude boundry (SplitPolygonOn360) due to ";
2189  msg += "geos illegal argument [" + IString(geosIll->what()) + "]";
2190  delete geosIll;
2191  throw IException(IException::Unknown, msg, _FILEINFO_);
2192  }
2193  catch(geos::util::GEOSException *geosExc) {
2194  std::string msg = "Unable to split polygon on longitude boundry (SplitPolygonOn360) due to ";
2195  msg += "geos exception [" + IString(geosExc->what()) + "]";
2196  delete geosExc;
2197  throw IException(IException::Unknown, msg, _FILEINFO_);
2198  }
2199  catch(IException &e) {
2200  std::string msg = "Unable to split polygon on longitude boundry (SplitPolygonOn360) due to ";
2201  msg += "isis operation exception [" + IString(e.what()) + "]";
2202  throw IException(e, IException::Unknown, msg, _FILEINFO_);
2203  }
2204  catch(...) {
2205  std::string msg = "Unable to split polygon on longitude boundry (SplitPolygonOn360) due to ";
2206  msg += "unknown exception";
2207  throw IException(IException::Unknown, msg, _FILEINFO_);
2208  }
2209  }
2210 } // end namespace isis
2211