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