geo: upgrade geographic algorithms

lib/geodesic.js

@@ -0,0 +1,111 @@
+
+/**
+  geodesic functions, javascript versions of aviation formulary:
+  http://williams.best.vwh.net/avform.htm
+
+  note: all point values are in radians, not degrees
+**/
+
+/**
+  distance between point A and point B (in radians)
+**/
+function distance( a, b ){
+  return Math.acos( Math.sin( a.lat ) * Math.sin( b.lat ) +
+    Math.cos( a.lat ) * Math.cos( b.lat ) * Math.cos( a.lon - b.lon ));
+}
+
+/**
+  distance between point A and point B (in radians)
+
+  note: for very short distances this version is less susceptible to rounding error
+**/
+function distance2( a, b ){
+  return 2 * Math.asin( Math.sqrt( Math.pow( Math.sin(( a.lat - b.lat ) / 2 ), 2 ) +
+    Math.cos( a.lat ) * Math.cos( b.lat ) * Math.pow( Math.sin(( a.lon - b.lon ) / 2), 2)));
+}
+
+/**
+  course from point A and point B (in radians)
+**/
+function course( a, b, d ){
+  return Math.acos(
+    ( Math.sin( b.lat ) - Math.sin( a.lat ) * Math.cos( d )) /
+    ( Math.sin( d ) * Math.cos( a.lat ) )
+  );
+}
+
+/**
+  cross track error (distance off course)
+
+  (positive XTD means right of course, negative means left)
+**/
+function crossTrack( d, crs1, crs2, A, B, D ){
+
+  var calc = ( crs1 - crs2 );
+
+  // north pole / south pole
+  if( A.lat === +90 ){ calc = ( D.lon - B.lon ); }
+  if( A.lat === -90 ){ calc = ( B.lon - D.lon ); }
+
+  return Math.asin( Math.sin( d ) * Math.sin( calc ) );
+}
+
+/**
+  along track distance (the distance from A along the course towards B to the point abeam D)
+**/
+function alongTrack( d, xtd ){
+  return Math.acos( Math.cos( d ) / Math.cos( xtd ) );
+}
+
+/**
+  along track distance (the distance from A along the course towards B to the point abeam D)
+
+  note: for very short distances this version is less susceptible to rounding error
+**/
+function alongTrack2( d, xtd ){
+  return Math.asin( Math.sqrt( Math.pow( Math.sin( d ), 2 ) - Math.pow( Math.sin( xtd ), 2 ) ) / Math.cos( xtd ) );
+}
+
+/**
+  interpolate f (percent) of the distance d along path A-B
+**/
+function interpolate( d, f, a, b ){
+  var A = Math.sin( (1-f) * d ) / Math.sin( d );
+  var B = Math.sin( f * d ) / Math.sin( d );
+  var X = A * Math.cos( a.lat ) * Math.cos( a.lon ) + B * Math.cos( b.lat ) * Math.cos( b.lon );
+  var Y = A * Math.cos( a.lat ) * Math.sin( a.lon ) + B * Math.cos( b.lat ) * Math.sin( b.lon );
+  var Z = A * Math.sin( a.lat ) + B * Math.sin( b.lat );
+  return {
+    lat: Math.atan2( Z, Math.sqrt( Math.pow( X, 2 ) + Math.pow( Y, 2 ) ) ),
+    lon: Math.atan2( Y, X )
+  }
+}
+
+/**
+  calculate the distance a using b and C
+
+       c
+  A -------B
+   \       |
+    \      |
+     \b    |a
+      \    |
+       \   |
+        \  |
+         \C|
+          \|
+
+  Napier's rules: tan(a) = tan(b)*cos(C)
+**/
+function project( b, C ){
+  return Math.atan( Math.tan( b ) * Math.cos( C ) );
+}
+
+module.exports.distance = distance;
+module.exports.distance2 = distance2;
+module.exports.course = course;
+module.exports.crossTrack = crossTrack;
+module.exports.alongTrack = alongTrack;
+module.exports.alongTrack2 = alongTrack2;
+module.exports.interpolate = interpolate;
+module.exports.project = project;

lib/interpolate.js

@@ -0,0 +1,34 @@
+
+function interpolate( addrPoints, vertexDistance ){
+
+  // cycle through calculated addrPoints and interpolate a fractional housenumber
+  // value which would sit at this vertexDistance.
+  for( var x=0; x<addrPoints.length-1; x++ ){
+
+    var thisAddr = addrPoints[x],
+        nextAddr = addrPoints[x+1];
+
+    // the vertex vertexDistance is less that the lowest housenumber
+    // @extrapolation
+    if( vertexDistance < thisAddr.dist ){ return null; }
+
+    // vertex vertexDistance is between two house number vertexDistance
+    if( nextAddr.dist > vertexDistance ){
+      var ratio = ((vertexDistance - thisAddr.dist) / (nextAddr.dist - thisAddr.dist));
+
+      // invert ratio if the street was drawn from high house number to low
+      if( thisAddr.housenumber > nextAddr.housenumber ){ ratio = 1 - ratio; }
+
+      if( ratio >= 1 || ratio <= 0 ){ break; } // will result in a duplicate value
+      var minHouseNumber = Math.min( thisAddr.housenumber, nextAddr.housenumber );
+      var maxHouseNumber = Math.max( thisAddr.housenumber, nextAddr.housenumber );
+      return minHouseNumber + (( maxHouseNumber - minHouseNumber ) * ratio);
+    }
+
+    // else the vertex is greater than the highest housenumber
+    // @extrapolation
+  }
+  return null;
+}
+
+module.exports = interpolate;

lib/project.js

@@ -1,4 +1,6 @@
 
+var geodesic = require('./geodesic');
+
 /**
   project point p on to edge v, u
 
@@ -10,11 +12,16 @@
 **/
 
 function pointOnEdge( v, u, p ){
+
+  var lon_scale = Math.cos( p[1] * ( Math.PI / 180 ) );
+
   var bx = v[0] - u[0];
   var by = v[1] - u[1];
-  var sq = bx*bx + by*by;
 
-  var scale = ( sq > 0 ) ? (( (p[0] - u[0])*bx + (p[1] - u[1])*by ) / sq) : 0.0;
+  var bx2 = bx * lon_scale;
+  var sq = bx2*bx2 + by*by;
+
+  var scale = ( sq > 0 ) ? (( (p[0] - u[0])*lon_scale*bx2 + (p[1] - u[1])*by ) / sq) : 0.0;
 
   if( scale <= 0.0 ){
     bx = u[0];
@@ -75,14 +82,17 @@ function pointOnLine( linestring, p ){
 }
 
 /**
-  Calculate the distance between two points (in the same units as the points)
+  Calculate the distance between two points (in degrees)
 
   ◯ < - - - > ◯ ?
 
 **/
 
 function distance( a, b ){
-  return Math.hypot( a[0]-b[0], a[1]-b[1] );
+  return toDeg( geodesic.distance2(
+    { lon: toRad( a[0] ), lat: toRad( a[1] ) },
+    { lon: toRad( b[0] ), lat: toRad( b[1] ) }
+  ));
 }
 
 /*
@@ -214,6 +224,16 @@ function bearing( p1, p2 ){
   return toDeg(Math.atan2(a, b));
 }
 
+// deduplicate an array or coordinates in geojson [ [ lon, lat ] ... ] format.
+function dedupe( coordinates ){
+  return coordinates.filter( function( coord, i ){
+    if( 0 === i ){ return true; }
+    if( coord[0] !== coordinates[i-1][0] ){ return true; }
+    if( coord[1] !== coordinates[i-1][1] ){ return true; }
+    return false;
+  });
+}
+
 function toRad(degree) { return degree * Math.PI / 180; }
 function toDeg(radian) { return radian * 180 / Math.PI; }
 
@@ -226,3 +246,6 @@ module.exports.lineDistance = lineDistance;
 module.exports.sliceLineAtProjection = sliceLineAtProjection;
 module.exports.parity = parity;
 module.exports.bearing = bearing;
+module.exports.dedupe = dedupe;
+module.exports.toRad = toRad;
+module.exports.toDeg = toDeg;

stream/oa/augment.js

@@ -2,7 +2,8 @@
 var through = require('through2'),
     polyline = require('polyline'),
     project = require('../../lib/project'),
-    analyze = require('../../lib/analyze');
+    analyze = require('../../lib/analyze'),
+    interpolate = require('../../lib/interpolate');
 
 // polyline precision
 var PRECISION = 6;
@@ -22,7 +23,7 @@ function streamFactory(db, done){
 
     // decode polylines
     lookup.streets.forEach( function( street, i ){
-      street.coordinates = polyline.toGeoJSON(street.line, PRECISION).coordinates;
+      street.coordinates = project.dedupe( polyline.toGeoJSON(street.line, PRECISION).coordinates );
       distances[i] = []; // init array
     });
 
@@ -87,9 +88,11 @@ function streamFactory(db, done){
     }, this);
 
     // ensure distances are sorted by distance ascending
-    distances = distances.map( function( d ){
-      return d.sort( function( a, b ){
-        return a.dist > b.dist;
+    // this is important because now the distances and coordinates
+    // arrays will run from the start of the street to the end.
+    distances.forEach( function( d ){
+      d.sort( function( a, b ){
+        return ( a.dist > b.dist ) ? 1 : -1;
       });
     });
 
@@ -136,7 +139,7 @@ function streamFactory(db, done){
     lookup.streets.forEach( function( street, si ){
 
       // distance travelled along the line string
-      var distance = 0;
+      var vertexDistance = 0;
 
       // insert each point on linestring in table
       // note: this allows us to ignore the linestring and simply linearly
@@ -146,38 +149,32 @@ function streamFactory(db, done){
         // not a line, just a single point;
         if( 0 === i ){ return; }
 
-        // ignore successive duplicate points in linestring
-        var previousVertex = street.coordinates[i-1];
-        if( previousVertex && vertex[0] === previousVertex[0] && vertex[1] === previousVertex[1] ){
-          return;
-        }
-
         // distance along line to this vertex
         var edge = street.coordinates.slice(i-1, i+1);
         if( edge.length == 2 ){
-          distance += project.lineDistance( edge );
-        }
+          vertexDistance += project.lineDistance( edge );
+        } // else should not have else!
 
         // projected fractional housenumber(s)
         var housenumbers = [];
 
         // zigzag interpolation
         // (one vertex interpolation produced)
         if( street.scheme === 'zigzag' ){
-          housenumbers.push( interpolate( distances[si], distance ) );
+          housenumbers.push( interpolate( distances[si], vertexDistance ) );
         }
         // updown interpolation
         // (two vertex interpolations produced)
         else {
           // left side
           housenumbers.push( interpolate( distances[si].filter( function( d ){
             return d.parity === 'L';
-          }), distance ) );
+          }), vertexDistance ) );
 
           // right side
           housenumbers.push( interpolate( distances[si].filter( function( d ){
             return d.parity === 'R';
-          }), distance ) );
+          }), vertexDistance ) );
         }
 
         // insert point values in db
@@ -202,34 +199,4 @@ function streamFactory(db, done){
   });
 }
 
-function interpolate( distances, distance ){
-
-  // cycle through calculated distances and interpolate a fractional housenumber
-  // value which would sit at this vertex.
-  for( var x=0; x<distances.length-1; x++ ){
-
-    var thisDist = distances[x],
-        nextDist = distances[x+1];
-
-    // the vertex distance is less that the lowest housenumber
-    // @extrapolation
-    if( distance < thisDist.dist ){
-      break;
-    }
-
-    // vertex distance is between two house number distance
-    if( nextDist.dist > distance ){
-      var ratio = 1 - ((distance - thisDist.dist) / (nextDist.dist - thisDist.dist));
-      if( ratio >= 1 || ratio <= 0 ){ break; } // will result in a duplicate value
-      var minHouseNumber = Math.min( thisDist.housenumber, nextDist.housenumber );
-      var maxHouseNumber = Math.max( thisDist.housenumber, nextDist.housenumber );
-      return minHouseNumber + (( maxHouseNumber - minHouseNumber ) * ratio);
-    }
-
-    // else the vertex is greater than the highest housenumber
-    // @extrapolation
-  }
-  return null;
-}
-
 module.exports = streamFactory;

test/_unit.js

@@ -5,7 +5,8 @@ var common = {};
 var tests = [
   require('./interface.js'),
   require('./lib/analyze.js'),
-  require('./lib/project.js')
+  require('./lib/project.js'),
+  require('./lib/interpolate.js')
 ];
 
 tests.map(function(t) {