index.js 9.9 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336
  1. 'use strict';
  2. var invariant = require('@turf/invariant');
  3. var D2R = Math.PI / 180;
  4. var R2D = 180 / Math.PI;
  5. var Coord = function (lon, lat) {
  6. this.lon = lon;
  7. this.lat = lat;
  8. this.x = D2R * lon;
  9. this.y = D2R * lat;
  10. };
  11. Coord.prototype.view = function () {
  12. return String(this.lon).slice(0, 4) + "," + String(this.lat).slice(0, 4);
  13. };
  14. Coord.prototype.antipode = function () {
  15. var anti_lat = -1 * this.lat;
  16. var anti_lon = this.lon < 0 ? 180 + this.lon : (180 - this.lon) * -1;
  17. return new Coord(anti_lon, anti_lat);
  18. };
  19. var LineString = function () {
  20. this.coords = [];
  21. this.length = 0;
  22. };
  23. LineString.prototype.move_to = function (coord) {
  24. this.length++;
  25. this.coords.push(coord);
  26. };
  27. var Arc = function (properties) {
  28. this.properties = properties || {};
  29. this.geometries = [];
  30. };
  31. Arc.prototype.json = function () {
  32. if (this.geometries.length <= 0) {
  33. return {
  34. geometry: { type: "LineString", coordinates: null },
  35. type: "Feature",
  36. properties: this.properties,
  37. };
  38. } else if (this.geometries.length === 1) {
  39. return {
  40. geometry: { type: "LineString", coordinates: this.geometries[0].coords },
  41. type: "Feature",
  42. properties: this.properties,
  43. };
  44. } else {
  45. var multiline = [];
  46. for (var i = 0; i < this.geometries.length; i++) {
  47. multiline.push(this.geometries[i].coords);
  48. }
  49. return {
  50. geometry: { type: "MultiLineString", coordinates: multiline },
  51. type: "Feature",
  52. properties: this.properties,
  53. };
  54. }
  55. };
  56. // TODO - output proper multilinestring
  57. Arc.prototype.wkt = function () {
  58. var wkt_string = "";
  59. var wkt = "LINESTRING(";
  60. var collect = function (c) {
  61. wkt += c[0] + " " + c[1] + ",";
  62. };
  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. * http://en.wikipedia.org/wiki/Great-circle_distance
  76. *
  77. */
  78. var GreatCircle = function (start, end, properties) {
  79. if (!start || start.x === undefined || start.y === undefined) {
  80. throw new Error(
  81. "GreatCircle constructor expects two args: start and end objects with x and y properties"
  82. );
  83. }
  84. if (!end || end.x === undefined || end.y === undefined) {
  85. throw new Error(
  86. "GreatCircle constructor expects two args: start and end objects with x and y properties"
  87. );
  88. }
  89. this.start = new Coord(start.x, start.y);
  90. this.end = new Coord(end.x, end.y);
  91. this.properties = properties || {};
  92. var w = this.start.x - this.end.x;
  93. var h = this.start.y - this.end.y;
  94. var z =
  95. Math.pow(Math.sin(h / 2.0), 2) +
  96. Math.cos(this.start.y) *
  97. Math.cos(this.end.y) *
  98. Math.pow(Math.sin(w / 2.0), 2);
  99. this.g = 2.0 * Math.asin(Math.sqrt(z));
  100. if (this.g === Math.PI) {
  101. throw new Error(
  102. "it appears " +
  103. start.view() +
  104. " and " +
  105. end.view() +
  106. " are 'antipodal', e.g diametrically opposite, thus there is no single route but rather infinite"
  107. );
  108. } else if (isNaN(this.g)) {
  109. throw new Error(
  110. "could not calculate great circle between " + start + " and " + end
  111. );
  112. }
  113. };
  114. /*
  115. * http://williams.best.vwh.net/avform.htm#Intermediate
  116. */
  117. GreatCircle.prototype.interpolate = function (f) {
  118. var A = Math.sin((1 - f) * this.g) / Math.sin(this.g);
  119. var B = Math.sin(f * this.g) / Math.sin(this.g);
  120. var x =
  121. A * Math.cos(this.start.y) * Math.cos(this.start.x) +
  122. B * Math.cos(this.end.y) * Math.cos(this.end.x);
  123. var y =
  124. A * Math.cos(this.start.y) * Math.sin(this.start.x) +
  125. B * Math.cos(this.end.y) * Math.sin(this.end.x);
  126. var z = A * Math.sin(this.start.y) + B * Math.sin(this.end.y);
  127. var lat = R2D * Math.atan2(z, Math.sqrt(Math.pow(x, 2) + Math.pow(y, 2)));
  128. var lon = R2D * Math.atan2(y, x);
  129. return [lon, lat];
  130. };
  131. /*
  132. * Generate points along the great circle
  133. */
  134. GreatCircle.prototype.Arc = function (npoints, options) {
  135. var first_pass = [];
  136. if (!npoints || npoints <= 2) {
  137. first_pass.push([this.start.lon, this.start.lat]);
  138. first_pass.push([this.end.lon, this.end.lat]);
  139. } else {
  140. var delta = 1.0 / (npoints - 1);
  141. for (var i = 0; i < npoints; ++i) {
  142. var step = delta * i;
  143. var pair = this.interpolate(step);
  144. first_pass.push(pair);
  145. }
  146. }
  147. /* partial port of dateline handling from:
  148. gdal/ogr/ogrgeometryfactory.cpp
  149. TODO - does not handle all wrapping scenarios yet
  150. */
  151. var bHasBigDiff = false;
  152. var dfMaxSmallDiffLong = 0;
  153. // from http://www.gdal.org/ogr2ogr.html
  154. // -datelineoffset:
  155. // (starting with GDAL 1.10) offset from dateline in degrees (default long. = +/- 10deg, geometries within 170deg to -170deg will be splited)
  156. var dfDateLineOffset = options && options.offset ? options.offset : 10;
  157. var dfLeftBorderX = 180 - dfDateLineOffset;
  158. var dfRightBorderX = -180 + dfDateLineOffset;
  159. var dfDiffSpace = 360 - dfDateLineOffset;
  160. // https://github.com/OSGeo/gdal/blob/7bfb9c452a59aac958bff0c8386b891edf8154ca/gdal/ogr/ogrgeometryfactory.cpp#L2342
  161. for (var j = 1; j < first_pass.length; ++j) {
  162. var dfPrevX = first_pass[j - 1][0];
  163. var dfX = first_pass[j][0];
  164. var dfDiffLong = Math.abs(dfX - dfPrevX);
  165. if (
  166. dfDiffLong > dfDiffSpace &&
  167. ((dfX > dfLeftBorderX && dfPrevX < dfRightBorderX) ||
  168. (dfPrevX > dfLeftBorderX && dfX < dfRightBorderX))
  169. ) {
  170. bHasBigDiff = true;
  171. } else if (dfDiffLong > dfMaxSmallDiffLong) {
  172. dfMaxSmallDiffLong = dfDiffLong;
  173. }
  174. }
  175. var poMulti = [];
  176. if (bHasBigDiff && dfMaxSmallDiffLong < dfDateLineOffset) {
  177. var poNewLS = [];
  178. poMulti.push(poNewLS);
  179. for (var k = 0; k < first_pass.length; ++k) {
  180. var dfX0 = parseFloat(first_pass[k][0]);
  181. if (k > 0 && Math.abs(dfX0 - first_pass[k - 1][0]) > dfDiffSpace) {
  182. var dfX1 = parseFloat(first_pass[k - 1][0]);
  183. var dfY1 = parseFloat(first_pass[k - 1][1]);
  184. var dfX2 = parseFloat(first_pass[k][0]);
  185. var dfY2 = parseFloat(first_pass[k][1]);
  186. if (
  187. dfX1 > -180 &&
  188. dfX1 < dfRightBorderX &&
  189. dfX2 === 180 &&
  190. k + 1 < first_pass.length &&
  191. first_pass[k - 1][0] > -180 &&
  192. first_pass[k - 1][0] < dfRightBorderX
  193. ) {
  194. poNewLS.push([-180, first_pass[k][1]]);
  195. k++;
  196. poNewLS.push([first_pass[k][0], first_pass[k][1]]);
  197. continue;
  198. } else if (
  199. dfX1 > dfLeftBorderX &&
  200. dfX1 < 180 &&
  201. dfX2 === -180 &&
  202. k + 1 < first_pass.length &&
  203. first_pass[k - 1][0] > dfLeftBorderX &&
  204. first_pass[k - 1][0] < 180
  205. ) {
  206. poNewLS.push([180, first_pass[k][1]]);
  207. k++;
  208. poNewLS.push([first_pass[k][0], first_pass[k][1]]);
  209. continue;
  210. }
  211. if (dfX1 < dfRightBorderX && dfX2 > dfLeftBorderX) {
  212. // swap dfX1, dfX2
  213. var tmpX = dfX1;
  214. dfX1 = dfX2;
  215. dfX2 = tmpX;
  216. // swap dfY1, dfY2
  217. var tmpY = dfY1;
  218. dfY1 = dfY2;
  219. dfY2 = tmpY;
  220. }
  221. if (dfX1 > dfLeftBorderX && dfX2 < dfRightBorderX) {
  222. dfX2 += 360;
  223. }
  224. if (dfX1 <= 180 && dfX2 >= 180 && dfX1 < dfX2) {
  225. var dfRatio = (180 - dfX1) / (dfX2 - dfX1);
  226. var dfY = dfRatio * dfY2 + (1 - dfRatio) * dfY1;
  227. poNewLS.push([
  228. first_pass[k - 1][0] > dfLeftBorderX ? 180 : -180,
  229. dfY,
  230. ]);
  231. poNewLS = [];
  232. poNewLS.push([
  233. first_pass[k - 1][0] > dfLeftBorderX ? -180 : 180,
  234. dfY,
  235. ]);
  236. poMulti.push(poNewLS);
  237. } else {
  238. poNewLS = [];
  239. poMulti.push(poNewLS);
  240. }
  241. poNewLS.push([dfX0, first_pass[k][1]]);
  242. } else {
  243. poNewLS.push([first_pass[k][0], first_pass[k][1]]);
  244. }
  245. }
  246. } else {
  247. // add normally
  248. var poNewLS0 = [];
  249. poMulti.push(poNewLS0);
  250. for (var l = 0; l < first_pass.length; ++l) {
  251. poNewLS0.push([first_pass[l][0], first_pass[l][1]]);
  252. }
  253. }
  254. var arc = new Arc(this.properties);
  255. for (var m = 0; m < poMulti.length; ++m) {
  256. var line = new LineString();
  257. arc.geometries.push(line);
  258. var points = poMulti[m];
  259. for (var j0 = 0; j0 < points.length; ++j0) {
  260. line.move_to(points[j0]);
  261. }
  262. }
  263. return arc;
  264. };
  265. /**
  266. * Calculate great circles routes as {@link LineString} or {@link MultiLineString}.
  267. * If the `start` and `end` points span the antimeridian, the resulting feature will
  268. * be split into a `MultiLineString`.
  269. *
  270. * @name greatCircle
  271. * @param {Coord} start source point feature
  272. * @param {Coord} end destination point feature
  273. * @param {Object} [options={}] Optional parameters
  274. * @param {Object} [options.properties={}] line feature properties
  275. * @param {number} [options.npoints=100] number of points
  276. * @param {number} [options.offset=10] offset controls the likelyhood that lines will
  277. * be split which cross the dateline. The higher the number the more likely.
  278. * @returns {Feature<LineString | MultiLineString>} great circle line feature
  279. * @example
  280. * var start = turf.point([-122, 48]);
  281. * var end = turf.point([-77, 39]);
  282. *
  283. * var greatCircle = turf.greatCircle(start, end, {properties: {name: 'Seattle to DC'}});
  284. *
  285. * //addToMap
  286. * var addToMap = [start, end, greatCircle]
  287. */
  288. function greatCircle(start, end, options) {
  289. // Optional parameters
  290. options = options || {};
  291. if (typeof options !== "object") throw new Error("options is invalid");
  292. var properties = options.properties;
  293. var npoints = options.npoints;
  294. var offset = options.offset;
  295. start = invariant.getCoord(start);
  296. end = invariant.getCoord(end);
  297. properties = properties || {};
  298. npoints = npoints || 100;
  299. offset = offset || 10;
  300. var generator = new GreatCircle(
  301. { x: start[0], y: start[1] },
  302. { x: end[0], y: end[1] },
  303. properties
  304. );
  305. var line = generator.Arc(npoints, { offset: offset });
  306. return line.json();
  307. }
  308. module.exports = greatCircle;
  309. module.exports.default = greatCircle;