From 9b0211f4b8c0cbd834197c480299ae6b878b025d Mon Sep 17 00:00:00 2001
From: Thomas Rasch <trasch@users.noreply.github.com>
Date: Mon, 29 Jul 2024 16:57:42 +0200
Subject: [PATCH] Bugfix: Holes in polygons where not subtracted from the outer
 ring area (#61)

---
 Sources/GISTools/Algorithms/Area.swift        | 33 ++++++++++---------
 .../GISToolsTests/Algorithms/AreaTests.swift  | 32 ++++++++++++++++++
 2 files changed, 49 insertions(+), 16 deletions(-)
 create mode 100644 Tests/GISToolsTests/Algorithms/AreaTests.swift

diff --git a/Sources/GISTools/Algorithms/Area.swift b/Sources/GISTools/Algorithms/Area.swift
index a35d59a..2f9f5a9 100644
--- a/Sources/GISTools/Algorithms/Area.swift
+++ b/Sources/GISTools/Algorithms/Area.swift
@@ -1,5 +1,5 @@
 #if !os(Linux)
-import CoreLocation
+    import CoreLocation
 #endif
 import Foundation
 
@@ -11,12 +11,12 @@ extension Polygon {
     ///
     /// - Returns: The area of the outer ring minus the areas of the inner rings, in square meters.
     public var area: Double {
-        guard let outerRing = outerRing else { return 0.0 }
+        guard let outerRing else { return 0.0 }
 
         var area: Double = abs(outerRing.area)
 
-        if let innerRings = innerRings {
-            area += innerRings.reduce(0.0, { $0 + abs($1.area) })
+        if let innerRings {
+            area -= innerRings.reduce(0.0, { $0 + abs($1.area) })
         }
 
         return area
@@ -47,33 +47,34 @@ extension Ring {
     /// Laboratory, Pasadena, CA, June 2007 https://trs.jpl.nasa.gov/handle/2014/41271
     public var area: Double {
         let coordinates = projection == .epsg4326
-            ? self.coordinates
-            : self.coordinates.map({ $0.projected(to: .epsg4326) })
+            ? coordinates
+            : coordinates.map({ $0.projected(to: .epsg4326) })
 
-        var area: Double = 0.0
+        var area = 0.0
         let coordinatesCount = coordinates.count
 
         if coordinatesCount > 2 {
             for index in 0 ..< coordinatesCount {
-                let controlPoints: (Coordinate3D, Coordinate3D, Coordinate3D)
-
-                if index == coordinatesCount - 2 {
-                    controlPoints = (
+                let controlPoints: (Coordinate3D, Coordinate3D, Coordinate3D) = if index == coordinatesCount - 2 {
+                    (
                         coordinates[coordinatesCount - 2],
                         coordinates[coordinatesCount - 1],
-                        coordinates[0])
+                        coordinates[0]
+                    )
                 }
                 else if index == coordinatesCount - 1 {
-                    controlPoints = (
+                    (
                         coordinates[coordinatesCount - 1],
                         coordinates[0],
-                        coordinates[1])
+                        coordinates[1]
+                    )
                 }
                 else {
-                    controlPoints = (
+                    (
                         coordinates[index],
                         coordinates[index + 1],
-                        coordinates[index + 2])
+                        coordinates[index + 2]
+                    )
                 }
 
                 area += (controlPoints.2.longitude.degreesToRadians - controlPoints.0.longitude.degreesToRadians) * sin(controlPoints.1.latitude.degreesToRadians)
diff --git a/Tests/GISToolsTests/Algorithms/AreaTests.swift b/Tests/GISToolsTests/Algorithms/AreaTests.swift
new file mode 100644
index 0000000..39c40e7
--- /dev/null
+++ b/Tests/GISToolsTests/Algorithms/AreaTests.swift
@@ -0,0 +1,32 @@
+@testable import GISTools
+import XCTest
+
+final class AreaTests: XCTestCase {
+
+    func testArea() throws {
+        let polygon1 = try XCTUnwrap(Polygon([[
+            Coordinate3D(x: 0.0, y: 0.0),
+            Coordinate3D(x: 100.0, y: 0.0),
+            Coordinate3D(x: 100.0, y: 100.0),
+            Coordinate3D(x: 0.0, y: 100.0),
+            Coordinate3D(x: 0.0, y: 0.0)
+        ]]))
+        let polygon2 = try XCTUnwrap(Polygon([[
+            Coordinate3D(x: 0.0, y: 0.0),
+            Coordinate3D(x: 100.0, y: 0.0),
+            Coordinate3D(x: 100.0, y: 100.0),
+            Coordinate3D(x: 0.0, y: 100.0),
+            Coordinate3D(x: 0.0, y: 0.0)
+        ], [
+            Coordinate3D(x: 25.0, y: 25.0),
+            Coordinate3D(x: 75.0, y: 25.0),
+            Coordinate3D(x: 75.0, y: 75.0),
+            Coordinate3D(x: 25.0, y: 75.0),
+            Coordinate3D(x: 25.0, y: 25.0)
+        ]]))
+
+        XCTAssertEqual(polygon1.area, 10000.0, accuracy: 0.1)
+        XCTAssertEqual(polygon2.area, 7500.0, accuracy: 0.1)
+    }
+
+}