Bienvenue sur PostGIS.fr

Bienvenue sur PostGIS.fr , le site de la communauté des utilisateurs francophones de PostGIS.

PostGIS ajoute le support d'objets géographique à la base de données PostgreSQL. En effet, PostGIS "spatialise" le serverur PostgreSQL, ce qui permet de l'utiliser comme une base de données SIG.

Maintenu à jour, en fonction de nos disponibilités et des diverses sorties des outils que nous testons, nous vous proposons l'ensemble de nos travaux publiés en langue française.

source: trunk/workshop-routing-foss4g/web/proj4js/lib/projCode/laea.js @ 76

Revision 76, 11.5 KB checked in by djay, 12 years ago (diff)

Ajout du répertoire web

  • Property svn:executable set to *
Line 
1/*******************************************************************************
2NAME                  LAMBERT AZIMUTHAL EQUAL-AREA
3 
4PURPOSE:        Transforms input longitude and latitude to Easting and
5                Northing for the Lambert Azimuthal Equal-Area projection.  The
6                longitude and latitude must be in radians.  The Easting
7                and Northing values will be returned in meters.
8
9PROGRAMMER              DATE           
10----------              ----           
11D. Steinwand, EROS      March, 1991   
12
13This function was adapted from the Lambert Azimuthal Equal Area projection
14code (FORTRAN) in the General Cartographic Transformation Package software
15which is available from the U.S. Geological Survey National Mapping Division.
16 
17ALGORITHM REFERENCES
18
191.  "New Equal-Area Map Projections for Noncircular Regions", John P. Snyder,
20    The American Cartographer, Vol 15, No. 4, October 1988, pp. 341-355.
21
222.  Snyder, John P., "Map Projections--A Working Manual", U.S. Geological
23    Survey Professional Paper 1395 (Supersedes USGS Bulletin 1532), United
24    State Government Printing Office, Washington D.C., 1987.
25
263.  "Software Documentation for GCTP General Cartographic Transformation
27    Package", U.S. Geological Survey National Mapping Division, May 1982.
28*******************************************************************************/
29
30Proj4js.Proj.laea = {
31  S_POLE: 1,
32  N_POLE: 2,
33  EQUIT: 3,
34  OBLIQ: 4,
35
36
37/* Initialize the Lambert Azimuthal Equal Area projection
38  ------------------------------------------------------*/
39  init: function() {
40    var t = Math.abs(this.lat0);
41    if (Math.abs(t - Proj4js.common.HALF_PI) < Proj4js.common.EPSLN) {
42      this.mode = this.lat0 < 0. ? this.S_POLE : this.N_POLE;
43    } else if (Math.abs(t) < Proj4js.common.EPSLN) {
44      this.mode = this.EQUIT;
45    } else {
46      this.mode = this.OBLIQ;
47    }
48    if (this.es > 0) {
49      var sinphi;
50 
51      this.qp = Proj4js.common.qsfnz(this.e, 1.0);
52      this.mmf = .5 / (1. - this.es);
53      this.apa = this.authset(this.es);
54      switch (this.mode) {
55        case this.N_POLE:
56        case this.S_POLE:
57          this.dd = 1.;
58          break;
59        case this.EQUIT:
60          this.rq = Math.sqrt(.5 * this.qp);
61          this.dd = 1. / this.rq;
62          this.xmf = 1.;
63          this.ymf = .5 * this.qp;
64          break;
65        case this.OBLIQ:
66          this.rq = Math.sqrt(.5 * this.qp);
67          sinphi = Math.sin(this.lat0);
68          this.sinb1 = Proj4js.common.qsfnz(this.e, sinphi) / this.qp;
69          this.cosb1 = Math.sqrt(1. - this.sinb1 * this.sinb1);
70          this.dd = Math.cos(this.lat0) / (Math.sqrt(1. - this.es * sinphi * sinphi) * this.rq * this.cosb1);
71          this.ymf = (this.xmf = this.rq) / this.dd;
72          this.xmf *= this.dd;
73          break;
74      }
75    } else {
76      if (this.mode == this.OBLIQ) {
77        this.sinph0 = Math.sin(this.lat0);
78        this.cosph0 = Math.cos(this.lat0);
79      }
80    }
81  },
82
83/* Lambert Azimuthal Equal Area forward equations--mapping lat,long to x,y
84  -----------------------------------------------------------------------*/
85  forward: function(p) {
86
87    /* Forward equations
88      -----------------*/
89    var x,y;
90    var lam=p.x;
91    var phi=p.y;
92    lam = Proj4js.common.adjust_lon(lam - this.long0);
93   
94    if (this.sphere) {
95        var coslam, cosphi, sinphi;
96     
97        sinphi = Math.sin(phi);
98        cosphi = Math.cos(phi);
99        coslam = Math.cos(lam);
100        switch (this.mode) {
101          case this.OBLIQ:
102          case this.EQUIT:
103            y = (this.mode == this.EQUIT) ? 1. + cosphi * coslam : 1. + this.sinph0 * sinphi + this.cosph0 * cosphi * coslam;
104            if (y <= Proj4js.common.EPSLN) {
105              Proj4js.reportError("laea:fwd:y less than eps");
106              return null;
107            }
108            y = Math.sqrt(2. / y);
109            x = y * cosphi * Math.sin(lam);
110            y *= (this.mode == this.EQUIT) ? sinphi : this.cosph0 * sinphi - this.sinph0 * cosphi * coslam;
111            break;
112          case this.N_POLE:
113            coslam = -coslam;
114          case this.S_POLE:
115            if (Math.abs(phi + this.phi0) < Proj4js.common.EPSLN) {
116              Proj4js.reportError("laea:fwd:phi < eps");
117              return null;
118            }
119            y = Proj4js.common.FORTPI - phi * .5;
120            y = 2. * ((this.mode == this.S_POLE) ? Math.cos(y) : Math.sin(y));
121            x = y * Math.sin(lam);
122            y *= coslam;
123            break;
124        }
125    } else {
126        var coslam, sinlam, sinphi, q, sinb=0.0, cosb=0.0, b=0.0;
127     
128        coslam = Math.cos(lam);
129        sinlam = Math.sin(lam);
130        sinphi = Math.sin(phi);
131        q = Proj4js.common.qsfnz(this.e, sinphi);
132        if (this.mode == this.OBLIQ || this.mode == this.EQUIT) {
133          sinb = q / this.qp;
134          cosb = Math.sqrt(1. - sinb * sinb);
135        }
136        switch (this.mode) {
137          case this.OBLIQ:
138            b = 1. + this.sinb1 * sinb + this.cosb1 * cosb * coslam;
139            break;
140          case this.EQUIT:
141            b = 1. + cosb * coslam;
142            break;
143          case this.N_POLE:
144            b = Proj4js.common.HALF_PI + phi;
145            q = this.qp - q;
146            break;
147          case this.S_POLE:
148            b = phi - Proj4js.common.HALF_PI;
149            q = this.qp + q;
150            break;
151        }
152        if (Math.abs(b) < Proj4js.common.EPSLN) {
153            Proj4js.reportError("laea:fwd:b < eps");
154            return null;
155        }
156        switch (this.mode) {
157          case this.OBLIQ:
158          case this.EQUIT:
159            b = Math.sqrt(2. / b);
160            if (this.mode == this.OBLIQ) {
161              y = this.ymf * b * (this.cosb1 * sinb - this.sinb1 * cosb * coslam);
162            } else {
163              y = (b = Math.sqrt(2. / (1. + cosb * coslam))) * sinb * this.ymf;
164            }
165            x = this.xmf * b * cosb * sinlam;
166            break;
167          case this.N_POLE:
168          case this.S_POLE:
169            if (q >= 0.) {
170              x = (b = Math.sqrt(q)) * sinlam;
171              y = coslam * ((this.mode == this.S_POLE) ? b : -b);
172            } else {
173              x = y = 0.;
174            }
175            break;
176        }
177    }
178
179    //v 1.0
180    /*
181    var sin_lat=Math.sin(lat);
182    var cos_lat=Math.cos(lat);
183
184    var sin_delta_lon=Math.sin(delta_lon);
185    var cos_delta_lon=Math.cos(delta_lon);
186
187    var g =this.sin_lat_o * sin_lat +this.cos_lat_o * cos_lat * cos_delta_lon;
188    if (g == -1.0) {
189      Proj4js.reportError("laea:fwd:Point projects to a circle of radius "+ 2.0 * R);
190      return null;
191    }
192    var ksp = this.a * Math.sqrt(2.0 / (1.0 + g));
193    var x = ksp * cos_lat * sin_delta_lon + this.x0;
194    var y = ksp * (this.cos_lat_o * sin_lat - this.sin_lat_o * cos_lat * cos_delta_lon) + this.y0;
195    */
196    p.x = this.a*x + this.x0;
197    p.y = this.a*y + this.y0;
198    return p;
199  },//lamazFwd()
200
201/* Inverse equations
202  -----------------*/
203  inverse: function(p) {
204    p.x -= this.x0;
205    p.y -= this.y0;
206    var x = p.x/this.a;
207    var y = p.y/this.a;
208   
209    if (this.sphere) {
210        var  cosz=0.0, rh, sinz=0.0;
211     
212        rh = Math.sqrt(x*x + y*y);
213        var phi = rh * .5;
214        if (phi > 1.) {
215          Proj4js.reportError("laea:Inv:DataError");
216          return null;
217        }
218        phi = 2. * Math.asin(phi);
219        if (this.mode == this.OBLIQ || this.mode == this.EQUIT) {
220          sinz = Math.sin(phi);
221          cosz = Math.cos(phi);
222        }
223        switch (this.mode) {
224        case this.EQUIT:
225          phi = (Math.abs(rh) <= Proj4js.common.EPSLN) ? 0. : Math.asin(y * sinz / rh);
226          x *= sinz;
227          y = cosz * rh;
228          break;
229        case this.OBLIQ:
230          phi = (Math.abs(rh) <= Proj4js.common.EPSLN) ? this.phi0 : Math.asin(cosz * sinph0 + y * sinz * cosph0 / rh);
231          x *= sinz * cosph0;
232          y = (cosz - Math.sin(phi) * sinph0) * rh;
233          break;
234        case this.N_POLE:
235          y = -y;
236          phi = Proj4js.common.HALF_PI - phi;
237          break;
238        case this.S_POLE:
239          phi -= Proj4js.common.HALF_PI;
240          break;
241        }
242        lam = (y == 0. && (this.mode == this.EQUIT || this.mode == this.OBLIQ)) ? 0. : Math.atan2(x, y);
243    } else {
244        var cCe, sCe, q, rho, ab=0.0;
245     
246        switch (this.mode) {
247          case this.EQUIT:
248          case this.OBLIQ:
249            x /= this.dd;
250            y *=  this.dd;
251            rho = Math.sqrt(x*x + y*y);
252            if (rho < Proj4js.common.EPSLN) {
253              p.x = 0.;
254              p.y = this.phi0;
255              return p;
256            }
257            sCe = 2. * Math.asin(.5 * rho / this.rq);
258            cCe = Math.cos(sCe);
259            x *= (sCe = Math.sin(sCe));
260            if (this.mode == this.OBLIQ) {
261              ab = cCe * this.sinb1 + y * sCe * this.cosb1 / rho
262              q = this.qp * ab;
263              y = rho * this.cosb1 * cCe - y * this.sinb1 * sCe;
264            } else {
265              ab = y * sCe / rho;
266              q = this.qp * ab;
267              y = rho * cCe;
268            }
269            break;
270          case this.N_POLE:
271            y = -y;
272          case this.S_POLE:
273            q = (x * x + y * y);
274            if (!q ) {
275              p.x = 0.;
276              p.y = this.phi0;
277              return p;
278            }
279            /*
280            q = this.qp - q;
281            */
282            ab = 1. - q / this.qp;
283            if (this.mode == this.S_POLE) {
284              ab = - ab;
285            }
286            break;
287        }
288        lam = Math.atan2(x, y);
289        phi = this.authlat(Math.asin(ab), this.apa);
290    }
291
292    /*
293    var Rh = Math.Math.sqrt(p.x *p.x +p.y * p.y);
294    var temp = Rh / (2.0 * this.a);
295
296    if (temp > 1) {
297      Proj4js.reportError("laea:Inv:DataError");
298      return null;
299    }
300
301    var z = 2.0 * Proj4js.common.asinz(temp);
302    var sin_z=Math.sin(z);
303    var cos_z=Math.cos(z);
304
305    var lon =this.long0;
306    if (Math.abs(Rh) > Proj4js.common.EPSLN) {
307       var lat = Proj4js.common.asinz(this.sin_lat_o * cos_z +this. cos_lat_o * sin_z *p.y / Rh);
308       var temp =Math.abs(this.lat0) - Proj4js.common.HALF_PI;
309       if (Math.abs(temp) > Proj4js.common.EPSLN) {
310          temp = cos_z -this.sin_lat_o * Math.sin(lat);
311          if(temp!=0.0) lon=Proj4js.common.adjust_lon(this.long0+Math.atan2(p.x*sin_z*this.cos_lat_o,temp*Rh));
312       } else if (this.lat0 < 0.0) {
313          lon = Proj4js.common.adjust_lon(this.long0 - Math.atan2(-p.x,p.y));
314       } else {
315          lon = Proj4js.common.adjust_lon(this.long0 + Math.atan2(p.x, -p.y));
316       }
317    } else {
318      lat = this.lat0;
319    }
320    */
321    //return(OK);
322    p.x = Proj4js.common.adjust_lon(this.long0+lam);
323    p.y = phi;
324    return p;
325  },//lamazInv()
326 
327/* determine latitude from authalic latitude */
328  P00: .33333333333333333333,
329  P01: .17222222222222222222,
330  P02: .10257936507936507936,
331  P10: .06388888888888888888,
332  P11: .06640211640211640211,
333  P20: .01641501294219154443,
334 
335  authset: function(es) {
336    var t;
337    var APA = new Array();
338    APA[0] = es * this.P00;
339    t = es * es;
340    APA[0] += t * this.P01;
341    APA[1] = t * this.P10;
342    t *= es;
343    APA[0] += t * this.P02;
344    APA[1] += t * this.P11;
345    APA[2] = t * this.P20;
346    return APA;
347  },
348 
349  authlat: function(beta, APA) {
350    var t = beta+beta;
351    return(beta + APA[0] * Math.sin(t) + APA[1] * Math.sin(t+t) + APA[2] * Math.sin(t+t+t));
352  }
353 
354};
355
356
357
Note: See TracBrowser for help on using the repository browser.