Merge pull request phayes#43 from Polzme/haversine

Adding Haversine method to calculate degrees between two coordinates.
imclab · Jul 20, 2012 · 2e9f112 · 2e9f112
2 parents 3a2a914 + 4096a6b
commit 2e9f112
Show file tree

Hide file tree

Showing 5 changed files with 67 additions and 22 deletions.
diff --git a/lib/geometry/Collection.class.php b/lib/geometry/Collection.class.php
@@ -165,12 +165,20 @@ public function length() {
 
   public function greatCircleLength($radius = 6378137) {
     $length = 0;
-    foreach ($this->components as $delta => $component) {
+    foreach ($this->components as $component) {
       $length += $component->greatCircleLength($radius);
     }
     return $length;
   }
 
+  public function haversineLength() {
+    $length = 0;
+    foreach ($this->components as $component) {
+      $length += $component->haversineLength();
+    }
+    return $length;
+  }
+
   public function dimension() {
     $dimension = 0;
     foreach ($this->components as $component) {

diff --git a/lib/geometry/LineString.class.php b/lib/geometry/LineString.class.php
@@ -77,31 +77,53 @@ public function length() {
 
   public function greatCircleLength($radius = 6378137) {
     $length = 0;
-    foreach ($this->getPoints() as $delta => $point) {
-      $previous_point = $this->geometryN($delta);
-      if ($previous_point) {
-        // Great circle method
-        $lat1 = deg2rad($point->getY());
-        $lat2 = deg2rad($previous_point->getY());
-        $lon1 = deg2rad($point->getX());
-        $lon2 = deg2rad($previous_point->getX());
-        $dlon = $lon2 - $lon1;
-        $length +=
-          $radius *
-            atan2(
-              sqrt(
-                pow(cos($lat2) * sin($dlon), 2) +
-                  pow(cos($lat1) * sin($lat2) - sin($lat1) * cos($lat2) * cos($dlon), 2)
-              )
-              ,
-              sin($lat1) * sin($lat2) +
-                cos($lat1) * cos($lat2) * cos($dlon)
-            );
-      }
+    $points = $this->getPoints();
+    for($i=0; $i<$this->numPoints()-1; $i++) {
+      $point = $points[$i];
+      $next_point = $points[$i+1];
+      if (!is_object($next_point)) {continue;}
+      // Great circle method
+      $lat1 = deg2rad($point->getY());
+      $lat2 = deg2rad($next_point->getY());
+      $lon1 = deg2rad($point->getX());
+      $lon2 = deg2rad($next_point->getX());
+      $dlon = $lon2 - $lon1;
+      $length +=
+        $radius *
+          atan2(
+            sqrt(
+              pow(cos($lat2) * sin($dlon), 2) +
+                pow(cos($lat1) * sin($lat2) - sin($lat1) * cos($lat2) * cos($dlon), 2)
+            )
+            ,
+            sin($lat1) * sin($lat2) +
+              cos($lat1) * cos($lat2) * cos($dlon)
+          );
     }
+    // Returns length in meters.
     return $length;
   }
 
+  public function haversineLength() {
+    $degrees = 0;
+    $points = $this->getPoints();
+    for($i=0; $i<$this->numPoints()-1; $i++) {
+      $point = $points[$i];
+      $next_point = $points[$i+1];
+      if (!is_object($next_point)) {continue;}
+      $degree = rad2deg(
+        acos(
+          sin(deg2rad($point->getY())) * sin(deg2rad($next_point->getY())) +
+            cos(deg2rad($point->getY())) * cos(deg2rad($next_point->getY())) *
+              cos(deg2rad(abs($point->getX() - $next_point->getX())))
+        )
+      );
+      $degrees += $degree;
+    }
+    // Returns degrees
+    return $degrees;
+  }
+
   public function explode() {
     $parts = array();
     $points = $this->getPoints();

diff --git a/lib/geometry/Point.class.php b/lib/geometry/Point.class.php
@@ -105,6 +105,10 @@ public function greatCircleLength() {
     return 0;
   }
 
+  public function haversineLength() {
+    return 0;
+  }
+
   // The boundary of a point is itself
   public function boundary() {
     return $this;

diff --git a/tests/input/opposite.gpx b/tests/input/opposite.gpx
@@ -0,0 +1,10 @@
+<?xml version="1.0" encoding="UTF-8" standalone="no" ?>
+<gpx xmlns="http://www.topografix.com/GPX/1/1" creator="" version="1.1" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsi:schemaLocation="http://www.topografix.com/GPX/1/1 http://www.topografix.com/GPX/1/1/gpx.xsd">
+ <trk>
+  <name>Mons and its opposite</name>
+  <trkseg>
+   <trkpt lat="50.45409" lon="3.95626"></trkpt>
+   <trkpt lat="-50.45409" lon="-176.04374"></trkpt>
+  </trkseg>
+ </trk>
+</gpx>
diff --git a/tests/test.php b/tests/test.php
@@ -40,6 +40,7 @@ function test_geometry($geometry) {
   $geometry->centroid();
   $geometry->length();
   $geometry->greatCircleLength();
+  $geometry->haversineLength();
   $geometry->y();
   $geometry->x();
   $geometry->numGeometries();