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"
25#include "SpecialPixel.h"
26#include "PolygonTools.h"
27#include "TProjection.h"
28#include "ProjectionFactory.h"
29#include "UniversalGroundMap.h"
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";
55 if(lonLatPolygon.isEmpty()) {
56 return globalFactory->createMultiPolygon().release();
59 vector<const geos::geom::Geometry *> *xyPolys =
new vector<const geos::geom::Geometry *>;
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));
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();
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(),
81 std::unique_ptr<geos::geom::LinearRing> hole = globalFactory->createLinearRing(*xycoords);
83 if(hole->isValid() && !hole->isEmpty()) {
84 holes.push_back(hole->clone());
92 geos::geom::CoordinateSequence *xycoords =
new geos::geom::CoordinateSequence();
93 geos::geom::CoordinateSequence *llcoords =
94 poly->getExteriorRing()->getCoordinates().release();
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(),
104 geos::geom::Polygon *newPoly = globalFactory->createPolygon(
105 globalFactory->createLinearRing(*xycoords), std::move(holes)).release();
107 if(!newPoly->isEmpty() && newPoly->getArea() > 1.0e-14) {
108 xyPolys->push_back(newPoly);
116 geos::geom::MultiPolygon *spikedPoly = globalFactory->createMultiPolygon(*xyPolys).release();
118 if(spikedPoly->isValid() && !spikedPoly->isEmpty()) {
123 geos::geom::MultiPolygon *despikedPoly =
Despike(spikedPoly);
131 IString msg =
"Unable to convert polygon from Lat/Lon to X/Y";
152 const geos::geom::MultiPolygon &xYPolygon,
TProjection *projection) {
154 if(projection == NULL) {
155 string msg =
"Unable to convert X/Y polygon to Lon/Lat. ";
156 msg +=
"No projection was supplied";
161 if(xYPolygon.isEmpty()) {
162 return globalFactory->createMultiPolygon().release();
165 vector<const geos::geom::Geometry *> *llPolys =
new vector<const geos::geom::Geometry *>;
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));
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();
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(),
186 holes.push_back(globalFactory->createLinearRing(*llcoords));
190 geos::geom::CoordinateSequence *llcoords =
new geos::geom::CoordinateSequence();
191 geos::geom::CoordinateSequence *xycoords =
192 poly->getExteriorRing()->getCoordinates().release();
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(),
202 llPolys->push_back(globalFactory->createPolygon(
203 globalFactory->createLinearRing(*llcoords), std::move(holes)).release());
208 geos::geom::MultiPolygon *spikedPoly = globalFactory->createMultiPolygon(*llPolys).release();
210 if(spikedPoly->isValid() && !spikedPoly->isEmpty()) {
215 geos::geom::MultiPolygon *despikedPoly =
Despike(spikedPoly);
223 IString msg =
"Unable to convert polygon from X/Y to Lat/Lon";
247 string msg =
"Unable to convert Lon/Lat polygon to Sample/Line. ";
248 msg +=
"No UniversalGroundMap was supplied";
253 if(lonLatPolygon.isEmpty()) {
254 return globalFactory->createMultiPolygon().release();
257 vector<const geos::geom::Geometry *> slPolys;
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));
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();
272 for(
unsigned int cord = 0; cord < llcoords->getSize(); ++cord) {
274 llcoords->getAt(cord).x);
275 slcoords->add(geos::geom::Coordinate(ugm->
Sample(),
278 holes.push_back(globalFactory->createLinearRing(*slcoords));
284 geos::geom::CoordinateSequence *slcoords =
new geos::geom::CoordinateSequence();
285 geos::geom::CoordinateSequence *llcoords =
286 poly->getExteriorRing()->getCoordinates().release();
289 for(
unsigned int cord = 0; cord < llcoords->getSize(); ++cord) {
291 llcoords->getAt(cord).x)) {
292 slcoords->add(geos::geom::Coordinate(ugm->
Sample(),
298 if (slcoords->getSize() > 0 && !slcoords->front().equals(slcoords->back())) {
299 slcoords->add(slcoords->front());
303 slPolys.push_back(globalFactory->createPolygon(
304 globalFactory->createLinearRing(*slcoords), std::move(holes)).release());
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()),
317 geos::geom::MultiPolygon *spikedPoly = globalFactory->createMultiPolygon(slPolys).release();
319 if(spikedPoly->isValid() && !spikedPoly->isEmpty()) {
324 geos::geom::MultiPolygon *despikedPoly =
Despike(spikedPoly);
332 IString msg =
"Unable to convert polygon from Lat/Lon to Sample/Line";
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());
358 return globalFactory->createMultiPolygon(*polys).release();
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());
380 return globalFactory->createMultiPolygon(*polys).release();
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>";
408 setprecision(15) << mpolygon->getEnvelopeInternal()->getMinX() <<
411 setprecision(15) << mpolygon->getEnvelopeInternal()->getMinY() <<
413 os <<
"</gml:coord>" << endl;
414 os <<
" <gml:coord>";
416 setprecision(15) << mpolygon->getEnvelopeInternal()->getMaxX() <<
419 setprecision(15) << mpolygon->getEnvelopeInternal()->getMaxY() <<
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>";
431 for(
unsigned int polyN = 0; polyN < mpolygon->getNumGeometries(); polyN++) {
432 geos::geom::CoordinateSequence *pts = mpolygon->getGeometryN(polyN)->getCoordinates().release();
434 for(
unsigned int i = 0; i < pts->getSize(); i++) {
435 double lon = pts->getAt(i).x;
436 double lat = pts->getAt(i).y;
438 os << setprecision(15) << lon <<
"," << setprecision(15) << lat <<
" ";
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>";
449 return os.str().c_str();
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>";
514 return os.str().c_str();
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));
540 unique_ptr<geos::geom::LinearRing> leftOf180Geom =
541 globalFactory->createLinearRing(leftOf180Pts);
543 geos::geom::Polygon *leftOf180Poly =
544 globalFactory->createPolygon(std::move(leftOf180Geom)).release();
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));
553 unique_ptr<geos::geom::LinearRing> rightOf180Geom =
554 globalFactory->createLinearRing(rightOf180Pts);
556 geos::geom::Polygon *rightOf180Poly =
557 globalFactory->createPolygon(std::move(rightOf180Geom)).release();
559 geos::geom::Geometry *preserved =
Intersect(leftOf180Poly, poly360);
560 geos::geom::Geometry *moving =
Intersect(rightOf180Poly, poly360);
562 geos::geom::CoordinateSequence *movingPts = moving->getCoordinates().release();
563 geos::geom::CoordinateSequence movedPts;
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));
570 if(movedPts.getSize()) {
571 movedPts.add(geos::geom::Coordinate(movedPts.getAt(0).x,
572 movedPts.getAt(0).y));
575 geos::geom::Geometry *moved = globalFactory->createPolygon(
576 globalFactory->createLinearRing(movedPts)).release();
578 std::vector<const geos::geom::Geometry *> geomsForCollection;
579 geomsForCollection.push_back(preserved);
580 geomsForCollection.push_back(moved);
582 geos::geom::GeometryCollection *the180Polys =
583 Isis::globalFactory->createGeometryCollection(geomsForCollection).release();
590 geos::geom::MultiPolygon *fixedResult = FixSeam(result);
596 catch(geos::util::GEOSException *exc) {
597 IString msg =
"Conversion to 180 failed. The reason given was [" +
603 IString msg =
"Conversion to 180 failed. Could not determine the reason";
620 const geos::geom::Envelope *envelope = mpolygon->getEnvelopeInternal();
622 double x = fabs(envelope->getMaxX() - envelope->getMinX());
623 double y = fabs(envelope->getMaxY() - envelope->getMinY());
624 double extent = max(x, y);
626 return mpolygon->getArea() / (extent * extent);
643 geos::geom::MultiPolygon *despiked =
Despike(spiked);
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));
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;
680 geos::geom::LinearRing *fixed =
FixGeometry(lr.release());
681 lr.reset(fixed->clone().release());
686 holes.push_back(lr->clone());
694 unique_ptr<geos::geom::LinearRing> ls = poly->getExteriorRing()->clone();
695 unique_ptr<geos::geom::LinearRing> lr;
697 lr =
Despike(ls.release())->clone();
701 geos::geom::LinearRing *fixed =
FixGeometry(lr.release());
702 lr.reset(fixed->clone().release());
708 if (ls && ls->isValid() && ls->getGeometryTypeId() == geos::geom::GEOS_LINEARRING) {
709 lr.reset(ls.release());
718 geos::geom::Polygon *tp = Isis::globalFactory->createPolygon(std::move(lr), std::move(holes)).release();
720 if(tp->isEmpty() || !tp->isValid()) {
722 newPolys->push_back(poly->clone().release());
725 newPolys->push_back(tp);
731 geos::geom::MultiPolygon *mp = Isis::globalFactory->createMultiPolygon(*newPolys).release();
733 if(!mp->isValid() || mp->isEmpty()) {
736 IString msg =
"Despike failed to correct the polygon";
741 if(fabs((mp->getArea() / multiPoly->getArea()) - 1.0) > 0.50) {
742 IString msg =
"Despike failed to correct the polygon " + mp->toString();
766 geos::geom::CoordinateSequence *vertices =
767 new geos::geom::CoordinateSequence(*(lineString->getCoordinates()));
770 if(vertices->getSize() < 4) {
772 return Isis::globalFactory->createLinearRing(geos::geom::CoordinateSequence()).release();
779 std::vector<geos::geom::Coordinate> coords;
780 vertices->toVector(coords);
782 vertices->setPoints(coords);
789 for(
long index = 0; index < (long)vertices->getSize(); index++) {
791 if(vertices->getSize() < 3) {
797 long testCoords[3] = {
805 for(
int j = 0; j < 3; j++) {
806 while(testCoords[j] < 0) {
807 testCoords[j] += vertices->getSize();
810 while(testCoords[j] >= (
long)vertices->getSize()) {
811 testCoords[j] -= vertices->getSize();
816 if(
IsSpiked(vertices->getAt(testCoords[0]),
817 vertices->getAt(testCoords[1]),
818 vertices->getAt(testCoords[2]))) {
820 std::vector<geos::geom::Coordinate> coords;
821 vertices->toVector(coords);
822 coords.erase(coords.begin()+testCoords[1]);
823 vertices->setPoints(coords);
830 if(vertices->getSize() < 3) {
834 return Isis::globalFactory->createLinearRing(geos::geom::CoordinateSequence()).release();
838 vertices->add(vertices->getAt(0));
839 return Isis::globalFactory->createLinearRing(*vertices).release();
853 geos::geom::Coordinate middle,
854 geos::geom::Coordinate last) {
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();
891 geos::geom::CoordinateSequence coords;
894 geos::geom::LineString *line = Isis::globalFactory->createLineString(coords).release();
900 double tolerance = line->getLength() / 100.0;
904 double distanceLastMiddle = geos::operation::distance::DistanceOp::distance(lastPt, middlePt);
905 double distanceLastLine = geos::operation::distance::DistanceOp::distance(lastPt, line);
907 if(distanceLastMiddle == 0.0)
return true;
911 if(distanceLastLine / distanceLastMiddle >= .05) {
916 if(spiked && distanceLastLine > tolerance) {
921 geos::geom::CoordinateSequence coords;
928 geos::geom::LinearRing *shell = Isis::globalFactory->createLinearRing(coords).release();
929 std::vector<geos::geom::LinearRing *> *empty =
new std::vector<geos::geom::LinearRing *>;
932 geos::geom::Polygon *poly = Isis::globalFactory->createPolygon(*shell, *empty);
937 if(poly->getArea() < 1.0e-10) {
967 const geos::geom::Geometry *geom2) {
969 return Operate(geom1, geom2, (
unsigned int)geos::operation::overlayng::OverlayNG::INTERSECTION);
971 catch(geos::util::GEOSException *exc) {
972 IString msg =
"Intersect operation failed. The reason given was [" +
IString(exc->what()) +
"]";
977 IString msg =
"Intersect operation failed";
981 IString msg =
"Intersect operation failed for an unknown reason";
987 geos::geom::Geometry *PolygonTools::Operate(
const geos::geom::Geometry *geom1,
988 const geos::geom::Geometry *geom2,
989 unsigned int opcode) {
994 geos::geom::Geometry *result = NULL;
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()) {
1006 geomFirst = geomSnapped;
1009 unsigned int precision = 15;
1010 unsigned int minPrecision = 13;
1014 std::unique_ptr< geos::geom::Geometry > resultAuto(
1015 geos::operation::overlay::snap::SnapOverlayOp::overlayOp(*geomFirst, *geomSecond, opcode));
1017 result = resultAuto->clone().release();
1019 catch(geos::util::GEOSException *exc) {
1021 if(!failed || precision == minPrecision)
throw;
1026 if(precision == minPrecision) {
1027 IString msg =
"An unknown geos error occurred";
1032 IString msg =
"An unknown geos error occurred when attempting to clone a geometry";
1056 if(result && !result->isValid()) {
1058 geos::geom::Geometry *newResult =
FixGeometry(result);
1060 if(fabs(newResult->getArea() / result->getArea() - 1.0) > 0.50) {
1064 IString msg =
"Operation [" + IString((
int)opcode) +
"] failed";
1071 catch(IException &e) {
1072 IString msg =
"Operation [" + IString((
int)opcode) +
"] failed";
1077 if(result == NULL) {
1078 IString msg =
"Operation [" + IString((
int)opcode) +
" failed";
1096 if(geom->getGeometryTypeId() == geos::geom::GEOS_MULTIPOLYGON) {
1097 return FixGeometry(
dynamic_cast<const geos::geom::MultiPolygon *
>(geom));
1099 if(geom->getGeometryTypeId() == geos::geom::GEOS_LINEARRING) {
1100 return FixGeometry(
dynamic_cast<const geos::geom::LinearRing *
>(geom));
1102 if(geom->getGeometryTypeId() == geos::geom::GEOS_POLYGON) {
1103 return FixGeometry(
dynamic_cast<const geos::geom::Polygon *
>(geom));
1105 if(geom->getGeometryTypeId() == geos::geom::GEOS_GEOMETRYCOLLECTION) {
1125 vector<const geos::geom::Geometry *> *newPolys =
new vector<const geos::geom::Geometry *>;
1128 for(
unsigned int geomIndex = 0; geomIndex < poly->getNumGeometries(); geomIndex ++) {
1130 dynamic_cast<const geos::geom::Polygon *
>(
1131 poly->getGeometryN(geomIndex)));
1132 if(fixedpoly->isValid()) {
1133 newPolys->push_back(fixedpoly);
1141 geos::geom::MultiPolygon *mp = Isis::globalFactory->createMultiPolygon(*newPolys).release();
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);
1154 if(thisHole->getGeometryTypeId() != geos::geom::GEOS_LINEARRING) {
1155 holes.push_back(thisHole->clone());
1161 geos::geom::LinearRing *newHole =
FixGeometry((geos::geom::LinearRing *)thisHole);
1162 holes.push_back(newHole->clone());
1164 catch (IException &e) {
1165 IString msg =
"Failed when attempting to fix interior ring of multipolygon";
1171 const geos::geom::LineString *exterior = poly->getExteriorRing();
1174 unique_ptr<geos::geom::LinearRing> newExterior;
1176 if(exterior->getGeometryTypeId() == geos::geom::GEOS_LINEARRING) {
1177 newExterior.reset(
FixGeometry((geos::geom::LinearRing *)exterior));
1180 IString msg =
"Failed when attempting to fix exterior ring of polygon. The exterior "
1181 "ring is not simple and closed";
1185 return globalFactory->createPolygon(std::move(newExterior), std::move(holes)).release();
1187 catch (IException &e) {
1188 IString msg =
"Failed when attempting to fix exterior ring of polygon";
1211 geos::geom::CoordinateSequence *coords = ring->getCoordinates().release();
1214 if(coords->getSize() < 4) {
1215 return globalFactory->createLinearRing(geos::geom::CoordinateSequence()).release();
1218 geos::geom::CoordinateSequence *newCoords =
new geos::geom::CoordinateSequence();
1219 const geos::geom::Coordinate *lastCoordinate = &coords->getAt(0);
1220 newCoords->add(*lastCoordinate);
1223 for(
unsigned int coordIndex = 1; coordIndex < coords->getSize() - 1; coordIndex ++) {
1224 const geos::geom::Coordinate *thisCoordinate = &coords->getAt(coordIndex);
1229 double difference[2] = {
1230 lastCoordinate->x - thisCoordinate->x,
1231 lastCoordinate->y - thisCoordinate->y,
1241 if(difference[0] == 0.0 && difference[1] != 0.0) {
1245 else if(difference[1] == 0.0 && difference[0] != 0.0) {
1249 else if(difference[0] == 0.0 && difference[1] == 0.0) {
1256 newCoords->add(*thisCoordinate);
1257 lastCoordinate = thisCoordinate;
1261 newCoords->add(geos::geom::Coordinate(newCoords->getAt(0).x, newCoords->getAt(0).y));
1262 geos::geom::LinearRing *newRing = NULL;
1266 if(newCoords->getSize() > 3) {
1267 newRing = globalFactory->createLinearRing(*newCoords).release();
1274 catch(geos::util::GEOSException *exc) {
1278 IString msg =
"Error when attempting to fix linear ring";
1282 if(newRing && !newRing->isValid() && ring->isValid()) {
1283 IString msg =
"Failed when attempting to fix linear ring";
1286 else if(!newRing || !newRing->isValid()) {
1291 newRing =
dynamic_cast<geos::geom::LinearRing *
>(ring->clone().release());
1311 if(num == 0.0)
return 0;
1315 int decimalPlace = 1;
1326 return decimalPlace;
1338 const geos::geom::Geometry *geom2) {
1340 return Operate(geom1, geom2, (
unsigned int)geos::operation::overlayng::OverlayNG::DIFFERENCE);
1342 catch(geos::util::GEOSException *exc) {
1343 IString msg =
"Difference operation failed. The reason given was [" +
1349 IString msg =
"Difference operation failed";
1353 IString msg =
"Difference operation failed for an unknown reason";
1374 if(geom->isEmpty()) {
1375 return Isis::globalFactory->createMultiPolygon().release();
1378 else if(geom->getArea() - DBL_EPSILON <= DBL_EPSILON) {
1379 return Isis::globalFactory->createMultiPolygon().release();
1382 else if(geom->getGeometryTypeId() == geos::geom::GEOS_MULTIPOLYGON) {
1383 return dynamic_cast<geos::geom::MultiPolygon *
>(geom->clone().release());
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();
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 =
1401 for(
unsigned int subPoly = 0;
1402 subPoly < subMultiPoly->getNumGeometries();
1404 const geos::geom::Polygon *poly =
1405 dynamic_cast<const geos::geom::Polygon *
>(
1406 subMultiPoly->getGeometryN(subPoly));
1407 polys.push_back(poly->clone().release());
1411 geos::geom::MultiPolygon *mp =
1412 Isis::globalFactory->createMultiPolygon(polys).release();
1413 if(mp->getArea() - DBL_EPSILON <= DBL_EPSILON) {
1415 mp = Isis::globalFactory->createMultiPolygon().release();
1423 return Isis::globalFactory->createMultiPolygon().release();
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();
1433 unsigned int aIntersectionBegin = 0;
1434 unsigned int aIntersectionEnd = 0;
1435 unsigned int bIntersectionBegin = 0;
1436 unsigned int bIntersectionEnd = 0;
1438 bool intersectionStarted =
false;
1439 bool intersectionEnded =
false;
1441 unsigned int lastBMatch = 0;
1442 for (
unsigned int i = 0;
1443 !intersectionEnded && i < polyAPoints->getSize();
1446 bool foundEquivalent =
false;
1448 geos::geom::Coordinate coordA = polyAPoints->getAt(i);
1451 for (
unsigned int j = 0;
1452 !foundEquivalent && j < polyBPoints->getSize();
1454 geos::geom::Coordinate coordB = polyBPoints->getAt(j);
1457 foundEquivalent = coordA.equals(coordB);
1459 if (foundEquivalent) lastBMatch = j;
1461 if (foundEquivalent && !intersectionStarted) {
1462 intersectionStarted =
true;
1463 aIntersectionBegin = i;
1464 bIntersectionBegin = j;
1468 if (!foundEquivalent && intersectionStarted && !intersectionEnded) {
1469 intersectionEnded =
true;
1470 aIntersectionEnd = i;
1471 bIntersectionEnd = lastBMatch;
1475 geos::geom::MultiPolygon * result = NULL;
1476 if (intersectionStarted && intersectionEnded) {
1477 geos::geom::CoordinateSequence merged;
1480 for (i = 0; i < aIntersectionBegin; i ++) {
1481 merged.add(polyAPoints->getAt(i));
1484 i = bIntersectionBegin;
1485 while (i != bIntersectionEnd) {
1486 merged.add(polyBPoints->getAt(i));
1488 if (i >= polyBPoints->getSize()) i = 0;
1491 for (i = aIntersectionEnd; i < polyAPoints->getSize() - 1; i++) {
1492 merged.add(polyAPoints->getAt(i));
1495 merged.add(merged.getAt(0));
1497 globalFactory->createLinearRing(merged)).release());
1504 geos::geom::MultiPolygon *PolygonTools::FixSeam(
1505 const geos::geom::MultiPolygon *poly) {
1507 std::vector<const geos::geom::Geometry *> polys;
1510 for(
unsigned int copyIndex = 0;
1511 copyIndex < poly->getNumGeometries();
1513 polys.push_back(poly->getGeometryN(copyIndex)->clone().release());
1516 unsigned int outerPolyIndex = 0;
1518 while(outerPolyIndex + 1 < polys.size()) {
1519 unsigned int innerPolyIndex = outerPolyIndex + 1;
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)));
1526 if(fixedPair != NULL) {
1527 const geos::geom::Geometry *oldInnerPoly = polys.at(innerPolyIndex);
1528 const geos::geom::Geometry *oldOuterPoly = polys.at(outerPolyIndex);
1530 polys.erase(polys.begin() + innerPolyIndex);
1531 polys[outerPolyIndex] = fixedPair->getGeometryN(0)->clone().release();
1532 innerPolyIndex = outerPolyIndex + 1;
1537 delete oldInnerPoly;
1538 oldInnerPoly = NULL;
1540 delete oldOuterPoly;
1541 oldOuterPoly = NULL;
1551 return globalFactory->createMultiPolygon(polys).release();
1565 unsigned int precision) {
1566 if(geom->getGeometryTypeId() == geos::geom::GEOS_MULTIPOLYGON) {
1568 dynamic_cast<const geos::geom::MultiPolygon *
>(geom), precision);
1570 if(geom->getGeometryTypeId() == geos::geom::GEOS_LINEARRING) {
1572 dynamic_cast<const geos::geom::LinearRing *
>(geom), precision);
1574 if(geom->getGeometryTypeId() == geos::geom::GEOS_POLYGON) {
1576 dynamic_cast<const geos::geom::Polygon *
>(geom), precision);
1578 if(geom->getGeometryTypeId() == geos::geom::GEOS_GEOMETRYCOLLECTION) {
1582 IString msg =
"PolygonTools::ReducePrecision does not support [" +
1599 unsigned int precision) {
1601 vector<const geos::geom::Geometry *> newPolys;
1604 for(
unsigned int geomIndex = 0; geomIndex < poly->getNumGeometries(); geomIndex ++) {
1606 dynamic_cast<const geos::geom::Polygon *
>(
1607 poly->getGeometryN(geomIndex)),
1610 if(!lowerPrecision->isEmpty()) {
1611 newPolys.push_back(lowerPrecision);
1614 delete lowerPrecision;
1618 geos::geom::MultiPolygon *mp = Isis::globalFactory->createMultiPolygon(newPolys).release();
1633 unsigned int precision) {
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);
1640 if(thisHole->getGeometryTypeId() != geos::geom::GEOS_LINEARRING) {
1641 holes.push_back(thisHole->clone());
1646 geos::geom::LinearRing *newHole =
ReducePrecision((geos::geom::LinearRing *)thisHole,
1649 if(!newHole->isEmpty()) {
1650 holes.push_back(newHole->clone());
1658 IString msg =
"Failed when attempting to fix interior ring of multipolygon";
1664 const geos::geom::LineString *exterior = poly->getExteriorRing();
1667 unique_ptr<geos::geom::LinearRing> newExterior;
1669 if(exterior->getGeometryTypeId() == geos::geom::GEOS_LINEARRING) {
1670 newExterior.reset(
ReducePrecision((geos::geom::LinearRing *)exterior, precision));
1673 IString msg =
"Failed when attempting to fix exterior ring of polygon. The exterior "
1674 "ring is not simple and closed";
1678 return globalFactory->createPolygon(std::move(newExterior), std::move(holes)).release();
1681 IString msg =
"Failed when attempting to fix exterior ring of polygon";
1697 unsigned int precision) {
1698 geos::geom::CoordinateSequence *coords = ring->getCoordinates().release();
1701 if(coords->getSize() <= 0) {
1702 return dynamic_cast<geos::geom::LinearRing *
>(ring->clone().release());
1705 geos::geom::CoordinateSequence newCoords;
1706 geos::geom::Coordinate *coord =
ReducePrecision(&coords->getAt(0), precision);
1707 newCoords.add(*coord);
1712 for(
unsigned int coordIndex = 1; coordIndex < coords->getSize() - 1; coordIndex ++) {
1713 const geos::geom::Coordinate *thisCoordinate = &coords->getAt(coordIndex);
1715 newCoords.add(*coord);
1720 newCoords.add(geos::geom::Coordinate(newCoords.getAt(0).x, newCoords.getAt(0).y));
1721 geos::geom::LinearRing *newRing = NULL;
1725 newRing = globalFactory->createLinearRing(newCoords).release();
1727 catch(geos::util::GEOSException *exc) {
1731 IString msg =
"Error when attempting to reduce precision of linear ring";
1737 geos::geom::LinearRing *tmp =
Despike(newRing);
1744 if(!newRing->isValid()) {
1745 IString msg =
"Failed when attempting to reduce precision of linear ring";
1763 unsigned int precision) {
1764 return new geos::geom::Coordinate(
1782 double result = num;
1787 double factor = pow(10.0, (
int)decimalPlace);
1790 double reducedNum = num / factor;
1792 double cutoff = pow(10.0, (
int)precision);
1793 double round_offset = (num < 0) ? -0.5 : 0.5;
1796 reducedNum = ((
long long)(reducedNum * cutoff + round_offset)) / cutoff;
1797 result = reducedNum * factor;
1813 switch(geom->getGeometryTypeId()) {
1814 case geos::geom::GEOS_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:
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";
1836 bool PolygonTools::Equal(
const geos::geom::MultiPolygon *poly1,
1837 const geos::geom::MultiPolygon *poly2) {
1839 vector<const geos::geom::Polygon *> polys1;
1840 vector<const geos::geom::Polygon *> polys2;
1842 if(poly1->getNumGeometries() != poly2->getNumGeometries())
return false;
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)));
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])) {
1856 polys1[p1] = polys1[polys1.size()-1];
1857 polys1.resize(polys1.size() - 1);
1859 polys2[p2] = polys2[polys2.size()-1];
1860 polys2.resize(polys2.size() - 1);
1865 return (polys1.size() == 0) && (polys2.size() == 0);
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;
1873 if(poly1->getNumInteriorRing() != poly2->getNumInteriorRing())
return false;
1875 if(!Equal(poly1->getExteriorRing(), poly2->getExteriorRing()))
return false;
1878 for(
unsigned int holeIndex = 0; holeIndex < poly1->getNumInteriorRing(); holeIndex ++) {
1881 if(poly1->getInteriorRingN(holeIndex)->getGeometryTypeId() == geos::geom::GEOS_LINESTRING) {
1882 holes1.push_back(poly1->getInteriorRingN(holeIndex));
1885 if(poly2->getInteriorRingN(holeIndex)->getGeometryTypeId() == geos::geom::GEOS_LINESTRING) {
1886 holes2.push_back(poly2->getInteriorRingN(holeIndex));
1891 if(holes1.size() != holes2.size())
return false;
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])) {
1897 holes1[h1] = holes1[holes1.size()-1];
1898 holes1.resize(holes1.size() - 1);
1900 holes2[h2] = holes2[holes2.size()-1];
1901 holes2.resize(holes2.size() - 1);
1906 return (holes1.size() == 0) && (holes2.size() == 0);
1910 bool PolygonTools::Equal(
const geos::geom::LineString *lineString1,
1911 const geos::geom::LineString *lineString2) {
1913 geos::geom::CoordinateSequence *coords1 = lineString1->getCoordinates().release();
1914 geos::geom::CoordinateSequence *coords2 = lineString2->getCoordinates().release();
1915 bool isEqual =
true;
1917 if(coords1->getSize() != coords2->getSize()) isEqual =
false;
1919 unsigned int index1 = 0;
1920 unsigned int index2 = 0;
1923 for(; index2 < coords2->getSize() - 1 && isEqual; index2 ++) {
1924 if(Equal(coords1->getAt(index1), coords2->getAt(index2)))
break;
1927 if(index2 == coords2->getSize() - 1) isEqual =
false;
1929 for(; index1 < coords1->getSize() - 1 && isEqual; index1 ++, index2 ++) {
1930 if(!Equal(coords1->getAt(index1), coords2->getAt(index2 % (coords2->getSize() - 1)))) {
1941 bool PolygonTools::Equal(
const geos::geom::Coordinate &coord1,
1942 const geos::geom::Coordinate &coord2) {
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;
1952 bool PolygonTools::Equal(
const double d1,
const double d2) {
1953 const double cutoff = 1e15;
1958 double factor = pow(10.0, (
int)decimalPlace);
1961 double reducedNum = d1 / factor;
1963 double round_offset = (d1 < 0) ? -0.5 : 0.5;
1966 long long num1 = ((
long long)(reducedNum * cutoff + round_offset));
1968 factor = pow(10.0, (
int)decimalPlace);
1971 reducedNum = d2 / factor;
1973 round_offset = (d2 < 0) ? -0.5 : 0.5;
1976 long long num2 = ((
long long)(reducedNum * cutoff + round_offset));
1979 return (num1 == num2);
1998 bool convertLon =
false;
1999 bool negAdjust =
false;
2000 bool newCoords =
false;
2001 geos::geom::CoordinateSequence *newLonLatPts =
new geos::geom::CoordinateSequence();
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;
2008 newLonLatPts->add(geos::geom::Coordinate(prevLon, prevLat));
2010 for (
unsigned int i = 1; i < inPolyCoords->getSize(); i++) {
2011 lon = inPolyCoords->getAt(i).x;
2012 lat = inPolyCoords->getAt(i).y;
2014 if (abs(lon - prevLon) > 180 && (prevLat != 90 && prevLat != -90)) {
2022 if ((lon - prevLon) > 0) {
2026 else if ((lon - prevLon) < 0) {
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));
2042 newLonLatPts->add(geos::geom::Coordinate(lon + lonOffset, lat));
2049 delete inPolyCoords;
2053 geos::geom::Polygon *newPoly = globalFactory->createPolygon
2054 (globalFactory->createLinearRing(*newLonLatPts)).release();
2056 delete newLonLatPts;
2057 return multi_polygon;
2062 geos::geom::Polygon *newPoly = globalFactory->createPolygon
2063 (globalFactory->createLinearRing(*newLonLatPts)).release();
2065 geos::geom::CoordinateSequence *pts =
new geos::geom::CoordinateSequence();
2066 geos::geom::CoordinateSequence *pts2 =
new geos::geom::CoordinateSequence();
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));
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));
2089 pts2->add(geos::geom::Coordinate(0., 90.));
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));
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));
2107 pts2->add(geos::geom::Coordinate(360., 90.));
2110 geos::geom::Polygon *boundaryPoly = globalFactory->createPolygon
2111 (globalFactory->createLinearRing(*pts)).release();
2112 geos::geom::Polygon *boundaryPoly2 = globalFactory->createPolygon
2113 (globalFactory->createLinearRing(*pts2)).release();
2121 delete intersection;
2125 delete intersection;
2132 vector<const geos::geom::Geometry *> finalpolys;
2133 geos::geom::Geometry *newGeom = NULL;
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();
2141 for (
size_t k = 0; k < pts3->getSize() ; k++) {
2142 double lon = pts3->getAt(k).x;
2143 double lat = pts3->getAt(k).y;
2150 newLonLatPts->add(geos::geom::Coordinate(lon, lat), k);
2153 finalpolys.push_back(globalFactory->createPolygon
2154 (globalFactory->createLinearRing(*newLonLatPts)).release());
2158 for (
unsigned int i = 0; i < convertPoly2->getNumGeometries(); i++) {
2159 newGeom = (convertPoly2->getGeometryN(i))->clone().release();
2160 finalpolys.push_back(newGeom);
2163 geos::geom::MultiPolygon *multi_polygon = globalFactory->createMultiPolygon(finalpolys).release();
2165 delete newLonLatPts;
2168 return multi_polygon;
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()) +
"]";
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()) +
"]";
2183 std::string msg =
"Unable to split polygon on longitude boundry (SplitPolygonOn360) due to ";
2184 msg +=
"isis operation exception [" +
IString(e.
what()) +
"]";
2188 std::string msg =
"Unable to split polygon on longitude boundry (SplitPolygonOn360) due to ";
2189 msg +=
"unknown exception";
@ Unknown
A type of error that cannot be classified as any of the other error types.
@ Programmer
This error is for when a programmer made an API call that was illegal.
const char * what() const
Returns a string representation of this exception in its current state.
Adds specific functionality to C++ strings.
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...
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.
Namespace for the standard library.