|
1
|
/*
|
|
2
|
Author: Mike Adair mike.adairATccrs.nrcan.gc.ca
|
|
3
|
Richard Greenwood rich@greenwoodmap.com
|
|
4
|
License: LGPL as per: http://www.gnu.org/copyleft/lesser.html
|
|
5
|
$Id$
|
|
6
|
*/
|
|
7
|
|
|
8
|
/**
|
|
9
|
* Provides latitude/longitude to map projection (and vice versa) transformation methods.
|
|
10
|
* Initialized with EPSG codes. Has properties for units and title strings.
|
|
11
|
* All coordinates are handled as points which is a 2 element array where x is
|
|
12
|
* the first element and y is the second.
|
|
13
|
* For the Forward() method pass in lat/lon and it returns map XY.
|
|
14
|
* For the Inverse() method pass in map XY and it returns lat/long.
|
|
15
|
*
|
|
16
|
* TBD: retrieve initialization params (and conversion code?) from a web service
|
|
17
|
*
|
|
18
|
* @constructor
|
|
19
|
*
|
|
20
|
* @param srs The SRS (EPSG code) of the projection
|
|
21
|
*/
|
|
22
|
|
|
23
|
function Proj(srs) {
|
|
24
|
this.srs = srs.toUpperCase();
|
|
25
|
switch(this.srs) {
|
|
26
|
case "EPSG:4326": //lat/lon projection WGS_84
|
|
27
|
case "EPSG:4269": //lat/lon projection WGS_84
|
|
28
|
case "CRS:84": //lat/lon projection WGS_84
|
|
29
|
case "EPSG:4965": //lat/lon projection RGF93G IGN-F FD 2005
|
|
30
|
case new String("http://www.opengis.net/gml/srs/epsg.xml#4326").toUpperCase():
|
|
31
|
this.Forward = identity;
|
|
32
|
this.Inverse = identity;
|
|
33
|
this.units = "degrees";
|
|
34
|
this.title = "Lat/Long";
|
|
35
|
break;
|
|
36
|
case "EPSG:42101": //North American LCC WGS_84
|
|
37
|
this.Init = lccinit;
|
|
38
|
this.Forward = ll2lcc;
|
|
39
|
this.Inverse = lcc2ll;
|
|
40
|
this.Init(new Array(6378137.0,6356752.314245,49.0,77.0,-95.0, 0.0, 0.0, -8000000));
|
|
41
|
this.units = "meters";
|
|
42
|
this.title = "Lambert Conformal Conic";
|
|
43
|
break;
|
|
44
|
case "EPSG:42304": //North American LCC NAD_83
|
|
45
|
this.Init = lccinit;
|
|
46
|
this.Forward = ll2lcc;
|
|
47
|
this.Inverse = lcc2ll;
|
|
48
|
this.Init(new Array(6378137.0,6356752.314,49.0,77.0,-95.0, 49.0, 0.0, 0));
|
|
49
|
this.units = "meters";
|
|
50
|
this.title = "Lambert Conformal Conic";
|
|
51
|
break;
|
|
52
|
case "EPSG:26986": //NAD83 / Massachusetts Mainland
|
|
53
|
this.Init = lccinit;
|
|
54
|
this.Forward = ll2lcc;
|
|
55
|
this.Inverse = lcc2ll;
|
|
56
|
this.Init(new Array(6378137.0,6356752.314,42.68333333333333,41.71666666666667,-71.5, 41.0, 200000, 750000));
|
|
57
|
this.units = "meters";
|
|
58
|
this.title = "Massachusetts Mainland (LCC)";
|
|
59
|
break;
|
|
60
|
case "EPSG:32761": //Polar Stereographic
|
|
61
|
case "EPSG:32661":
|
|
62
|
this.Init = psinit;
|
|
63
|
this.Forward = ll2ps;
|
|
64
|
this.Inverse = ps2ll;
|
|
65
|
this.Init(new Array(6378137.0, 6356752.314245, 0.0, -90.0, 2000000, 2000000));
|
|
66
|
this.units = "meters";
|
|
67
|
this.title = "Polar Stereographic";
|
|
68
|
break;
|
|
69
|
case "EPSG:27561"://LAMB1 IGN-F modification FD 2005
|
|
70
|
this.Init = lccinit;
|
|
71
|
this.Forward = ll2lcc;
|
|
72
|
this.Inverse = lcc2ll;
|
|
73
|
this.Init(new Array(6378249.2, 6356515.0, 49.50, 49.50, 2.33722916655, 49.50, 600000.0, 200000.0));
|
|
74
|
this.units = "meters";
|
|
75
|
this.title = "Lambert Conformal Conic";
|
|
76
|
break;
|
|
77
|
case "EPSG:27562"://LAMB2 IGN-F modification FD 2005
|
|
78
|
this.Init = lccinit;
|
|
79
|
this.Forward = ll2lcc;
|
|
80
|
this.Inverse = lcc2ll;
|
|
81
|
this.Init(new Array(6378249.2, 6356515.0, 46.80, 46.80, 2.33722916655, 46.8, 600000.0, 200000.0));
|
|
82
|
this.units = "meters";
|
|
83
|
this.title = "Lambert Conformal Conic";
|
|
84
|
break;
|
|
85
|
case "EPSG:27572"://LAMBE IGN-F modification FD 2005
|
|
86
|
case "EPSG:27582"://LAMB2 IGN-F modification (deprecated) FD 2005
|
|
87
|
this.Init = lccinit;
|
|
88
|
this.Forward = ll2lcc;
|
|
89
|
this.Inverse = lcc2ll;
|
|
90
|
this.Init(new Array(6378249.2, 6356515.0, 46.80, 46.80, 2.33722916655, 46.8, 600000.0, 2200000.0));
|
|
91
|
this.units = "meters";
|
|
92
|
this.title = "Lambert Conformal Conic";
|
|
93
|
break;
|
|
94
|
case "EPSG:27563"://LAMB3 IGN-F modification FD 2005
|
|
95
|
this.Init = lccinit;
|
|
96
|
this.Forward = ll2lcc;
|
|
97
|
this.Inverse = lcc2ll;
|
|
98
|
this.Init(new Array(6378249.2, 6356515.0, 44.10, 44.10, 2.33722916655, 44.10, 600000.0, 200000.0));
|
|
99
|
this.units = "meters";
|
|
100
|
this.title = "Lambert Conformal Conic";
|
|
101
|
break;
|
|
102
|
case "EPSG:27564"://LAMB4 IGN-F modification FD 2005
|
|
103
|
this.Init = lccinit;
|
|
104
|
this.Forward = ll2lcc;
|
|
105
|
this.Inverse = lcc2ll;
|
|
106
|
this.Init(new Array(6378249.2, 6356515.0, 42.17, 42.17, 2.33722916655, 42.17, 234.358, 185861.369));
|
|
107
|
this.units = "meters";
|
|
108
|
this.title = "Lambert Conformal Conic";
|
|
109
|
break;
|
|
110
|
case "EPSG:2154" ://LAMB93 IGN-F modification FD 2005
|
|
111
|
this.Init = lccinit;
|
|
112
|
this.Forward = ll2lcc;
|
|
113
|
this.Inverse = lcc2ll;
|
|
114
|
this.Init(new Array(6378137.0, 6356752.3141, 44.00, 49.00, 3.00000000001, 46.50, 700000.0, 6600000.0));
|
|
115
|
this.units = "meters";
|
|
116
|
this.title = "Lambert Conformal Conic";
|
|
117
|
break;
|
|
118
|
case "EPSG:4326":case "EPSG:4269":case "CRS:84":case "EPSG:4965":case new String("http://www.opengis.net/gml/srs/epsg.xml#4326").toUpperCase():
|
|
119
|
this.Forward=identity;
|
|
120
|
this.Inverse=identity;
|
|
121
|
this.units="degrees";
|
|
122
|
this.title="Lat/Long";
|
|
123
|
break;
|
|
124
|
case "EPSG:102758":this.title="NAD 1983 StatePlane Wyoming West FIPS 4904 US Survey Feet";
|
|
125
|
this.Init=tminit;
|
|
126
|
this.Forward=ll2tm
|
|
127
|
this.Inverse=tm2ll;
|
|
128
|
this.Init(new Array(grs80[0], grs80[1], 0.9999375, -110.0833333333333, 40.5, 800000, 100000));
|
|
129
|
this.units="usfeet";
|
|
130
|
break;
|
|
131
|
case "EPSG:32158":this.title="NAD 1983 StatePlane Wyoming West meters";
|
|
132
|
this.Init=tminit;
|
|
133
|
this.Forward=ll2tm
|
|
134
|
this.Inverse=tm2ll;
|
|
135
|
this.Init(new Array(grs80[0], grs80[1], 0.9999375, -110.0833333333333, 40.5, 800000, 100000));
|
|
136
|
this.units="meters";
|
|
137
|
break;
|
|
138
|
// UTM NAD83 Zones 3 thru 23
|
|
139
|
case"EPSG:26903":case"EPSG:26904":case"EPSG:26905":case"EPSG:26906":case"EPSG:26907":case"EPSG:26908":case"EPSG:26909":
|
|
140
|
case"EPSG:26910":case"EPSG:26911":case"EPSG:26912":case"EPSG:26913":case"EPSG:26914":case"EPSG:26915":case"EPSG:26916":
|
|
141
|
case"EPSG:26917":case"EPSG:26918":case"EPSG:26919":case"EPSG:26920":case"EPSG:26921":case"EPSG:26922":case"EPSG:26923":
|
|
142
|
this.title="NAD83 / UTM zone "+ srs.substr(8,2)+"N";
|
|
143
|
this.Init=utminit;
|
|
144
|
this.Forward=ll2tm;
|
|
145
|
this.Inverse=tm2ll;
|
|
146
|
this.Init(new Array(grs80[0], grs80[1], 0.9996, srs.substr(8,2)));
|
|
147
|
this.units="meters";
|
|
148
|
break;
|
|
149
|
// UTM WGS84 Zones 1 thru 60 North
|
|
150
|
case"EPSG:32601":case"EPSG:32602":
|
|
151
|
case"EPSG:32603":case"EPSG:32604":case"EPSG:32605":case"EPSG:32606":case"EPSG:32607":case"EPSG:32608":case"EPSG:32609":
|
|
152
|
case"EPSG:32610":case"EPSG:32611":case"EPSG:32612":case"EPSG:32613":case"EPSG:32614":case"EPSG:32615":case"EPSG:32616":
|
|
153
|
case"EPSG:32617":case"EPSG:32618":case"EPSG:32619":case"EPSG:32620":case"EPSG:32621":case"EPSG:32622":case"EPSG:32623":
|
|
154
|
case"EPSG:32624":case"EPSG:32625":case"EPSG:32626":case"EPSG:32627":case"EPSG:32628":case"EPSG:32629":
|
|
155
|
case"EPSG:32630":case"EPSG:32631":case"EPSG:32632":case"EPSG:32633":case"EPSG:32634":case"EPSG:32635":case"EPSG:32636":
|
|
156
|
case"EPSG:32637":case"EPSG:32638":case"EPSG:32639":case"EPSG:32640":case"EPSG:32641":case"EPSG:32642":
|
|
157
|
case"EPSG:32643":case"EPSG:32644":case"EPSG:32645":case"EPSG:32646":case"EPSG:32647":case"EPSG:32648":case"EPSG:32649":
|
|
158
|
case"EPSG:32650":case"EPSG:32651":case"EPSG:32652":case"EPSG:32653":case"EPSG:32654":case"EPSG:32655":case"EPSG:32656":
|
|
159
|
case"EPSG:32657":case"EPSG:32658":case"EPSG:32659":case"EPSG:32660":
|
|
160
|
this.title="WGS84 / UTM zone " + srs.substr(8,2) + "N";
|
|
161
|
this.Init=utminit;
|
|
162
|
this.Forward=ll2tm;
|
|
163
|
this.Inverse=tm2ll;
|
|
164
|
this.Init(new Array(wgs84[0], wgs84[1], 0.9996, srs.substr(8,2)));
|
|
165
|
this.units="meters";
|
|
166
|
break;
|
|
167
|
// UTM WGS84 Zones 1 thru 60 South
|
|
168
|
case"EPSG:32701":case"EPSG:32702":
|
|
169
|
case"EPSG:32703":case"EPSG:32704":case"EPSG:32705":case"EPSG:32706":case"EPSG:32707":case"EPSG:32708":case"EPSG:32709":
|
|
170
|
case"EPSG:32710":case"EPSG:32711":case"EPSG:32712":case"EPSG:32713":case"EPSG:32714":case"EPSG:32715":case"EPSG:32716":
|
|
171
|
case"EPSG:32717":case"EPSG:32718":case"EPSG:32719":case"EPSG:32720":case"EPSG:32721":case"EPSG:32722":case"EPSG:32723":
|
|
172
|
case"EPSG:32724":case"EPSG:32725":case"EPSG:32726":case"EPSG:32727":case"EPSG:32728":case"EPSG:32729":
|
|
173
|
case"EPSG:32730":case"EPSG:32731":case"EPSG:32732":case"EPSG:32733":case"EPSG:32734":case"EPSG:32735":case"EPSG:32736":
|
|
174
|
case"EPSG:32737":case"EPSG:32738":case"EPSG:32739":case"EPSG:32740":case"EPSG:32741":case"EPSG:32742":
|
|
175
|
case"EPSG:32743":case"EPSG:32744":case"EPSG:32745":case"EPSG:32746":case"EPSG:32747":case"EPSG:32748":case"EPSG:32749":
|
|
176
|
case"EPSG:32750":case"EPSG:32751":case"EPSG:32752":case"EPSG:32753":case"EPSG:32754":case"EPSG:32755":case"EPSG:32756":
|
|
177
|
case"EPSG:32757":case"EPSG:32758":case"EPSG:32759":case"EPSG:32760":
|
|
178
|
this.title="WGS84 / UTM zone " + srs.substr(8,2) + "S";
|
|
179
|
this.Init=utminit;
|
|
180
|
this.Forward=ll2tm;
|
|
181
|
this.Inverse=tm2ll;
|
|
182
|
this.Init(new Array(wgs84[0], wgs84[1], 0.9996, "-"+srs.substr(8,2)));
|
|
183
|
this.units="meters";
|
|
184
|
break;
|
|
185
|
case "EPSG:26591":this.title="Monte Mario (Rome) / Italy zone 1";
|
|
186
|
this.Init=tminit;
|
|
187
|
this.Forward=ll2tm
|
|
188
|
this.Inverse=tm2ll;
|
|
189
|
this.Init(new Array(6378388.0, 6356911.94612795,0.9996, 9, 0.0, 1500000.0, 0.0));
|
|
190
|
this.units="meters";
|
|
191
|
break;
|
|
192
|
case "SCENE": //this is really a pixel projection with bilinear interpolation of the corners to get ll
|
|
193
|
this.Init = sceneInit;
|
|
194
|
this.Forward = ll2scene;
|
|
195
|
this.Inverse = scene2ll;
|
|
196
|
this.GetXYCoords = identity; //override to get line 0 at top left
|
|
197
|
this.GetPLCoords = identity; //
|
|
198
|
break;
|
|
199
|
case "PIXEL":
|
|
200
|
this.Forward = ll2pixel;
|
|
201
|
this.Inverse = pixel2ll;
|
|
202
|
this.units = "pixels";
|
|
203
|
this.GetXYCoords = identity; //override to get line 0 at top left
|
|
204
|
this.GetPLCoords = identity; //
|
|
205
|
break;
|
|
206
|
default:
|
|
207
|
//or retrieve parameters from web service based on SRS lookup
|
|
208
|
alert("unsupported map projection: "+this.srs);
|
|
209
|
}
|
|
210
|
|
|
211
|
this.matchSrs = function(otherSrs) {
|
|
212
|
if (this.srs == otherSrs.toUpperCase() ) return true;
|
|
213
|
return false;
|
|
214
|
}
|
|
215
|
|
|
216
|
}
|
|
217
|
|
|
218
|
function identity(coords) {
|
|
219
|
return coords;
|
|
220
|
}
|
|
221
|
|
|
222
|
/**
|
|
223
|
* Scene projection forward transformation.
|
|
224
|
* Forward trasnformation need reverse bilinear interpolation or orbit modelling
|
|
225
|
* (to be implemented)
|
|
226
|
* @param coords Lat/Long coords passed in
|
|
227
|
* @return map coordinates
|
|
228
|
*/
|
|
229
|
function ll2scene(coords) {
|
|
230
|
alert("ll2scene not defined");
|
|
231
|
//return new Array(124, 15+256); //for testing only,
|
|
232
|
return null;
|
|
233
|
}
|
|
234
|
|
|
235
|
/**
|
|
236
|
* Scene projection Inverse transformation.
|
|
237
|
* This is really a pixel representation with bi-linear interpolation of the corner coords.
|
|
238
|
* @param coords map coordinates passed in
|
|
239
|
* @return Lat/Long coords
|
|
240
|
*/
|
|
241
|
function scene2ll(coords) {
|
|
242
|
var xpct = (coords[0]-this.ul[0])/(this.lr[0]-this.ul[0]);
|
|
243
|
var ypct = (coords[1]-this.ul[1])/(this.lr[1]-this.ul[1]);
|
|
244
|
// alert("pct:"+xpct+":"+ypct);
|
|
245
|
var lon = bilinterp(xpct, ypct, this.cul[0], this.cur[0], this.cll[0], this.clr[0])
|
|
246
|
var lat = bilinterp(xpct, ypct, this.cul[1], this.cur[1], this.cll[1], this.clr[1])
|
|
247
|
return new Array(lon, lat);
|
|
248
|
}
|
|
249
|
|
|
250
|
/**
|
|
251
|
* Scene projection initialization function
|
|
252
|
* @param param array of the corner coordinates (which are in turn 2D point arrays)
|
|
253
|
* in the order upper-left, upper-right, lower-left, lower-right
|
|
254
|
*/
|
|
255
|
function sceneInit(param) {
|
|
256
|
this.cul = param[0];
|
|
257
|
this.cur = param[1];
|
|
258
|
this.cll = param[2];
|
|
259
|
this.clr = param[3];
|
|
260
|
}
|
|
261
|
|
|
262
|
/**
|
|
263
|
* Bilinear interpolation function to return an interpolated value for either
|
|
264
|
* x or y for some point located within a box where the XY of the corners are known.
|
|
265
|
* This should be applied to the x and y coordinates separately,
|
|
266
|
* ie. first interpolate for the x value, then the y.
|
|
267
|
* the a,b,c,d params are thus the x or y values at the corners.
|
|
268
|
* @param x distance from the left side of the box as a percentage
|
|
269
|
* @param y distance from the top of the box as a percentage
|
|
270
|
* @param a either x or y value at the UL corner of the box
|
|
271
|
* @param b either x or y value at the UR corner of the box
|
|
272
|
* @param c either x or y value at the LL corner of the box
|
|
273
|
* @param d either x or y value at the LR corner of the box
|
|
274
|
*/
|
|
275
|
function bilinterp(x, y, a, b, c, d) {
|
|
276
|
var top = x*(b-a) + a;
|
|
277
|
var bot = x*(d-c) + c;
|
|
278
|
//alert("top:"+top+" bot:"+bot);
|
|
279
|
return y*(bot-top) + top;
|
|
280
|
}
|
|
281
|
|
|
282
|
// pixel projection object definition
|
|
283
|
// a pixel representation
|
|
284
|
// forward transformation
|
|
285
|
function ll2pixel(coords) {
|
|
286
|
alert("ll2pixel not defined");
|
|
287
|
return null;
|
|
288
|
}
|
|
289
|
|
|
290
|
// inverse transformation
|
|
291
|
function pixel2ll(coords) {
|
|
292
|
alert("pixel2ll not defined");
|
|
293
|
// return new Array(lon, lat);
|
|
294
|
return null;
|
|
295
|
}
|
|
296
|
|
|
297
|
|
|
298
|
/****************************************************************************
|
|
299
|
The following code is a direct port of GCTPC coordinate transformation code
|
|
300
|
from C to Javascript. For more information go to http://edcftp.cr.usgs.gov/pub//software/gctpc/
|
|
301
|
Currently suppported projections include: Lambert Conformal Conic (LCC),
|
|
302
|
Lat/Long, Polar Stereographic, Transverse Mercator (including UTM).
|
|
303
|
Porting C to Javascript is fairly straightforward so other support for more
|
|
304
|
projections is easy to add.
|
|
305
|
*/
|
|
306
|
|
|
307
|
var PI = Math.PI;
|
|
308
|
var HALF_PI = PI*0.5;
|
|
309
|
var TWO_PI = PI*2.0;
|
|
310
|
var EPSLN = 1.0e-10;
|
|
311
|
var R2D = 57.2957795131;
|
|
312
|
var D2R =0.0174532925199;
|
|
313
|
var R = 6370997.0; // Radius of the earth (sphere)
|
|
314
|
|
|
315
|
|
|
316
|
// Initialize the Lambert Conformal conic projection
|
|
317
|
// -----------------------------------------------------------------
|
|
318
|
function lccinit(param) {
|
|
319
|
// array of: r_maj,r_min,lat1,lat2,c_lon,c_lat,false_east,false_north
|
|
320
|
//double c_lat; /* center latitude */
|
|
321
|
//double c_lon; /* center longitude */
|
|
322
|
//double lat1; /* first standard parallel */
|
|
323
|
//double lat2; /* second standard parallel */
|
|
324
|
//double r_maj; /* major axis */
|
|
325
|
//double r_min; /* minor axis */
|
|
326
|
//double false_east; /* x offset in meters */
|
|
327
|
//double false_north; /* y offset in meters */
|
|
328
|
|
|
329
|
this.r_major = param[0];
|
|
330
|
this.r_minor = param[1];
|
|
331
|
var lat1 = param[2] * D2R;
|
|
332
|
var lat2 = param[3] * D2R;
|
|
333
|
this.center_lon = param[4] * D2R;
|
|
334
|
this.center_lat = param[5] * D2R;
|
|
335
|
this.false_easting = param[6];
|
|
336
|
this.false_northing = param[7];
|
|
337
|
|
|
338
|
// Standard Parallels cannot be equal and on opposite sides of the equator
|
|
339
|
if (Math.abs(lat1+lat2) < EPSLN) {
|
|
340
|
alert("Equal Latitiudes for St. Parallels on opposite sides of equator - lccinit");
|
|
341
|
return;
|
|
342
|
}
|
|
343
|
|
|
344
|
var temp = this.r_minor / this.r_major;
|
|
345
|
this.e = Math.sqrt(1.0 - temp*temp);
|
|
346
|
|
|
347
|
var sin1 = Math.sin(lat1);
|
|
348
|
var cos1 = Math.cos(lat1);
|
|
349
|
var ms1 = msfnz(this.e, sin1, cos1);
|
|
350
|
var ts1 = tsfnz(this.e, lat1, sin1);
|
|
351
|
|
|
352
|
var sin2 = Math.sin(lat2);
|
|
353
|
var cos2 = Math.cos(lat2);
|
|
354
|
var ms2 = msfnz(this.e, sin2, cos2);
|
|
355
|
var ts2 = tsfnz(this.e, lat2, sin2);
|
|
356
|
|
|
357
|
var ts0 = tsfnz(this.e, this.center_lat, Math.sin(this.center_lat));
|
|
358
|
|
|
359
|
if (Math.abs(lat1 - lat2) > EPSLN) {
|
|
360
|
this.ns = Math.log(ms1/ms2)/Math.log(ts1/ts2);
|
|
361
|
} else {
|
|
362
|
this.ns = sin1;
|
|
363
|
}
|
|
364
|
this.f0 = ms1 / (this.ns * Math.pow(ts1, this.ns));
|
|
365
|
this.rh = this.r_major * this.f0 * Math.pow(ts0, this.ns);
|
|
366
|
}
|
|
367
|
|
|
368
|
|
|
369
|
// Lambert Conformal conic forward equations--mapping lat,long to x,y
|
|
370
|
// -----------------------------------------------------------------
|
|
371
|
function ll2lcc(coords) {
|
|
372
|
|
|
373
|
var lon = coords[0];
|
|
374
|
var lat = coords[1];
|
|
375
|
|
|
376
|
// convert to radians
|
|
377
|
if ( lat <= 90.0 && lat >= -90.0 && lon <= 180.0 && lon >= -180.0) {
|
|
378
|
lat *= D2R;
|
|
379
|
lon *= D2R;
|
|
380
|
} else {
|
|
381
|
alert("*** Input out of range ***: lon: "+lon+" - lat: "+lat);
|
|
382
|
return null;
|
|
383
|
}
|
|
384
|
|
|
385
|
var con = Math.abs( Math.abs(lat) - HALF_PI);
|
|
386
|
var ts;
|
|
387
|
if (con > EPSLN) {
|
|
388
|
ts = tsfnz(this.e, lat, Math.sin(lat) );
|
|
389
|
rh1 = this.r_major * this.f0 * Math.pow(ts, this.ns);
|
|
390
|
} else {
|
|
391
|
con = lat * this.ns;
|
|
392
|
if (con <= 0) {
|
|
393
|
alert("Point can not be projected - ll2lcc");
|
|
394
|
return null;
|
|
395
|
}
|
|
396
|
rh1 = 0;
|
|
397
|
}
|
|
398
|
var theta = this.ns * adjust_lon(lon - this.center_lon);
|
|
399
|
var x = rh1 * Math.sin(theta) + this.false_easting;
|
|
400
|
var y = this.rh - rh1 * Math.cos(theta) + this.false_northing;
|
|
401
|
|
|
402
|
return new Array(x,y);
|
|
403
|
}
|
|
404
|
|
|
405
|
// Lambert Conformal Conic inverse equations--mapping x,y to lat/long
|
|
406
|
// -----------------------------------------------------------------
|
|
407
|
function lcc2ll(coords) {
|
|
408
|
|
|
409
|
var rh1, con, ts;
|
|
410
|
var lat, lon;
|
|
411
|
x = coords[0] - this.false_easting;
|
|
412
|
y = this.rh - coords[1] + this.false_northing;
|
|
413
|
if (this.ns > 0) {
|
|
414
|
rh1 = Math.sqrt (x * x + y * y);
|
|
415
|
con = 1.0;
|
|
416
|
} else {
|
|
417
|
rh1 = -Math.sqrt (x * x + y * y);
|
|
418
|
con = -1.0;
|
|
419
|
}
|
|
420
|
var theta = 0.0;
|
|
421
|
if (rh1 != 0) {
|
|
422
|
theta = Math.atan2((con * x),(con * y));
|
|
423
|
}
|
|
424
|
if ((rh1 != 0) || (this.ns > 0.0)) {
|
|
425
|
con = 1.0/this.ns;
|
|
426
|
ts = Math.pow((rh1/(this.r_major * this.f0)), con);
|
|
427
|
lat = phi2z(this.e, ts);
|
|
428
|
if (lat == -9999) return null;
|
|
429
|
} else {
|
|
430
|
lat = -HALF_PI;
|
|
431
|
}
|
|
432
|
lon = adjust_lon(theta/this.ns + this.center_lon);
|
|
433
|
return new Array(R2D*lon, R2D*lat);
|
|
434
|
}
|
|
435
|
|
|
436
|
//*******************************************************************************
|
|
437
|
//NAME POLAR STEREOGRAPHIC
|
|
438
|
// converted from the GCTPC package
|
|
439
|
//* Variables common to all subroutines in this code file
|
|
440
|
// static double r_major; /* major axis */
|
|
441
|
// static double r_minor; /* minor axis */
|
|
442
|
// static double e; /* eccentricity */
|
|
443
|
// static double e4; /* e4 calculated from eccentricity*/
|
|
444
|
// static double center_lon; /* center longitude */
|
|
445
|
// static double center_lat; /* center latitude */
|
|
446
|
// static double fac; /* sign variable */
|
|
447
|
// static double ind; /* flag variable */
|
|
448
|
// static double mcs; /* small m */
|
|
449
|
// static double tcs; /* small t */
|
|
450
|
// static double false_northing; /* y offset in meters */
|
|
451
|
// static double false_easting; /* x offset in meters */
|
|
452
|
|
|
453
|
|
|
454
|
//* Initialize the Polar Stereographic projection
|
|
455
|
function psinit(param) {
|
|
456
|
// array consisting of: r_maj,r_min,c_lon,c_lat,false_east,false_north)
|
|
457
|
//double c_lon; /* center longitude in degrees */
|
|
458
|
//double c_lat; /* center latitude in degrees */
|
|
459
|
//double r_maj; /* major axis */
|
|
460
|
//double r_min; /* minor axis */
|
|
461
|
//double false_east; /* x offset in meters */
|
|
462
|
//double false_north; /* y offset in meters */
|
|
463
|
|
|
464
|
this.r_major = param[0];
|
|
465
|
this.r_minor = param[1];
|
|
466
|
this.center_lon = param[2] * D2R;
|
|
467
|
this.center_lat = param[3] * D2R;
|
|
468
|
this.false_easting = param[4];
|
|
469
|
this.false_northing = param[5];
|
|
470
|
|
|
471
|
var temp = this.r_minor / this.r_major;
|
|
472
|
this.e = 1.0 - temp*temp;
|
|
473
|
this.e = Math.sqrt(this.e);
|
|
474
|
var con = 1.0 + this.e;
|
|
475
|
var com = 1.0 - this.e;
|
|
476
|
this.e4 = Math.sqrt( Math.pow(con,con) * Math.pow(com,com) );
|
|
477
|
this.fac = (this.center_lat < 0) ? -1.0 : 1.0;
|
|
478
|
this.ind = 0;
|
|
479
|
if (Math.abs(Math.abs(this.center_lat) - HALF_PI) > EPSLN) {
|
|
480
|
this.ind = 1;
|
|
481
|
var con1 = this.fac * this.center_lat;
|
|
482
|
var sinphi = Math.sin(con1);
|
|
483
|
this.mcs = msfnz(this.e, sinphi, Math.cos(con1));
|
|
484
|
this.tcs = tsfnz(this.e, con1, sinphi);
|
|
485
|
}
|
|
486
|
}
|
|
487
|
|
|
488
|
//* Polar Stereographic forward equations--mapping lat,long to x,y
|
|
489
|
// --------------------------------------------------------------*/
|
|
490
|
function ll2ps(coords) {
|
|
491
|
|
|
492
|
var lon = coords[0];
|
|
493
|
var lat = coords[1];
|
|
494
|
|
|
495
|
var con1 = this.fac * adjust_lon(lon - this.center_lon);
|
|
496
|
var con2 = this.fac * lat;
|
|
497
|
var sinphi = Math.sin(con2);
|
|
498
|
var ts = tsfnz(this.e, con2, sinphi);
|
|
499
|
if (this.ind != 0) {
|
|
500
|
rh = this.r_major * this.mcs * ts / this.tcs;
|
|
501
|
} else {
|
|
502
|
rh = 2.0 * this.r_major * ts / this.e4;
|
|
503
|
}
|
|
504
|
var x = this.fac * rh * Math.sin(con1) + this.false_easting;
|
|
505
|
var y = -this.fac * rh * Math.cos(con1) + this.false_northing;;
|
|
506
|
|
|
507
|
return new Array(x,y);
|
|
508
|
}
|
|
509
|
|
|
510
|
//* Polar Stereographic inverse equations--mapping x,y to lat/long
|
|
511
|
// --------------------------------------------------------------*/
|
|
512
|
function ps2ll(coords) {
|
|
513
|
|
|
514
|
x = (coords[0] - this.false_easting) * this.fac;
|
|
515
|
y = (coords[1] - this.false_northing) * this.fac;
|
|
516
|
var rh = Math.sqrt(x * x + y * y);
|
|
517
|
if (this.ind != 0) {
|
|
518
|
ts = rh * this.tcs/(this.r_major * this.mcs);
|
|
519
|
} else {
|
|
520
|
ts = rh * this.e4/(this.r_major * 2.0);
|
|
521
|
}
|
|
522
|
var lat = this.fac * phi2z(this.e, ts);
|
|
523
|
if (lat == -9999) return null;
|
|
524
|
var lon = 0;
|
|
525
|
if (rh == 0) {
|
|
526
|
lon = this.fac * this.center_lon;
|
|
527
|
} else {
|
|
528
|
lon = adjust_lon(this.fac * Math.atan2(x, -y) + this.center_lon);
|
|
529
|
}
|
|
530
|
return new Array(R2D*lon, R2D*lat);
|
|
531
|
}
|
|
532
|
|
|
533
|
// following constants added by Greenwood
|
|
534
|
function semi_minor(a, f) { return a-(a*(1/f)); }
|
|
535
|
|
|
536
|
var grs80 = new Array(6378137.0, 6356752.31414036); // r_maj, r_min
|
|
537
|
var wgs84 = new Array(6378137.0, 6356752.31424518);
|
|
538
|
var wgs72 = new Array(6378135.0, 6356750.52001609);
|
|
539
|
var intl = new Array(6378388.0, 6356911.94612795); // (f=297) from ESRI
|
|
540
|
|
|
541
|
var usfeet = 1200/3937; // US Survey foot
|
|
542
|
var feet = 0.3048; // International foot
|
|
543
|
|
|
544
|
// following functions from gctpc cproj.c for transverse mercator projections
|
|
545
|
function e0fn(x){return(1.0-0.25*x*(1.0+x/16.0*(3.0+1.25*x)));}
|
|
546
|
function e1fn(x){return(0.375*x*(1.0+0.25*x*(1.0+0.46875*x)));}
|
|
547
|
function e2fn(x){return(0.05859375*x*x*(1.0+0.75*x));}
|
|
548
|
function e3fn(x){return(x*x*x*(35.0/3072.0));}
|
|
549
|
function mlfn(e0,e1,e2,e3,phi){return(e0*phi-e1*Math.sin(2.0*phi)+e2*Math.sin(4.0*phi)-e3*Math.sin(6.0*phi));}
|
|
550
|
|
|
551
|
/**
|
|
552
|
Initialize Transverse Mercator projection
|
|
553
|
*/
|
|
554
|
function tminit(param) {
|
|
555
|
this.r_maj=param[0];
|
|
556
|
this.r_min=param[1];
|
|
557
|
this.scale_fact=param[2];
|
|
558
|
this.lon_center=param[3]*D2R;
|
|
559
|
this.lat_origin=param[4]*D2R;
|
|
560
|
this.false_easting=param[5];
|
|
561
|
this.false_northing=param[6];
|
|
562
|
var temp = this.r_min / this.r_maj;
|
|
563
|
this.es = 1.0 - Math.pow(temp,2);
|
|
564
|
// this.e = Math.sqrt(this.es);
|
|
565
|
this.e0 = e0fn(this.es);
|
|
566
|
this.e1 = e1fn(this.es);
|
|
567
|
this.e2 = e2fn(this.es);
|
|
568
|
this.e3 = e3fn(this.es);
|
|
569
|
this.ml0 = this.r_maj * mlfn(this.e0, this.e1, this.e2, this.e3, this.lat_origin);
|
|
570
|
this.esp = this.es / (1.0 - this.es);
|
|
571
|
this.ind = (this.es < .00001) ? 1 : 0;
|
|
572
|
}
|
|
573
|
|
|
574
|
/**
|
|
575
|
Initialize UTM projection
|
|
576
|
*/
|
|
577
|
function utminit(param) {
|
|
578
|
this.r_maj=param[0];
|
|
579
|
this.r_min=param[1];
|
|
580
|
this.scale_fact=param[2];
|
|
581
|
var zone=param[3];
|
|
582
|
this.lat_origin = 0.0;
|
|
583
|
this.lon_center = ((6 * Math.abs(zone)) - 183) * D2R;
|
|
584
|
this.false_easting = 500000.0;
|
|
585
|
this.false_northing = (zone < 0) ? 10000000.0 : 0.0;
|
|
586
|
var temp = this.r_min / this.r_maj;
|
|
587
|
this.es = 1.0 - Math.pow(temp,2);
|
|
588
|
// this.e = Math.sqrt(this.es);
|
|
589
|
this.e0 = e0fn(this.es);
|
|
590
|
this.e1 = e1fn(this.es);
|
|
591
|
this.e2 = e2fn(this.es);
|
|
592
|
this.e3 = e3fn(this.es);
|
|
593
|
this.ml0 = this.r_maj * mlfn(this.e0, this.e1, this.e2, this.e3, this.lat_origin);
|
|
594
|
this.esp = this.es / (1.0 - this.es);
|
|
595
|
this.ind = (this.es < .00001) ? 1 : 0;
|
|
596
|
} // utminit()
|
|
597
|
|
|
598
|
/**
|
|
599
|
Transverse Mercator Forward - long/lat to x/y
|
|
600
|
*/
|
|
601
|
function ll2tm(coords) {
|
|
602
|
var lon=coords[0]*D2R;
|
|
603
|
var lat=coords[1]*D2R;
|
|
604
|
var delta_lon = adjust_lon(lon - this.lon_center); /* Delta longitude (Given longitude - center */
|
|
605
|
var con; /* cone constant */
|
|
606
|
var x, y;
|
|
607
|
var sin_phi=Math.sin(lat);
|
|
608
|
var cos_phi=Math.cos(lat);
|
|
609
|
|
|
610
|
if (this.ind != 0) {
|
|
611
|
var b = cos_phi * Math.sin(delta_lon);
|
|
612
|
if ((Math.abs(Math.abs(b) - 1.0)) < .0000000001) {
|
|
613
|
alert("Error in ll2tm(): Point projects into infinity");
|
|
614
|
return(93);
|
|
615
|
} else {
|
|
616
|
x = .5 * this.r_maj * this.scale_fact * Math.log((1.0 + b)/(1.0 - b));
|
|
617
|
con = Math.acos(cos_phi * Math.cos(delta_lon)/Math.sqrt(1.0 - b*b));
|
|
618
|
if (lat < 0)
|
|
619
|
con = - con;
|
|
620
|
y = this.r_maj * this.scale_fact * (con - this.lat_origin);
|
|
621
|
}
|
|
622
|
} else {
|
|
623
|
var al = cos_phi * delta_lon;
|
|
624
|
var als = Math.pow(al,2);
|
|
625
|
var c = this.esp * Math.pow(cos_phi,2);
|
|
626
|
var tq = Math.tan(lat);
|
|
627
|
var t = Math.pow(tq,2);
|
|
628
|
con = 1.0 - this.es * Math.pow(sin_phi,2);
|
|
629
|
var n = this.r_maj / Math.sqrt(con);
|
|
630
|
var ml = this.r_maj * mlfn(this.e0, this.e1, this.e2, this.e3, lat);
|
|
631
|
|
|
632
|
x = this.scale_fact * n * al * (1.0 + als / 6.0 * (1.0 - t + c + als / 20.0 * (5.0 - 18.0 * t + Math.pow(t,2) + 72.0 * c - 58.0 * this.esp))) + this.false_easting;
|
|
633
|
y = this.scale_fact * (ml - this.ml0 + n * tq * (als * (0.5 + als / 24.0 * (5.0 - t + 9.0 * c + 4.0 * Math.pow(c,2) + als / 30.0 * (61.0 - 58.0 * t + Math.pow(t,2) + 600.0 * c - 330.0 * this.esp))))) + this.false_northing;
|
|
634
|
|
|
635
|
switch(this.units){
|
|
636
|
case "usfeet":
|
|
637
|
x /= usfeet;
|
|
638
|
y /= usfeet
|
|
639
|
break;
|
|
640
|
case "feet":
|
|
641
|
x = x / feet;
|
|
642
|
y = y / feet;
|
|
643
|
break;
|
|
644
|
} // switch()
|
|
645
|
}
|
|
646
|
return new Array(x,y);
|
|
647
|
} // ll2tm()
|
|
648
|
|
|
649
|
/**
|
|
650
|
Transverse Mercator Inverse - x/y to long/lat
|
|
651
|
*/
|
|
652
|
function tm2ll(coords) {
|
|
653
|
var x=coords[0];
|
|
654
|
var y=coords[1];
|
|
655
|
var con, phi; /* temporary angles */
|
|
656
|
var delta_phi; /* difference between longitudes */
|
|
657
|
var i;
|
|
658
|
var max_iter = 6; /* maximun number of iterations */
|
|
659
|
var lat, lon;
|
|
660
|
|
|
661
|
if (this.ind != 0) { /* spherical form */
|
|
662
|
var f = exp(x/(this.r_maj * this.scale_fact));
|
|
663
|
var g = .5 * (f - 1/f);
|
|
664
|
var temp = this.lat_origin + y/(this.r_maj * this.scale_fact);
|
|
665
|
var h = cos(temp);
|
|
666
|
con = sqrt((1.0 - h * h)/(1.0 + g * g));
|
|
667
|
lat = asinz(con);
|
|
668
|
if (temp < 0)
|
|
669
|
lat = -lat;
|
|
670
|
if ((g == 0) && (h == 0)) {
|
|
671
|
lon = this.lon_center;
|
|
672
|
} else {
|
|
673
|
lon = adjust_lon(atan2(g,h) + this.lon_center);
|
|
674
|
}
|
|
675
|
} else { // ellipsoidal form
|
|
676
|
x = x - this.false_easting;
|
|
677
|
y = y - this.false_northing;
|
|
678
|
|
|
679
|
con = (this.ml0 + y / this.scale_fact) / this.r_maj;
|
|
680
|
phi = con;
|
|
681
|
for (i=0;;i++) {
|
|
682
|
delta_phi=((con + this.e1 * Math.sin(2.0*phi) - this.e2 * Math.sin(4.0*phi) + this.e3 * Math.sin(6.0*phi)) / this.e0) - phi;
|
|
683
|
phi += delta_phi;
|
|
684
|
if (Math.abs(delta_phi) <= EPSLN) break;
|
|
685
|
if (i >= max_iter) {
|
|
686
|
alert ("Error in tm2ll(): Latitude failed to converge");
|
|
687
|
return(95);
|
|
688
|
}
|
|
689
|
} // for()
|
|
690
|
if (Math.abs(phi) < HALF_PI) {
|
|
691
|
// sincos(phi, &sin_phi, &cos_phi);
|
|
692
|
var sin_phi=Math.sin(phi);
|
|
693
|
var cos_phi=Math.cos(phi);
|
|
694
|
var tan_phi = Math.tan(phi);
|
|
695
|
var c = this.esp * Math.pow(cos_phi,2);
|
|
696
|
var cs = Math.pow(c,2);
|
|
697
|
var t = Math.pow(tan_phi,2);
|
|
698
|
var ts = Math.pow(t,2);
|
|
699
|
con = 1.0 - this.es * Math.pow(sin_phi,2);
|
|
700
|
var n = this.r_maj / Math.sqrt(con);
|
|
701
|
var r = n * (1.0 - this.es) / con;
|
|
702
|
var d = x / (n * this.scale_fact);
|
|
703
|
var ds = Math.pow(d,2);
|
|
704
|
lat = phi - (n * tan_phi * ds / r) * (0.5 - ds / 24.0 * (5.0 + 3.0 * t + 10.0 * c - 4.0 * cs - 9.0 * this.esp - ds / 30.0 * (61.0 + 90.0 * t + 298.0 * c + 45.0 * ts - 252.0 * this.esp - 3.0 * cs)));
|
|
705
|
lon = adjust_lon(this.lon_center + (d * (1.0 - ds / 6.0 * (1.0 + 2.0 * t + c - ds / 20.0 * (5.0 - 2.0 * c + 28.0 * t - 3.0 * cs + 8.0 * this.esp + 24.0 * ts))) / cos_phi));
|
|
706
|
} else {
|
|
707
|
lat = HALF_PI * sign(y);
|
|
708
|
lon = this.lon_center;
|
|
709
|
}
|
|
710
|
}
|
|
711
|
return new Array(lon*R2D, lat*R2D);
|
|
712
|
} // tm2ll()
|
|
713
|
|
|
714
|
|
|
715
|
// Function to compute the constant small m which is the radius of
|
|
716
|
// a parallel of latitude, phi, divided by the semimajor axis.
|
|
717
|
// -----------------------------------------------------------------
|
|
718
|
function msfnz(eccent, sinphi, cosphi) {
|
|
719
|
var con = eccent * sinphi;
|
|
720
|
return cosphi/(Math.sqrt(1.0 - con * con));
|
|
721
|
}
|
|
722
|
|
|
723
|
// Function to compute the constant small t for use in the forward
|
|
724
|
// computations in the Lambert Conformal Conic and the Polar
|
|
725
|
// Stereographic projections.
|
|
726
|
// -----------------------------------------------------------------
|
|
727
|
function tsfnz(eccent, phi, sinphi) {
|
|
728
|
var con = eccent * sinphi;
|
|
729
|
var com = .5 * eccent;
|
|
730
|
con = Math.pow(((1.0 - con) / (1.0 + con)), com);
|
|
731
|
return (Math.tan(.5 * (HALF_PI - phi))/con);
|
|
732
|
}
|
|
733
|
|
|
734
|
|
|
735
|
// Function to compute the latitude angle, phi2, for the inverse of the
|
|
736
|
// Lambert Conformal Conic and Polar Stereographic projections.
|
|
737
|
// ----------------------------------------------------------------
|
|
738
|
function phi2z(eccent, ts) {
|
|
739
|
var eccnth = .5 * eccent;
|
|
740
|
var con, dphi;
|
|
741
|
var phi = HALF_PI - 2 * Math.atan(ts);
|
|
742
|
for (i = 0; i <= 15; i++) {
|
|
743
|
con = eccent * Math.sin(phi);
|
|
744
|
dphi = HALF_PI - 2 * Math.atan(ts *(Math.pow(((1.0 - con)/(1.0 + con)),eccnth))) - phi;
|
|
745
|
phi += dphi;
|
|
746
|
if (Math.abs(dphi) <= .0000000001) return phi;
|
|
747
|
}
|
|
748
|
alert("Convergence error - phi2z");
|
|
749
|
return -9999;
|
|
750
|
}
|
|
751
|
|
|
752
|
// Function to return the sign of an argument
|
|
753
|
function sign(x) { if (x < 0.0) return(-1); else return(1);}
|
|
754
|
|
|
755
|
// Function to adjust longitude to -180 to 180; input in radians
|
|
756
|
function adjust_lon(x) {x=(Math.abs(x)<PI)?x:(x-(sign(x)*TWO_PI));return(x);}
|
|
757
|
|