Isis 3 Programmer Reference
GisGeometry.cpp
1 
6 /* SPDX-License-Identifier: CC0-1.0 */
7 #include "GisGeometry.h"
8 
9 // Qt library
10 #include <QDebug>
11 #include <QScopedPointer>
12 #include <QString>
13 
14 // geos library
15 #include <geos_c.h>
16 
17 // other ISIS
18 #include "Cube.h"
19 #include "IException.h"
20 #include "ImagePolygon.h"
21 #include "GisTopology.h"
22 #include "SpecialPixel.h"
23 
24 //#define DISABLE_PREPARED_GEOMETRY 1
25 
26 using namespace std;
27 
28 namespace Isis {
29 
33  GisGeometry::GisGeometry() : m_type(None), m_geom(0), m_preparedGeom(0) {
34  // Must ensure GEOS is initialized
35  GisTopology::instance();
36  }
37 
38 
50  GisGeometry::GisGeometry(const double xlongitude, const double ylatitude) :
51  m_type(GeosGis), m_geom(0), m_preparedGeom(0){
52 
53  GisTopology::instance();
54  m_geom = makePoint(xlongitude, ylatitude);
55  m_preparedGeom = makePrepared(this->m_geom);
56  return;
57  }
58 
59 
68  GisGeometry::GisGeometry(Cube &cube) : m_type(IsisCube), m_geom(0), m_preparedGeom(0) {
69 
70  GisTopology::instance();
71  m_geom = fromCube(cube);
72  m_preparedGeom = makePrepared(this->m_geom);
73  return;
74  }
75 
76 
85  GisGeometry::GisGeometry(const QString &gisSource, const GisGeometry::Type t) :
86  m_type(t), m_geom(0), m_preparedGeom(0) {
87  GisTopology *gis(GisTopology::instance());
88  if (WKT == t) {
89  m_geom = gis->geomFromWKT(gisSource);
90  }
91  else if (WKB == t) {
92  m_geom = gis->geomFromWKB(gisSource);
93  }
94  else if (IsisCube == t) {
95  Cube cube;
96  cube.open(gisSource);
97  m_geom = fromCube(cube);
98  }
99  else {
100  throw IException(IException::Programmer,
101  QString("Unknown GIS type given [%1]").arg(typeToString(t)),
102  _FILEINFO_);
103  }
104 
105  // Get the prepared geometry
106  m_preparedGeom = makePrepared(this->m_geom);
107  m_type = t;
108  return;
109  }
110 
111 
120  GisGeometry::GisGeometry(const GisGeometry &geom) : m_type(geom.m_type),
121  m_geom(0),
122  m_preparedGeom(0) {
123 
124  GisTopology *gis(GisTopology::instance());
125  if ( geom.isDefined() ) {
126  m_geom = gis->clone(geom.m_geom);
127  m_preparedGeom = makePrepared(this->m_geom);
128  }
129  }
130 
131 
140  GisGeometry::GisGeometry(GEOSGeometry *geom) : m_type(GeosGis), m_geom(geom),
141  m_preparedGeom(0) {
142  GisTopology::instance();
143  m_preparedGeom = makePrepared(this->m_geom);
144  return;
145  }
146 
147 
151  GisGeometry::~GisGeometry() {
152  destroy();
153  }
154 
155 
164  GisGeometry &GisGeometry::operator=(GisGeometry const &geom) {
165  if ( this != &geom ) {
166  destroy();
167  m_type = geom.m_type;
168  if ( geom.isDefined() ) {
169  GisTopology *gis(GisTopology::instance());
170  m_geom = gis->clone(geom.m_geom);
171  m_preparedGeom = makePrepared(this->m_geom);
172  }
173  }
174  return (*this);
175  }
176 
177 
187  void GisGeometry::setGeometry(GEOSGeometry *geom) {
188  destroy();
189  m_geom = geom;
190  m_preparedGeom = makePrepared(this->m_geom);
191  return;
192  }
193 
194 
200  bool GisGeometry::isDefined() const {
201  return (m_geom != 0);
202  }
203 
204 
213  bool GisGeometry::isValid() const {
214  if (!isDefined()) {
215  return (false);
216  }
217 
218  return (1 == GEOSisValid(this->m_geom));
219  }
220 
221 
227  QString GisGeometry::isValidReason() const {
228  QString result = "Not defined!";
229  if ( isDefined() ) {
230  GisTopology *gis(GisTopology::instance());
231  char *reason = GEOSisValidReason(this->m_geom);
232  result = reason;
233  gis->destroy(reason);
234  }
235  return (result);
236  }
237 
238 
244  GisGeometry::Type GisGeometry::type() const {
245  return (m_type);
246  }
247 
248 
256  GisGeometry::Type GisGeometry::type(const QString &gstrType) {
257  QString gtype = gstrType.toLower();
258  if ("wkt" == gtype) return (WKT);
259  if ("wkb" == gtype) return (WKB);
260  if ("cube" == gtype) return (IsisCube);
261  if ("isiscube" == gtype) return (IsisCube);
262  if ("geometry" == gtype) return (GeosGis);
263  if ("geosgis" == gtype) return (GeosGis);
264  if ("gis" == gtype) return (GeosGis);
265  if ("geos" == gtype) return (GeosGis);
266  return (None);
267  }
268 
269 
277  QString GisGeometry::typeToString(const GisGeometry::Type &type) {
278  if (WKT == type) return "WKT";
279  if (WKB == type) return "WKB";
280  if (IsisCube == type) return "IsisCube";
281  if (GeosGis == type) return "GeosGis";
282  return "None";
283  }
284 
285 
291  const GEOSGeometry *GisGeometry::geometry() const {
292  return (m_geom);
293  }
294 
295 
305  const GEOSPreparedGeometry *GisGeometry::preparedGeometry() const {
306  return (m_preparedGeom);
307  }
308 
309 
318  GisGeometry *GisGeometry::clone() const {
319  if (!isDefined()) {
320  return (new GisGeometry());
321  }
322 
323  GisTopology *gis(GisTopology::instance());
324  QScopedPointer<GisGeometry> geom(new GisGeometry());
325 
326  geom->m_type = m_type;
327  geom->m_geom = gis->clone(m_geom);
328  geom->m_preparedGeom = makePrepared(geom->m_geom);
329  return (geom.take());
330  }
331 
332 
338  bool GisGeometry::isEmpty() const {
339  if ( !isValid() ) {
340  return (true);
341  }
342  return (1 == GEOSisEmpty(this->m_geom));
343  }
344 
345 
359  double GisGeometry::area( ) const {
360  if ( !isValid() ) {
361  return (0.0);
362  }
363 
364  int result = 0;
365  double gisArea = 0.0;
366  result = GEOSArea(this->m_geom, &gisArea);
367  if (1 != result) {
368  gisArea = 0.0;
369  }
370  return (gisArea);
371  }
372 
373 
382  double GisGeometry::length( ) const {
383  if ( !isValid() ) {
384  return (0.0);
385  }
386 
387  int result = 0;
388  double gisLength = 0.0;
389  result = GEOSLength(this->m_geom, &gisLength);
390  if (1 != result) {
391  gisLength = 0.0;
392  }
393  return (gisLength);
394  }
395 
396 
406  double GisGeometry::distance(const GisGeometry &target) const {
407  if ( !isValid() ) {
408  return (false);
409  }
410  if ( !target.isValid() ) {
411  return (false);
412  }
413 
414  int result = 0;
415  double dist = Null;
416  result = GEOSDistance(this->m_geom, target.geometry(), &dist);
417  if ( 1 != result ) {
418  dist = Null;
419  }
420  return (dist);
421  }
422 
423 
429  int GisGeometry::points() const {
430  if (!isValid() ) { return (0); }
431 
432  int ngeoms = GEOSGetNumGeometries(m_geom);
433  int npoints = 0;
434  for (int i = 0 ; i < ngeoms ; i++) {
435  npoints += GEOSGetNumCoordinates(GEOSGetGeometryN(m_geom, i));
436  }
437 
438  return ( npoints );
439  }
440 
441 
451  bool GisGeometry::intersects(const GisGeometry &target) const {
452  if ( !isValid() ) {
453  return (false);
454  }
455  if ( !target.isValid() ) {
456  return (false);
457  }
458 
459  int result = 0;
460  if ( 0 != this->m_preparedGeom) {
461  result = GEOSPreparedIntersects(this->m_preparedGeom, target.geometry());
462  }
463  else {
464  result = GEOSIntersects(this->m_geom, target.geometry());
465  }
466 
467  return (1 == result);
468  }
469 
470 
478  bool GisGeometry::contains(const GisGeometry &target) const {
479  if ( !isValid() ) {
480  return (false);
481  }
482  if ( !target.isValid() ) {
483  return (false);
484  }
485 
486  int result = 0;
487  if ( 0 != this->m_preparedGeom) {
488  result = GEOSPreparedContains(this->m_preparedGeom, target.geometry());
489  }
490  else {
491  result = GEOSContains(this->m_geom, target.geometry());
492  }
493 
494  return (1 == result);
495  }
496 
497 
504  bool GisGeometry::disjoint(const GisGeometry &target) const {
505  if ( !isValid() ) {
506  return (false);
507  }
508  if ( !target.isValid() ) {
509  return (false);
510  }
511 
512  int result = 0;
513  if ( 0 != m_preparedGeom) {
514  result = GEOSPreparedDisjoint(m_preparedGeom, target.geometry());
515  }
516  else {
517  result = GEOSDisjoint(m_geom, target.geometry());
518  }
519 
520  return (1 == result);
521  }
522 
523 
531  bool GisGeometry::overlaps(const GisGeometry &target) const {
532  if ( !isValid() ) {
533  return (false);
534  }
535  if ( !target.isValid() ) {
536  return (false);
537  }
538 
539  int result = 0;
540  if ( 0 != m_preparedGeom) {
541  result = GEOSPreparedOverlaps(m_preparedGeom, target.geometry());
542  }
543  else {
544  result = GEOSOverlaps(m_geom, target.geometry());
545  }
546 
547  return (1 == result);
548  }
549 
550 
557  bool GisGeometry::equals(const GisGeometry &target) const {
558  if ( !isValid() ) {
559  return (false);
560  }
561 
562  if ( !target.isValid() ) {
563  return (false);
564  }
565 
566  int result = GEOSEquals(this->m_geom, target.geometry());
567  return ( 1 == result );
568  }
569 
570 
582  double GisGeometry::intersectRatio(const GisGeometry &target) const {
583  if ( !isValid() ) {
584  return (0.0);
585  }
586  if ( !target.isValid() ) {
587  return (0.0);
588  }
589 
590  // Check for any intersection at all
591  // if ( !intersects(target) ) {
592  // return (0.0);
593  // }
594 
595  // Prevent dividing by 0
596  if (this->area() == 0) {
597  return (0.0);
598  }
599 
600  QScopedPointer<GisGeometry> inCommon(intersection(target));
601  double ratio = inCommon->area() / this->area();
602  return (ratio);
603  }
604 
605 
619  GisGeometry *GisGeometry::envelope() const {
620  if ( !isValid() ) {
621  return (new GisGeometry());
622  }
623 
624  GEOSGeometry *geom = GEOSEnvelope(m_geom);
625  return (new GisGeometry(geom));
626  }
627 
628 
636  GisGeometry *GisGeometry::convexHull() const {
637  if ( !isValid() ) {
638  return (new GisGeometry());
639  }
640 
641  GEOSGeometry *geom = GEOSConvexHull(m_geom);
642  return (new GisGeometry(geom));
643  }
644 
645 
659  GisGeometry *GisGeometry::simplify(const double &tolerance) const {
660  if ( !isValid() ) {
661  return (0);
662  }
663  GEOSGeometry *geom = GEOSTopologyPreserveSimplify(m_geom, tolerance);
664  if ( 0 == geom ) {
665  return (0);
666  }
667 
668  return (new GisGeometry(geom));
669  }
670 
671 
681  GisGeometry *GisGeometry::intersection(const GisGeometry &target) const {
682  // Non-valid geometries return empty geometries
683  if ( !isValid() ) {
684  return (new GisGeometry());
685  }
686  if ( !target.isValid() ) {
687  return (new GisGeometry());
688  }
689 
690  GEOSGeometry *geom = GEOSIntersection(m_geom, target.geometry());
691  return (new GisGeometry(geom));
692  }
693 
694 
701  GisGeometry *GisGeometry::g_union(const GisGeometry &target) const {
702  if ( !isValid() ) {
703  return (new GisGeometry());
704  }
705  if ( !target.isValid() ) {
706  return (new GisGeometry());
707  }
708 
709  GEOSGeometry *geom = GEOSUnion(m_geom, target.geometry());
710  return (new GisGeometry(geom));
711  }
712 
713 
723  bool GisGeometry::centroid(double &xlongitude, double &ylatitude) const {
724  xlongitude = ylatitude = Null;
725  if ( !isValid() ) {
726  return (false);
727  }
728 
729  GEOSGeometry *center = GEOSGetCentroid(m_geom);
730  if ( 0 != center ) {
731  GEOSGeomGetX(center, &xlongitude);
732  GEOSGeomGetY(center, &ylatitude);
733  GisTopology::instance()->destroy(center);
734  return (true);
735  }
736 
737  return (false);
738  }
739 
740 
746  GisGeometry *GisGeometry::centroid() const {
747  if ( !isValid() ) {
748  return (new GisGeometry());
749  }
750 
751  GEOSGeometry *center = GEOSGetCentroid(m_geom);
752  return (new GisGeometry(center));
753  }
754 
755 
762  GEOSPreparedGeometry const *GisGeometry::makePrepared(const GEOSGeometry *geom)
763  const {
764  #if defined(DISABLE_PREPARED_GEOMETRY)
765  return (0);
766  #else
767  GisTopology *gis(GisTopology::instance());
768  return (gis->preparedGeometry(geom));
769  #endif
770  }
771 
772 
779  GEOSGeometry *GisGeometry::makePoint(const double x, const double y) const {
780 
781  GEOSCoordSequence *point = GEOSCoordSeq_create(1, 2);
782  GEOSCoordSeq_setX(point, 0, x);
783  GEOSCoordSeq_setY(point, 0, y);
784 
785  return (GEOSGeom_createPoint(point));
786  }
787 
788 
795  GEOSGeometry *GisGeometry::fromCube(Cube &cube) const {
796  ImagePolygon myGis = cube.readFootprint();
797  GisTopology *gis(GisTopology::instance());
798  QString polyStr = QString::fromStdString(myGis.polyStr());
799  return (gis->geomFromWKT(polyStr));
800  }
801 
802 
806  void GisGeometry::destroy() {
807  GisTopology *gis(GisTopology::instance());
808  gis->destroy(m_geom);
809  gis->destroy(m_preparedGeom);
810  m_geom = 0;
811  m_preparedGeom = 0;
812  return;
813  }
814 
815 } //namespace Isis
Isis::GisGeometry::points
int points() const
Get number of points in geometry.
Definition: GisGeometry.cpp:429
Isis::GisGeometry::IsisCube
@ IsisCube
An ISIS Cube is used to create the geometry.
Definition: GisGeometry.h:60
Isis::GisTopology::preparedGeometry
const GEOSPreparedGeometry * preparedGeometry(const GEOSGeometry *geom) const
Gets a GEOSPreparedGeometry from the given GEOSGeometry.
Definition: GisTopology.cpp:163
Isis::GisGeometry::overlaps
bool overlaps(const GisGeometry &target) const
Test for overlapping geometries.
Definition: GisGeometry.cpp:531
Isis::GisGeometry::envelope
GisGeometry * envelope() const
Computes the envelope or bounding box of this geometry.
Definition: GisGeometry.cpp:619
Isis::GisTopology
This class models GIS topology.
Definition: GisTopology.h:34
Isis::GisGeometry::area
double area() const
Computes the area of a geometry.
Definition: GisGeometry.cpp:359
Isis::GisGeometry::centroid
GisGeometry * centroid() const
Computes the centroid of the geometry and returns it as a new geometry.
Definition: GisGeometry.cpp:746
Isis::GisGeometry::m_type
Type m_type
Geometry type of GIS source.
Definition: GisGeometry.h:118
Isis::GisGeometry::GisGeometry
GisGeometry()
Fundamental constructor of an empty object.
Definition: GisGeometry.cpp:33
Isis::GisGeometry::WKT
@ WKT
The GEOS library WKT reader is used to create the geometry.
Definition: GisGeometry.h:58
Isis::Cube::readFootprint
ImagePolygon readFootprint() const
Read the footprint polygon for the Cube.
Definition: Cube.cpp:866
Isis::ImagePolygon::polyStr
std::string polyStr() const
Return a geos Multipolygon.
Definition: ImagePolygon.h:215
Isis::GisGeometry::isValidReason
QString isValidReason() const
Returns a string describing reason for invalid geometry.
Definition: GisGeometry.cpp:227
Isis::GisGeometry::intersection
GisGeometry * intersection(const GisGeometry &geom) const
Computes the intersection of two geometries.
Definition: GisGeometry.cpp:681
Isis::GisGeometry::type
Type type() const
Returns the type (origin) of the geometry.
Definition: GisGeometry.cpp:244
Isis::GisGeometry::geometry
const GEOSGeometry * geometry() const
Returns the GEOSGeometry object to extend functionality.
Definition: GisGeometry.cpp:291
Isis::GisGeometry::isEmpty
bool isEmpty() const
Tests for a defined but empty geometry.
Definition: GisGeometry.cpp:338
Isis::GisGeometry::makePrepared
GEOSPreparedGeometry const * makePrepared(const GEOSGeometry *geom) const
Creates a prepared geometry of current geometry.
Definition: GisGeometry.cpp:762
Isis::ImagePolygon
Create cube polygons, read/write polygons to blobs.
Definition: ImagePolygon.h:153
Isis::GisGeometry::setGeometry
void setGeometry(GEOSGeometry *geom)
Set the geometry directly taking ownership.
Definition: GisGeometry.cpp:187
Isis::GisTopology::clone
GEOSGeometry * clone(const GEOSGeometry *geom) const
Clones the given GEOSGeometry pointer.
Definition: GisTopology.cpp:148
Isis::GisGeometry
Encapsulation class provides support for GEOS-C API.
Definition: GisGeometry.h:50
Isis::GisGeometry::GeosGis
@ GeosGis
GEOS GIS. A geometry object cannot be created with this geometry type.
Definition: GisGeometry.h:61
Isis::GisGeometry::isDefined
bool isDefined() const
Determines if the current geometry is valid.
Definition: GisGeometry.cpp:200
Isis::GisGeometry::clone
GisGeometry * clone() const
Clones the contents of this geometry to a new instance.
Definition: GisGeometry.cpp:318
Isis::GisGeometry::equals
bool equals(const GisGeometry &target) const
Test if target and this geometry are equal.
Definition: GisGeometry.cpp:557
Isis::GisTopology::geomFromWKT
GEOSGeometry * geomFromWKT(const QString &wkt)
Reads in the geometry from the given well-known text formatted string.
Definition: GisTopology.cpp:128
Isis::GisTopology::instance
static GisTopology * instance()
Gets the singleton instance of this class.
Definition: GisTopology.cpp:83
Isis::GisTopology::geomFromWKB
GEOSGeometry * geomFromWKB(const QString &wkb)
Reads in the geometry from the given well-known binary formatted string.
Definition: GisTopology.cpp:100
Isis::GisGeometry::contains
bool contains(const GisGeometry &target) const
Test if the target geometry is contained within this geometry.
Definition: GisGeometry.cpp:478
Isis::GisGeometry::preparedGeometry
const GEOSPreparedGeometry * preparedGeometry() const
Returns special GEOS prepared geometry if it exists.
Definition: GisGeometry.cpp:305
Isis::GisTopology::destroy
void destroy(GEOSGeometry *geom) const
Destroys the given GEOS geometry.
Definition: GisTopology.cpp:232
Isis::GisGeometry::convexHull
GisGeometry * convexHull() const
Computes the convex hull of the geometry.
Definition: GisGeometry.cpp:636
Isis::GisGeometry::g_union
GisGeometry * g_union(const GisGeometry &geom) const
Computes the union of two geometries.
Definition: GisGeometry.cpp:701
Isis::GisGeometry::simplify
GisGeometry * simplify(const double &tolerance) const
Simplify complex or overdetermined geoemtry.
Definition: GisGeometry.cpp:659
Isis::GisGeometry::None
@ None
No geometry. A geometry object cannot be created with this geometry type.
Definition: GisGeometry.h:57
Isis::Cube
IO Handler for Isis Cubes.
Definition: Cube.h:167
Isis::GisGeometry::m_preparedGeom
GEOSPreparedGeometry const * m_preparedGeom
A prepared geometry from the GEOS library.
Definition: GisGeometry.h:120
Isis::GisGeometry::makePoint
GEOSGeometry * makePoint(const double x, const double y) const
Create a point geometry.
Definition: GisGeometry.cpp:779
Isis::GisGeometry::WKB
@ WKB
The GEOS library WKB reader is used to create the geometry.
Definition: GisGeometry.h:59
Isis::IException
Isis exception class.
Definition: IException.h:91
Isis::Null
const double Null
Value for an Isis Null pixel.
Definition: SpecialPixel.h:95
Isis::GisGeometry::intersects
bool intersects(const GisGeometry &target) const
Computes a new geometry from the intersection of the two geomtries.
Definition: GisGeometry.cpp:451
Isis::GisGeometry::typeToString
static QString typeToString(const Type &type)
Returns the type of the Geometry as a QString.
Definition: GisGeometry.cpp:277
Isis::GisGeometry::fromCube
GEOSGeometry * fromCube(Cube &cube) const
Reads Polygon from ISIS Cube and returns geometry from contents.
Definition: GisGeometry.cpp:795
Isis::IException::Programmer
@ Programmer
This error is for when a programmer made an API call that was illegal.
Definition: IException.h:146
std
Namespace for the standard library.
Isis::GisGeometry::intersectRatio
double intersectRatio(const GisGeometry &geom) const
Computes intersect ratio between two geometries.
Definition: GisGeometry.cpp:582
Isis::GisGeometry::m_geom
GEOSGeometry * m_geom
Pointer to GEOS-C opaque structure.
Definition: GisGeometry.h:119
Isis::GisGeometry::length
double length() const
Computes the length of a geometry.
Definition: GisGeometry.cpp:382
Isis::GisGeometry::destroy
void destroy()
Destroys the GEOS elements of this geometry object.
Definition: GisGeometry.cpp:806
Isis::GisGeometry::distance
double distance(const GisGeometry &target) const
Computes the distance between two geometries.
Definition: GisGeometry.cpp:406
Isis::GisGeometry::disjoint
bool disjoint(const GisGeometry &target) const
Tests for disjoint geometries.
Definition: GisGeometry.cpp:504
Isis::GisGeometry::isValid
bool isValid() const
Determines validity of the geometry contained in this object.
Definition: GisGeometry.cpp:213
Isis::Cube::open
void open(const QString &cfile, QString access="r")
This method will open an isis cube for reading or reading/writing.
Definition: Cube.cpp:627
Isis::GisGeometry::Type
Type
Source type of the geometry.
Definition: GisGeometry.h:56
Isis::GisGeometry::operator=
GisGeometry & operator=(GisGeometry const &geom)
Assignment operator for GISGeomtries.
Definition: GisGeometry.cpp:164
Isis::GisGeometry::~GisGeometry
virtual ~GisGeometry()
Destructor.
Definition: GisGeometry.cpp:151
Isis
This is free and unencumbered software released into the public domain.
Definition: Apollo.h:16