https://www.apostaseis.gr/loc_ap/scripts/arc.js

1'use strict';
2
3var D2R = Math.PI / 180;
4var R2D = 180 / Math.PI;
5
6var Coord = function(lon,lat) {
7    this.lon = lon;
8    this.lat = lat;
9    this.x = D2R * lon;
10    this.y = D2R * lat;
11};
12
13Coord.prototype.view = function() {
14    return String(this.lon).slice(0, 4) + ',' + String(this.lat).slice(0, 4);
15};
16
17Coord.prototype.antipode = function() {
18    var anti_lat = -1 * this.lat;
19    var anti_lon = (this.lon < 0) ? 180 + this.lon : (180 - this.lon) * -1;
20    return new Coord(anti_lon, anti_lat);
21};
22
23var LineString = function() {
24    this.coords = [];
25    this.length = 0;
26};
27
28LineString.prototype.move_to = function(coord) {
29    this.length++;
30    this.coords.push(coord);
31};
32
33var Arc = function(properties) {
34    this.properties = properties || {};
35    this.geometries = [];
36};
37
38Arc.prototype.json = function() {
39    if (this.geometries.length <= 0) {
40        return {'geometry': { 'type': 'LineString', 'coordinates': null },
41                'type': 'Feature', 'properties': this.properties
42               };
43    } else if (this.geometries.length == 1) {
44        return {'geometry': { 'type': 'LineString', 'coordinates': this.geometries[0].coords },
45                'type': 'Feature', 'properties': this.properties
46               };
47    } else {
48        var multiline = [];
49        for (var i = 0; i < this.geometries.length; i++) {
50            multiline.push(this.geometries[i].coords);
51        }
52        return {'geometry': { 'type': 'MultiLineString', 'coordinates': multiline },
53                'type': 'Feature', 'properties': this.properties
54               };
55    }
56};
57
58// TODO - output proper multilinestring
59Arc.prototype.wkt = function() {
60    var wkt_string = '';
61    var wkt = 'LINESTRING(';
62    var collect = function(c) { wkt += c[0] + ' ' + c[1] + ','; };
63    for (var i = 0; i < this.geometries.length; i++) {
64        if (this.geometries[i].coords.length === 0) {
65            return 'LINESTRING(empty)';
66        } else {
67            var coords = this.geometries[i].coords;
68            coords.forEach(collect);
69            wkt_string += wkt.substring(0, wkt.length - 1) + ')';
70        }
71    }
72    return wkt_string;
73};
74
75/*
76 * http://en.wikipedia.org/wiki/Great-circle_distance
77 *
78 */
79var GreatCircle = function(start,end,properties) {
80    if (!start || start.x === undefined || start.y === undefined) {
81        throw new Error("GreatCircle constructor expects two args: start and end objects with x and y properties");
82    }
83    if (!end || end.x === undefined || end.y === undefined) {
84        throw new Error("GreatCircle constructor expects two args: start and end objects with x and y properties");
85    }
86    this.start = new Coord(start.x,start.y);
87    this.end = new Coord(end.x,end.y);
88    this.properties = properties || {};
89
90    var w = this.start.x - this.end.x;
91    var h = this.start.y - this.end.y;
92    var z = Math.pow(Math.sin(h / 2.0), 2) +
93                Math.cos(this.start.y) *
94                   Math.cos(this.end.y) *
95                     Math.pow(Math.sin(w / 2.0), 2);
96    this.g = 2.0 * Math.asin(Math.sqrt(z));
97
98    if (this.g == Math.PI) {
99        throw new Error('it appears ' + start.view() + ' and ' + end.view() + " are 'antipodal', e.g diametrically opposite, thus there is no single route but rather infinite");
100    } else if (isNaN(this.g)) {
101        throw new Error('could not calculate great circle between ' + start + ' and ' + end);
102    }
103};
104
105/*
106 * http://williams.best.vwh.net/avform.htm#Intermediate
107 */
108GreatCircle.prototype.interpolate = function(f) {
109    var A = Math.sin((1 - f) * this.g) / Math.sin(this.g);
110    var B = Math.sin(f * this.g) / Math.sin(this.g);
111    var x = A * Math.cos(this.start.y) * Math.cos(this.start.x) + B * Math.cos(this.end.y) * Math.cos(this.end.x);
112    var y = A * Math.cos(this.start.y) * Math.sin(this.start.x) + B * Math.cos(this.end.y) * Math.sin(this.end.x);
113    var z = A * Math.sin(this.start.y) + B * Math.sin(this.end.y);
114    var lat = R2D * Math.atan2(z, Math.sqrt(Math.pow(x, 2) + Math.pow(y, 2)));
115    var lon = R2D * Math.atan2(y, x);
116    return [lon, lat];
117};
118
119
120
121/*
122 * Generate points along the great circle
123 */
124GreatCircle.prototype.Arc = function(npoints,options) {
125    var first_pass = [];
126    if (!npoints || npoints <= 2) {
127        first_pass.push([this.start.lon, this.start.lat]);
128        first_pass.push([this.end.lon, this.end.lat]);
129    } else {
130        var delta = 1.0 / (npoints - 1);
131        for (var i = 0; i < npoints; ++i) {
132            var step = delta * i;
133            var pair = this.interpolate(step);
134            first_pass.push(pair);
135        }
136    }
137    /* partial port of dateline handling from:
138      gdal/ogr/ogrgeometryfactory.cpp
139
140      TODO - does not handle all wrapping scenarios yet
141    */
142    var bHasBigDiff = false;

143    var dfMaxSmallDiffLong = 0;
144    // from http://www.gdal.org/ogr2ogr.html
145    // -datelineoffset:
146    // (starting with GDAL 1.10) offset from dateline in degrees (default long. = +/- 10deg, geometries within 170deg to -170deg will be splited)
147    var dfDateLineOffset = options && options.offset ? options.offset : 10;
148    var dfLeftBorderX = 180 - dfDateLineOffset;
149    var dfRightBorderX = -180 + dfDateLineOffset;
150    var dfDiffSpace = 360 - dfDateLineOffset;
151
152    // https://github.com/OSGeo/gdal/blob/7bfb9c452a59aac958bff0c8386b891edf8154ca/gdal/ogr/ogrgeometryfactory.cpp#L2342
153    for (var j = 1; j < first_pass.length; ++j) {
154        var dfPrevX = first_pass[j-1][0];
155        var dfX = first_pass[j][0];
156        var dfDiffLong = Math.abs(dfX - dfPrevX);
157        if (dfDiffLong > dfDiffSpace &&
158            ((dfX > dfLeftBorderX && dfPrevX < dfRightBorderX) || (dfPrevX > dfLeftBorderX && dfX < dfRightBorderX))) {
159            bHasBigDiff = true;
160        } else if (dfDiffLong > dfMaxSmallDiffLong) {
161            dfMaxSmallDiffLong = dfDiffLong;
162        }
163    }
164
165    var poMulti = [];
166    if (bHasBigDiff && dfMaxSmallDiffLong < dfDateLineOffset) {
167        var poNewLS = [];
168        poMulti.push(poNewLS);
169        for (var k = 0; k < first_pass.length; ++k) {
170            var dfX0 = parseFloat(first_pass[k][0]);
171            if (k > 0 &&  Math.abs(dfX0 - first_pass[k-1][0]) > dfDiffSpace) {
172                var dfX1 = parseFloat(first_pass[k-1][0]);
173                var dfY1 = parseFloat(first_pass[k-1][1]);
174                var dfX2 = parseFloat(first_pass[k][0]);
175                var dfY2 = parseFloat(first_pass[k][1]);
176                if (dfX1 > -180 && dfX1 < dfRightBorderX && dfX2 == 180 &&
177                    k+1 < first_pass.length &&
178                   first_pass[k-1][0] > -180 && first_pass[k-1][0] < dfRightBorderX)
179                {
180                     poNewLS.push([-180, first_pass[k][1]]);
181                     k++;
182                     poNewLS.push([first_pass[k][0], first_pass[k][1]]);
183                     continue;
184                } else if (dfX1 > dfLeftBorderX && dfX1 < 180 && dfX2 == -180 &&
185                     k+1 < first_pass.length &&
186                     first_pass[k-1][0] > dfLeftBorderX && first_pass[k-1][0] < 180)
187                {
188                     poNewLS.push([180, first_pass[k][1]]);
189                     k++;
190                     poNewLS.push([first_pass[k][0], first_pass[k][1]]);
191                     continue;
192                }
193
194                if (dfX1 < dfRightBorderX && dfX2 > dfLeftBorderX)
195                {
196                    // swap dfX1, dfX2
197                    var tmpX = dfX1;
198                    dfX1 = dfX2;
199                    dfX2 = tmpX;
200                    // swap dfY1, dfY2
201                    var tmpY = dfY1;
202                    dfY1 = dfY2;
203                    dfY2 = tmpY;
204                }
205                if (dfX1 > dfLeftBorderX && dfX2 < dfRightBorderX) {
206                    dfX2 += 360;
207                }
208
209                if (dfX1 <= 180 && dfX2 >= 180 && dfX1 < dfX2)
210                {
211                    var dfRatio = (180 - dfX1) / (dfX2 - dfX1);
212                    var dfY = dfRatio * dfY2 + (1 - dfRatio) * dfY1;
213                    poNewLS.push([first_pass[k-1][0] > dfLeftBorderX ? 180 : -180, dfY]);
214                    poNewLS = [];
215                    poNewLS.push([first_pass[k-1][0] > dfLeftBorderX ? -180 : 180, dfY]);
216                    poMulti.push(poNewLS);
217                }
218                else
219                {
220                    poNewLS = [];
221                    poMulti.push(poNewLS);
222                }
223                poNewLS.push([dfX0, first_pass[k][1]]);
224            } else {
225                poNewLS.push([first_pass[k][0], first_pass[k][1]]);
226            }
227        }
228    } else {
229        // add normally
230        var poNewLS0 = [];
231        poMulti.push(poNewLS0);
232        for (var l = 0; l < first_pass.length; ++l) {
233            poNewLS0.push([first_pass[l][0],first_pass[l][1]]);
234        }
235    }
236
237    var arc = new Arc(this.properties);
238    for (var m = 0; m < poMulti.length; ++m) {
239        var line = new LineString();
240        arc.geometries.push(line);
241        var points = poMulti[m];
242        for (var j0 = 0; j0 < points.length; ++j0) {
243            line.move_to(points[j0]);
244        }
245    }
246    return arc;
247};
248
249if (typeof module !== 'undefined' && typeof module.exports !== 'undefined') {
250  // nodejs
251  module.exports.Coord = Coord;
252  module.exports.Arc = Arc;
253  module.exports.GreatCircle = GreatCircle;
254
255} else {
256  // browser
257  var arc = {};
258  arc.Coord = Coord;
259  arc.Arc = Arc;
260  arc.GreatCircle = GreatCircle;
261}