Isis 3 Programmer Reference
PolygonTools.cpp
1 
6 /* SPDX-License-Identifier: CC0-1.0 */
7 
8 #include <string>
9 #include <iostream>
10 #include <sstream>
11 #include <iomanip>
12 #include <vector>
13 #include <cmath>
14 
15 #include <QDebug>
16 #include "geos/geom/BinaryOp.h"
17 #include "geos/geom/CoordinateArraySequence.h"
18 #include "geos/geom/CoordinateSequence.h"
19 #include "geos/geom/LinearRing.h"
20 #include "geos/geom/Point.h"
21 #include "geos/geom/Polygon.h"
22 #include "geos/operation/distance/DistanceOp.h"
23 #include "geos/opOverlay.h"
24 #include "geos/operation/overlay/snap/GeometrySnapper.h"
25 
26 #include "SpecialPixel.h"
27 #include "PolygonTools.h"
28 #include "TProjection.h"
29 #include "ProjectionFactory.h"
30 #include "UniversalGroundMap.h"
31 
32 using namespace std;
33 namespace Isis {
34 
47  geos::geom::MultiPolygon *PolygonTools::LatLonToXY(
48  const geos::geom::MultiPolygon &lonLatPolygon, TProjection *projection) {
49  if (projection == NULL) {
50  string msg = "Unable to convert Lon/Lat polygon to X/Y. ";
51  msg += "No projection has was supplied";
52  throw IException(IException::Programmer, msg, _FILEINFO_);
53  }
54 
55  // Convert the Lat/Lon poly coordinates to X/Y coordinates
56  if(lonLatPolygon.isEmpty()) {
57  return globalFactory->createMultiPolygon();
58  }
59  else {
60  vector<geos::geom::Geometry *> *xyPolys = new vector<geos::geom::Geometry *>;
61  // Convert each polygon in this multi-polygon
62  for(unsigned int g = 0; g < lonLatPolygon.getNumGeometries(); ++g) {
63  const geos::geom::Polygon *poly =
64  dynamic_cast<const geos::geom::Polygon *>(
65  lonLatPolygon.getGeometryN(g));
66 
67  // Convert each hole inside this polygon
68  vector<geos::geom::Geometry *> *holes = new vector<geos::geom::Geometry *>;
69  for(unsigned int h = 0; h < poly->getNumInteriorRing(); ++h) {
70  geos::geom::CoordinateSequence *xycoords = new geos::geom::CoordinateArraySequence();
71  geos::geom::CoordinateSequence *llcoords =
72  poly->getInteriorRingN(h)->getCoordinates();
73 
74  // Convert each coordinate in this hole
75  for(unsigned int cord = 0; cord < llcoords->getSize(); ++cord) {
76  projection->SetGround(llcoords->getAt(cord).y,
77  llcoords->getAt(cord).x);
78  xycoords->add(geos::geom::Coordinate(projection->XCoord(),
79  projection->YCoord()));
80  } // end num coords in hole loop
81 
82  geos::geom::LinearRing *hole = globalFactory->createLinearRing(xycoords);
83 
84  if(hole->isValid() && !hole->isEmpty()) {
85  holes->push_back(hole);
86  }
87  else {
88  delete hole;
89  }
90  } // end num holes in polygon loop
91 
92  // Convert the exterior ring of this polygon
93  geos::geom::CoordinateSequence *xycoords = new geos::geom::CoordinateArraySequence();
94  geos::geom::CoordinateSequence *llcoords =
95  poly->getExteriorRing()->getCoordinates();
96 
97  // Convert each coordinate in the exterior ring of this polygon
98  for (unsigned int cord = 0; cord < llcoords->getSize(); ++cord) {
99  projection->SetGround(llcoords->getAt(cord).y,
100  llcoords->getAt(cord).x);
101  xycoords->add(geos::geom::Coordinate(projection->XCoord(),
102  projection->YCoord()));
103  } // end exterior ring coordinate loop
104 
105  geos::geom::Polygon *newPoly = globalFactory->createPolygon(
106  globalFactory->createLinearRing(xycoords), holes);
107 
108  if(newPoly->isValid() && !newPoly->isEmpty() && newPoly->getArea() > 1.0e-14) {
109  xyPolys->push_back(newPoly);
110  }
111  else {
112  delete newPoly;
113  }
114  } // end num geometry in multi-poly
115 
116  // Create a new multipoly from all the new X/Y polygon(s)
117  geos::geom::MultiPolygon *spikedPoly = globalFactory->createMultiPolygon(xyPolys);
118 
119  if(spikedPoly->isValid() && !spikedPoly->isEmpty()) {
120  return spikedPoly;
121  }
122  else {
123  try {
124  geos::geom::MultiPolygon *despikedPoly = Despike(spikedPoly);
125 
126  delete spikedPoly;
127  spikedPoly = NULL;
128 
129  return despikedPoly;
130  }
131  catch(IException &e) {
132  IString msg = "Unable to convert polygon from Lat/Lon to X/Y";
133  throw IException(IException::Programmer, msg, _FILEINFO_);
134  }
135  }
136 
137  } // end else
138  }
139 
140 
153  geos::geom::MultiPolygon *PolygonTools::XYToLatLon(
154  const geos::geom::MultiPolygon &xYPolygon, TProjection *projection) {
155 
156  if(projection == NULL) {
157  string msg = "Unable to convert X/Y polygon to Lon/Lat. ";
158  msg += "No projection was supplied";
159  throw IException(IException::Programmer, msg, _FILEINFO_);
160  }
161 
162  // Convert the X/Y poly coordinates to Lat/Lon coordinates
163  if(xYPolygon.isEmpty()) {
164  return globalFactory->createMultiPolygon();
165  }
166  else {
167  vector<geos::geom::Geometry *> *llPolys = new vector<geos::geom::Geometry *>;
168  // Convert each polygon in this multi-polygon
169  for(unsigned int g = 0; g < xYPolygon.getNumGeometries(); ++g) {
170  const geos::geom::Polygon *poly =
171  dynamic_cast<const geos::geom::Polygon *>(
172  xYPolygon.getGeometryN(g));
173 
174  // Convert each hole inside this polygon
175  vector<geos::geom::Geometry *> *holes = new vector<geos::geom::Geometry *>;
176  for(unsigned int h = 0; h < poly->getNumInteriorRing(); ++h) {
177  geos::geom::CoordinateSequence *llcoords = new geos::geom::CoordinateArraySequence();
178  geos::geom::CoordinateSequence *xycoords =
179  poly->getInteriorRingN(h)->getCoordinates();
180 
181  // Convert each coordinate in this hole
182  for(unsigned int cord = 0; cord < xycoords->getSize(); ++cord) {
183  projection->SetWorld(xycoords->getAt(cord).x,
184  xycoords->getAt(cord).y);
185  llcoords->add(geos::geom::Coordinate(projection->Longitude(),
186  projection->Latitude()));
187  } // end num coords in hole loop
188  holes->push_back(globalFactory->createLinearRing(llcoords));
189  } // end num holes in polygon loop
190 
191  // Convert the exterior ring of this polygon
192  geos::geom::CoordinateSequence *llcoords = new geos::geom::DefaultCoordinateSequence();
193  geos::geom::CoordinateSequence *xycoords =
194  poly->getExteriorRing()->getCoordinates();
195 
196  // Convert each coordinate in the exterior ring of this polygon
197  for(unsigned int cord = 0; cord < xycoords->getSize(); ++cord) {
198  projection->SetWorld(xycoords->getAt(cord).x,
199  xycoords->getAt(cord).y);
200  llcoords->add(geos::geom::Coordinate(projection->Longitude(),
201  projection->Latitude()));
202  } // end exterior ring coordinate loop
203 
204  llPolys->push_back(globalFactory->createPolygon(
205  globalFactory->createLinearRing(llcoords), holes));
206  } // end num geometry in multi-poly
207 
208 
209  // Create a new multipoly from all the new Lat/Lon polygon(s)
210  geos::geom::MultiPolygon *spikedPoly = globalFactory->createMultiPolygon(llPolys);
211 
212  if(spikedPoly->isValid() && !spikedPoly->isEmpty()) {
213  return spikedPoly;
214  }
215  else {
216  try {
217  geos::geom::MultiPolygon *despikedPoly = Despike(spikedPoly);
218 
219  delete spikedPoly;
220  spikedPoly = NULL;
221 
222  return despikedPoly;
223  }
224  catch(IException &e) {
225  IString msg = "Unable to convert polygon from X/Y to Lat/Lon";
226  throw IException(IException::Programmer, msg, _FILEINFO_);
227  }
228  }
229  } // end else
230  }
231 
232 
245  geos::geom::MultiPolygon *PolygonTools::LatLonToSampleLine(
246  const geos::geom::MultiPolygon &lonLatPolygon, UniversalGroundMap *ugm) {
247 
248  if(ugm == NULL) {
249  string msg = "Unable to convert Lon/Lat polygon to Sample/Line. ";
250  msg += "No UniversalGroundMap was supplied";
251  throw IException(IException::Programmer, msg, _FILEINFO_);
252  }
253 
254  // Convert the Lon/Lat poly coordinates to Sample/Line coordinates
255  if(lonLatPolygon.isEmpty()) {
256  return globalFactory->createMultiPolygon();
257  }
258  else {
259  vector<geos::geom::Geometry *> *slPolys = new vector<geos::geom::Geometry *>;
260  // Convert each polygon in this multi-polygon
261  for(unsigned int g = 0; g < lonLatPolygon.getNumGeometries(); g++) {
262  const geos::geom::Polygon *poly =
263  dynamic_cast<const geos::geom::Polygon *>(
264  lonLatPolygon.getGeometryN(g));
265 
266  // Convert each hole inside this polygon
267  vector<geos::geom::Geometry *> *holes = new vector<geos::geom::Geometry *>;
268  for(unsigned int h = 0; h < poly->getNumInteriorRing(); ++h) {
269  geos::geom::CoordinateSequence *slcoords = new geos::geom::DefaultCoordinateSequence();
270  geos::geom::CoordinateSequence *llcoords =
271  poly->getInteriorRingN(h)->getCoordinates();
272 
273  // Convert each coordinate in this hole
274  for(unsigned int cord = 0; cord < llcoords->getSize(); ++cord) {
275  ugm->SetUniversalGround(llcoords->getAt(cord).y,
276  llcoords->getAt(cord).x);
277  slcoords->add(geos::geom::Coordinate(ugm->Sample(),
278  ugm->Line()));
279  } // end num coords in hole loop
280  holes->push_back(globalFactory->createLinearRing(slcoords));
281  delete slcoords;
282  delete llcoords;
283  } // end num holes in polygon loop
284 
285  // Convert the exterior ring of this polygon
286  geos::geom::CoordinateSequence *slcoords = new geos::geom::CoordinateArraySequence();
287  geos::geom::CoordinateSequence *llcoords =
288  poly->getExteriorRing()->getCoordinates();
289 
290  // Convert each coordinate in the exterior ring of this polygon
291  for(unsigned int cord = 0; cord < llcoords->getSize(); ++cord) {
292  if (ugm->SetUniversalGround(llcoords->getAt(cord).y,
293  llcoords->getAt(cord).x)) {
294  slcoords->add(geos::geom::Coordinate(ugm->Sample(),
295  ugm->Line()));
296  }
297  } // end exterior ring coordinate loop
298 
299  // Make sure that the line string is closed.
300  if (slcoords->getSize() > 0 && !slcoords->front().equals(slcoords->back())) {
301  slcoords->add(slcoords->front());
302  }
303 
304  try {
305  slPolys->push_back(globalFactory->createPolygon(
306  globalFactory->createLinearRing(slcoords), holes));
307  }
308  catch (std::exception &e) {
309  throw IException(IException::Unknown,
310  QObject::tr("Unable to convert polygon from Lat/Lon to Sample/Line. The "
311  "error given was [%1].").arg(e.what()),
312  _FILEINFO_);
313  }
314 
315  delete llcoords;
316  } // end num geometry in multi-poly
317 
318  // Create a new multipoly from all the new Sample/Line polygon(s)
319  geos::geom::MultiPolygon *spikedPoly = globalFactory->createMultiPolygon(slPolys);
320 
321  if(spikedPoly->isValid() && !spikedPoly->isEmpty()) {
322  return spikedPoly;
323  }
324  else {
325  try {
326  geos::geom::MultiPolygon *despikedPoly = Despike(spikedPoly);
327 
328  delete spikedPoly;
329  spikedPoly = NULL;
330 
331  return despikedPoly;
332  }
333  catch (IException &e) {
334  IString msg = "Unable to convert polygon from Lat/Lon to Sample/Line";
335  throw IException(IException::Programmer, msg, _FILEINFO_);
336  }
337  }
338  } // end else
339  }
340 
341 
354  geos::geom::MultiPolygon *PolygonTools::CopyMultiPolygon(const geos::geom::MultiPolygon *mpolygon) {
355 
356  vector<geos::geom::Geometry *> *polys = new vector<geos::geom::Geometry *>;
357  for(unsigned int i = 0; i < mpolygon->getNumGeometries(); ++i) {
358  polys->push_back((mpolygon->getGeometryN(i))->clone());
359  }
360  return globalFactory->createMultiPolygon(polys);
361  }
362 
363 
376  geos::geom::MultiPolygon *PolygonTools::CopyMultiPolygon(const geos::geom::MultiPolygon &mpolygon) {
377 
378  vector<geos::geom::Geometry *> *polys = new vector<geos::geom::Geometry *>;
379  for(unsigned int i = 0; i < mpolygon.getNumGeometries(); ++i) {
380  polys->push_back((mpolygon.getGeometryN(i))->clone());
381  }
382  return globalFactory->createMultiPolygon(polys);
383  }
384 
385 
394  QString PolygonTools::ToGML(const geos::geom::MultiPolygon *mpolygon, QString idString,
395  QString schema) {
396 
397  ostringstream os;
398 
399  //Write out the GML header
400  os << "<?xml version=\"1.0\" encoding=\"utf-8\" ?>" << endl;
401  os << "<ogr:FeatureCollection" << endl;
402  os << " xmlns:xsi=\"http://www.w3.org/2001/XMLSchema-instance\"" << endl;
403  os << " xsi:schemaLocation=\"http://ogr.maptools.org/ " << schema << "\"" << endl;
404  os << " xmlns:ogr=\"http://ogr.maptools.org/\"" << endl;
405  os << " xmlns:gml=\"http://www.opengis.net/gml\">" << endl;
406  os << " <gml:boundedBy>" << endl;
407  os << " <gml:Box>" << endl;
408  os << " <gml:coord>";
409  os << "<gml:X>" <<
410  setprecision(15) << mpolygon->getEnvelopeInternal()->getMinX() <<
411  "</gml:X>";
412  os << "<gml:Y>" <<
413  setprecision(15) << mpolygon->getEnvelopeInternal()->getMinY() <<
414  "</gml:Y>";
415  os << "</gml:coord>" << endl;
416  os << " <gml:coord>";
417  os << "<gml:X>" <<
418  setprecision(15) << mpolygon->getEnvelopeInternal()->getMaxX() <<
419  "</gml:X>";
420  os << "<gml:Y>" <<
421  setprecision(15) << mpolygon->getEnvelopeInternal()->getMaxY() <<
422  "</gml:Y>";
423  os << "</gml:coord>" << endl;
424  os << " </gml:Box>" << endl;
425  os << " </gml:boundedBy>" << endl << endl;
426  os << " <gml:featureMember>" << endl;
427  os << " <ogr:multi_polygon fid=\"0\">" << endl;
428  os << " <ogr:ID>" << idString << "</ogr:ID>" << endl;
429  os << " <ogr:geometryProperty><gml:MultiPolygon><gml:polygonMember>" <<
430  "<gml:Polygon><gml:outerBoundaryIs><gml:LinearRing><gml:coordinates>";
431 
432 
433  for(unsigned int polyN = 0; polyN < mpolygon->getNumGeometries(); polyN++) {
434  geos::geom::CoordinateSequence *pts = mpolygon->getGeometryN(polyN)->getCoordinates();
435 
436  for(unsigned int i = 0; i < pts->getSize(); i++) {
437  double lon = pts->getAt(i).x;
438  double lat = pts->getAt(i).y;
439 
440  os << setprecision(15) << lon << "," << setprecision(15) << lat << " ";
441  }
442  }
443 
444  os << "</gml:coordinates></gml:LinearRing></gml:outerBoundaryIs>" <<
445  "</gml:Polygon></gml:polygonMember></gml:MultiPolygon>" <<
446  "</ogr:geometryProperty>" << endl;
447  os << " </ogr:multi_polygon>" << endl;
448  os << " </gml:featureMember>" << endl;
449  os << "</ogr:FeatureCollection>";
450 
451  return os.str().c_str();
452  }
453 
454 
461  QString PolygonTools::GMLSchema() {
462 
463  ostringstream os;
464 
465  os << "<?xml version=\"1.0\" encoding=\"UTF-8\"?>" << endl;
466  os << "<xs:schema targetNamespace=\"http://ogr.maptools.org/\" "
467  "xmlns:ogr=\"http://ogr.maptools.org/\" "
468  "xmlns:xs=\"http://www.w3.org/2001/XMLSchema\" "
469  "xmlns:gml=\"http://www.opengis.net/gml\" "
470  "elementFormDefault=\"qualified\" "
471  "version=\"1.0\">" << endl;
472  os << " <xs:import namespace=\"http://www.opengis.net/gml\" "
473  "schemaLocation=\"http://schemas.opengis.net/gml/2.1.2/feature.xsd\"/>" << endl;
474  os << " <xs:element name=\"FeatureCollection\" "
475  "type=\"ogr:FeatureCollectionType\" "
476  "substitutionGroup=\"gml:_FeatureCollection\"/>" << endl;
477  os << " <xs:complexType name=\"FeatureCollectionType\">" << endl;
478  os << " <xs:complexContent>" << endl;
479  os << " <xs:extension base=\"gml:AbstractFeatureCollectionType\">" << endl;
480  os << " <xs:attribute name=\"lockId\" type=\"xs:string\" use=\"optional\"/>" << endl;
481  os << " <xs:attribute name=\"scope\" type=\"xs:string\" use=\"optional\"/>" << endl;
482  os << " </xs:extension>" << endl;
483  os << " </xs:complexContent>" << endl;
484  os << " </xs:complexType>" << endl;
485  os << " <xs:element name=\"multi_polygon\" "
486  "type=\"ogr:multi_polygon_Type\" "
487  "substitutionGroup=\"gml:_Feature\"/>" << endl;
488  os << " <xs:complexType name=\"multi_polygon_Type\">" << endl;
489  os << " <xs:complexContent>" << endl;
490  os << " <xs:extension base=\"gml:AbstractFeatureType\">" << endl;
491  os << " <xs:sequence>" << endl;
492  os << " <xs:element name=\"geometryProperty\" "
493  "type=\"gml:MultiPolygonPropertyType\" "
494  "nillable=\"true\" minOccurs=\"0\" maxOccurs=\"1\"/>" << endl;
495  os << " <xs:element name=\"fid\" nillable=\"true\" minOccurs=\"0\" "
496  "maxOccurs=\"1\">" << endl;
497  os << " <xs:simpleType>" << endl;
498  os << " <xs:restriction base=\"xs:string\">" << endl;
499  os << " </xs:restriction>" << endl;
500  os << " </xs:simpleType>" << endl;
501  os << " </xs:element>" << endl;
502  os << " <xs:element name=\"ID\" nillable=\"true\" minOccurs=\"0\" "
503  "maxOccurs=\"1\">" << endl;
504  os << " <xs:simpleType>" << endl;
505  os << " <xs:restriction base=\"xs:integer\">" << endl;
506  os << " <xs:totalDigits value=\"16\"/>" << endl;
507  os << " </xs:restriction>" << endl;
508  os << " </xs:simpleType>" << endl;
509  os << " </xs:element>" << endl;
510  os << " </xs:sequence>" << endl;
511  os << " </xs:extension>" << endl;
512  os << " </xs:complexContent>" << endl;
513  os << " </xs:complexType>" << endl;
514  os << "</xs:schema>";
515 
516  return os.str().c_str();
517  }
518 
519 
527  geos::geom::MultiPolygon *PolygonTools::To180(geos::geom::MultiPolygon *poly360) {
528  try {
529  // Let's take the 360 pieces that sit between 180 and 360 and move them
530  // to -180 to 0. To accomplish this, make a poly that fits 180 -> 360
531  // degrees longitude and intersect (this gives us the pieces that sit
532  // >180). Move this intersection to the left. Then make a poly that fits
533  // 0 to 180 and intersect with the original. These two combined are the
534  // result.
535  geos::geom::CoordinateArraySequence *leftOf180Pts =
536  new geos::geom::CoordinateArraySequence();
537  leftOf180Pts->add(geos::geom::Coordinate(0, -90));
538  leftOf180Pts->add(geos::geom::Coordinate(0, 90));
539  leftOf180Pts->add(geos::geom::Coordinate(180, 90));
540  leftOf180Pts->add(geos::geom::Coordinate(180, -90));
541  leftOf180Pts->add(geos::geom::Coordinate(0, -90));
542 
543  geos::geom::LinearRing *leftOf180Geom =
544  globalFactory->createLinearRing(leftOf180Pts);
545 
546  geos::geom::Polygon *leftOf180Poly =
547  globalFactory->createPolygon(leftOf180Geom, NULL);
548 
549  geos::geom::CoordinateArraySequence *rightOf180Pts =
550  new geos::geom::CoordinateArraySequence();
551  rightOf180Pts->add(geos::geom::Coordinate(180, -90));
552  rightOf180Pts->add(geos::geom::Coordinate(180, 90));
553  rightOf180Pts->add(geos::geom::Coordinate(360, 90));
554  rightOf180Pts->add(geos::geom::Coordinate(360, -90));
555  rightOf180Pts->add(geos::geom::Coordinate(180, -90));
556 
557  geos::geom::LinearRing *rightOf180Geom =
558  globalFactory->createLinearRing(rightOf180Pts);
559 
560  geos::geom::Polygon *rightOf180Poly =
561  globalFactory->createPolygon(rightOf180Geom, NULL);
562 
563  geos::geom::Geometry *preserved = Intersect(leftOf180Poly, poly360);
564  geos::geom::Geometry *moving = Intersect(rightOf180Poly, poly360);
565 
566  geos::geom::CoordinateSequence *movingPts = moving->getCoordinates();
567  geos::geom::CoordinateSequence *movedPts =
568  new geos::geom::CoordinateArraySequence();
569 
570  for(unsigned int i = 0; i < movingPts->getSize(); i ++) {
571  movedPts->add(geos::geom::Coordinate(movingPts->getAt(i).x - 360.0,
572  movingPts->getAt(i).y));
573  }
574 
575  if(movedPts->getSize()) {
576  movedPts->add(geos::geom::Coordinate(movedPts->getAt(0).x,
577  movedPts->getAt(0).y));
578  }
579 
580  geos::geom::Geometry *moved = globalFactory->createPolygon(
581  globalFactory->createLinearRing(movedPts), NULL);
582 
583  std::vector<geos::geom::Geometry *> *geomsForCollection = new
584  std::vector<geos::geom::Geometry *>;
585  geomsForCollection->push_back(preserved);
586  geomsForCollection->push_back(moved);
587 
588  geos::geom::GeometryCollection *the180Polys =
589  Isis::globalFactory->createGeometryCollection(geomsForCollection);
590 
591  geos::geom::MultiPolygon *result = MakeMultiPolygon(the180Polys);
592  delete the180Polys;
593  the180Polys = NULL;
594 
595  geos::geom::MultiPolygon *fixedResult = FixSeam(result);
596  delete result;
597  result = NULL;
598 
599  return fixedResult;
600  }
601  catch(geos::util::GEOSException *exc) {
602  IString msg = "Conversion to 180 failed. The reason given was [" +
603  IString(exc->what()) + "]";
604  delete exc;
605  throw IException(IException::Programmer, msg, _FILEINFO_);
606  }
607  catch(...) {
608  IString msg = "Conversion to 180 failed. Could not determine the reason";
609  throw IException(IException::Programmer, msg, _FILEINFO_);
610  }
611  }
612 
613 
624  double PolygonTools::Thickness(const geos::geom::MultiPolygon *mpolygon) {
625  const geos::geom::Envelope *envelope = mpolygon->getEnvelopeInternal();
626 
627  double x = fabs(envelope->getMaxX() - envelope->getMinX());
628  double y = fabs(envelope->getMaxY() - envelope->getMinY());
629  double extent = max(x, y);
630 
631  return mpolygon->getArea() / (extent * extent);
632  }
633 
634 
646  geos::geom::MultiPolygon *PolygonTools::Despike(const geos::geom::Geometry *geom) {
647  geos::geom::MultiPolygon *spiked = MakeMultiPolygon(geom);
648  geos::geom::MultiPolygon *despiked = Despike(spiked);
649 
650  delete spiked;
651  spiked = NULL;
652 
653  return despiked;
654  }
655 
656 
667  geos::geom::MultiPolygon *PolygonTools::Despike(const geos::geom::MultiPolygon *multiPoly) {
668  // Despike each polygon in the multipolygon
669  vector<geos::geom::Geometry *> *newPolys = new vector<geos::geom::Geometry *>;
670  for(unsigned int g = 0; g < multiPoly->getNumGeometries(); ++g) {
671  const geos::geom::Polygon *poly =
672  dynamic_cast<const geos::geom::Polygon *>(multiPoly->getGeometryN(g));
673 
674  // Despike each hole inside this polygon
675  vector<geos::geom::Geometry *> *holes = new vector<geos::geom::Geometry *>;
676  for(unsigned int h = 0; h < poly->getNumInteriorRing(); ++h) {
677  const geos::geom::LineString *ls = poly->getInteriorRingN(h);
678  geos::geom::LinearRing *lr;
679 
680  // If the hole is not valid fix it
681  // If the hole is NOT valid despike it
682  lr = Despike(ls);
683 
684  if(!lr->isValid()) {
685  geos::geom::LinearRing *fixed = FixGeometry(lr);
686  delete lr;
687  lr = fixed;
688  }
689 
690  // Save this hole if it is not too small
691  if(!lr->isEmpty()) {
692  holes->push_back(lr);
693  }
694  else {
695  delete lr;
696  }
697  } // End holes loop
698 
699  // Work on the main polygon
700  const geos::geom::LineString *ls = poly->getExteriorRing();
701  geos::geom::LinearRing *lr;
702 
703  lr = Despike(ls);
704 
705  try {
706  if(!lr->isValid()) {
707  geos::geom::LinearRing *fixed = FixGeometry(lr);
708  delete lr;
709  lr = fixed;
710  }
711  }
712  catch(IException &e) {
713  // Sometimes despike and fix fail, but the input is really valid. We can just go
714  // with the non-despiked polygon.
715  if(ls->isValid() && ls->getGeometryTypeId() == geos::geom::GEOS_LINEARRING) {
716  lr = dynamic_cast<geos::geom::LinearRing *>(ls->clone());
717  }
718  else {
719  throw;
720  }
721  }
722 
723  // Create a new polygon with the holes and save it
724  if(!lr->isEmpty()) {
725  geos::geom::Polygon *tp = Isis::globalFactory->createPolygon(lr, holes);
726 
727  if(tp->isEmpty() || !tp->isValid()) {
728  delete tp;
729  newPolys->push_back(poly->clone());
730  }
731  else {
732  newPolys->push_back(tp);
733  }
734  }
735  } // End polygons loop
736 
737  // Create a new multipoly from the polygon(s)
738  geos::geom::MultiPolygon *mp = Isis::globalFactory->createMultiPolygon(newPolys);
739 
740  if(!mp->isValid() || mp->isEmpty()) {
741  delete mp;
742  mp = NULL;
743  IString msg = "Despike failed to correct the polygon";
744  throw IException(IException::Programmer, msg, _FILEINFO_);
745  }
746 
747  // if multipoly area changes more than 25% we did something bad to the multipolygon
748  if(fabs((mp->getArea() / multiPoly->getArea()) - 1.0) > 0.50) {
749  IString msg = "Despike failed to correct the polygon " + mp->toString();
750  delete mp;
751 
752  throw IException(IException::Programmer, msg, _FILEINFO_);
753  }
754 
755  return mp;
756  }
757 
758 
772  geos::geom::LinearRing *PolygonTools::Despike(const geos::geom::LineString *lineString) {
773  geos::geom::CoordinateSequence *vertices = lineString->getCoordinates();
774 
775  // We need a full polygon to despike = 3 points (+1 for end==first) = at least 4 points
776  if(vertices->getSize() < 4) {
777  delete vertices;
778  return Isis::globalFactory->createLinearRing(geos::geom::CoordinateArraySequence());
779  }
780 
781  // delete one of the duplicate first/end coordinates,
782  // spikes can occur here and the duplicate points throws off the test
783  vertices->deleteAt(vertices->getSize() - 1);
784 
785  // Index will become negative if our first points are spiked, so we need
786  // to make it big enough to encapsulate the (unsigned int) from getSize()
787  // and allow negative numbers.
788  for(long index = 0; index < (long)vertices->getSize(); index++) {
789  // If we're under 3 vertices, we've despiked all the points out of the polygon :(
790  if(vertices->getSize() < 3) {
791  break;
792  }
793 
794  // These are the array indices that we're going to test for a spike,
795  // the middle one is the one that's spiked if IsSpiked(...) is true.
796  long testCoords[3] = {
797  index - 1,
798  index,
799  index + 1
800  };
801 
802  // Make sure the index is inside of our coordinates array (we'll have both too small/too
803  // big of an index)
804  for(int j = 0; j < 3; j++) {
805  while(testCoords[j] < 0) {
806  testCoords[j] += vertices->getSize();
807  }
808 
809  while(testCoords[j] >= (long)vertices->getSize()) {
810  testCoords[j] -= vertices->getSize();
811  }
812  }
813 
814  // Test the middle point for a spike
815  if(IsSpiked(vertices->getAt(testCoords[0]),
816  vertices->getAt(testCoords[1]),
817  vertices->getAt(testCoords[2]))) {
818  // It's spiked, delete it
819  vertices->deleteAt(testCoords[1]);
820 
821  // Back up to the first test that is affected by this change
822  index -= 2;
823  }
824  }
825 
826  if(vertices->getSize() < 3) {
827  delete vertices;
828  vertices = NULL;
829 
830  return Isis::globalFactory->createLinearRing(geos::geom::CoordinateArraySequence());
831  }
832  else {
833  // Duplicate the first vertex as the last to close the polygon
834  vertices->add(vertices->getAt(0));
835  return Isis::globalFactory->createLinearRing(vertices);
836  }
837  }
838 
848  bool PolygonTools::IsSpiked(geos::geom::Coordinate first,
849  geos::geom::Coordinate middle,
850  geos::geom::Coordinate last) {
851  return TestSpiked(first, middle, last) || TestSpiked(last, middle, first);
852  }
853 
882  bool PolygonTools::TestSpiked(geos::geom::Coordinate first, geos::geom::Coordinate middle,
883  geos::geom::Coordinate last) {
884  geos::geom::Point *firstPt = Isis::globalFactory->createPoint(first);
885  geos::geom::Point *middlePt = Isis::globalFactory->createPoint(middle);
886  geos::geom::Point *lastPt = Isis::globalFactory->createPoint(last);
887 
888  geos::geom::CoordinateSequence *coords = new geos::geom::CoordinateArraySequence();
889  coords->add(first);
890  coords->add(middle);
891  geos::geom::LineString *line = Isis::globalFactory->createLineString(coords); // line takes ownership
892 
893  // The lower the tolerance, the less this algorithm removes and thus
894  // the better chance of success in findimageoverlaps. However, if you
895  // lower the tolerance then there is also a greater chance of programs
896  // such as autoseed failing. 1% is the current tolerance.
897  double tolerance = line->getLength() / 100.0;
898 
899  bool spiked = true;
900 
901  double distanceLastMiddle = geos::operation::distance::DistanceOp::distance(lastPt, middlePt);
902  double distanceLastLine = geos::operation::distance::DistanceOp::distance(lastPt, line);
903 
904  if(distanceLastMiddle == 0.0) return true; // get rid of same point
905 
906  // Checks the ratio of the distance between the last point and the line, and the last point
907  // and the middle point if the ratio is very small then there is a spike
908  if(distanceLastLine / distanceLastMiddle >= .05) {
909  spiked = false;
910  }
911 
912  // If the point is away from the line, keep it
913  if(spiked && distanceLastLine > tolerance) {
914  spiked = false;
915  }
916 
917  if(!spiked) {
918  geos::geom::CoordinateSequence *coords = new geos::geom::CoordinateArraySequence();
919  coords->add(first);
920  coords->add(middle);
921  coords->add(last);
922  coords->add(first);
923 
924  // shell takes ownership of coords
925  geos::geom::LinearRing *shell = Isis::globalFactory->createLinearRing(coords);
926  std::vector<geos::geom::Geometry *> *empty = new std::vector<geos::geom::Geometry *>;
927 
928  // poly takes ownership of shell and empty
929  geos::geom::Polygon *poly = Isis::globalFactory->createPolygon(shell, empty);
930 
931 
932  // if these 3 points define a straight line then the middle is worthless (defines nothing)
933  // or problematic
934  if(poly->getArea() < 1.0e-10) {
935  spiked = true;
936  }
937 
938  delete poly;
939  }
940 
941 
942  delete firstPt;
943  delete middlePt;
944  delete lastPt;
945  delete line;
946 
947  return spiked;
948  }
949 
950 
964  geos::geom::Geometry *PolygonTools::Intersect(const geos::geom::Geometry *geom1,
965  const geos::geom::Geometry *geom2) {
966  try {
967  return Operate(geom1, geom2, (unsigned int)geos::operation::overlay::OverlayOp::opINTERSECTION);
968  }
969  catch(geos::util::GEOSException *exc) {
970  IString msg = "Intersect operation failed. The reason given was [" + IString(exc->what()) + "]";
971  delete exc;
972  throw IException(IException::Programmer, msg, _FILEINFO_);
973  }
974  catch(IException &e) {
975  IString msg = "Intersect operation failed";
976  throw IException(e, IException::Programmer, msg, _FILEINFO_);
977  }
978  catch(...) {
979  IString msg = "Intersect operation failed for an unknown reason";
980  throw IException(IException::Programmer, msg, _FILEINFO_);
981  }
982  }
983 
984 
985  geos::geom::Geometry *PolygonTools::Operate(const geos::geom::Geometry *geom1,
986  const geos::geom::Geometry *geom2,
987  unsigned int opcode) {
988 
989  geos::operation::overlay::OverlayOp::OpCode code =
990  (geos::operation::overlay::OverlayOp::OpCode)opcode;
991 
992  geos::geom::Geometry *result = NULL;
993  bool failed = true;
994  geos::geom::Geometry *geomFirst = MakeMultiPolygon(geom1);
995  geos::geom::Geometry *geomSecond = MakeMultiPolygon(geom2);
996 
997  geos::operation::overlay::snap::GeometrySnapper snap(*geomFirst);
998  geos::geom::Geometry *geomSnapped = snap.snapTo(*geomSecond, 1.0e-10)->clone();
999  if(!geomSnapped->isValid()) {
1000  delete geomSnapped;
1001  }
1002  else {
1003  delete geomFirst;
1004  geomFirst = geomSnapped;
1005  }
1006 
1007  unsigned int precision = 15;
1008  unsigned int minPrecision = 13;
1009  while(failed) {
1010  try {
1011  // C++11: the geos BinaryOp returns an auto_ptr, we use release() to create a unique_ptr.
1012  std::unique_ptr< geos::geom::Geometry > resultAuto(
1013  BinaryOp(geomFirst, geomSecond, geos::operation::overlay::overlayOp(code)).release());
1014  failed = false;
1015  result = resultAuto->clone();
1016  }
1017  catch(geos::util::GEOSException *exc) {
1018  // Just in case the clone failed....
1019  if(!failed || precision == minPrecision) throw;
1020 
1021  delete exc;
1022  }
1023  catch(...) {
1024  if(precision == minPrecision) {
1025  IString msg = "An unknown geos error occurred";
1026  throw IException(IException::Programmer, msg, _FILEINFO_);
1027  }
1028 
1029  if(!failed) {
1030  IString msg = "An unknown geos error occurred when attempting to clone a geometry";
1031  throw IException(IException::Programmer, msg, _FILEINFO_);
1032  }
1033  }
1034 
1035  // try reducing the precision
1036  if(failed) {
1037  precision --;
1038  geos::geom::Geometry *tmp = ReducePrecision(geomFirst, precision);
1039  delete geomFirst;
1040  geomFirst = tmp;
1041 
1042  tmp = ReducePrecision(geomSecond, precision);
1043  delete geomSecond;
1044  geomSecond = tmp;
1045  }
1046  }
1047 
1048  delete geomFirst;
1049  geomFirst = NULL;
1050 
1051  delete geomSecond;
1052  geomSecond = NULL;
1053 
1054  if(result && !result->isValid()) {
1055  try {
1056  geos::geom::Geometry *newResult = FixGeometry(result);
1057 
1058  if(fabs(newResult->getArea() / result->getArea() - 1.0) > 0.50) {
1059  delete newResult;
1060  delete result;
1061 
1062  IString msg = "Operation [" + IString((int)opcode) + "] failed";
1063  throw IException(IException::Programmer, msg, _FILEINFO_);
1064  }
1065 
1066  delete result;
1067  result = newResult;
1068  }
1069  catch(IException &e) {
1070  IString msg = "Operation [" + IString((int)opcode) + "] failed";
1071  throw IException(IException::Programmer, msg, _FILEINFO_);
1072  }
1073  }
1074 
1075  if(result == NULL) {
1076  IString msg = "Operation [" + IString((int)opcode) + " failed";
1077  throw IException(IException::Programmer, msg, _FILEINFO_);
1078  }
1079 
1080  return result;
1081  }
1082 
1093  geos::geom::Geometry *PolygonTools::FixGeometry(const geos::geom::Geometry *geom) {
1094  if(geom->getGeometryTypeId() == geos::geom::GEOS_MULTIPOLYGON) {
1095  return FixGeometry(dynamic_cast<const geos::geom::MultiPolygon *>(geom));
1096  }
1097  if(geom->getGeometryTypeId() == geos::geom::GEOS_LINEARRING) {
1098  return FixGeometry(dynamic_cast<const geos::geom::LinearRing *>(geom));
1099  }
1100  if(geom->getGeometryTypeId() == geos::geom::GEOS_POLYGON) {
1101  return FixGeometry(dynamic_cast<const geos::geom::Polygon *>(geom));
1102  }
1103  if(geom->getGeometryTypeId() == geos::geom::GEOS_GEOMETRYCOLLECTION) {
1104  return FixGeometry(MakeMultiPolygon(geom));
1105  }
1106  else {
1107  IString msg = "PolygonTools::FixGeometry does not support [" + GetGeometryName(geom) + "]";
1108  throw IException(IException::Programmer, msg, _FILEINFO_);
1109  }
1110  }
1111 
1121  geos::geom::MultiPolygon *PolygonTools::FixGeometry(const geos::geom::MultiPolygon *poly) {
1122  // Maybe two points are on top of each other
1123  vector<geos::geom::Geometry *> *newPolys = new vector<geos::geom::Geometry *>;
1124 
1125  // Convert each polygon in this multi-polygon
1126  for(unsigned int geomIndex = 0; geomIndex < poly->getNumGeometries(); geomIndex ++) {
1127  geos::geom::Polygon *fixedpoly = FixGeometry(
1128  dynamic_cast<const geos::geom::Polygon *>(
1129  poly->getGeometryN(geomIndex)));
1130  if(fixedpoly->isValid()) {
1131  newPolys->push_back(fixedpoly);
1132  }
1133  else {
1134  delete fixedpoly;
1135  }
1136  fixedpoly = NULL;
1137  }
1138 
1139  geos::geom::MultiPolygon *mp = Isis::globalFactory->createMultiPolygon(newPolys);
1140  return mp;
1141  }
1142 
1143 
1144  geos::geom::Polygon *PolygonTools::FixGeometry(const geos::geom::Polygon *poly) {
1145 
1146  // Convert each hole inside this polygon
1147  vector<geos::geom::Geometry *> *holes = new vector<geos::geom::Geometry *>;
1148  for(unsigned int holeIndex = 0; holeIndex < poly->getNumInteriorRing(); holeIndex ++) {
1149  const geos::geom::LineString *thisHole = poly->getInteriorRingN(holeIndex);
1150 
1151  // We hope they are all linear rings (closed/simple), but if not just leave it be
1152  if(thisHole->getGeometryTypeId() != geos::geom::GEOS_LINEARRING) {
1153  holes->push_back(thisHole->clone());
1154  continue;
1155  }
1156 
1157  try {
1158  geos::geom::LinearRing *newHole = FixGeometry((geos::geom::LinearRing *)thisHole);
1159  holes->push_back(newHole);
1160  }
1161  catch (IException &e) {
1162  IString msg = "Failed when attempting to fix interior ring of multipolygon";
1163  throw IException(IException::Programmer, msg, _FILEINFO_);
1164  }
1165  } // end num holes in polygon loop
1166 
1167 
1168  const geos::geom::LineString *exterior = poly->getExteriorRing();
1169 
1170  try {
1171  geos::geom::LinearRing *newExterior = NULL;
1172 
1173  if(exterior->getGeometryTypeId() == geos::geom::GEOS_LINEARRING) {
1174  newExterior = FixGeometry((geos::geom::LinearRing *)exterior);
1175  }
1176  else {
1177  IString msg = "Failed when attempting to fix exterior ring of polygon. The exterior "
1178  "ring is not simple and closed";
1179  throw IException(IException::Programmer, msg, _FILEINFO_);
1180  }
1181 
1182  return globalFactory->createPolygon(newExterior, holes);
1183  }
1184  catch (IException &e) {
1185  IString msg = "Failed when attempting to fix exterior ring of polygon";
1186  throw IException(e, IException::Programmer, msg, _FILEINFO_);
1187  }
1188  }
1189 
1190 
1206  geos::geom::LinearRing *PolygonTools::FixGeometry(const geos::geom::LinearRing *ring) {
1207 
1208  geos::geom::CoordinateSequence *coords = ring->getCoordinates();
1209 
1210  // Check this, just in case
1211  if(coords->getSize() < 4) {
1212  return globalFactory->createLinearRing(new geos::geom::DefaultCoordinateSequence());
1213  }
1214 
1215  geos::geom::CoordinateSequence *newCoords = new geos::geom::DefaultCoordinateSequence();
1216  const geos::geom::Coordinate *lastCoordinate = &coords->getAt(0);
1217  newCoords->add(*lastCoordinate);
1218 
1219  // Convert each coordinate in this hole
1220  for(unsigned int coordIndex = 1; coordIndex < coords->getSize() - 1; coordIndex ++) {
1221  const geos::geom::Coordinate *thisCoordinate = &coords->getAt(coordIndex);
1222 
1223  // we're going to compare the decimal place of the current point to the decimal place
1224  // of the difference, if they are drastically different then geos might not be seeing them
1225  // correctly.
1226  double difference[2] = {
1227  lastCoordinate->x - thisCoordinate->x,
1228  lastCoordinate->y - thisCoordinate->y,
1229  };
1230 
1231  // geos isnt differentiating between points this close
1232  double minDiff = fabs(DecimalPlace(thisCoordinate->x) - DecimalPlace(difference[0]));
1233 
1234  minDiff = min(minDiff, fabs(DecimalPlace(thisCoordinate->y) - DecimalPlace(difference[1])));
1235 
1236  // Cases where the difference in one direction is exactly zero, and the other direction is
1237  // next to zero appear often enough (especially in despike).
1238  if(difference[0] == 0.0 && difference[1] != 0.0) {
1239  // subtracted the two points, got deltaX = 0.0, use the y difference decimal place
1240  minDiff = fabs(DecimalPlace(thisCoordinate->y) - DecimalPlace(difference[1]));
1241  }
1242  else if(difference[1] == 0.0 && difference[0] != 0.0) {
1243  // subtracted the two points, got deltaY = 0.0, use the x difference decimal place
1244  minDiff = fabs(DecimalPlace(thisCoordinate->x) - DecimalPlace(difference[0]));
1245  }
1246  else if(difference[0] == 0.0 && difference[1] == 0.0) {
1247  // subtracted the two points, got 0.0, so it's same point... make sure it gets ignored!
1248  minDiff = 1E99;
1249  }
1250 
1251  // geos has a hard time differentiating when points get too close...
1252  if(minDiff <= 15) {
1253  newCoords->add(*thisCoordinate);
1254  lastCoordinate = thisCoordinate;
1255  }
1256  } // end num coords in hole loop
1257 
1258  newCoords->add(geos::geom::Coordinate(newCoords->getAt(0).x, newCoords->getAt(0).y));
1259  geos::geom::LinearRing *newRing = NULL;
1260 
1261  // Now that we've weeded out any bad coordinates, let's rebuild the geometry
1262  try {
1263  if(newCoords->getSize() > 3) {
1264  newRing = globalFactory->createLinearRing(newCoords);
1265  }
1266  else {
1267  delete newCoords;
1268  newCoords = NULL;
1269  }
1270  }
1271  catch(geos::util::GEOSException *exc) {
1272  delete exc;
1273  exc = NULL;
1274 
1275  IString msg = "Error when attempting to fix linear ring";
1276  throw IException(IException::Programmer, msg, _FILEINFO_);
1277  }
1278 
1279  if(newRing && !newRing->isValid() && ring->isValid()) {
1280  IString msg = "Failed when attempting to fix linear ring";
1281  throw IException(IException::Programmer, msg, _FILEINFO_);
1282  }
1283  else if(!newRing || !newRing->isValid()) {
1284  if(newRing) {
1285  delete newRing;
1286  }
1287 
1288  newRing = dynamic_cast<geos::geom::LinearRing *>(ring->clone());
1289  }
1290 
1291  return newRing;
1292  }
1293 
1294 
1306  int PolygonTools::DecimalPlace(double num) {
1307  // 0.0 = decimal place 0
1308  if(num == 0.0) return 0;
1309 
1310  num = fabs(num);
1311 
1312  int decimalPlace = 1;
1313  while(num < 1.0) {
1314  num *= 10.0;
1315  decimalPlace --;
1316  }
1317 
1318  while(num > 10.0) {
1319  num /= 10.0;
1320  decimalPlace ++;
1321  }
1322 
1323  return decimalPlace;
1324  }
1325 
1334  geos::geom::Geometry *PolygonTools::Difference(const geos::geom::Geometry *geom1,
1335  const geos::geom::Geometry *geom2) {
1336  try {
1337  return Operate(geom1, geom2, (unsigned int)geos::operation::overlay::OverlayOp::opDIFFERENCE);
1338  }
1339  catch(geos::util::GEOSException *exc) {
1340  IString msg = "Difference operation failed. The reason given was [" +
1341  IString(exc->what()) + "]";
1342  delete exc;
1343  throw IException(IException::Programmer, msg, _FILEINFO_);
1344  }
1345  catch(IException &e) {
1346  IString msg = "Difference operation failed";
1347  throw IException(e, IException::Programmer, msg, _FILEINFO_);
1348  }
1349  catch(...) {
1350  IString msg = "Difference operation failed for an unknown reason";
1351  throw IException(IException::Programmer, msg, _FILEINFO_);
1352  }
1353  }
1354 
1355 
1369  geos::geom::MultiPolygon *PolygonTools::MakeMultiPolygon(const geos::geom::Geometry *geom) {
1370  // The area of the geometry is too small, so just ignore it.
1371  if(geom->isEmpty()) {
1372  return Isis::globalFactory->createMultiPolygon();
1373  }
1374 
1375  else if(geom->getArea() - DBL_EPSILON <= DBL_EPSILON) {
1376  return Isis::globalFactory->createMultiPolygon();
1377  }
1378 
1379  else if(geom->getGeometryTypeId() == geos::geom::GEOS_MULTIPOLYGON) {
1380  return dynamic_cast<geos::geom::MultiPolygon *>(geom->clone());
1381  }
1382 
1383  else if(geom->getGeometryTypeId() == geos::geom::GEOS_POLYGON) {
1384  vector<geos::geom::Geometry *> *polys = new vector<geos::geom::Geometry *>;
1385  polys->push_back(geom->clone());
1386  geos::geom::MultiPolygon *mp = Isis::globalFactory->createMultiPolygon(polys);
1387  return mp;
1388  }
1389 
1390  else if(geom->getGeometryTypeId() == geos::geom::GEOS_GEOMETRYCOLLECTION) {
1391  vector<geos::geom::Geometry *> * polys =
1392  new vector<geos::geom::Geometry *>;
1393  const geos::geom::GeometryCollection *gc =
1394  dynamic_cast<const geos::geom::GeometryCollection *>(geom);
1395  for(unsigned int i = 0; i < gc->getNumGeometries(); ++i) {
1396  geos::geom::MultiPolygon *subMultiPoly =
1397  MakeMultiPolygon(gc->getGeometryN(i));
1398 
1399  for(unsigned int subPoly = 0;
1400  subPoly < subMultiPoly->getNumGeometries();
1401  subPoly ++) {
1402  const geos::geom::Polygon *poly =
1403  dynamic_cast<const geos::geom::Polygon *>(
1404  subMultiPoly->getGeometryN(subPoly));
1405  polys->push_back(dynamic_cast<geos::geom::Polygon *>(poly->clone()));
1406  }
1407  }
1408 
1409  geos::geom::MultiPolygon *mp =
1410  Isis::globalFactory->createMultiPolygon(polys);
1411  if(mp->getArea() - DBL_EPSILON <= DBL_EPSILON) {
1412  delete mp;
1413  mp = Isis::globalFactory->createMultiPolygon();
1414  }
1415 
1416  return mp;
1417  }
1418 
1419  // All other geometry types are invalid so ignore them
1420  else {
1421  return Isis::globalFactory->createMultiPolygon();
1422  }
1423  }
1424 
1425 
1426  geos::geom::MultiPolygon *PolygonTools::FixSeam(
1427  const geos::geom::Polygon *polyA, const geos::geom::Polygon *polyB) {
1428  geos::geom::CoordinateSequence *polyAPoints = polyA->getCoordinates();
1429  geos::geom::CoordinateSequence *polyBPoints = polyB->getCoordinates();
1430 
1431  unsigned int aIntersectionBegin = 0;
1432  unsigned int aIntersectionEnd = 0;
1433  unsigned int bIntersectionBegin = 0;
1434  unsigned int bIntersectionEnd = 0;
1435 
1436  bool intersectionStarted = false;
1437  bool intersectionEnded = false;
1438 
1439  unsigned int lastBMatch = 0;
1440  for (unsigned int i = 0;
1441  !intersectionEnded && i < polyAPoints->getSize();
1442  i++) {
1443 
1444  bool foundEquivalent = false;
1445 
1446  geos::geom::Coordinate coordA = polyAPoints->getAt(i);
1447  coordA = *ReducePrecision(&coordA, 13);
1448 
1449  for (unsigned int j = 0;
1450  !foundEquivalent && j < polyBPoints->getSize();
1451  j++) {
1452  geos::geom::Coordinate coordB = polyBPoints->getAt(j);
1453  coordB = *ReducePrecision(&coordB, 13);
1454 
1455  foundEquivalent = coordA.equals(coordB);
1456 
1457  if (foundEquivalent) lastBMatch = j;
1458 
1459  if (foundEquivalent && !intersectionStarted) {
1460  intersectionStarted = true;
1461  aIntersectionBegin = i;
1462  bIntersectionBegin = j;
1463  }
1464  }
1465 
1466  if (!foundEquivalent && intersectionStarted && !intersectionEnded) {
1467  intersectionEnded = true;
1468  aIntersectionEnd = i;
1469  bIntersectionEnd = lastBMatch;
1470  }
1471  }
1472 
1473  geos::geom::MultiPolygon * result = NULL;
1474  if (intersectionStarted && intersectionEnded) {
1475  geos::geom::CoordinateSequence *merged =
1476  new geos::geom::CoordinateArraySequence;
1477 
1478  unsigned int i = 0;
1479  for (i = 0; i < aIntersectionBegin; i ++) {
1480  merged->add(polyAPoints->getAt(i));
1481  }
1482 
1483  i = bIntersectionBegin;
1484  while (i != bIntersectionEnd) {
1485  merged->add(polyBPoints->getAt(i));
1486  i++;
1487  if (i >= polyBPoints->getSize()) i = 0;
1488  }
1489 
1490  for (i = aIntersectionEnd; i < polyAPoints->getSize() - 1; i++) {
1491  merged->add(polyAPoints->getAt(i));
1492  }
1493 
1494  merged->add(merged->getAt(0));
1495  result = MakeMultiPolygon(globalFactory->createPolygon(
1496  globalFactory->createLinearRing(merged), NULL));
1497  }
1498 
1499  return result;
1500  }
1501 
1502 
1503  geos::geom::MultiPolygon *PolygonTools::FixSeam(
1504  const geos::geom::MultiPolygon *poly) {
1505 
1506  std::vector<geos::geom::Geometry *> *polys =
1507  new std::vector<geos::geom::Geometry *>;
1508 
1509 
1510  for(unsigned int copyIndex = 0;
1511  copyIndex < poly->getNumGeometries();
1512  copyIndex ++) {
1513  polys->push_back(poly->getGeometryN(copyIndex)->clone());
1514  }
1515 
1516  unsigned int outerPolyIndex = 0;
1517 
1518  while(outerPolyIndex + 1 < polys->size()) {
1519  unsigned int innerPolyIndex = outerPolyIndex + 1;
1520 
1521  while(innerPolyIndex < polys->size()) {
1522  geos::geom::MultiPolygon *fixedPair = FixSeam(
1523  dynamic_cast<geos::geom::Polygon *>(polys->at(outerPolyIndex)),
1524  dynamic_cast<geos::geom::Polygon *>(polys->at(innerPolyIndex)));
1525 
1526  if(fixedPair != NULL) {
1527  geos::geom::Geometry *oldInnerPoly = polys->at(innerPolyIndex);
1528  geos::geom::Geometry *oldOuterPoly = polys->at(outerPolyIndex);
1529 
1530  polys->erase(polys->begin() + innerPolyIndex);
1531  (*polys)[outerPolyIndex] = fixedPair->getGeometryN(0)->clone();
1532  innerPolyIndex = outerPolyIndex + 1;
1533 
1534  delete fixedPair;
1535  fixedPair = NULL;
1536 
1537  delete oldInnerPoly;
1538  oldInnerPoly = NULL;
1539 
1540  delete oldOuterPoly;
1541  oldOuterPoly = NULL;
1542  }
1543  else {
1544  innerPolyIndex ++;
1545  }
1546  }
1547 
1548  outerPolyIndex ++;
1549  }
1550 
1551  return globalFactory->createMultiPolygon(polys);
1552  }
1553 
1554 
1564  geos::geom::Geometry *PolygonTools::ReducePrecision(const geos::geom::Geometry *geom,
1565  unsigned int precision) {
1566  if(geom->getGeometryTypeId() == geos::geom::GEOS_MULTIPOLYGON) {
1567  return ReducePrecision(
1568  dynamic_cast<const geos::geom::MultiPolygon *>(geom), precision);
1569  }
1570  if(geom->getGeometryTypeId() == geos::geom::GEOS_LINEARRING) {
1571  return ReducePrecision(
1572  dynamic_cast<const geos::geom::LinearRing *>(geom), precision);
1573  }
1574  if(geom->getGeometryTypeId() == geos::geom::GEOS_POLYGON) {
1575  return ReducePrecision(
1576  dynamic_cast<const geos::geom::Polygon *>(geom), precision);
1577  }
1578  if(geom->getGeometryTypeId() == geos::geom::GEOS_GEOMETRYCOLLECTION) {
1579  return ReducePrecision(MakeMultiPolygon(geom), precision);
1580  }
1581  else {
1582  IString msg = "PolygonTools::ReducePrecision does not support [" +
1583  GetGeometryName(geom) + "]";
1584  throw IException(IException::Programmer, msg, _FILEINFO_);
1585  }
1586  }
1587 
1588 
1598  geos::geom::MultiPolygon *PolygonTools::ReducePrecision(const geos::geom::MultiPolygon *poly,
1599  unsigned int precision) {
1600  // Maybe two points are on top of each other
1601  vector<geos::geom::Geometry *> *newPolys = new vector<geos::geom::Geometry *>;
1602 
1603  // Convert each polygon in this multi-polygon
1604  for(unsigned int geomIndex = 0; geomIndex < poly->getNumGeometries(); geomIndex ++) {
1605  geos::geom::Geometry *lowerPrecision = ReducePrecision(
1606  dynamic_cast<const geos::geom::Polygon *>(
1607  poly->getGeometryN(geomIndex)),
1608  precision);
1609 
1610  if(!lowerPrecision->isEmpty()) {
1611  newPolys->push_back(lowerPrecision);
1612  }
1613  else {
1614  delete lowerPrecision;
1615  }
1616  }
1617 
1618  geos::geom::MultiPolygon *mp = Isis::globalFactory->createMultiPolygon(newPolys);
1619  return mp;
1620  }
1621 
1622 
1632  geos::geom::Polygon *PolygonTools::ReducePrecision(const geos::geom::Polygon *poly,
1633  unsigned int precision) {
1634  // Convert each hole inside this polygon
1635  vector<geos::geom::Geometry *> *holes = new vector<geos::geom::Geometry *>;
1636  for(unsigned int holeIndex = 0; holeIndex < poly->getNumInteriorRing(); holeIndex ++) {
1637  const geos::geom::LineString *thisHole = poly->getInteriorRingN(holeIndex);
1638 
1639  // We hope they are all linear rings (closed/simple), but if not just leave it be
1640  if(thisHole->getGeometryTypeId() != geos::geom::GEOS_LINEARRING) {
1641  holes->push_back(thisHole->clone());
1642  continue;
1643  }
1644 
1645  try {
1646  geos::geom::LinearRing *newHole = ReducePrecision((geos::geom::LinearRing *)thisHole,
1647  precision);
1648 
1649  if(!newHole->isEmpty()) {
1650  holes->push_back(newHole);
1651  }
1652  else {
1653  delete newHole;
1654  }
1655 
1656  }
1657  catch(IException &e) {
1658  IString msg = "Failed when attempting to fix interior ring of multipolygon";
1659  throw IException(e, IException::Programmer, msg, _FILEINFO_);
1660  }
1661  } // end num holes in polygon loop
1662 
1663 
1664  const geos::geom::LineString *exterior = poly->getExteriorRing();
1665 
1666  try {
1667  geos::geom::LinearRing *newExterior = NULL;
1668 
1669  if(exterior->getGeometryTypeId() == geos::geom::GEOS_LINEARRING) {
1670  newExterior = ReducePrecision((geos::geom::LinearRing *)exterior, precision);
1671  }
1672  else {
1673  IString msg = "Failed when attempting to fix exterior ring of polygon. The exterior "
1674  "ring is not simple and closed";
1675  throw IException(IException::Programmer, msg, _FILEINFO_);
1676  }
1677 
1678  return globalFactory->createPolygon(newExterior, holes);
1679  }
1680  catch(IException &e) {
1681  IString msg = "Failed when attempting to fix exterior ring of polygon";
1682  throw IException(e, IException::Programmer, msg, _FILEINFO_);
1683  }
1684  }
1685 
1686 
1696  geos::geom::LinearRing *PolygonTools::ReducePrecision(const geos::geom::LinearRing *ring,
1697  unsigned int precision) {
1698  geos::geom::CoordinateSequence *coords = ring->getCoordinates();
1699 
1700  // Check this, just in case
1701  if(coords->getSize() <= 0) return dynamic_cast<geos::geom::LinearRing *>(
1702  ring->clone());
1703 
1704  geos::geom::CoordinateSequence *newCoords = new geos::geom::DefaultCoordinateSequence();
1705  geos::geom::Coordinate *coord = ReducePrecision(&coords->getAt(0), precision);
1706  newCoords->add(*coord);
1707  delete coord;
1708  coord = NULL;
1709 
1710  // Convert each coordinate in this ring
1711  for(unsigned int coordIndex = 1; coordIndex < coords->getSize() - 1; coordIndex ++) {
1712  const geos::geom::Coordinate *thisCoordinate = &coords->getAt(coordIndex);
1713  coord = ReducePrecision(thisCoordinate, precision);
1714  newCoords->add(*coord);
1715  delete coord;
1716  coord = NULL;
1717  }
1718 
1719  newCoords->add(geos::geom::Coordinate(newCoords->getAt(0).x, newCoords->getAt(0).y));
1720  geos::geom::LinearRing *newRing = NULL;
1721 
1722  // Now that we've weeded out any bad coordinates, let's rebuild the geometry
1723  try {
1724  newRing = globalFactory->createLinearRing(newCoords);
1725  }
1726  catch(geos::util::GEOSException *exc) {
1727  delete exc;
1728  exc = NULL;
1729 
1730  IString msg = "Error when attempting to reduce precision of linear ring";
1731  throw IException(IException::Programmer, msg, _FILEINFO_);
1732  }
1733 
1734  // try to despike
1735  try {
1736  geos::geom::LinearRing *tmp = Despike(newRing);
1737  delete newRing;
1738  newRing = tmp;
1739  }
1740  catch(IException &e) {
1741  }
1742 
1743  if(!newRing->isValid()) {
1744  IString msg = "Failed when attempting to reduce precision of linear ring";
1745  throw IException(IException::Programmer, msg, _FILEINFO_);
1746  }
1747 
1748  return newRing;
1749  }
1750 
1751 
1761  geos::geom::Coordinate *PolygonTools::ReducePrecision(const geos::geom::Coordinate *coord,
1762  unsigned int precision) {
1763  return new geos::geom::Coordinate(
1764  ReducePrecision(coord->x, precision),
1765  ReducePrecision(coord->y, precision),
1766  ReducePrecision(coord->z, precision)
1767  );
1768  }
1769 
1770 
1780  double PolygonTools::ReducePrecision(double num, unsigned int precision) {
1781  double result = num;
1782 
1783  // If num == nan then this test will fail
1784  if (num == num) {
1785  int decimalPlace = DecimalPlace(num);
1786  double factor = pow(10.0, (int)decimalPlace);
1787 
1788  // reduced num is in the form 0.nnnnnnnnnn...
1789  double reducedNum = num / factor;
1790 
1791  double cutoff = pow(10.0, (int)precision);
1792  double round_offset = (num < 0) ? -0.5 : 0.5;
1793 
1794  // cast off the digits past the precision's place
1795  reducedNum = ((long long)(reducedNum * cutoff + round_offset)) / cutoff;
1796  result = reducedNum * factor;
1797  }
1798 
1799  return result;
1800  }
1801 
1802 
1811  QString PolygonTools::GetGeometryName(const geos::geom::Geometry *geom) {
1812  switch(geom->getGeometryTypeId()) {
1813  case geos::geom::GEOS_POINT:
1814  return "Point";
1815  case geos::geom::GEOS_LINESTRING:
1816  return "Line String";
1817  case geos::geom::GEOS_LINEARRING:
1818  return "Linear Ring";
1819  case geos::geom::GEOS_POLYGON:
1820  return "Polygon";
1821  case geos::geom::GEOS_MULTIPOINT:
1822  return "Multi Point";
1823  case geos::geom::GEOS_MULTILINESTRING:
1824  return "Multi Line String";
1825  case geos::geom::GEOS_MULTIPOLYGON:
1826  return "Multi Polygon";
1827  case geos::geom::GEOS_GEOMETRYCOLLECTION:
1828  return "Geometry Collection";
1829  default:
1830  return "UNKNOWN";
1831  }
1832  }
1833 
1834 
1835  bool PolygonTools::Equal(const geos::geom::MultiPolygon *poly1,
1836  const geos::geom::MultiPolygon *poly2) {
1837 
1838  vector<const geos::geom::Polygon *> polys1;
1839  vector<const geos::geom::Polygon *> polys2;
1840 
1841  if(poly1->getNumGeometries() != poly2->getNumGeometries()) return false;
1842 
1843  // Convert each polygon in this multi-polygon
1844  for(unsigned int geomIndex = 0; geomIndex < poly1->getNumGeometries(); geomIndex ++) {
1845  polys1.push_back(dynamic_cast<const geos::geom::Polygon *>(
1846  poly1->getGeometryN(geomIndex)));
1847  polys2.push_back(dynamic_cast<const geos::geom::Polygon *>(
1848  poly2->getGeometryN(geomIndex)));
1849  }
1850 
1851  for(int p1 = polys1.size() - 1; (p1 >= 0) && polys1.size(); p1 --) {
1852  for(int p2 = polys2.size() - 1; (p2 >= 0) && polys2.size(); p2 --) {
1853  if(Equal(polys1[p1], polys2[p2])) {
1854  // Delete polys1[p1] by replacing it with the last Polygon in polys1
1855  polys1[p1] = polys1[polys1.size()-1];
1856  polys1.resize(polys1.size() - 1);
1857  // Delete polys2[p2] by replacing it with the last Polygon in polys2
1858  polys2[p2] = polys2[polys2.size()-1];
1859  polys2.resize(polys2.size() - 1);
1860  }
1861  }
1862  }
1863 
1864  return (polys1.size() == 0) && (polys2.size() == 0);
1865  }
1866 
1867 
1868  bool PolygonTools::Equal(const geos::geom::Polygon *poly1, const geos::geom::Polygon *poly2) {
1869  vector<const geos::geom::LineString *> holes1;
1870  vector<const geos::geom::LineString *> holes2;
1871 
1872  if(poly1->getNumInteriorRing() != poly2->getNumInteriorRing()) return false;
1873 
1874  if(!Equal(poly1->getExteriorRing(), poly2->getExteriorRing())) return false;
1875 
1876  // Convert each hole inside this polygon
1877  for(unsigned int holeIndex = 0; holeIndex < poly1->getNumInteriorRing(); holeIndex ++) {
1878 
1879  // We hope they are all linear rings (closed/simple), but if not just leave it be
1880  if(poly1->getInteriorRingN(holeIndex)->getGeometryTypeId() == geos::geom::GEOS_LINESTRING) {
1881  holes1.push_back(poly1->getInteriorRingN(holeIndex));
1882  }
1883 
1884  if(poly2->getInteriorRingN(holeIndex)->getGeometryTypeId() == geos::geom::GEOS_LINESTRING) {
1885  holes2.push_back(poly2->getInteriorRingN(holeIndex));
1886  }
1887 
1888  }
1889 
1890  if(holes1.size() != holes2.size()) return false;
1891 
1892  for(int h1 = holes1.size() - 1; (h1 >= 0) && holes1.size(); h1 --) {
1893  for(int h2 = holes2.size() - 1; (h2 >= 0) && holes2.size(); h2 --) {
1894  if(Equal(holes1[h1], holes2[h2])) {
1895  // Delete holes1[h1] by replacing it with the last Polygon in holes1
1896  holes1[h1] = holes1[holes1.size()-1];
1897  holes1.resize(holes1.size() - 1);
1898  // Delete holes2[h2] by replacing it with the last Polygon in holes2
1899  holes2[h2] = holes2[holes2.size()-1];
1900  holes2.resize(holes2.size() - 1);
1901  }
1902  }
1903  }
1904 
1905  return (holes1.size() == 0) && (holes2.size() == 0);
1906  }
1907 
1908 
1909  bool PolygonTools::Equal(const geos::geom::LineString *lineString1,
1910  const geos::geom::LineString *lineString2) {
1911 
1912  geos::geom::CoordinateSequence *coords1 = lineString1->getCoordinates();
1913  geos::geom::CoordinateSequence *coords2 = lineString2->getCoordinates();
1914  bool isEqual = true;
1915 
1916  if(coords1->getSize() != coords2->getSize()) isEqual = false;
1917 
1918  unsigned int index1 = 0;
1919  unsigned int index2 = 0;
1920 
1921  // -1 extra for dupicate start/end coordinates
1922  for(; index2 < coords2->getSize() - 1 && isEqual; index2 ++) {
1923  if(Equal(coords1->getAt(index1), coords2->getAt(index2))) break;
1924  }
1925 
1926  if(index2 == coords2->getSize() - 1) isEqual = false;
1927 
1928  for(; index1 < coords1->getSize() - 1 && isEqual; index1 ++, index2 ++) {
1929  if(!Equal(coords1->getAt(index1), coords2->getAt(index2 % (coords2->getSize() - 1)))) {
1930  isEqual = false;
1931  }
1932  }
1933 
1934  delete coords1;
1935  delete coords2;
1936  return isEqual;
1937  }
1938 
1939 
1940  bool PolygonTools::Equal(const geos::geom::Coordinate &coord1,
1941  const geos::geom::Coordinate &coord2) {
1942 
1943  if(!Equal(coord1.x, coord2.x)) return false;
1944  if(!Equal(coord1.y, coord2.y)) return false;
1945  if(!Equal(coord1.y, coord2.y)) return false;
1946 
1947  return true;
1948  }
1949 
1950 
1951  bool PolygonTools::Equal(const double d1, const double d2) {
1952  const double cutoff = 1e15;
1953 
1954  if(DecimalPlace(d1) != DecimalPlace(d2)) return false;
1955 
1956  int decimalPlace = DecimalPlace(d1);
1957  double factor = pow(10.0, (int)decimalPlace);
1958 
1959  // reduced num is in the form 0.nnnnnnnnnn...
1960  double reducedNum = d1 / factor;
1961 
1962  double round_offset = (d1 < 0) ? -0.5 : 0.5;
1963 
1964  // cast off the digits past the precision's place
1965  long long num1 = ((long long)(reducedNum * cutoff + round_offset));
1966 
1967  factor = pow(10.0, (int)decimalPlace);
1968 
1969  // reduced num is in the form 0.nnnnnnnnnn...
1970  reducedNum = d2 / factor;
1971 
1972  round_offset = (d2 < 0) ? -0.5 : 0.5;
1973 
1974  // cast off the digits past the precision's place
1975  long long num2 = ((long long)(reducedNum * cutoff + round_offset));
1976 
1977 
1978  return (num1 == num2);
1979  }
1980 
1981 
1996  geos::geom::MultiPolygon *PolygonTools::SplitPolygonOn360(const geos::geom::Polygon *inPoly) {
1997  bool convertLon = false;
1998  bool negAdjust = false;
1999  bool newCoords = false; // coordinates have been adjusted
2000  geos::geom::CoordinateSequence *newLonLatPts = new geos::geom::CoordinateArraySequence();
2001  double lon, lat;
2002  double lonOffset = 0;
2003  geos::geom::CoordinateSequence *inPolyCoords = inPoly->getCoordinates();
2004  double prevLon = inPolyCoords->getAt(0).x;
2005  double prevLat = inPolyCoords->getAt(0).y;
2006 
2007  newLonLatPts->add(geos::geom::Coordinate(prevLon, prevLat));
2008  double dist = 0.0;
2009  for (unsigned int i = 1; i < inPolyCoords->getSize(); i++) {
2010  lon = inPolyCoords->getAt(i).x;
2011  lat = inPolyCoords->getAt(i).y;
2012  // check to see if you just crossed the Meridian
2013  if (abs(lon - prevLon) > 180 && (prevLat != 90 && prevLat != -90)) {
2014  newCoords = true;
2015  // if you were already converting then stop (crossed Meridian even number of times)
2016  if (convertLon) {
2017  convertLon = false;
2018  lonOffset = 0;
2019  }
2020  else { // Need to start converting again, deside how to adjust coordinates
2021  if ((lon - prevLon) > 0) {
2022  lonOffset = -360.;
2023  negAdjust = true;
2024  }
2025  else if ((lon - prevLon) < 0) {
2026  lonOffset = 360.;
2027  negAdjust = false;
2028  }
2029  convertLon = true;
2030  }
2031  }
2032 
2033  // Change to a minimum calculation
2034  if (newCoords && dist == 0.0) {
2035  double longitude = (lon + lonOffset) - prevLon;
2036  double latitude = lat - prevLat;
2037  dist = std::sqrt((longitude * longitude) + (latitude * latitude));
2038  }
2039 
2040  // add coord
2041  newLonLatPts->add(geos::geom::Coordinate(lon + lonOffset, lat));
2042 
2043  // set current to old
2044  prevLon = lon;
2045  prevLat = lat;
2046  }
2047 
2048  delete inPolyCoords;
2049 
2050  // Nothing was done so return
2051  if (!newCoords) {
2052  geos::geom::Polygon *newPoly = globalFactory->createPolygon
2053  (globalFactory->createLinearRing(newLonLatPts), NULL);
2054  geos::geom::MultiPolygon *multi_polygon = PolygonTools::MakeMultiPolygon(newPoly);
2055  delete newLonLatPts;
2056  return multi_polygon;
2057  }
2058 
2059  // bisect into seperate polygons
2060  try {
2061  geos::geom::Polygon *newPoly = globalFactory->createPolygon
2062  (globalFactory->createLinearRing(newLonLatPts), NULL);
2063 
2064  geos::geom::CoordinateSequence *pts = new geos::geom::CoordinateArraySequence();
2065  geos::geom::CoordinateSequence *pts2 = new geos::geom::CoordinateArraySequence();
2066 
2067  // Depending on direction of compensation bound accordingly
2068  //***************************************************
2069 
2070  // please verify correct if you change these values
2071  //***************************************************
2072  if (negAdjust) {
2073  pts->add(geos::geom::Coordinate(0., 90.));
2074  pts->add(geos::geom::Coordinate(-360., 90.));
2075  pts->add(geos::geom::Coordinate(-360., -90.));
2076  pts->add(geos::geom::Coordinate(0., -90.));
2077  for (double lat = -90.0 + dist; lat < 90.0; lat += dist) {
2078  pts->add(geos::geom::Coordinate(0.0, lat));
2079  }
2080  pts->add(geos::geom::Coordinate(0., 90.));
2081  pts2->add(geos::geom::Coordinate(0., 90.));
2082  pts2->add(geos::geom::Coordinate(360., 90.));
2083  pts2->add(geos::geom::Coordinate(360., -90.));
2084  pts2->add(geos::geom::Coordinate(0., -90.));
2085  for (double lat = -90.0 + dist; lat < 90.0; lat += dist) {
2086  pts2->add(geos::geom::Coordinate(0.0, lat));
2087  }
2088  pts2->add(geos::geom::Coordinate(0., 90.));
2089  }
2090  else {
2091  pts->add(geos::geom::Coordinate(360., 90.));
2092  pts->add(geos::geom::Coordinate(720., 90.));
2093  pts->add(geos::geom::Coordinate(720., -90.));
2094  pts->add(geos::geom::Coordinate(360., -90.));
2095  for (double lat = -90.0 + dist; lat < 90.0; lat += dist) {
2096  pts->add(geos::geom::Coordinate(360.0, lat));
2097  }
2098  pts->add(geos::geom::Coordinate(360., 90.));
2099  pts2->add(geos::geom::Coordinate(360., 90.));
2100  pts2->add(geos::geom::Coordinate(0., 90.));
2101  pts2->add(geos::geom::Coordinate(0., -90.));
2102  pts2->add(geos::geom::Coordinate(360., -90.));
2103  for (double lat = -90.0 + dist; lat < 90.0; lat += dist) {
2104  pts2->add(geos::geom::Coordinate(360.0, lat));
2105  }
2106  pts2->add(geos::geom::Coordinate(360., 90.));
2107  }
2108 
2109  geos::geom::Polygon *boundaryPoly = globalFactory->createPolygon
2110  (globalFactory->createLinearRing(pts), NULL);
2111  geos::geom::Polygon *boundaryPoly2 = globalFactory->createPolygon
2112  (globalFactory->createLinearRing(pts2), NULL);
2113  /*------------------------------------------------------------------------
2114  / Intersecting the original polygon (converted coordinates) with the
2115  / boundary polygons will create the multi polygons with the converted coordinates.
2116  / These will need to be converted back to the original coordinates.
2117  /-----------------------------------------------------------------------*/
2118  geos::geom::Geometry *intersection = PolygonTools::Intersect(newPoly, boundaryPoly);
2119  geos::geom::MultiPolygon *convertPoly = PolygonTools::MakeMultiPolygon(intersection);
2120  delete intersection;
2121 
2122  intersection = PolygonTools::Intersect(newPoly, boundaryPoly2);
2123  geos::geom::MultiPolygon *convertPoly2 = PolygonTools::MakeMultiPolygon(intersection);
2124  delete intersection;
2125 
2126  /*------------------------------------------------------------------------
2127  / Adjust points created in the negative space or >360 space to be back in
2128  / the 0-360 world. This will always only need to be done on convertPoly.
2129  / Then add geometries to finalpolys.
2130  /-----------------------------------------------------------------------*/
2131  vector<geos::geom::Geometry *> *finalpolys = new vector<geos::geom::Geometry *>;
2132  geos::geom::Geometry *newGeom = NULL;
2133 
2134  for (unsigned int i = 0; i < convertPoly->getNumGeometries(); i++) {
2135  newGeom = (convertPoly->getGeometryN(i))->clone();
2136  pts = convertPoly->getGeometryN(i)->getCoordinates();
2137  geos::geom::CoordinateSequence *newLonLatPts = new geos::geom::CoordinateArraySequence();
2138  // fix the points
2139 
2140  for (unsigned int k = 0; k < pts->getSize() ; k++) {
2141  double lon = pts->getAt(k).x;
2142  double lat = pts->getAt(k).y;
2143  if (negAdjust) {
2144  lon = lon + 360;
2145  }
2146  else {
2147  lon = lon - 360;
2148  }
2149  newLonLatPts->add(geos::geom::Coordinate(lon, lat), k);
2150 
2151  }
2152  // Add the points to polys
2153  finalpolys->push_back(globalFactory->createPolygon
2154  (globalFactory->createLinearRing(newLonLatPts), NULL));
2155  }
2156 
2157  // This loop is over polygons that will always be in 0-360 space no need to convert
2158  for (unsigned int i = 0; i < convertPoly2->getNumGeometries(); i++) {
2159  newGeom = (convertPoly2->getGeometryN(i))->clone();
2160  finalpolys->push_back(newGeom);
2161  }
2162 
2163  geos::geom::MultiPolygon *multi_polygon = globalFactory->createMultiPolygon(*finalpolys);
2164 
2165  delete finalpolys;
2166  delete newGeom;
2167  delete newLonLatPts;
2168  delete pts;
2169  delete pts2;
2170  return multi_polygon;
2171  }
2172  catch(geos::util::IllegalArgumentException *geosIll) {
2173  std::string msg = "Unable to split polygon on longitude boundry (SplitPolygonOn360) due to ";
2174  msg += "geos illegal argument [" + IString(geosIll->what()) + "]";
2175  delete geosIll;
2176  throw IException(IException::Unknown, msg, _FILEINFO_);
2177  }
2178  catch(geos::util::GEOSException *geosExc) {
2179  std::string msg = "Unable to split polygon on longitude boundry (SplitPolygonOn360) due to ";
2180  msg += "geos exception [" + IString(geosExc->what()) + "]";
2181  delete geosExc;
2182  throw IException(IException::Unknown, msg, _FILEINFO_);
2183  }
2184  catch(IException &e) {
2185  std::string msg = "Unable to split polygon on longitude boundry (SplitPolygonOn360) due to ";
2186  msg += "isis operation exception [" + IString(e.what()) + "]";
2187  throw IException(e, IException::Unknown, msg, _FILEINFO_);
2188  }
2189  catch(...) {
2190  std::string msg = "Unable to split polygon on longitude boundry (SplitPolygonOn360) due to ";
2191  msg += "unknown exception";
2192  throw IException(IException::Unknown, msg, _FILEINFO_);
2193  }
2194  }
2195 } // end namespace isis
2196 
Isis::UniversalGroundMap
Universal Ground Map.
Definition: UniversalGroundMap.h:69
Isis::UniversalGroundMap::Sample
double Sample() const
Returns the current line value of the camera model or projection.
Definition: UniversalGroundMap.cpp:200
Isis::TProjection::Longitude
virtual double Longitude() const
This returns a longitude with correct longitude direction and domain as specified in the label object...
Definition: TProjection.cpp:823
Isis::IException::what
const char * what() const
Returns a string representation of this exception in its current state.
Definition: IException.cpp:375
Isis::TProjection::Latitude
virtual double Latitude() const
This returns a latitude with correct latitude type as specified in the label object.
Definition: TProjection.cpp:811
Isis::TProjection
Base class for Map TProjections.
Definition: TProjection.h:166
Isis::IException
Isis exception class.
Definition: IException.h:91
Isis::Projection::SetWorld
virtual bool SetWorld(const double x, const double y)
This method is used to set a world coordinate.
Definition: Projection.cpp:497
Isis::TProjection::SetGround
virtual bool SetGround(const double lat, const double lon)
This method is used to set the latitude/longitude (assumed to be of the correct LatitudeType,...
Definition: TProjection.cpp:760
std
Namespace for the standard library.
Isis::UniversalGroundMap::SetUniversalGround
bool SetUniversalGround(double lat, double lon)
Returns whether the lat/lon position was set successfully in the camera model or projection.
Definition: UniversalGroundMap.cpp:102
Isis::IString
Adds specific functionality to C++ strings.
Definition: IString.h:165
Isis::Projection::YCoord
double YCoord() const
This returns the projection Y provided SetGround, SetCoordinate, SetUniversalGround,...
Definition: Projection.cpp:400
Isis::UniversalGroundMap::Line
double Line() const
Returns the current line value of the camera model or projection.
Definition: UniversalGroundMap.cpp:214
Isis::Projection::XCoord
double XCoord() const
This returns the projection X provided SetGround, SetCoordinate, SetUniversalGround,...
Definition: Projection.cpp:387
Isis
This is free and unencumbered software released into the public domain.
Definition: Apollo.h:16