WGS84 Compliant Great Circle Distance Calculation

Tue, 08 Mar 2011 17:43:11 GMT

Hi,

I am trying to create a system that shows the distance between two coordinates. So far I've managed to get part way there using the Haversine formula to calculate great circle distance based on a spherical earth model. However, the requirements have just been tightened up and I now need to try and find a way of using a WGS84 compliant method (oblate spheroid model) which is much more complex.

I've found a suitable formula (the Vincenty algorithm http://www.movable-type.co.uk/scripts/latlong-vincenty.html ) but have yet to figure out how to get this into QlikView! The database the information is coming out of is clunky, in-house written (and those that wrote it no longer work for the company) and not Oracle based so ideally I need to figure this out without involving the source database. I have a a list of journeys in a qvd file with to and from coordinates in them and I need to calculate the distance ... any ideas? 🙂

Ever hopeful!

Emma

WGS84 Compliant Great Circle Distance Calculation

Wed, 09 Mar 2011 10:19:08 GMT

hi,

Take a look at the source of the linked HTML page == there is Javascript with the formula 🙂 !

Makea QlikView file with a module written in Jscript. Paste the code there.


/**
 * Calculates geodetic distance between two points specified by latitude/longitude using
 * Vincenty inverse formula for ellipsoids
 *
 * @param {Number} lat1, lon1: first point in decimal degrees
 * @param {Number} lat2, lon2: second point in decimal degrees
 * @returns (Number} distance in metres between points
 */
function distVincenty(lat1, lon1, lat2, lon2) {
 var a = 6378137, b = 6356752.3142, f = 1/298.257223563; // WGS-84 ellipsoid params
 var L = (lon2-lon1).toRad();
 var U1 = Math.atan((1-f) * Math.tan(lat1.toRad()));
 var U2 = Math.atan((1-f) * Math.tan(lat2.toRad()));
 var sinU1 = Math.sin(U1), cosU1 = Math.cos(U1);
 var sinU2 = Math.sin(U2), cosU2 = Math.cos(U2);

 var lambda = L, lambdaP, iterLimit = 100;
 do {
 var sinLambda = Math.sin(lambda), cosLambda = Math.cos(lambda);
 var sinSigma = Math.sqrt((cosU2*sinLambda) * (cosU2*sinLambda) +
 (cosU1*sinU2-sinU1*cosU2*cosLambda) * (cosU1*sinU2-sinU1*cosU2*cosLambda));
 if (sinSigma==0) return 0; // co-incident points
 var cosSigma = sinU1*sinU2 + cosU1*cosU2*cosLambda;
 var sigma = Math.atan2(sinSigma, cosSigma);
 var sinAlpha = cosU1 * cosU2 * sinLambda / sinSigma;
 var cosSqAlpha = 1 - sinAlpha*sinAlpha;
 var cos2SigmaM = cosSigma - 2*sinU1*sinU2/cosSqAlpha;
 if (isNaN(cos2SigmaM)) cos2SigmaM = 0; // equatorial line: cosSqAlpha=0 (§6)
 var C = f/16*cosSqAlpha*(4+f*(4-3*cosSqAlpha));
 lambdaP = lambda;
 lambda = L + (1-C) * f * sinAlpha *
 (sigma + C*sinSigma*(cos2SigmaM+C*cosSigma*(-1+2*cos2SigmaM*cos2SigmaM)));
 } while (Math.abs(lambda-lambdaP) > 1e-12 && --iterLimit>0);
 if (iterLimit==0) return NaN // formula failed to converge

 var uSq = cosSqAlpha * (a*a - b*b) / (b*b);
 var A = 1 + uSq/16384*(4096+uSq*(-768+uSq*(320-175*uSq)));
 var B = uSq/1024 * (256+uSq*(-128+uSq*(74-47*uSq)));
 var deltaSigma = B*sinSigma*(cos2SigmaM+B/4*(cosSigma*(-1+2*cos2SigmaM*cos2SigmaM)-
 B/6*cos2SigmaM*(-3+4*sinSigma*sinSigma)*(-3+4*cos2SigmaM*cos2SigmaM)));
 var s = b*A*(sigma-deltaSigma);

 s = s.toFixed(3); // round to 1mm precision
 return s;

 // note: to return initial/final bearings in addition to distance, use something like:
 var fwdAz = Math.atan2(cosU2*sinLambda, cosU1*sinU2-sinU1*cosU2*cosLambda);
 var revAz = Math.atan2(cosU1*sinLambda, -sinU1*cosU2+cosU1*sinU2*cosLambda);
 return { distance: s, initialBearing: fwdAz.toDeg(), finalBearing: revAz.toDeg() };
}

/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

/** extend Number object with methods for converting degrees/radians */
toRad = function() { // convert degrees to radians
 return this * Math.PI / 180;
}
toDeg = function() { // convert radians to degrees
 return this * 180 / Math.PI;
}

/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/* Latitude/longitude spherical geodesy formulae & scripts (c) Chris Veness 2002-2010 */
/* - www.movable-type.co.uk/scripts/latlong.html */
/* */
/* Sample usage: */
/* var p1 = new LatLon(51.5136, -0.0983); */
/* var p2 = new LatLon(51.4778, -0.0015); */
/* var dist = p1.distanceTo(p2); // in km */
/* var brng = p1.bearingTo(p2); // in degrees clockwise from north */
/* ... etc */
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/* Note that minimal error checking is performed in this example code! */
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */


/**
 * Creates a point on the earth's surface at the supplied latitude / longitude
 *
 * @constructor
 * @param {Number} lat: latitude in numeric degrees
 * @param {Number} lon: longitude in numeric degrees
 * @param {Number} [rad=6371]: radius of earth if different value is required from standard 6,371km
 */
function LatLon(lat, lon, rad) {
 if (typeof(rad) == 'undefined') rad = 6371; // earth's mean radius in km
 // only accept numbers or valid numeric strings
 this._lat = typeof(lat)=='number' ? lat : typeof(lat)=='string' && lat.trim()!='' ? +lat : NaN;
 this._lon = typeof(lat)=='number' ? lon : typeof(lon)=='string' && lon.trim()!='' ? +lon : NaN;
 this._radius = typeof(rad)=='number' ? rad : typeof(rad)=='string' && trim(lon)!='' ? +rad : NaN;
}


/**
 * Returns the distance from this point to the supplied point, in km
 * (using Haversine formula)
 *
 * from: Haversine formula - R. W. Sinnott, "Virtues of the Haversine",
 * Sky and Telescope, vol 68, no 2, 1984
 *
 * @param {LatLon} point: Latitude/longitude of destination point
 * @param {Number} [precision=4]: no of significant digits to use for returned value
 * @returns {Number} Distance in km between this point and destination point
 */
distanceTo = function(point, precision) {
 // default 4 sig figs reflects typical 0.3% accuracy of spherical model
 if (typeof precision == 'undefined') precision = 4;

 var R = this._radius;
 var lat1 = this._lat.toRad(), lon1 = this._lon.toRad();
 var lat2 = point._lat.toRad(), lon2 = point._lon.toRad();
 var dLat = lat2 - lat1;
 var dLon = lon2 - lon1;

 var a = Math.sin(dLat/2) * Math.sin(dLat/2) +
 Math.cos(lat1) * Math.cos(lat2) *
 Math.sin(dLon/2) * Math.sin(dLon/2);
 var c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1-a));
 var d = R * c;
 return d.toPrecisionFixed(precision);
}


/**
 * Returns the (initial) bearing from this point to the supplied point, in degrees
 * see http://williams.best.vwh.net/avform.htm#Crs
 *
 * @param {LatLon} point: Latitude/longitude of destination point
 * @returns {Number} Initial bearing in degrees from North
 */
bearingTo = function(point) {
 var lat1 = this._lat.toRad(), lat2 = point._lat.toRad();
 var dLon = (point._lon-this._lon).toRad();

 var y = Math.sin(dLon) * Math.cos(lat2);
 var x = Math.cos(lat1)*Math.sin(lat2) -
 Math.sin(lat1)*Math.cos(lat2)*Math.cos(dLon);
 var brng = Math.atan2(y, x);

 return (brng.toDeg()+360) % 360;
}


/**
 * Returns final bearing arriving at supplied destination point from this point; the final bearing
 * will differ from the initial bearing by varying degrees according to distance and latitude
 *
 * @param {LatLon} point: Latitude/longitude of destination point
 * @returns {Number} Final bearing in degrees from North
 */
finalBearingTo = function(point) {
 // get initial bearing from supplied point back to this point...
 var lat1 = point._lat.toRad(), lat2 = this._lat.toRad();
 var dLon = (this._lon-point._lon).toRad();

 var y = Math.sin(dLon) * Math.cos(lat2);
 var x = Math.cos(lat1)*Math.sin(lat2) -
 Math.sin(lat1)*Math.cos(lat2)*Math.cos(dLon);
 var brng = Math.atan2(y, x);

 // ... & reverse it by adding 180°
 return (brng.toDeg()+180) % 360;
}


/**
 * Returns the midpoint between this point and the supplied point.
 * see http://mathforum.org/library/drmath/view/51822.html for derivation
 *
 * @param {LatLon} point: Latitude/longitude of destination point
 * @returns {LatLon} Midpoint between this point and the supplied point
 */
midpointTo = function(point) {
 lat1 = this._lat.toRad(), lon1 = this._lon.toRad();
 lat2 = point._lat.toRad();
 var dLon = (point._lon-this._lon).toRad();

 var Bx = Math.cos(lat2) * Math.cos(dLon);
 var By = Math.cos(lat2) * Math.sin(dLon);

 lat3 = Math.atan2(Math.sin(lat1)+Math.sin(lat2),
 Math.sqrt( (Math.cos(lat1)+Bx)*(Math.cos(lat1)+Bx) + By*By) );
 lon3 = lon1 + Math.atan2(By, Math.cos(lat1) + Bx);

 return new LatLon(lat3.toDeg(), lon3.toDeg());
}


/**
 * Returns the destination point from this point having travelled the given distance (in km) on the
 * given initial bearing (bearing may vary before destination is reached)
 *
 * see http://williams.best.vwh.net/avform.htm#LL
 *
 * @param {Number} brng: Initial bearing in degrees
 * @param {Number} dist: Distance in km
 * @returns {LatLon} Destination point
 */
destinationPoint = function(brng, dist) {
 dist = typeof(dist)=='number' ? dist : typeof(dist)=='string' && dist.trim()!='' ? +dist : NaN;
 dist = dist/this._radius; // convert dist to angular distance in radians
 brng = brng.toRad(); //
 var lat1 = this._lat.toRad(), lon1 = this._lon.toRad();

 var lat2 = Math.asin( Math.sin(lat1)*Math.cos(dist) +
 Math.cos(lat1)*Math.sin(dist)*Math.cos(brng) );
 var lon2 = lon1 + Math.atan2(Math.sin(brng)*Math.sin(dist)*Math.cos(lat1),
 Math.cos(dist)-Math.sin(lat1)*Math.sin(lat2));
 lon2 = (lon2+3*Math.PI)%(2*Math.PI) - Math.PI; // normalise to -180...+180

 return new LatLon(lat2.toDeg(), lon2.toDeg());
}


/**
 * Returns the point of intersection of two paths defined by point and bearing
 *
 * see http://williams.best.vwh.net/avform.htm#Intersection
 *
 * @param {LatLon} p1: First point
 * @param {Number} brng1: Initial bearing from first point
 * @param {LatLon} p2: Second point
 * @param {Number} brng2: Initial bearing from second point
 * @returns {LatLon} Destination point (null if no unique intersection defined)
 */
intersection = function(p1, brng1, p2, brng2) {
 brng1 = typeof brng1 == 'number' ? brng1 : typeof brng1 == 'string' && trim(brng1)!='' ? +brng1 : NaN;
 brng2 = typeof brng2 == 'number' ? brng2 : typeof brng2 == 'string' && trim(brng2)!='' ? +brng2 : NaN;
 lat1 = p1._lat.toRad(), lon1 = p1._lon.toRad();
 lat2 = p2._lat.toRad(), lon2 = p2._lon.toRad();
 brng13 = brng1.toRad(), brng23 = brng2.toRad();
 dLat = lat2-lat1, dLon = lon2-lon1;

 dist12 = 2*Math.asin( Math.sqrt( Math.sin(dLat/2)*Math.sin(dLat/2) +
 Math.cos(lat1)*Math.cos(lat2)*Math.sin(dLon/2)*Math.sin(dLon/2) ) );
 if (dist12 == 0) return null;

 // initial/final bearings between points
 brngA = Math.acos( ( Math.sin(lat2) - Math.sin(lat1)*Math.cos(dist12) ) /
 ( Math.sin(dist12)*Math.cos(lat1) ) );
 if (isNaN(brngA)) brngA = 0; // protect against rounding
 brngB = Math.acos( ( Math.sin(lat1) - Math.sin(lat2)*Math.cos(dist12) ) /
 ( Math.sin(dist12)*Math.cos(lat2) ) );

 if (Math.sin(lon2-lon1) > 0) {
 brng12 = brngA;
 brng21 = 2*Math.PI - brngB;
 } else {
 brng12 = 2*Math.PI - brngA;
 brng21 = brngB;
 }

 alpha1 = (brng13 - brng12 + Math.PI) % (2*Math.PI) - Math.PI; // angle 2-1-3
 alpha2 = (brng21 - brng23 + Math.PI) % (2*Math.PI) - Math.PI; // angle 1-2-3

 if (Math.sin(alpha1)==0 && Math.sin(alpha2)==0) return null; // infinite intersections
 if (Math.sin(alpha1)*Math.sin(alpha2) < 0) return null; // ambiguous intersection

 //alpha1 = Math.abs(alpha1);
 //alpha2 = Math.abs(alpha2);
 // ... Ed Williams takes abs of alpha1/alpha2, but seems to break calculation?

 alpha3 = Math.acos( -Math.cos(alpha1)*Math.cos(alpha2) +
 Math.sin(alpha1)*Math.sin(alpha2)*Math.cos(dist12) );
 dist13 = Math.atan2( Math.sin(dist12)*Math.sin(alpha1)*Math.sin(alpha2),
 Math.cos(alpha2)+Math.cos(alpha1)*Math.cos(alpha3) )
 lat3 = Math.asin( Math.sin(lat1)*Math.cos(dist13) +
 Math.cos(lat1)*Math.sin(dist13)*Math.cos(brng13) );
 dLon13 = Math.atan2( Math.sin(brng13)*Math.sin(dist13)*Math.cos(lat1),
 Math.cos(dist13)-Math.sin(lat1)*Math.sin(lat3) );
 lon3 = lon1+dLon13;
 lon3 = (lon3+Math.PI) % (2*Math.PI) - Math.PI; // normalise to -180..180º

 return new LatLon(lat3.toDeg(), lon3.toDeg());
}


/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

/**
 * Returns the distance from this point to the supplied point, in km, travelling along a rhumb line
 *
 * see http://williams.best.vwh.net/avform.htm#Rhumb
 *
 * @param {LatLon} point: Latitude/longitude of destination point
 * @returns {Number} Distance in km between this point and destination point
 */
rhumbDistanceTo = function(point) {
 var R = this._radius;
 var lat1 = this._lat.toRad(), lat2 = point._lat.toRad();
 var dLat = (point._lat-this._lat).toRad();
 var dLon = Math.abs(point._lon-this._lon).toRad();

 var dPhi = Math.log(Math.tan(lat2/2+Math.PI/4)/Math.tan(lat1/2+Math.PI/4));
 var q = (!isNaN(dLat/dPhi)) ? dLat/dPhi : Math.cos(lat1); // E-W line gives dPhi=0
 // if dLon over 180° take shorter rhumb across 180° meridian:
 if (dLon > Math.PI) dLon = 2*Math.PI - dLon;
 var dist = Math.sqrt(dLat*dLat + q*q*dLon*dLon) * R;

 return dist.toPrecisionFixed(4); // 4 sig figs reflects typical 0.3% accuracy of spherical model
}

/**
 * Returns the bearing from this point to the supplied point along a rhumb line, in degrees
 *
 * @param {LatLon} point: Latitude/longitude of destination point
 * @returns {Number} Bearing in degrees from North
 */
LatLon.prototype.rhumbBearingTo = function(point) {
 var lat1 = this._lat.toRad(), lat2 = point._lat.toRad();
 var dLon = (point._lon-this._lon).toRad();

 var dPhi = Math.log(Math.tan(lat2/2+Math.PI/4)/Math.tan(lat1/2+Math.PI/4));
 if (Math.abs(dLon) > Math.PI) dLon = dLon>0 ? -(2*Math.PI-dLon) : (2*Math.PI+dLon);
 var brng = Math.atan2(dLon, dPhi);

 return (brng.toDeg()+360) % 360;
}

/**
 * Returns the destination point from this point having travelled the given distance (in km) on the
 * given bearing along a rhumb line
 *
 * @param {Number} brng: Bearing in degrees from North
 * @param {Number} dist: Distance in km
 * @returns {LatLon} Destination point
 */
rhumbDestinationPoint = function(brng, dist) {
 var R = this._radius;
 var d = parseFloat(dist)/R; // d = angular distance covered on earth's surface
 var lat1 = this._lat.toRad(), lon1 = this._lon.toRad();
 brng = brng.toRad();

 var lat2 = lat1 + d*Math.cos(brng);
 var dLat = lat2-lat1;
 var dPhi = Math.log(Math.tan(lat2/2+Math.PI/4)/Math.tan(lat1/2+Math.PI/4));
 var q = (!isNaN(dLat/dPhi)) ? dLat/dPhi : Math.cos(lat1); // E-W line gives dPhi=0
 var dLon = d*Math.sin(brng)/q;
 // check for some daft bugger going past the pole
 if (Math.abs(lat2) > Math.PI/2) lat2 = lat2>0 ? Math.PI-lat2 : -(Math.PI-lat2);
 lon2 = (lon1+dLon+3*Math.PI)%(2*Math.PI) - Math.PI;

 return new LatLon(lat2.toDeg(), lon2.toDeg());
}

/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */


/**
 * Returns the latitude of this point; signed numeric degrees if no format, otherwise format & dp
 * as per Geo.toLat()
 *
 * @param {String} [format]: Return value as 'd', 'dm', 'dms'
 * @param {Number} [dp=0|2|4]: No of decimal places to display
 * @returns {Number|String} Numeric degrees if no format specified, otherwise deg/min/sec
 *
 * @requires Geo
 */
lat = function(format, dp) {
 if (typeof format == 'undefined') return this._lat;

 return Geo.toLat(this._lat, format, dp);
}

/**
 * Returns the longitude of this point; signed numeric degrees if no format, otherwise format & dp
 * as per Geo.toLon()
 *
 * @param {String} [format]: Return value as 'd', 'dm', 'dms'
 * @param {Number} [dp=0|2|4]: No of decimal places to display
 * @returns {Number|String} Numeric degrees if no format specified, otherwise deg/min/sec
 *
 * @requires Geo
 */
lon = function(format, dp) {
 if (typeof format == 'undefined') return this._lon;

 return Geo.toLon(this._lon, format, dp);
}

/**
 * Returns a string representation of this point; format and dp as per lat()/lon()
 *
 * @param {String} [format]: Return value as 'd', 'dm', 'dms'
 * @param {Number} [dp=0|2|4]: No of decimal places to display
 * @returns {String} Comma-separated latitude/longitude
 *
 * @requires Geo
 */
toString = function(format, dp) {
 if (typeof format == 'undefined') format = 'dms';

 if (isNaN(this._lat) || isNaN(this._lon)) return '-,-';

 return Geo.toLat(this._lat, format, dp) + ', ' + Geo.toLon(this._lon, format, dp);
}

/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

// ---- extend Number object with methods for converting degrees/radians

/** Converts numeric degrees to radians */
if (typeof(Number.prototype.toRad) === "undefined") {
 Number.prototype.toRad = function() {
 return this * Math.PI / 180;
 }
}

/** Converts radians to numeric (signed) degrees */
if (typeof(Number.prototype.toDeg) === "undefined") {
 Number.prototype.toDeg = function() {
 return this * 180 / Math.PI;
 }
}

/**
 * Formats the significant digits of a number, using only fixed-point notation (no exponential)
 *
 * @param {Number} precision: Number of significant digits to appear in the returned string
 * @returns {String} A string representation of number which contains precision significant digits
 */
if (typeof(Number.prototype.toPrecisionFixed) === "undefined") {
 Number.prototype.toPrecisionFixed = function(precision) {
 if (isNaN(this)) return 'NaN';
 var numb = this < 0 ? -this : this; // can't take log of -ve number...
 var sign = this < 0 ? '-' : '';

 if (numb == 0) { n = '0.'; while (precision--) n += '0'; return n }; // can't take log of zero

 var scale = Math.ceil(Math.log(numb)*Math.LOG10E); // no of digits before decimal
 var n = String(Math.round(numb * Math.pow(10, precision-scale)));
 if (scale > 0) { // add trailing zeros & insert decimal as required
 l = scale - n.length;
 while (l-- > 0) n = n + '0';
 if (scale < n.length) n = n.slice(0,scale) + '.' + n.slice(scale);
 } else { // prefix decimal and leading zeros if required
 while (scale++ < 0) n = '0' + n;
 n = '0.' + n;
 }
 return sign + n;
 }
}

/** Trims whitespace from string (q.v. blog.stevenlevithan.com/archives/faster-trim-javascript) */
if (typeof(String.prototype.trim) === "undefined") {
 String.prototype.trim = function() {
 return String(this).replace(/^\s\s*/, '').replace(/\s\s*$/, '');
 }
}


/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/* Geodesy representation conversion functions (c) Chris Veness 2002-2010 */
/* - www.movable-type.co.uk/scripts/latlong.html */
/* */
/* Sample usage: */
/* var lat = Geo.parseDMS('51° 28′ 40.12″ N'); */
/* var lon = Geo.parseDMS('000° 00′ 05.31″ W'); */
/* var p1 = new LatLon(lat, lon); */
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

var Geo = {}; // Geo namespace, representing static class

/**
 * Parses string representing degrees/minutes/seconds into numeric degrees
 *
 * This is very flexible on formats, allowing signed decimal degrees, or deg-min-sec optionally
 * suffixed by compass direction (NSEW). A variety of separators are accepted (eg 3º 37' 09"W)
 * or fixed-width format without separators (eg 0033709W). Seconds and minutes may be omitted.
 * (Note minimal validation is done).
 *
 * @param {String|Number} dmsStr: Degrees or deg/min/sec in variety of formats
 * @returns {Number} Degrees as decimal number
 * @throws {TypeError} dmsStr is an object, perhaps DOM object without .value?
 */
parseDMS = function(dmsStr) {
 if (typeof deg == 'object') throw new TypeError('Geo.parseDMS - dmsStr is [DOM?] object');

 // check for signed decimal degrees without NSEW, if so return it directly
 if (typeof dmsStr === 'number' && isFinite(dmsStr)) return Number(dmsStr);

 // strip off any sign or compass dir'n & split out separate d/m/s
 var dms = String(dmsStr).trim().replace(/^-/,'').replace(/[NSEW]$/i,'').split(/[^0-9.,]+/);
 if (dms[dms.length-1]=='') dms.splice(dms.length-1); // from trailing symbol

 if (dms == '') return NaN;

 // and convert to decimal degrees...
 switch (dms.length) {
 case 3: // interpret 3-part result as d/m/s
 var deg = dms[0]/1 + dms[1]/60 + dms[2]/3600;
 break;
 case 2: // interpret 2-part result as d/m
 var deg = dms[0]/1 + dms[1]/60;
 break;
 case 1: // just d (possibly decimal) or non-separated dddmmss
 var deg = dms[0];
 // check for fixed-width unseparated format eg 0033709W
 if (/[NS]/i.test(dmsStr)) deg = '0' + deg; // - normalise N/S to 3-digit degrees
 if (/[0-9]{7}/.test(deg)) deg = deg.slice(0,3)/1 + deg.slice(3,5)/60 + deg.slice(5)/3600;
 break;
 default:
 return NaN;
 }
 if (/^-|[WS]$/i.test(dmsStr.trim())) deg = -deg; // take '-', west and south as -ve
 return Number(deg);
}

/**
 * Convert decimal degrees to deg/min/sec format
 * - degree, prime, double-prime symbols are added, but sign is discarded, though no compass
 * direction is added
 *
 * @private
 * @param {Number} deg: Degrees
 * @param {String} [format=dms]: Return value as 'd', 'dm', 'dms'
 * @param {Number} [dp=0|2|4]: No of decimal places to use - default 0 for dms, 2 for dm, 4 for d
 * @returns {String} deg formatted as deg/min/secs according to specified format
 * @throws {TypeError} deg is an object, perhaps DOM object without .value?
 */
toDMS = function(deg, format, dp) {
 if (typeof deg == 'object') throw new TypeError('Geo.toDMS - deg is [DOM?] object');
 if (isNaN(deg)) return 'NaN'; // give up here if we can't make a number from deg

 // default values
 if (typeof format == 'undefined') format = 'dms';
 if (typeof dp == 'undefined') {
 switch (format) {
 case 'd': dp = 4; break;
 case 'dm': dp = 2; break;
 case 'dms': dp = 0; break;
 default: format = 'dms'; dp = 0; // be forgiving on invalid format
 }
 }

 deg = Math.abs(deg); // (unsigned result ready for appending compass dir'n)

 switch (format) {
 case 'd':
 d = deg.toFixed(dp); // round degrees
 if (d<100) d = '0' + d; // pad with leading zeros
 if (d<10) d = '0' + d;
 dms = d + '\u00B0'; // add º symbol
 break;
 case 'dm':
 var min = (deg*60).toFixed(dp); // convert degrees to minutes & round
 var d = Math.floor(min / 60); // get component deg/min
 var m = (min % 60).toFixed(dp); // pad with trailing zeros
 if (d<100) d = '0' + d; // pad with leading zeros
 if (d<10) d = '0' + d;
 if (m<10) m = '0' + m;
 dms = d + '\u00B0' + m + '\u2032'; // add º, ' symbols
 break;
 case 'dms':
 var sec = (deg*3600).toFixed(dp); // convert degrees to seconds & round
 var d = Math.floor(sec / 3600); // get component deg/min/sec
 var m = Math.floor(sec/60) % 60;
 var s = (sec % 60).toFixed(dp); // pad with trailing zeros
 if (d<100) d = '0' + d; // pad with leading zeros
 if (d<10) d = '0' + d;
 if (m<10) m = '0' + m;
 if (s<10) s = '0' + s;
 dms = d + '\u00B0' + m + '\u2032' + s + '\u2033'; // add º, ', " symbols
 break;
 }

 return dms;
}

/**
 * Convert numeric degrees to deg/min/sec latitude (suffixed with N/S)
 *
 * @param {Number} deg: Degrees
 * @param {String} [format=dms]: Return value as 'd', 'dm', 'dms'
 * @param {Number} [dp=0|2|4]: No of decimal places to use - default 0 for dms, 2 for dm, 4 for d
 * @returns {String} Deg/min/seconds
 */
toLat = function(deg, format, dp) {
 var lat = Geo.toDMS(deg, format, dp);
 return lat=='' ? '' : lat.slice(1) + (deg<0 ? 'S' : 'N'); // knock off initial '0' for lat!
}

/**
 * Convert numeric degrees to deg/min/sec longitude (suffixed with E/W)
 *
 * @param {Number} deg: Degrees
 * @param {String} [format=dms]: Return value as 'd', 'dm', 'dms'
 * @param {Number} [dp=0|2|4]: No of decimal places to use - default 0 for dms, 2 for dm, 4 for d
 * @returns {String} Deg/min/seconds
 */
toLon = function(deg, format, dp) {
 var lon = Geo.toDMS(deg, format, dp);
 return lon=='' ? '' : lon + (deg<0 ? 'W' : 'E');
}

/**
 * Convert numeric degrees to deg/min/sec as a bearing (0º..360º)
 *
 * @param {Number} deg: Degrees
 * @param {String} [format=dms]: Return value as 'd', 'dm', 'dms'
 * @param {Number} [dp=0|2|4]: No of decimal places to use - default 0 for dms, 2 for dm, 4 for d
 * @returns {String} Deg/min/seconds
 */
toBrng = function(deg, format, dp) {
 deg = (Number(deg)+360) % 360; // normalise -ve values to 180º..360º
 var brng = Geo.toDMS(deg, format, dp);
 return brng.replace('360', '0'); // just in case rounding took us up to 360º!
}

Load some example data:


GEO:
LOAD * INLINE [
 ROUTE, LAT_1, LONG_1, LAT_2, LONG_2
 R1, 53 09 02N, 001 50 40W, 52 12 19N, 000 08 33W
];


DIST:
LOAD
 *,
 distVincenty(parseDMS(LAT_1), parseDMS(LONG_1), parseDMS(LAT_2), parseDMS(LONG_2)) as DIST
resident GEO;

drop table GEO;

-Alex

WGS84 Compliant Great Circle Distance Calculation

Thu, 10 Mar 2011 11:20:13 GMT

Hi Alex,

Many thanks again for the script and instructions, I had a small (well ... quite large) panic yesterday when I compared the results I got using this to the results I had using the Haversine (spherical model) formula. Luckily, after a bit of investigation I found that all it was, was that the Lat and Long format in my file has the NSEW indicator at the front not the end so although the formula will still run (and give valid results on N and E coordinates which was confusing) it didn't seem to spot the indicator and do anything with it for the S and W directions.

I'm just really putting this on here for completeness in case anyone else uses the script ...

I have found someone here with more knowledge of Javascipt than me (not difficult!!) and there is one line of code we've amended to give a little more flexibility and allow for either a leading or trailing compass letter:

if (/^-|[WS]$/i.test(dmsStr.trim())) deg = -deg; // take '-', west and south as -ve

we have changed to

if (/^-|[WS]|[WS]$/i.test(dmsStr.trim())) deg = -deg; // take '-', west and south as -ve

A small change but means that I can plug it straight into our data without needing to convert to degrees first or try and move the compass letter from the beginning to the end of the coordinates!

Many thanks once again.

Cheers,

Emma

topic WGS84 Compliant Great Circle Distance Calculation in QlikView

WGS84 Compliant Great Circle Distance Calculation

WGS84 Compliant Great Circle Distance Calculation

WGS84 Compliant Great Circle Distance Calculation

WGS84 Compliant Great Circle Distance Calculation

WGS84 Compliant Great Circle Distance Calculation

WGS84 Compliant Great Circle Distance Calculation