index.js 9.9 KB

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