Curved geometries: Implement getArea

include/geos/algorithm/Area.h

@@ -23,12 +23,19 @@
 #include <geos/geom/CoordinateSequence.h>
 
 namespace geos {
+
+namespace geom {
+class Curve;
+}
+
 namespace algorithm { // geos::algorithm
 
 
 class GEOS_DLL Area {
 public:
 
+    static double ofClosedCurve(const geom::Curve& ring);
+
     /**
     * Computes the area for a ring.
     *

include/geos/geom/CircularArc.h

@@ -108,6 +108,16 @@ class GEOS_DLL CircularArc {
         return angle()*radius();
     }
 
+    double getArea() const {
+        if (isLinear()) {
+            return 0;
+        }
+
+        auto R = radius();
+        auto theta = angle();
+        return R*R/2*(theta - std::sin(theta));
+    }
+
     double theta0() const {
         return std::atan2(p0.y - center().y, p0.x - center().x);
     }

include/geos/geom/CircularString.h

@@ -41,6 +41,10 @@ class GEOS_DLL CircularString : public SimpleCurve {
         return true;
     }
 
+    bool isCurved() const override {
+        return true;
+    }
+
     std::unique_ptr<CircularString> reverse() const
     {
         return std::unique_ptr<CircularString>(reverseImpl());

include/geos/geom/LineString.h

@@ -92,6 +92,10 @@ class GEOS_DLL LineString: public SimpleCurve {
 
     double getLength() const override;
 
+    bool isCurved() const override {
+        return false;
+    }
+
     /**
      * Creates a LineString whose coordinates are in the reverse
      * order of this object's

include/geos/geom/SimpleCurve.h

@@ -94,6 +94,8 @@ class GEOS_DLL SimpleCurve : public Curve {
 
     virtual bool isCoordinate(CoordinateXY& pt) const;
 
+    virtual bool isCurved() const = 0;
+
     bool isEmpty() const override;
 
     /** \brief

src/algorithm/Area.cpp

@@ -20,6 +20,10 @@
 #include <vector>
 
 #include <geos/algorithm/Area.h>
+#include <geos/geom/CircularArc.h>
+#include <geos/geom/Curve.h>
+#include <geos/geom/SimpleCurve.h>
+#include <geos/util/IllegalArgumentException.h>
 
 using geos::geom::CoordinateXY;
 
@@ -93,7 +97,58 @@ Area::ofRingSigned(const geom::CoordinateSequence* ring)
     return sum / 2.0;
 }
 
+double
+Area::ofClosedCurve(const geom::Curve& ring) {
+    if (!ring.isClosed()) {
+        throw util::IllegalArgumentException("Argument is not closed");
+    }
+
+    double sum = 0;
+
+    for (std::size_t i = 0; i < ring.getNumCurves(); i++) {
+        const geom::SimpleCurve& section = *ring.getCurveN(i);
+
+        if (section.isEmpty()) {
+            continue;
+        }
+
+        const geom::CoordinateSequence& coords = *section.getCoordinatesRO();
+
+        if (section.isCurved()) {
+            for (std::size_t j = 2; j < coords.size(); j += 2) {
+                const CoordinateXY& p0 = coords.getAt<CoordinateXY>(j-2);
+                const CoordinateXY& p1 = coords.getAt<CoordinateXY>(j-1);
+                const CoordinateXY& p2 = coords.getAt<CoordinateXY>(j);
 
+                double triangleArea = 0.5*(p0.x*p2.y - p2.x*p0.y);
+                sum += triangleArea;
+
+                geom::CircularArc arc(p0, p1, p2);
+                if (arc.isLinear()) {
+                    continue;
+                }
+
+                double circularSegmentArea = arc.getArea();
+
+                if (algorithm::Orientation::index(p0, p2, p1) == algorithm::Orientation::CLOCKWISE) {
+                    sum += circularSegmentArea;
+                } else {
+                    sum -= circularSegmentArea;
+                }
+            }
+        } else {
+            for (std::size_t j = 1; j < coords.size(); j++) {
+                const CoordinateXY& p0 = coords.getAt<CoordinateXY>(j-1);
+                const CoordinateXY& p1 = coords.getAt<CoordinateXY>(j);
+
+                double triangleArea = 0.5*(p0.x*p1.y - p1.x*p0.y);
+                sum += triangleArea;
+            }
+        }
+    }
+
+    return std::abs(sum);
+}
 
 } // namespace geos.algorithm
 } //namespace geos

src/geom/CurvePolygon.cpp

@@ -12,6 +12,7 @@
  *
  **********************************************************************/
 
+#include <geos/algorithm/Area.h>
 #include <geos/geom/Curve.h>
 #include <geos/geom/CurvePolygon.h>
 #include <geos/geom/CoordinateSequence.h>
@@ -53,7 +54,11 @@ namespace geom {
     }
 
     double CurvePolygon::getArea() const {
-        throw util::UnsupportedOperationException();
+        double sum = algorithm::Area::ofClosedCurve(*shell);
+        for (const auto& hole : holes) {
+            sum -= algorithm::Area::ofClosedCurve(*hole);
+        }
+        return sum;
     }
 
     bool CurvePolygon::hasCurvedComponents() const {

tests/unit/algorithm/AreaTest.cpp

@@ -19,6 +19,7 @@
 
 using namespace geos;
 using namespace geos::geom;
+using geos::algorithm::Area;
 
 namespace tut {
 //
@@ -33,7 +34,7 @@ struct test_area_data {
     geos::io::WKTReader reader_;
     test_area_data():
         geom_(nullptr),
-        pm_(1),
+        pm_(),
         factory_(GeometryFactory::create(&pm_, 0)), reader_(factory_.get())
     {
         assert(nullptr == geom_);
@@ -48,19 +49,23 @@ struct test_area_data {
     void
     checkAreaOfRing(std::string wkt, double expectedArea)
     {
-        std::unique_ptr<Geometry> lineGeom(reader_.read(wkt));
-        std::unique_ptr<LineString> line(dynamic_cast<LineString*>(lineGeom.release()));
-        ensure(nullptr != line.get());
-        const CoordinateSequence* ringSeq = line->getCoordinatesRO();
+        auto ringGeom = reader_.read<Curve>(wkt);
 
-        std::vector<Coordinate> ringCoords;
-        ringSeq->toVector(ringCoords);
+        if (const LineString* line = dynamic_cast<const LineString*>(ringGeom.get())) {
+            const CoordinateSequence* ringSeq = line->getCoordinatesRO();
 
-        double actual1 = algorithm::Area::ofRing(ringCoords);
-        double actual2 = algorithm::Area::ofRing(ringSeq);
+            std::vector<Coordinate> ringCoords;
+            ringSeq->toVector(ringCoords);
 
-        ensure_equals(actual1, expectedArea);
-        ensure_equals(actual2, expectedArea);
+            double actual1 = algorithm::Area::ofRing(ringCoords);
+            double actual2 = algorithm::Area::ofRing(ringSeq);
+
+            ensure_equals(actual1, expectedArea);
+            ensure_equals(actual2, expectedArea);
+        }
+
+        double actual3 = algorithm::Area::ofClosedCurve(*ringGeom);
+        ensure_equals("Area::ofClosedCurve", actual3, expectedArea, 1e-6);
     }
 
     void
@@ -115,6 +120,41 @@ void object::test<3>
     checkAreaOfRingSigned("LINESTRING (100 200, 100 100, 200 100, 200 200, 100 200)", -10000.0);
 }
 
+template<>
+template<>
+void object::test<4>
+()
+{
+    checkAreaOfRing("CIRCULARSTRING (0 0, 2 2, 4 0, 2 -2, 0 0)", 4*MATH_PI);
+}
+
+template<>
+template<>
+void object::test<5>
+()
+{
+    checkAreaOfRing("COMPOUNDCURVE (CIRCULARSTRING (0 0, 2 2, 4 0), (4 0, 0 0))", 2*MATH_PI);
+}
+
+template<>
+template<>
+void object::test<6>
+()
+{
+    // expected ares from PostGIS after ST_CurveToLine(geom, 1e-13, 1)
+    checkAreaOfRing("CIRCULARSTRING (0 0, 2 2, 4 0, 2 1, 0 0)", 3.48759);
+}
+
+template<>
+template<>
+void object::test<7>
+()
+{
+    // expected ares from PostGIS after ST_CurveToLine(geom, 1e-13, 1)
+    checkAreaOfRing("COMPOUNDCURVE (CIRCULARSTRING (0 0, 2 0, 2 1, 2 3, 4 3, 4 5, 1 4, 0.5 0.8, 0 0))", 11.243342);
+    checkAreaOfRing("COMPOUNDCURVE (CIRCULARSTRING (0 0, 2 0, 2 1, 2 3, 4 3), (4 3, 4 5, 1 4, 0 0))", 9.321903);
+}
+
 
 } // namespace tut
 

tests/unit/geom/CircularArcTest.cpp

@@ -73,4 +73,19 @@ void object::test<2>()
     checkLength({1.6, 0.4}, {1.6, 0.5}, {1.7, 1}, 0.6122445326877711);
 }
 
+
+// test getArea()
+template<>
+template<>
+void object::test<3>()
+{
+    ensure_equals("half circle, R=2", CircularArc({-2, 0}, {0, 2}, {2, 0}).getArea(), MATH_PI*2);
+
+    ensure_equals("full circle, R=3", CircularArc({-3, 0}, {3, 0}, {-3, 0}).getArea(), MATH_PI*3*3);
+
+    ensure_equals("3/4, mouth up, R=2", CircularArc({-std::sqrt(2), std::sqrt(2)}, {0, -2}, {std::sqrt(2), std::sqrt(2)}).getArea(), MATH_PI*4 - 2*(MATH_PI/2-1), 1e-8);
+
+    ensure_equals("1/4, pointing up, R=2", CircularArc({-std::sqrt(2), std::sqrt(2)}, {0, 2}, {std::sqrt(2), std::sqrt(2)}).getArea(), 2*(MATH_PI/2-1), 1e-8);
+}
+
 }

tests/unit/geom/CurvePolygonTest.cpp

@@ -74,7 +74,7 @@ void object::test<1>()
     ensure("getCoordinates", cp->getCoordinates()->isEmpty());
     ensure("getCoordinate", cp->getCoordinate() == nullptr);
 
-    ensure_THROW(cp->getArea(), geos::util::UnsupportedOperationException);
+    ensure_equals("getArea", cp->getArea(), 0.0);
     ensure_equals("getLength", cp->getLength(), 0.0);
 }
 
@@ -90,7 +90,7 @@ void object::test<2>()
 
     // Geometry size functions
     ensure("isEmpty", !cp_->isEmpty());
-    ensure_THROW(cp_->getArea(), geos::util::UnsupportedOperationException);
+    ensure_equals("getArea", cp_->getArea(), 9.0526564962674, 1e-8); // expected value from PostGIS with ST_CurveToLine(geom, 1e-13, 1)
     ensure_equals("getLength", cp_->getLength(), 19.236489581872586, 1e-8);
     ensure_equals("getNumGeometries", cp_->getNumGeometries(), 1u);
     ensure_equals("getNumPoints", cp_->getNumPoints(), 14u);

tests/unit/geom/MultiSurfaceTest.cpp

@@ -85,7 +85,7 @@ void object::test<2>()
 
     // Geometry size functions
     ensure("isEmpty", !ms_->isEmpty());
-    ensure_THROW(ms_->getArea(), geos::util::UnsupportedOperationException);
+    ensure_equals("getArea", ms_->getArea(), 4.141592653589132, 1e-6); // expected value from PostGIS with ST_CurveToLine(geom, 1e-13, 1)
     ensure_equals("getLength", ms_->getLength(), 10.283185307179586);
     ensure_equals("getNumGeometries", ms_->getNumGeometries(), 2u);
     ensure_equals("getNumPoints", ms_->getNumPoints(), 10u);