aea.js 3.3 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129
  1. import msfnz from '../common/msfnz';
  2. import qsfnz from '../common/qsfnz';
  3. import adjust_lon from '../common/adjust_lon';
  4. import asinz from '../common/asinz';
  5. import {EPSLN} from '../constants/values';
  6. export function init() {
  7. if (Math.abs(this.lat1 + this.lat2) < EPSLN) {
  8. return;
  9. }
  10. this.temp = this.b / this.a;
  11. this.es = 1 - Math.pow(this.temp, 2);
  12. this.e3 = Math.sqrt(this.es);
  13. this.sin_po = Math.sin(this.lat1);
  14. this.cos_po = Math.cos(this.lat1);
  15. this.t1 = this.sin_po;
  16. this.con = this.sin_po;
  17. this.ms1 = msfnz(this.e3, this.sin_po, this.cos_po);
  18. this.qs1 = qsfnz(this.e3, this.sin_po);
  19. this.sin_po = Math.sin(this.lat2);
  20. this.cos_po = Math.cos(this.lat2);
  21. this.t2 = this.sin_po;
  22. this.ms2 = msfnz(this.e3, this.sin_po, this.cos_po);
  23. this.qs2 = qsfnz(this.e3, this.sin_po);
  24. this.sin_po = Math.sin(this.lat0);
  25. this.cos_po = Math.cos(this.lat0);
  26. this.t3 = this.sin_po;
  27. this.qs0 = qsfnz(this.e3, this.sin_po);
  28. if (Math.abs(this.lat1 - this.lat2) > EPSLN) {
  29. this.ns0 = (this.ms1 * this.ms1 - this.ms2 * this.ms2) / (this.qs2 - this.qs1);
  30. }
  31. else {
  32. this.ns0 = this.con;
  33. }
  34. this.c = this.ms1 * this.ms1 + this.ns0 * this.qs1;
  35. this.rh = this.a * Math.sqrt(this.c - this.ns0 * this.qs0) / this.ns0;
  36. }
  37. /* Albers Conical Equal Area forward equations--mapping lat,long to x,y
  38. -------------------------------------------------------------------*/
  39. export function forward(p) {
  40. var lon = p.x;
  41. var lat = p.y;
  42. this.sin_phi = Math.sin(lat);
  43. this.cos_phi = Math.cos(lat);
  44. var qs = qsfnz(this.e3, this.sin_phi);
  45. var rh1 = this.a * Math.sqrt(this.c - this.ns0 * qs) / this.ns0;
  46. var theta = this.ns0 * adjust_lon(lon - this.long0);
  47. var x = rh1 * Math.sin(theta) + this.x0;
  48. var y = this.rh - rh1 * Math.cos(theta) + this.y0;
  49. p.x = x;
  50. p.y = y;
  51. return p;
  52. }
  53. export function inverse(p) {
  54. var rh1, qs, con, theta, lon, lat;
  55. p.x -= this.x0;
  56. p.y = this.rh - p.y + this.y0;
  57. if (this.ns0 >= 0) {
  58. rh1 = Math.sqrt(p.x * p.x + p.y * p.y);
  59. con = 1;
  60. }
  61. else {
  62. rh1 = -Math.sqrt(p.x * p.x + p.y * p.y);
  63. con = -1;
  64. }
  65. theta = 0;
  66. if (rh1 !== 0) {
  67. theta = Math.atan2(con * p.x, con * p.y);
  68. }
  69. con = rh1 * this.ns0 / this.a;
  70. if (this.sphere) {
  71. lat = Math.asin((this.c - con * con) / (2 * this.ns0));
  72. }
  73. else {
  74. qs = (this.c - con * con) / this.ns0;
  75. lat = this.phi1z(this.e3, qs);
  76. }
  77. lon = adjust_lon(theta / this.ns0 + this.long0);
  78. p.x = lon;
  79. p.y = lat;
  80. return p;
  81. }
  82. /* Function to compute phi1, the latitude for the inverse of the
  83. Albers Conical Equal-Area projection.
  84. -------------------------------------------*/
  85. export function phi1z(eccent, qs) {
  86. var sinphi, cosphi, con, com, dphi;
  87. var phi = asinz(0.5 * qs);
  88. if (eccent < EPSLN) {
  89. return phi;
  90. }
  91. var eccnts = eccent * eccent;
  92. for (var i = 1; i <= 25; i++) {
  93. sinphi = Math.sin(phi);
  94. cosphi = Math.cos(phi);
  95. con = eccent * sinphi;
  96. com = 1 - con * con;
  97. dphi = 0.5 * com * com / cosphi * (qs / (1 - eccnts) - sinphi / com + 0.5 / eccent * Math.log((1 - con) / (1 + con)));
  98. phi = phi + dphi;
  99. if (Math.abs(dphi) <= 1e-7) {
  100. return phi;
  101. }
  102. }
  103. return null;
  104. }
  105. export var names = ["Albers_Conic_Equal_Area", "Albers", "aea"];
  106. export default {
  107. init: init,
  108. forward: forward,
  109. inverse: inverse,
  110. names: names,
  111. phi1z: phi1z
  112. };