omerc.js 6.3 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241
  1. import tsfnz from '../common/tsfnz';
  2. import adjust_lon from '../common/adjust_lon';
  3. import phi2z from '../common/phi2z';
  4. import { D2R, EPSLN, HALF_PI, TWO_PI, FORTPI } from '../constants/values';
  5. var TOL = 1e-7;
  6. function isTypeA(P) {
  7. var typeAProjections = ['Hotine_Oblique_Mercator','Hotine_Oblique_Mercator_Azimuth_Natural_Origin'];
  8. var projectionName = typeof P.PROJECTION === "object" ? Object.keys(P.PROJECTION)[0] : P.PROJECTION;
  9. return 'no_uoff' in P || 'no_off' in P || typeAProjections.indexOf(projectionName) !== -1;
  10. }
  11. /* Initialize the Oblique Mercator projection
  12. ------------------------------------------*/
  13. export function init() {
  14. var con, com, cosph0, D, F, H, L, sinph0, p, J, gamma = 0,
  15. gamma0, lamc = 0, lam1 = 0, lam2 = 0, phi1 = 0, phi2 = 0, alpha_c = 0, AB;
  16. // only Type A uses the no_off or no_uoff property
  17. // https://github.com/OSGeo/proj.4/issues/104
  18. this.no_off = isTypeA(this);
  19. this.no_rot = 'no_rot' in this;
  20. var alp = false;
  21. if ("alpha" in this) {
  22. alp = true;
  23. }
  24. var gam = false;
  25. if ("rectified_grid_angle" in this) {
  26. gam = true;
  27. }
  28. if (alp) {
  29. alpha_c = this.alpha;
  30. }
  31. if (gam) {
  32. gamma = (this.rectified_grid_angle * D2R);
  33. }
  34. if (alp || gam) {
  35. lamc = this.longc;
  36. } else {
  37. lam1 = this.long1;
  38. phi1 = this.lat1;
  39. lam2 = this.long2;
  40. phi2 = this.lat2;
  41. if (Math.abs(phi1 - phi2) <= TOL || (con = Math.abs(phi1)) <= TOL ||
  42. Math.abs(con - HALF_PI) <= TOL || Math.abs(Math.abs(this.lat0) - HALF_PI) <= TOL ||
  43. Math.abs(Math.abs(phi2) - HALF_PI) <= TOL) {
  44. throw new Error();
  45. }
  46. }
  47. var one_es = 1.0 - this.es;
  48. com = Math.sqrt(one_es);
  49. if (Math.abs(this.lat0) > EPSLN) {
  50. sinph0 = Math.sin(this.lat0);
  51. cosph0 = Math.cos(this.lat0);
  52. con = 1 - this.es * sinph0 * sinph0;
  53. this.B = cosph0 * cosph0;
  54. this.B = Math.sqrt(1 + this.es * this.B * this.B / one_es);
  55. this.A = this.B * this.k0 * com / con;
  56. D = this.B * com / (cosph0 * Math.sqrt(con));
  57. F = D * D -1;
  58. if (F <= 0) {
  59. F = 0;
  60. } else {
  61. F = Math.sqrt(F);
  62. if (this.lat0 < 0) {
  63. F = -F;
  64. }
  65. }
  66. this.E = F += D;
  67. this.E *= Math.pow(tsfnz(this.e, this.lat0, sinph0), this.B);
  68. } else {
  69. this.B = 1 / com;
  70. this.A = this.k0;
  71. this.E = D = F = 1;
  72. }
  73. if (alp || gam) {
  74. if (alp) {
  75. gamma0 = Math.asin(Math.sin(alpha_c) / D);
  76. if (!gam) {
  77. gamma = alpha_c;
  78. }
  79. } else {
  80. gamma0 = gamma;
  81. alpha_c = Math.asin(D * Math.sin(gamma0));
  82. }
  83. this.lam0 = lamc - Math.asin(0.5 * (F - 1 / F) * Math.tan(gamma0)) / this.B;
  84. } else {
  85. H = Math.pow(tsfnz(this.e, phi1, Math.sin(phi1)), this.B);
  86. L = Math.pow(tsfnz(this.e, phi2, Math.sin(phi2)), this.B);
  87. F = this.E / H;
  88. p = (L - H) / (L + H);
  89. J = this.E * this.E;
  90. J = (J - L * H) / (J + L * H);
  91. con = lam1 - lam2;
  92. if (con < -Math.pi) {
  93. lam2 -=TWO_PI;
  94. } else if (con > Math.pi) {
  95. lam2 += TWO_PI;
  96. }
  97. this.lam0 = adjust_lon(0.5 * (lam1 + lam2) - Math.atan(J * Math.tan(0.5 * this.B * (lam1 - lam2)) / p) / this.B);
  98. gamma0 = Math.atan(2 * Math.sin(this.B * adjust_lon(lam1 - this.lam0)) / (F - 1 / F));
  99. gamma = alpha_c = Math.asin(D * Math.sin(gamma0));
  100. }
  101. this.singam = Math.sin(gamma0);
  102. this.cosgam = Math.cos(gamma0);
  103. this.sinrot = Math.sin(gamma);
  104. this.cosrot = Math.cos(gamma);
  105. this.rB = 1 / this.B;
  106. this.ArB = this.A * this.rB;
  107. this.BrA = 1 / this.ArB;
  108. AB = this.A * this.B;
  109. if (this.no_off) {
  110. this.u_0 = 0;
  111. } else {
  112. this.u_0 = Math.abs(this.ArB * Math.atan(Math.sqrt(D * D - 1) / Math.cos(alpha_c)));
  113. if (this.lat0 < 0) {
  114. this.u_0 = - this.u_0;
  115. }
  116. }
  117. F = 0.5 * gamma0;
  118. this.v_pole_n = this.ArB * Math.log(Math.tan(FORTPI - F));
  119. this.v_pole_s = this.ArB * Math.log(Math.tan(FORTPI + F));
  120. }
  121. /* Oblique Mercator forward equations--mapping lat,long to x,y
  122. ----------------------------------------------------------*/
  123. export function forward(p) {
  124. var coords = {};
  125. var S, T, U, V, W, temp, u, v;
  126. p.x = p.x - this.lam0;
  127. if (Math.abs(Math.abs(p.y) - HALF_PI) > EPSLN) {
  128. W = this.E / Math.pow(tsfnz(this.e, p.y, Math.sin(p.y)), this.B);
  129. temp = 1 / W;
  130. S = 0.5 * (W - temp);
  131. T = 0.5 * (W + temp);
  132. V = Math.sin(this.B * p.x);
  133. U = (S * this.singam - V * this.cosgam) / T;
  134. if (Math.abs(Math.abs(U) - 1.0) < EPSLN) {
  135. throw new Error();
  136. }
  137. v = 0.5 * this.ArB * Math.log((1 - U)/(1 + U));
  138. temp = Math.cos(this.B * p.x);
  139. if (Math.abs(temp) < TOL) {
  140. u = this.A * p.x;
  141. } else {
  142. u = this.ArB * Math.atan2((S * this.cosgam + V * this.singam), temp);
  143. }
  144. } else {
  145. v = p.y > 0 ? this.v_pole_n : this.v_pole_s;
  146. u = this.ArB * p.y;
  147. }
  148. if (this.no_rot) {
  149. coords.x = u;
  150. coords.y = v;
  151. } else {
  152. u -= this.u_0;
  153. coords.x = v * this.cosrot + u * this.sinrot;
  154. coords.y = u * this.cosrot - v * this.sinrot;
  155. }
  156. coords.x = (this.a * coords.x + this.x0);
  157. coords.y = (this.a * coords.y + this.y0);
  158. return coords;
  159. }
  160. export function inverse(p) {
  161. var u, v, Qp, Sp, Tp, Vp, Up;
  162. var coords = {};
  163. p.x = (p.x - this.x0) * (1.0 / this.a);
  164. p.y = (p.y - this.y0) * (1.0 / this.a);
  165. if (this.no_rot) {
  166. v = p.y;
  167. u = p.x;
  168. } else {
  169. v = p.x * this.cosrot - p.y * this.sinrot;
  170. u = p.y * this.cosrot + p.x * this.sinrot + this.u_0;
  171. }
  172. Qp = Math.exp(-this.BrA * v);
  173. Sp = 0.5 * (Qp - 1 / Qp);
  174. Tp = 0.5 * (Qp + 1 / Qp);
  175. Vp = Math.sin(this.BrA * u);
  176. Up = (Vp * this.cosgam + Sp * this.singam) / Tp;
  177. if (Math.abs(Math.abs(Up) - 1) < EPSLN) {
  178. coords.x = 0;
  179. coords.y = Up < 0 ? -HALF_PI : HALF_PI;
  180. } else {
  181. coords.y = this.E / Math.sqrt((1 + Up) / (1 - Up));
  182. coords.y = phi2z(this.e, Math.pow(coords.y, 1 / this.B));
  183. if (coords.y === Infinity) {
  184. throw new Error();
  185. }
  186. coords.x = -this.rB * Math.atan2((Sp * this.cosgam - Vp * this.singam), Math.cos(this.BrA * u));
  187. }
  188. coords.x += this.lam0;
  189. return coords;
  190. }
  191. export var names = ["Hotine_Oblique_Mercator", "Hotine Oblique Mercator", "Hotine_Oblique_Mercator_Azimuth_Natural_Origin", "Hotine_Oblique_Mercator_Two_Point_Natural_Origin", "Hotine_Oblique_Mercator_Azimuth_Center", "Oblique_Mercator", "omerc"];
  192. export default {
  193. init: init,
  194. forward: forward,
  195. inverse: inverse,
  196. names: names
  197. };