Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add polygon area calculation #21

Merged
merged 15 commits into from
Jan 8, 2018
Merged
6 changes: 6 additions & 0 deletions Turf.xcodeproj/project.pbxproj
Original file line number Diff line number Diff line change
Expand Up @@ -7,6 +7,8 @@
objects = {

/* Begin PBXBuildFile section */
052521D81FC47A1300DD266A /* polygon.geojson in Resources */ = {isa = PBXBuildFile; fileRef = 052521D71FC47A0300DD266A /* polygon.geojson */; };
052521D91FC47A1900DD266A /* polygon.geojson in Resources */ = {isa = PBXBuildFile; fileRef = 052521D71FC47A0300DD266A /* polygon.geojson */; };
353E9B101F3E093A007CFA23 /* Turf.framework in Frameworks */ = {isa = PBXBuildFile; fileRef = 353E9B071F3E093A007CFA23 /* Turf.framework */; };
353E9B1E1F3E09D3007CFA23 /* CoreLocation.swift in Sources */ = {isa = PBXBuildFile; fileRef = 356D24531F17948C003BBB9D /* CoreLocation.swift */; };
353E9B1F1F3E09D8007CFA23 /* Turf.swift in Sources */ = {isa = PBXBuildFile; fileRef = 35650B0A1F150DDB00B5C158 /* Turf.swift */; };
Expand Down Expand Up @@ -41,6 +43,7 @@
/* End PBXContainerItemProxy section */

/* Begin PBXFileReference section */
052521D71FC47A0300DD266A /* polygon.geojson */ = {isa = PBXFileReference; lastKnownFileType = text; path = polygon.geojson; sourceTree = "<group>"; };
353E9B071F3E093A007CFA23 /* Turf.framework */ = {isa = PBXFileReference; explicitFileType = wrapper.framework; includeInIndex = 0; path = Turf.framework; sourceTree = BUILT_PRODUCTS_DIR; };
353E9B0F1F3E093A007CFA23 /* TurfMacTests.xctest */ = {isa = PBXFileReference; explicitFileType = wrapper.cfbundle; includeInIndex = 0; path = TurfMacTests.xctest; sourceTree = BUILT_PRODUCTS_DIR; };
35650AF01F150DC500B5C158 /* Turf.framework */ = {isa = PBXFileReference; explicitFileType = wrapper.framework; includeInIndex = 0; path = Turf.framework; sourceTree = BUILT_PRODUCTS_DIR; };
Expand Down Expand Up @@ -135,6 +138,7 @@
isa = PBXGroup;
children = (
356D24581F179B72003BBB9D /* dc-line.geojson */,
052521D71FC47A0300DD266A /* polygon.geojson */,
);
path = Fixtures;
sourceTree = "<group>";
Expand Down Expand Up @@ -291,6 +295,7 @@
buildActionMask = 2147483647;
files = (
35CB7F6F1F798A51008A18C8 /* dc-line.geojson in Resources */,
052521D91FC47A1900DD266A /* polygon.geojson in Resources */,
);
runOnlyForDeploymentPostprocessing = 0;
};
Expand All @@ -306,6 +311,7 @@
buildActionMask = 2147483647;
files = (
356D24591F179B72003BBB9D /* dc-line.geojson in Resources */,
052521D81FC47A1300DD266A /* polygon.geojson in Resources */,
);
runOnlyForDeploymentPostprocessing = 0;
};
Expand Down
61 changes: 60 additions & 1 deletion Turf/Turf.swift
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,9 @@ public typealias LocationRadians = Double
public typealias RadianDistance = Double
public typealias RadianDirection = Double

let metersPerRadian = 6_373_000.0

let metersPerRadian: CLLocationDistance = 6_373_000.0
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Apparently this value of 6 373 000 m came from an old version of turf-distance that I originally ported to an internal Swift application, before it got copied into another internal Swift application, then copied into the navigation SDK, then moved here. It’s close to the value of 6 372 797.560 856 that Osmium describes as “Earth’s quadratic mean radius for WGS84”.

These days, Turf.js uses a spherical approximation of 6 671 008.8 m for its radian-to-meter conversion and Haversine formula. The Haversine formula was always meant to be used with a spherical meters-per-radian value. We should probably change this value to match Turf.js. 🙀

/ref Turfjs/turf#978 Turfjs/turf#1012 Turfjs/turf#1176
/cc @frederoni @bsudekum

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Since the polygon area calculation doesn’t depend on metersPerRadian, we can treat this change as tail work: #26.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If we're going to change metersPerRadian to 6 671 008.8 m, should this happen in a different PR since it would globally impact this library? Or are you ok with me making the change here?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yeah, let’s do it in a separate PR for #26, so that any fallout is easier to track down.

In the meantime, can you add some comments explaining why these two constants differ?

let equatorialRadius: CLLocationDistance = 6378137
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Nit: similar to in metersPerRadian, use underscores to group the digits.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

equatorialRadius is currently set to the equatorial radius according to WGS84 and the IUGG. This value is used in OSRM and mbgl for Web Mercator calculations.

Since different radii are used in different formulas and the errors can sometimes build up, let’s document which is which.


/**
A `RadianCoordinate2D` is a coordinate represented in radians as opposed to
Expand Down Expand Up @@ -285,3 +287,60 @@ public struct Polyline {
return closestCoordinate
}
}

struct Ring {
var coordinates: [CLLocationCoordinate2D]
Copy link
Contributor

@frederoni frederoni Nov 22, 2017

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Capturing from chat: The coordinates property on a GeoJSON Polygon should be represented as a [[CLLocationCoordate2D]]

Even better now encapsulated in a Ring.


/**
* Calculate the approximate area of the polygon were it projected onto the earth.
* Note that this area will be positive if ring is oriented clockwise, otherwise it will be negative.
*
* Reference:
* Robert. G. Chamberlain and William H. Duquette, "Some Algorithms for Polygons on a Sphere", JPL Publication 07-03, Jet Propulsion
* Laboratory, Pasadena, CA, June 2007 http://trs-new.jpl.nasa.gov/dspace/handle/2014/40409
*
*/
internal func area() -> Double {
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Turn this method into a computed property to reduce the friction of using it.

var area: Double = 0
let coordinatesCount: Int = coordinates.count

if coordinatesCount > 2 {
for index in 0..<coordinatesCount {

let controlPoints: (CLLocationCoordinate2D, CLLocationCoordinate2D, CLLocationCoordinate2D)

if index == coordinatesCount - 2 {
controlPoints = (coordinates[coordinatesCount - 2],
coordinates[coordinatesCount - 1],
coordinates[0])
} else if index == coordinatesCount - 1 {
controlPoints = (coordinates[coordinatesCount - 1],
coordinates[0],
coordinates[1])
} else {
controlPoints = (coordinates[index],
coordinates[index + 1],
coordinates[index + 2])
}

area += (controlPoints.2.longitude.toRadians() - controlPoints.0.longitude.toRadians()) * sin(controlPoints.1.latitude.toRadians())
}

area *= equatorialRadius * equatorialRadius / 2
}
return area
}
}

public struct Polygon {
var outerRing: Ring
var innerRings: [Ring]

// Ported from https://github.com/Turfjs/turf/blob/a94151418cb969868fdb42955a19a133512da0fd/packages/turf-area/index.js

var area: Double {
return abs(outerRing.area()) - innerRings
.map { abs($0.area()) }
.reduce(0, +)
}
}
19 changes: 19 additions & 0 deletions TurfTests/Fixtures/polygon.geojson
Original file line number Diff line number Diff line change
@@ -0,0 +1,19 @@
{
"type": "Feature",
"properties": {},
"geometry": {
"type": "Polygon",
"coordinates": [
[
[125, -15],
[113, -22],
[117, -37],
[130, -33],
[148, -39],
[154, -27],
[144, -15],
[125, -15]
]
]
}
}
14 changes: 14 additions & 0 deletions TurfTests/TurfTests.swift
Original file line number Diff line number Diff line change
Expand Up @@ -286,4 +286,18 @@ class TurfTests: XCTestCase {
let b = radian.toDegrees()
XCTAssertEqual(b, 229, accuracy: 1)
}

func testPolygonArea() {
let json = Fixture.JSONFromFileNamed(name: "polygon")
let geometry = json["geometry"] as! [String: Any]
let geoJSONCoordinates = geometry["coordinates"] as! [[[Double]]]
let allRings = geoJSONCoordinates.map {
$0.map { CLLocationCoordinate2D(latitude: $0[1], longitude: $0[0]) }
}
let outerRing = Ring(coordinates: allRings.first!)

let polygon = Polygon(outerRing: outerRing, innerRings: [])
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What about the inner rings in allRings?


XCTAssertEqual(polygon.area, 7766240997209, accuracy: 0.1)
}
}