1/* 2 proj4js.js -- Javascript reprojection library. 3 4 Authors: Mike Adair madairATdmsolutions.ca 5 Richard Greenwood richATgreenwoodmap.com 6 Didier Richard didier.richardATign.fr 7 Stephen Irons stephen.ironsATclear.net.nz 8 Olivier Terral oterralATgmail.com 9 10 License: 11 Copyright (c) 2012, Mike Adair, Richard Greenwood, Didier Richard, 12 Stephen Irons and Olivier Terral 13 14 Permission is hereby granted, free of charge, to any person obtaining a 15 copy of this software and associated documentation files (the "Software"), 16 to deal in the Software without restriction, including without limitation 17 the rights to use, copy, modify, merge, publish, distribute, sublicense, 18 and/or sell copies of the Software, and to permit persons to whom the 19 Software is furnished to do so, subject to the following conditions: 20 21 The above copyright notice and this permission notice shall be included 22 in all copies or substantial portions of the Software. 23 24 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS 25 OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, 26 FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL 27 THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER 28 LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING 29 FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER 30 DEALINGS IN THE SOFTWARE. 31 32 Note: This program is an almost direct port of the C library PROJ.4. 33*/ 34/* ====================================================================== 35 proj4js.js 36 ====================================================================== */ 37 38/* 39Author: Mike Adair madairATdmsolutions.ca 40 Richard Greenwood [email protected] 41License: LGPL as per: http://www.gnu.org/copyleft/lesser.html 42 43$Id: Proj.js 2956 2007-07-09 12:17:52Z steven $ 44*/ 45 46/** 47 * Namespace: Proj4js 48 * 49 * Proj4js is a JavaScript library to transform point coordinates from one 50 * coordinate system to another, including datum transformations. 51 * 52 * This library is a port of both the Proj.4 and GCTCP C libraries to JavaScript. 53 * Enabling these transformations in the browser allows geographic data stored 54 * in different projections to be combined in browser-based web mapping 55 * applications. 56 * 57 * Proj4js must have access to coordinate system initialization strings (which 58 * are the same as for PROJ.4 command line). Thes can be included in your 59 * application using a <script> tag or Proj4js can load CS initialization 60 * strings from a local directory or a web service such as spatialreference.org. 61 * 62 * Similarly, Proj4js must have access to projection transform code. These can 63 * be included individually using a <script> tag in your page, built into a 64 * custom build of Proj4js or loaded dynamically at run-time. Using the 65 * -combined and -compressed versions of Proj4js includes all projection class 66 * code by default. 67 * 68 * Note that dynamic loading of defs and code happens ascynchrously, check the 69 * Proj.readyToUse flag before using the Proj object. If the defs and code 70 * required by your application are loaded through script tags, dynamic loading 71 * is not required and the Proj object will be readyToUse on return from the 72 * constructor. 73 * 74 * All coordinates are handled as points which have a .x and a .y property 75 * which will be modified in place. 76 * 77 * Override Proj4js.reportError for output of alerts and warnings. 78 * 79 * See http://trac.osgeo.org/proj4js/wiki/UserGuide for full details. 80*/ 81 82/** 83 * Global namespace object for Proj4js library 84 */ 85var Proj4js = { 86 87 /** 88 * Property: defaultDatum 89 * The datum to use when no others a specified 90 */ 91 defaultDatum: 'WGS84', //default datum 92 93 /** 94 * Method: transform(source, dest, point) 95 * Transform a point coordinate from one map projection to another. This is 96 * really the only public method you should need to use. 97 * 98 * Parameters: 99 * source - {Proj4js.Proj} source map projection for the transformation 100 * dest - {Proj4js.Proj} destination map projection for the transformation 101 * point - {Object} point to transform, may be geodetic (long, lat) or 102 * projected Cartesian (x,y), but should always have x,y properties. 103 */ 104 transform: function(source, dest, point) { 105 if (!source.readyToUse) { 106 this.reportError("Proj4js initialization for:"+source.srsCode+" not yet complete"); 107 return point; 108 } 109 if (!dest.readyToUse) { 110 this.reportError("Proj4js initialization for:"+dest.srsCode+" not yet complete"); 111 return point; 112 } 113 114 // Workaround for datum shifts towgs84, if either source or destination projection is not wgs84 115 if (source.datum && dest.datum && ( 116 ((source.datum.datum_type == Proj4js.common.PJD_3PARAM || source.datum.datum_type == Proj4js.common.PJD_7PARAM) && dest.datumCode != "WGS84") || 117 ((dest.datum.datum_type == Proj4js.common.PJD_3PARAM || dest.datum.datum_type == Proj4js.common.PJD_7PARAM) && source.datumCode != "WGS84"))) { 118 var wgs84 = Proj4js.WGS84; 119 this.transform(source, wgs84, point); 120 source = wgs84; 121 } 122 123 // DGR, 2010/11/12 124 if (source.axis!="enu") { 125 this.adjust_axis(source,false,point); 126 } 127 128 // Transform source points to long/lat, if they aren't already. 129 if ( source.projName=="longlat") { 130 point.x *= Proj4js.common.D2R; // convert degrees to radians 131 point.y *= Proj4js.common.D2R; 132 } else { 133 if (source.to_meter) { 134 point.x *= source.to_meter; 135 point.y *= source.to_meter; 136 } 137 source.inverse(point); // Convert Cartesian to longlat 138 } 139 140 // Adjust for the prime meridian if necessary 141 if (source.from_greenwich) { 142 point.x += source.from_greenwich; 143 } 144 145 // Convert datums if needed, and if possible. 146 point = this.datum_transform( source.datum, dest.datum, point ); 147 148 // Adjust for the prime meridian if necessary 149 if (dest.from_greenwich) { 150 point.x -= dest.from_greenwich; 151 } 152 153 if( dest.projName=="longlat" ) { 154 // convert radians to decimal degrees 155 point.x *= Proj4js.common.R2D; 156 point.y *= Proj4js.common.R2D; 157 } else { // else project 158 dest.forward(point); 159 if (dest.to_meter) { 160 point.x /= dest.to_meter; 161 point.y /= dest.to_meter; 162 } 163 } 164 165 // DGR, 2010/11/12 166 if (dest.axis!="enu") { 167 this.adjust_axis(dest,true,point); 168 } 169 170 return point; 171 }, // transform() 172 173 /** datum_transform() 174 source coordinate system definition, 175 destination coordinate system definition, 176 point to transform in geodetic coordinates (long, lat, height) 177 */ 178 datum_transform : function( source, dest, point ) { 179 180 // Short cut if the datums are identical.
181 if( source.compare_datums( dest ) ) { 182 return point; // in this case, zero is sucess, 183 // whereas cs_compare_datums returns 1 to indicate TRUE 184 // confusing, should fix this 185 } 186 187 // Explicitly skip datum transform by setting 'datum=none' as parameter for either source or dest 188 if( source.datum_type == Proj4js.common.PJD_NODATUM 189 || dest.datum_type == Proj4js.common.PJD_NODATUM) { 190 return point; 191 } 192 193 // Do we need to go through geocentric coordinates? 194 if( source.es != dest.es || source.a != dest.a 195 || source.datum_type == Proj4js.common.PJD_3PARAM 196 || source.datum_type == Proj4js.common.PJD_7PARAM 197 || dest.datum_type == Proj4js.common.PJD_3PARAM 198 || dest.datum_type == Proj4js.common.PJD_7PARAM) 199 { 200 201 // Convert to geocentric coordinates. 202 source.geodetic_to_geocentric( point ); 203 // CHECK_RETURN; 204 205 // Convert between datums 206 if( source.datum_type == Proj4js.common.PJD_3PARAM || source.datum_type == Proj4js.common.PJD_7PARAM ) { 207 source.geocentric_to_wgs84(point); 208 // CHECK_RETURN; 209 } 210 211 if( dest.datum_type == Proj4js.common.PJD_3PARAM || dest.datum_type == Proj4js.common.PJD_7PARAM ) { 212 dest.geocentric_from_wgs84(point); 213 // CHECK_RETURN; 214 } 215 216 // Convert back to geodetic coordinates 217 dest.geocentric_to_geodetic( point ); 218 // CHECK_RETURN; 219 } 220 221 return point; 222 }, // cs_datum_transform 223 224 /** 225 * Function: adjust_axis 226 * Normalize or de-normalized the x/y/z axes. The normal form is "enu" 227 * (easting, northing, up). 228 * Parameters: 229 * crs {Proj4js.Proj} the coordinate reference system 230 * denorm {Boolean} when false, normalize 231 * point {Object} the coordinates to adjust 232 */ 233 adjust_axis: function(crs, denorm, point) { 234 var xin= point.x, yin= point.y, zin= point.z || 0.0; 235 var v, t; 236 for (var i= 0; i<3; i++) { 237 if (denorm && i==2 && point.z===undefined) { continue; } 238 if (i==0) { v= xin; t= 'x'; } 239 else if (i==1) { v= yin; t= 'y'; } 240 else { v= zin; t= 'z'; } 241 switch(crs.axis[i]) { 242 case 'e': 243 point[t]= v; 244 break; 245 case 'w': 246 point[t]= -v; 247 break; 248 case 'n': 249 point[t]= v; 250 break; 251 case 's': 252 point[t]= -v; 253 break; 254 case 'u': 255 if (point[t]!==undefined) { point.z= v; } 256 break; 257 case 'd': 258 if (point[t]!==undefined) { point.z= -v; } 259 break; 260 default : 261 alert("ERROR: unknow axis ("+crs.axis[i]+") - check definition of "+crs.projName); 262 return null; 263 } 264 } 265 return point; 266 }, 267 268 /** 269 * Function: reportError 270 * An internal method to report errors back to user. 271 * Override this in applications to report error messages or throw exceptions. 272 */ 273 reportError: function(msg) { 274 //console.log(msg); 275 }, 276 277/** 278 * 279 * Title: Private Methods 280 * The following properties and methods are intended for internal use only. 281 * 282 * This is a minimal implementation of JavaScript inheritance methods so that 283 * Proj4js can be used as a stand-alone library. 284 * These are copies of the equivalent OpenLayers methods at v2.7 285 */ 286 287/** 288 * Function: extend 289 * Copy all properties of a source object to a destination object. Modifies 290 * the passed in destination object. Any properties on the source object 291 * that are set to undefined will not be (re)set on the destination object. 292 * 293 * Parameters: 294 * destination - {Object} The object that will be modified 295 * source - {Object} The object with properties to be set on the destination 296 * 297 * Returns: 298 * {Object} The destination object. 299 */ 300 extend: function(destination, source) { 301 destination = destination || {}; 302 if(source) { 303 for(var property in source) { 304 var value = source[property]; 305 if(value !== undefined) { 306 destination[property] = value; 307 } 308 } 309 } 310 return destination; 311 }, 312 313/** 314 * Constructor: Class 315 * Base class used to construct all other classes. Includes support for 316 * multiple inheritance. 317 * 318 */ 319 Class: function() { 320 var Class = function() { 321 this.initialize.apply(this, arguments); 322 }; 323 324 var extended = {}; 325 var parent; 326 for(var i=0; i<arguments.length; ++i) { 327 if(typeof arguments[i] == "function") { 328 // get the prototype of the superclass 329 parent = arguments[i].prototype; 330 } else { 331 // in this case we're extending with the prototype 332 parent = arguments[i]; 333 } 334 Proj4js.extend(extended, parent); 335 } 336 Class.prototype = extended; 337 338 return Class; 339 }, 340 341 /** 342 * Function: bind 343 * Bind a function to an object. Method to easily create closures with 344 * 'this' altered. 345 * 346 * Parameters: 347 * func - {Function} Input function. 348 * object - {Object} The object to bind to the input function (as this). 349 * 350 * Returns: 351 * {Function} A closure with 'this' set to the passed in object. 352 */ 353 bind: function(func, object) { 354 // create a reference to all arguments past the second one 355 var args = Array.prototype.slice.apply(arguments, [2]); 356 return function() { 357 // Push on any additional arguments from the actual function call.
358 // These will come after those sent to the bind call. 359 var newArgs = args.concat( 360 Array.prototype.slice.apply(arguments, [0]) 361 ); 362 return func.apply(object, newArgs); 363 }; 364 }, 365 366/** 367 * The following properties and methods handle dynamic loading of JSON objects. 368 */ 369 370 /** 371 * Property: scriptName 372 * {String} The filename of this script without any path. 373 */ 374 scriptName: "proj4js-combined.js", 375 376 /** 377 * Property: defsLookupService 378 * AJAX service to retreive projection definition parameters from 379 */ 380 defsLookupService: 'http://spatialreference.org/ref', 381 382 /** 383 * Property: libPath 384 * internal: http server path to library code. 385 */ 386 libPath: null, 387 388 /** 389 * Function: getScriptLocation 390 * Return the path to this script. 391 * 392 * Returns: 393 * Path to this script 394 */ 395 getScriptLocation: function () { 396 if (this.libPath) return this.libPath; 397 var scriptName = this.scriptName; 398 var scriptNameLen = scriptName.length; 399 400 var scripts = document.getElementsByTagName('script'); 401 for (var i = 0; i < scripts.length; i++) { 402 var src = scripts[i].getAttribute('src'); 403 if (src) { 404 var index = src.lastIndexOf(scriptName); 405 // is it found, at the end of the URL? 406 if ((index > -1) && (index + scriptNameLen == src.length)) { 407 this.libPath = src.slice(0, -scriptNameLen); 408 break; 409 } 410 } 411 } 412 return this.libPath||""; 413 }, 414 415 /** 416 * Function: loadScript 417 * Load a JS file from a URL into a <script> tag in the page. 418 * 419 * Parameters: 420 * url - {String} The URL containing the script to load 421 * onload - {Function} A method to be executed when the script loads successfully 422 * onfail - {Function} A method to be executed when there is an error loading the script 423 * loadCheck - {Function} A boolean method that checks to see if the script 424 * has loaded. Typically this just checks for the existance of 425 * an object in the file just loaded. 426 */ 427 loadScript: function(url, onload, onfail, loadCheck) { 428 var script = document.createElement('script'); 429 script.defer = false; 430 script.type = "text/javascript"; 431 script.id = url; 432 script.src = url; 433 script.onload = onload; 434 script.onerror = onfail; 435 script.loadCheck = loadCheck; 436 if (/MSIE/.test(navigator.userAgent)) { 437 script.onreadystatechange = this.checkReadyState; 438 } 439 document.getElementsByTagName('head')[0].appendChild(script); 440 }, 441 442 /** 443 * Function: checkReadyState 444 * IE workaround since there is no onerror handler. Calls the user defined 445 * loadCheck method to determine if the script is loaded. 446 * 447 */ 448 checkReadyState: function() { 449 if (this.readyState == 'loaded') { 450 if (!this.loadCheck()) { 451 this.onerror(); 452 } else { 453 this.onload(); 454 } 455 } 456 } 457}; 458 459/** 460 * Class: Proj4js.Proj 461 * 462 * Proj objects provide transformation methods for point coordinates 463 * between geodetic latitude/longitude and a projected coordinate system. 464 * once they have been initialized with a projection code. 465 * 466 * Initialization of Proj objects is with a projection code, usually EPSG codes, 467 * which is the key that will be used with the Proj4js.defs array. 468 * 469 * The code passed in will be stripped of colons and converted to uppercase 470 * to locate projection definition files. 471 * 472 * A projection object has properties for units and title strings. 473 */ 474Proj4js.Proj = Proj4js.Class({ 475 476 /** 477 * Property: readyToUse 478 * Flag to indicate if initialization is complete for this Proj object 479 */ 480 readyToUse: false, 481 482 /** 483 * Property: title
484 * The title to describe the projection 485 */ 486 title: null, 487 488 /** 489 * Property: projName 490 * The projection class for this projection, e.g. lcc (lambert conformal conic, 491 * or merc for mercator). These are exactly equivalent to their Proj4 492 * counterparts. 493 */ 494 projName: null, 495 /** 496 * Property: units 497 * The units of the projection. Values include 'm' and 'degrees' 498 */ 499 units: null, 500 /** 501 * Property: datum 502 * The datum specified for the projection 503 */ 504 datum: null, 505 /** 506 * Property: x0 507 * The x coordinate origin 508 */ 509 x0: 0, 510 /** 511 * Property: y0 512 * The y coordinate origin 513 */ 514 y0: 0, 515 /** 516 * Property: localCS 517 * Flag to indicate if the projection is a local one in which no transforms 518 * are required. 519 */ 520 localCS: false, 521 522 /** 523 * Property: queue 524 * Buffer (FIFO) to hold callbacks waiting to be called when projection loaded. 525 */ 526 queue: null, 527 528 /** 529 * Constructor: initialize 530 * Constructor for Proj4js.Proj objects 531 * 532 * Parameters: 533 * srsCode - a code for map projection definition parameters. These are usually 534 * (but not always) EPSG codes. 535 */ 536 initialize: function(srsCode, callback) { 537 this.srsCodeInput = srsCode; 538 539 //Register callbacks prior to attempting to process definition 540 this.queue = []; 541 if( callback ){ 542 this.queue.push( callback ); 543 } 544 545 //check to see if this is a WKT string 546 if ((srsCode.indexOf('GEOGCS') >= 0) || 547 (srsCode.indexOf('GEOCCS') >= 0) || 548 (srsCode.indexOf('PROJCS') >= 0) || 549 (srsCode.indexOf('LOCAL_CS') >= 0)) { 550 this.parseWKT(srsCode); 551 this.deriveConstants(); 552 this.loadProjCode(this.projName); 553 return; 554 } 555 556 // DGR 2008-08-03 : support urn and url 557 if (srsCode.indexOf('urn:') == 0) { 558 //urn:ORIGINATOR:def:crs:CODESPACE:VERSION:ID 559 var urn = srsCode.split(':'); 560 if ((urn[1] == 'ogc' || urn[1] =='x-ogc') && 561 (urn[2] =='def') && 562 (urn[3] =='crs')) { 563 srsCode = urn[4]+':'+urn[urn.length-1]; 564 } 565 } else if (srsCode.indexOf('http://') == 0) { 566 //url#ID 567 var url = srsCode.split('#'); 568 if (url[0].match(/epsg.org/)) { 569 // http://www.epsg.org/# 570 srsCode = 'EPSG:'+url[1]; 571 } else if (url[0].match(/RIG.xml/)) { 572 //http://librairies.ign.fr/geoportail/resources/RIG.xml# 573 //http://interop.ign.fr/registers/ign/RIG.xml# 574 srsCode = 'IGNF:'+url[1]; 575 } 576 } 577 this.srsCode = srsCode.toUpperCase(); 578 if (this.srsCode.indexOf("EPSG") == 0) { 579 this.srsCode = this.srsCode; 580 this.srsAuth = 'epsg'; 581 this.srsProjNumber = this.srsCode.substring(5); 582 // DGR 2007-11-20 : authority IGNF 583 } else if (this.srsCode.indexOf("IGNF") == 0) { 584 this.srsCode = this.srsCode; 585 this.srsAuth = 'IGNF'; 586 this.srsProjNumber = this.srsCode.substring(5); 587 // DGR 2008-06-19 : pseudo-authority CRS for WMS 588 } else if (this.srsCode.indexOf("CRS") == 0) { 589 this.srsCode = this.srsCode; 590 this.srsAuth = 'CRS'; 591 this.srsProjNumber = this.srsCode.substring(4); 592 } else { 593 this.srsAuth = ''; 594 this.srsProjNumber = this.srsCode; 595 } 596 597 this.loadProjDefinition(); 598 }, 599 600/** 601 * Function: loadProjDefinition 602 * Loads the coordinate system initialization string if required. 603 * Note that dynamic loading happens asynchronously so an application must 604 * wait for the readyToUse property is set to true. 605 * To prevent dynamic loading, include the defs through a script tag in 606 * your application. 607 * 608 */ 609 loadProjDefinition: function() { 610 //check in memory 611 if (Proj4js.defs[this.srsCode]) { 612 this.defsLoaded(); 613 return; 614 } 615 616 //else check for def on the server 617 var url = Proj4js.getScriptLocation() + 'defs/' + this.srsAuth.toUpperCase() + this.srsProjNumber + '.js'; 618 Proj4js.loadScript(url, 619 Proj4js.bind(this.defsLoaded, this), 620 Proj4js.bind(this.loadFromService, this), 621 Proj4js.bind(this.checkDefsLoaded, this) ); 622 }, 623 624/**
625 * Function: loadFromService 626 * Creates the REST URL for loading the definition from a web service and 627 * loads it. 628 * 629 */ 630 loadFromService: function() { 631 //else load from web service 632 var url = Proj4js.defsLookupService +'/' + this.srsAuth +'/'+ this.srsProjNumber + '/proj4js/'; 633 Proj4js.loadScript(url, 634 Proj4js.bind(this.defsLoaded, this), 635 Proj4js.bind(this.defsFailed, this), 636 Proj4js.bind(this.checkDefsLoaded, this) ); 637 }, 638 639/** 640 * Function: defsLoaded 641 * Continues the Proj object initilization once the def file is loaded 642 * 643 */ 644 defsLoaded: function() { 645 this.parseDefs(); 646 this.loadProjCode(this.projName); 647 }, 648 649/** 650 * Function: checkDefsLoaded 651 * This is the loadCheck method to see if the def object exists 652 * 653 */ 654 checkDefsLoaded: function() { 655 if (Proj4js.defs[this.srsCode]) { 656 return true; 657 } else { 658 return false; 659 } 660 }, 661 662 /** 663 * Function: defsFailed 664 * Report an error in loading the defs file, but continue on using WGS84 665 * 666 */ 667 defsFailed: function() { 668 Proj4js.reportError('failed to load projection definition for: '+this.srsCode); 669 Proj4js.defs[this.srsCode] = Proj4js.defs['WGS84']; //set it to something so it can at least continue 670 this.defsLoaded(); 671 }, 672 673/** 674 * Function: loadProjCode 675 * Loads projection class code dynamically if required. 676 * Projection code may be included either through a script tag or in 677 * a built version of proj4js 678 * 679 */ 680 loadProjCode: function(projName) { 681 if (Proj4js.Proj[projName]) { 682 this.initTransforms(); 683 return; 684 } 685 686 //the URL for the projection code 687 var url = Proj4js.getScriptLocation() + 'projCode/' + projName + '.js'; 688 Proj4js.loadScript(url, 689 Proj4js.bind(this.loadProjCodeSuccess, this, projName), 690 Proj4js.bind(this.loadProjCodeFailure, this, projName), 691 Proj4js.bind(this.checkCodeLoaded, this, projName) ); 692 }, 693 694 /** 695 * Function: loadProjCodeSuccess 696 * Loads any proj dependencies or continue on to final initialization. 697 * 698 */ 699 loadProjCodeSuccess: function(projName) { 700 if (Proj4js.Proj[projName].dependsOn){ 701 this.loadProjCode(Proj4js.Proj[projName].dependsOn); 702 } else { 703 this.initTransforms(); 704 } 705 }, 706 707 /** 708 * Function: defsFailed 709 * Report an error in loading the proj file. Initialization of the Proj 710 * object has failed and the readyToUse flag will never be set. 711 * 712 */ 713 loadProjCodeFailure: function(projName) { 714 Proj4js.reportError("failed to find projection file for: " + projName); 715 //TBD initialize with identity transforms so proj will still work? 716 }, 717 718/** 719 * Function: checkCodeLoaded 720 * This is the loadCheck method to see if the projection code is loaded 721 * 722 */ 723 checkCodeLoaded: function(projName) { 724 if (Proj4js.Proj[projName]) { 725 return true; 726 } else { 727 return false; 728 } 729 }, 730 731/** 732 * Function: initTransforms 733 * Finalize the initialization of the Proj object 734 * 735 */ 736 initTransforms: function() { 737 Proj4js.extend(this, Proj4js.Proj[this.projName]); 738 this.init(); 739 this.readyToUse = true; 740 if( this.queue ) { 741 var item; 742 while( (item = this.queue.shift()) ) { 743 item.call( this, this ); 744 } 745 } 746 }, 747 748/** 749 * Function: parseWKT 750 * Parses a WKT string to get initialization parameters 751 * 752 */ 753 wktRE: /^(\w+)\[(.*)\]$/, 754 parseWKT: function(wkt) { 755 var wktMatch = wkt.match(this.wktRE); 756 if (!wktMatch) return; 757 var wktObject = wktMatch[1]; 758 var wktContent = wktMatch[2]; 759 var wktTemp = wktContent.split(","); 760 var wktName; 761 if (wktObject.toUpperCase() == "TOWGS84") { 762 wktName = wktObject; //no name supplied for the TOWGS84 array 763 } else { 764 wktName = wktTemp.shift(); 765 } 766 wktName = wktName.replace(/^\"/,""); 767 wktName = wktName.replace(/\"$/,""); 768 769 /* 770 wktContent = wktTemp.join(","); 771 var wktArray = wktContent.split("],"); 772 for (var i=0; i<wktArray.length-1; ++i) { 773 wktArray[i] += "]"; 774 } 775 */ 776 777 var wktArray = new Array(); 778 var bkCount = 0; 779 var obj = ""; 780 for (var i=0; i<wktTemp.length; ++i) { 781 var token = wktTemp[i]; 782 for (var j=0; j<token.length; ++j) { 783 if (token.charAt(j) == "[") ++bkCount; 784 if (token.charAt(j) == "]") --bkCount; 785 } 786 obj += token; 787 if (bkCount === 0) { 788 wktArray.push(obj); 789 obj = ""; 790 } else { 791 obj += ","; 792 } 793 } 794 795 //do something based on the type of the wktObject being parsed 796 //add in variations in the spelling as required 797 switch (wktObject) { 798 case 'LOCAL_CS': 799 this.projName = 'identity' 800 this.localCS = true; 801 this.srsCode = wktName; 802 break; 803 case 'GEOGCS': 804 this.projName = 'longlat' 805 this.geocsCode = wktName; 806 if (!this.srsCode) this.srsCode = wktName; 807 break; 808 case 'PROJCS': 809 this.srsCode = wktName; 810 break; 811 case 'GEOCCS': 812 break; 813 case 'PROJECTION': 814 this.projName = Proj4js.wktProjections[wktName] 815 break; 816 case 'DATUM': 817 this.datumName = wktName; 818 break; 819 case 'LOCAL_DATUM': 820 this.datumCode = 'none'; 821 break; 822 case 'SPHEROID': 823 this.ellps = wktName; 824 this.a = parseFloat(wktArray.shift()); 825 this.rf = parseFloat(wktArray.shift()); 826 break; 827 case 'PRIMEM': 828 this.from_greenwich = parseFloat(wktArray.shift()); //to radians? 829 break; 830 case 'UNIT': 831 this.units = wktName; 832 this.unitsPerMeter = parseFloat(wktArray.shift()); 833 break; 834 case 'PARAMETER': 835 var name = wktName.toLowerCase(); 836 var value = parseFloat(wktArray.shift()); 837 //there may be many variations on the wktName values, add in case 838 //statements as required 839 switch (name) { 840 case 'false_easting': 841 this.x0 = value; 842 break; 843 case 'false_northing': 844 this.y0 = value; 845 break; 846 case 'scale_factor': 847 this.k0 = value; 848 break; 849 case 'central_meridian': 850 this.long0 = value*Proj4js.common.D2R; 851 break; 852 case 'latitude_of_origin': 853 this.lat0 = value*Proj4js.common.D2R; 854 break; 855 case 'more_here': 856 break; 857 default: 858 break; 859 } 860 break; 861 case 'TOWGS84': 862 this.datum_params = wktArray; 863 break; 864 //DGR 2010-11-12: AXIS 865 case 'AXIS': 866 var name= wktName.toLowerCase(); 867 var value= wktArray.shift(); 868 switch (value) { 869 case 'EAST' : value= 'e'; break; 870 case 'WEST' : value= 'w'; break; 871 case 'NORTH': value= 'n'; break; 872 case 'SOUTH': value= 's'; break; 873 case 'UP' : value= 'u'; break; 874 case 'DOWN' : value= 'd'; break; 875 case 'OTHER': 876 default : value= ' '; break;//FIXME 877 } 878 if (!this.axis) { this.axis= "enu"; } 879 switch(name) { 880 case 'x': this.axis= value + this.axis.substr(1,2); break; 881 case 'y': this.axis= this.axis.substr(0,1) + value + this.axis.substr(2,1); break; 882 case 'z': this.axis= this.axis.substr(0,2) + value ; break; 883 default : break; 884 } 885 case 'MORE_HERE': 886 break; 887 default: 888 break; 889 } 890 for (var i=0; i<wktArray.length; ++i) { 891 this.parseWKT(wktArray[i]); 892 } 893 }, 894 895/**
896 * Function: parseDefs 897 * Parses the PROJ.4 initialization string and sets the associated properties. 898 * 899 */ 900 parseDefs: function() { 901 this.defData = Proj4js.defs[this.srsCode]; 902 var paramName, paramVal; 903 if (!this.defData) { 904 return; 905 } 906 var paramArray=this.defData.split("+"); 907 908 for (var prop=0; prop<paramArray.length; prop++) { 909 var property = paramArray[prop].split("="); 910 paramName = property[0].toLowerCase(); 911 paramVal = property[1]; 912 913 switch (paramName.replace(/\s/gi,"")) { // trim out spaces 914 case "": break; // throw away nameless parameter 915 case "title": this.title = paramVal; break; 916 case "proj": this.projName = paramVal.replace(/\s/gi,""); break; 917 case "units": this.units = paramVal.replace(/\s/gi,""); break; 918 case "datum": this.datumCode = paramVal.replace(/\s/gi,""); break; 919 case "nadgrids": this.nagrids = paramVal.replace(/\s/gi,""); break; 920 case "ellps": this.ellps = paramVal.replace(/\s/gi,""); break; 921 case "a": this.a = parseFloat(paramVal); break; // semi-major radius 922 case "b": this.b = parseFloat(paramVal); break; // semi-minor radius 923 // DGR 2007-11-20 924 case "rf": this.rf = parseFloat(paramVal); break; // inverse flattening rf= a/(a-b) 925 case "lat_0": this.lat0 = paramVal*Proj4js.common.D2R; break; // phi0, central latitude 926 case "lat_1": this.lat1 = paramVal*Proj4js.common.D2R; break; //standard parallel 1 927 case "lat_2": this.lat2 = paramVal*Proj4js.common.D2R; break; //standard parallel 2 928 case "lat_ts": this.lat_ts = paramVal*Proj4js.common.D2R; break; // used in merc and eqc 929 case "lon_0": this.long0 = paramVal*Proj4js.common.D2R; break; // lam0, central longitude 930 case "alpha": this.alpha = parseFloat(paramVal)*Proj4js.common.D2R; break; //for somerc projection 931 case "lonc": this.longc = paramVal*Proj4js.common.D2R; break; //for somerc projection 932 case "x_0": this.x0 = parseFloat(paramVal); break; // false easting 933 case "y_0": this.y0 = parseFloat(paramVal); break; // false northing 934 case "k_0": this.k0 = parseFloat(paramVal); break; // projection scale factor 935 case "k": this.k0 = parseFloat(paramVal); break; // both forms returned 936 case "r_a": this.R_A = true; break; // sphere--area of ellipsoid 937 case "zone": this.zone = parseInt(paramVal,10); break; // UTM Zone 938 case "south": this.utmSouth = true; break; // UTM north/south 939 case "towgs84":this.datum_params = paramVal.split(","); break; 940 case "to_meter": this.to_meter = parseFloat(paramVal); break; // cartesian scaling 941 case "from_greenwich": this.from_greenwich = paramVal*Proj4js.common.D2R; break; 942 // DGR 2008-07-09 : if pm is not a well-known prime meridian take 943 // the value instead of 0.0, then convert to radians 944 case "pm": paramVal = paramVal.replace(/\s/gi,""); 945 this.from_greenwich = Proj4js.PrimeMeridian[paramVal] ? 946 Proj4js.PrimeMeridian[paramVal] : parseFloat(paramVal); 947 this.from_greenwich *= Proj4js.common.D2R; 948 break; 949 // DGR 2010-11-12: axis 950 case "axis": paramVal = paramVal.replace(/\s/gi,""); 951 var legalAxis= "ewnsud"; 952 if (paramVal.length==3 && 953 legalAxis.indexOf(paramVal.substr(0,1))!=-1 && 954 legalAxis.indexOf(paramVal.substr(1,1))!=-1 && 955 legalAxis.indexOf(paramVal.substr(2,1))!=-1) { 956 this.axis= paramVal; 957 } //FIXME: be silent ? 958 break 959 case "no_defs": break; 960 default: //alert("Unrecognized parameter: " + paramName); 961 } // switch() 962 } // for paramArray 963 this.deriveConstants(); 964 }, 965 966/**
967 * Function: deriveConstants 968 * Sets several derived constant values and initialization of datum and ellipse 969 * parameters. 970 * 971 */ 972 deriveConstants: function() { 973 if (this.nagrids == '@null') this.datumCode = 'none'; 974 if (this.datumCode && this.datumCode != 'none') { 975 var datumDef = Proj4js.Datum[this.datumCode]; 976 if (datumDef) { 977 this.datum_params = datumDef.towgs84 ? datumDef.towgs84.split(',') : null; 978 this.ellps = datumDef.ellipse; 979 this.datumName = datumDef.datumName ? datumDef.datumName : this.datumCode; 980 } 981 } 982 if (!this.a) { // do we have an ellipsoid? 983 var ellipse = Proj4js.Ellipsoid[this.ellps] ? Proj4js.Ellipsoid[this.ellps] : Proj4js.Ellipsoid['WGS84']; 984 Proj4js.extend(this, ellipse); 985 } 986 if (this.rf && !this.b) this.b = (1.0 - 1.0/this.rf) * this.a; 987 if (this.rf === 0 || Math.abs(this.a - this.b)<Proj4js.common.EPSLN) { 988 this.sphere = true; 989 this.b= this.a; 990 } 991 this.a2 = this.a * this.a; // used in geocentric 992 this.b2 = this.b * this.b; // used in geocentric 993 this.es = (this.a2-this.b2)/this.a2; // e ^ 2 994 this.e = Math.sqrt(this.es); // eccentricity 995 if (this.R_A) { 996 this.a *= 1. - this.es * (Proj4js.common.SIXTH + this.es * (Proj4js.common.RA4 + this.es * Proj4js.common.RA6)); 997 this.a2 = this.a * this.a; 998 this.b2 = this.b * this.b; 999 this.es = 0.; 1000 } 1001 this.ep2=(this.a2-this.b2)/this.b2; // used in geocentric 1002 if (!this.k0) this.k0 = 1.0; //default value 1003 //DGR 2010-11-12: axis 1004 if (!this.axis) { this.axis= "enu"; } 1005 1006 this.datum = new Proj4js.datum(this); 1007 } 1008}); 1009 1010Proj4js.Proj.longlat = { 1011 init: function() { 1012 //no-op for longlat 1013 }, 1014 forward: function(pt) { 1015 //identity transform 1016 return pt; 1017 }, 1018 inverse: function(pt) { 1019 //identity transform 1020 return pt; 1021 } 1022}; 1023Proj4js.Proj.identity = Proj4js.Proj.longlat; 1024 1025/** 1026 Proj4js.defs is a collection of coordinate system definition objects in the 1027 PROJ.4 command line format. 1028 Generally a def is added by means of a separate .js file for example: 1029 1030 <SCRIPT type="text/javascript" src="defs/EPSG26912.js"></SCRIPT> 1031 1032 def is a CS definition in PROJ.4 WKT format, for example: 1033 +proj="tmerc" //longlat, etc. 1034 +a=majorRadius 1035 +b=minorRadius 1036 +lat0=somenumber 1037 +long=somenumber 1038*/ 1039Proj4js.defs = { 1040 // These are so widely used, we'll go ahead and throw them in 1041 // without requiring a separate .js file 1042 'WGS84': "+title=long/lat:WGS84 +proj=longlat +ellps=WGS84 +datum=WGS84 +units=degrees", 1043 'EPSG:4326': "+title=long/lat:WGS84 +proj=longlat +a=6378137.0 +b=6356752.31424518 +ellps=WGS84 +datum=WGS84 +units=degrees", 1044 'EPSG:4269': "+title=long/lat:NAD83 +proj=longlat +a=6378137.0 +b=6356752.31414036 +ellps=GRS80 +datum=NAD83 +units=degrees", 1045 'EPSG:3875': "+title= Google Mercator +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs" 1046}; 1047Proj4js.defs['EPSG:3785'] = Proj4js.defs['EPSG:3875']; //maintain backward compat, official code is 3875 1048Proj4js.defs['GOOGLE'] = Proj4js.defs['EPSG:3875']; 1049Proj4js.defs['EPSG:900913'] = Proj4js.defs['EPSG:3875']; 1050Proj4js.defs['EPSG:102113'] = Proj4js.defs['EPSG:3875']; 1051 1052Proj4js.common = { 1053 PI : 3.141592653589793238, //Math.PI, 1054 HALF_PI : 1.570796326794896619, //Math.PI*0.5, 1055 TWO_PI : 6.283185307179586477, //Math.PI*2, 1056 FORTPI : 0.78539816339744833, 1057 R2D : 57.29577951308232088, 1058 D2R : 0.01745329251994329577, 1059 SEC_TO_RAD : 4.84813681109535993589914102357e-6, /* SEC_TO_RAD = Pi/180/3600 */ 1060 EPSLN : 1.0e-10, 1061 MAX_ITER : 20, 1062 // following constants from geocent.c 1063 COS_67P5 : 0.38268343236508977, /* cosine of 67.5 degrees */ 1064 AD_C : 1.0026000, /* Toms region 1 constant */ 1065 1066 /* datum_type values */ 1067 PJD_UNKNOWN : 0, 1068 PJD_3PARAM : 1, 1069 PJD_7PARAM : 2, 1070 PJD_GRIDSHIFT: 3, 1071 PJD_WGS84 : 4, // WGS84 or equivalent 1072 PJD_NODATUM : 5, // WGS84 or equivalent 1073 SRS_WGS84_SEMIMAJOR : 6378137.0, // only used in grid shift transforms 1074 1075 // ellipoid pj_set_ell.c 1076 SIXTH : .1666666666666666667, /* 1/6 */ 1077 RA4 : .04722222222222222222, /* 17/360 */ 1078 RA6 : .02215608465608465608, /* 67/3024 */ 1079 RV4 : .06944444444444444444, /* 5/72 */ 1080 RV6 : .04243827160493827160, /* 55/1296 */ 1081 1082// Function to compute the constant small m which is the radius of 1083// a parallel of latitude, phi, divided by the semimajor axis. 1084// ----------------------------------------------------------------- 1085 msfnz : function(eccent, sinphi, cosphi) { 1086 var con = eccent * sinphi; 1087 return cosphi/(Math.sqrt(1.0 - con * con)); 1088 }, 1089 1090// Function to compute the constant small t for use in the forward 1091// computations in the Lambert Conformal Conic and the Polar 1092// Stereographic projections. 1093// ----------------------------------------------------------------- 1094 tsfnz : function(eccent, phi, sinphi) { 1095 var con = eccent * sinphi; 1096 var com = .5 * eccent; 1097 con = Math.pow(((1.0 - con) / (1.0 + con)), com); 1098 return (Math.tan(.5 * (this.HALF_PI - phi))/con); 1099 }, 1100 1101// Function to compute the latitude angle, phi2, for the inverse of the 1102// Lambert Conformal Conic and Polar Stereographic projections. 1103// ---------------------------------------------------------------- 1104 phi2z : function(eccent, ts) { 1105 var eccnth = .5 * eccent; 1106 var con, dphi; 1107 var phi = this.HALF_PI - 2 * Math.atan(ts); 1108 for (var i = 0; i <= 15; i++) { 1109 con = eccent * Math.sin(phi); 1110 dphi = this.HALF_PI - 2 * Math.atan(ts *(Math.pow(((1.0 - con)/(1.0 + con)),eccnth))) - phi; 1111 phi += dphi; 1112 if (Math.abs(dphi) <= .0000000001) return phi; 1113 } 1114 alert("phi2z has NoConvergence"); 1115 return (-9999); 1116 }, 1117 1118/* Function to compute constant small q which is the radius of a 1119 parallel of latitude, phi, divided by the semimajor axis. 1120------------------------------------------------------------*/ 1121 qsfnz : function(eccent,sinphi) { 1122 var con; 1123 if (eccent > 1.0e-7) { 1124 con = eccent * sinphi; 1125 return (( 1.0- eccent * eccent) * (sinphi /(1.0 - con * con) - (.5/eccent)*Math.log((1.0 - con)/(1.0 + con)))); 1126 } else { 1127 return(2.0 * sinphi); 1128 } 1129 }, 1130 1131/* Function to eliminate roundoff errors in asin 1132----------------------------------------------*/ 1133 asinz : function(x) { 1134 if (Math.abs(x)>1.0) { 1135 x=(x>1.0)?1.0:-1.0; 1136 } 1137 return Math.asin(x); 1138 }, 1139 1140// following functions from gctpc cproj.c for transverse mercator projections 1141 e0fn : function(x) {return(1.0-0.25*x*(1.0+x/16.0*(3.0+1.25*x)));}, 1142 e1fn : function(x) {return(0.375*x*(1.0+0.25*x*(1.0+0.46875*x)));}, 1143 e2fn : function(x) {return(0.05859375*x*x*(1.0+0.75*x));}, 1144 e3fn : function(x) {return(x*x*x*(35.0/3072.0));}, 1145 mlfn : function(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));}, 1146 1147 srat : function(esinp, exp) { 1148 return(Math.pow((1.0-esinp)/(1.0+esinp), exp)); 1149 }, 1150 1151// Function to return the sign of an argument 1152 sign : function(x) { if (x < 0.0) return(-1); else return(1);}, 1153 1154// Function to adjust longitude to -180 to 180; input in radians 1155 adjust_lon : function(x) { 1156 x = (Math.abs(x) < this.PI) ? x: (x - (this.sign(x)*this.TWO_PI) ); 1157 return x; 1158 }, 1159 1160// IGNF - DGR : algorithms used by IGN France 1161 1162// Function to adjust latitude to -90 to 90; input in radians 1163 adjust_lat : function(x) { 1164 x= (Math.abs(x) < this.HALF_PI) ? x: (x - (this.sign(x)*this.PI) ); 1165 return x; 1166 }, 1167 1168// Latitude Isometrique - close to tsfnz ... 1169 latiso : function(eccent, phi, sinphi) { 1170 if (Math.abs(phi) > this.HALF_PI) return +Number.NaN; 1171 if (phi==this.HALF_PI) return Number.POSITIVE_INFINITY; 1172 if (phi==-1.0*this.HALF_PI) return -1.0*Number.POSITIVE_INFINITY; 1173
1174 var con= eccent*sinphi; 1175 return Math.log(Math.tan((this.HALF_PI+phi)/2.0))+eccent*Math.log((1.0-con)/(1.0+con))/2.0; 1176 }, 1177 1178 fL : function(x,L) { 1179 return 2.0*Math.atan(x*Math.exp(L)) - this.HALF_PI; 1180 }, 1181 1182// Inverse Latitude Isometrique - close to ph2z 1183 invlatiso : function(eccent, ts) { 1184 var phi= this.fL(1.0,ts); 1185 var Iphi= 0.0; 1186 var con= 0.0; 1187 do { 1188 Iphi= phi; 1189 con= eccent*Math.sin(Iphi); 1190 phi= this.fL(Math.exp(eccent*Math.log((1.0+con)/(1.0-con))/2.0),ts) 1191 } while (Math.abs(phi-Iphi)>1.0e-12); 1192 return phi; 1193 }, 1194 1195// Needed for Gauss Schreiber 1196// Original: Denis Makarov ([email protected]) 1197// Web Site: http://www.binarythings.com 1198 sinh : function(x) 1199 { 1200 var r= Math.exp(x); 1201 r= (r-1.0/r)/2.0; 1202 return r; 1203 }, 1204 1205 cosh : function(x) 1206 { 1207 var r= Math.exp(x); 1208 r= (r+1.0/r)/2.0; 1209 return r; 1210 }, 1211 1212 tanh : function(x) 1213 { 1214 var r= Math.exp(x); 1215 r= (r-1.0/r)/(r+1.0/r); 1216 return r; 1217 }, 1218 1219 asinh : function(x) 1220 { 1221 var s= (x>= 0? 1.0:-1.0); 1222 return s*(Math.log( Math.abs(x) + Math.sqrt(x*x+1.0) )); 1223 }, 1224 1225 acosh : function(x) 1226 { 1227 return 2.0*Math.log(Math.sqrt((x+1.0)/2.0) + Math.sqrt((x-1.0)/2.0)); 1228 }, 1229 1230 atanh : function(x) 1231 { 1232 return Math.log((x-1.0)/(x+1.0))/2.0; 1233 }, 1234 1235// Grande Normale 1236 gN : function(a,e,sinphi) 1237 { 1238 var temp= e*sinphi; 1239 return a/Math.sqrt(1.0 - temp*temp); 1240 }, 1241 1242 //code from the PROJ.4 pj_mlfn.c file; this may be useful for other projections 1243 pj_enfn: function(es) { 1244 var en = new Array(); 1245 en[0] = this.C00 - es * (this.C02 + es * (this.C04 + es * (this.C06 + es * this.C08))); 1246 en[1] = es * (this.C22 - es * (this.C04 + es * (this.C06 + es * this.C08))); 1247 var t = es * es; 1248 en[2] = t * (this.C44 - es * (this.C46 + es * this.C48)); 1249 t *= es; 1250 en[3] = t * (this.C66 - es * this.C68); 1251 en[4] = t * es * this.C88; 1252 return en; 1253 }, 1254 1255 pj_mlfn: function(phi, sphi, cphi, en) { 1256 cphi *= sphi; 1257 sphi *= sphi; 1258 return(en[0] * phi - cphi * (en[1] + sphi*(en[2]+ sphi*(en[3] + sphi*en[4])))); 1259 }, 1260 1261 pj_inv_mlfn: function(arg, es, en) { 1262 var k = 1./(1.-es); 1263 var phi = arg; 1264 for (var i = Proj4js.common.MAX_ITER; i ; --i) { /* rarely goes over 2 iterations */ 1265 var s = Math.sin(phi); 1266 var t = 1. - es * s * s; 1267 //t = this.pj_mlfn(phi, s, Math.cos(phi), en) - arg; 1268 //phi -= t * (t * Math.sqrt(t)) * k; 1269 t = (this.pj_mlfn(phi, s, Math.cos(phi), en) - arg) * (t * Math.sqrt(t)) * k; 1270 phi -= t; 1271 if (Math.abs(t) < Proj4js.common.EPSLN) 1272 return phi; 1273 } 1274 Proj4js.reportError("cass:pj_inv_mlfn: Convergence error"); 1275 return phi; 1276 }, 1277 1278/* meridinal distance for ellipsoid and inverse 1279** 8th degree - accurate to < 1e-5 meters when used in conjuction 1280** with typical major axis values. 1281** Inverse determines phi to EPS (1e-11) radians, about 1e-6 seconds. 1282*/ 1283 C00: 1.0, 1284 C02: .25, 1285 C04: .046875, 1286 C06: .01953125, 1287 C08: .01068115234375, 1288 C22: .75, 1289 C44: .46875, 1290 C46: .01302083333333333333, 1291 C48: .00712076822916666666, 1292 C66: .36458333333333333333, 1293 C68: .00569661458333333333, 1294 C88: .3076171875 1295 1296}; 1297 1298/** datum object 1299*/ 1300Proj4js.datum = Proj4js.Class({ 1301 1302 initialize : function(proj) { 1303 this.datum_type = Proj4js.common.PJD_WGS84; //default setting 1304 if (proj.datumCode && proj.datumCode == 'none') { 1305 this.datum_type = Proj4js.common.PJD_NODATUM; 1306 } 1307 if (proj && proj.datum_params) { 1308 for (var i=0; i<proj.datum_params.length; i++) { 1309 proj.datum_params[i]=parseFloat(proj.datum_params[i]); 1310 } 1311 if (proj.datum_params[0] != 0 || proj.datum_params[1] != 0 || proj.datum_params[2] != 0 ) { 1312 this.datum_type = Proj4js.common.PJD_3PARAM; 1313 } 1314 if (proj.datum_params.length > 3) { 1315 if (proj.datum_params[3] != 0 || proj.datum_params[4] != 0 || 1316 proj.datum_params[5] != 0 || proj.datum_params[6] != 0 ) { 1317 this.datum_type = Proj4js.common.PJD_7PARAM; 1318 proj.datum_params[3] *= Proj4js.common.SEC_TO_RAD; 1319 proj.datum_params[4] *= Proj4js.common.SEC_TO_RAD; 1320 proj.datum_params[5] *= Proj4js.common.SEC_TO_RAD; 1321 proj.datum_params[6] = (proj.datum_params[6]/1000000.0) + 1.0; 1322 } 1323 } 1324 } 1325 if (proj) { 1326 this.a = proj.a; //datum object also uses these values 1327 this.b = proj.b; 1328 this.es = proj.es; 1329 this.ep2 = proj.ep2; 1330 this.datum_params = proj.datum_params; 1331 } 1332 }, 1333 1334 /****************************************************************/ 1335 // cs_compare_datums() 1336 // Returns TRUE if the two datums match, otherwise FALSE. 1337 compare_datums : function( dest ) { 1338 if( this.datum_type != dest.datum_type ) { 1339 return false; // false, datums are not equal 1340 } else if( this.a != dest.a || Math.abs(this.es-dest.es) > 0.000000000050 ) { 1341 // the tolerence for es is to ensure that GRS80 and WGS84 1342 // are considered identical 1343 return false;
1344 } else if( this.datum_type == Proj4js.common.PJD_3PARAM ) { 1345 return (this.datum_params[0] == dest.datum_params[0] 1346 && this.datum_params[1] == dest.datum_params[1] 1347 && this.datum_params[2] == dest.datum_params[2]); 1348 } else if( this.datum_type == Proj4js.common.PJD_7PARAM ) { 1349 return (this.datum_params[0] == dest.datum_params[0] 1350 && this.datum_params[1] == dest.datum_params[1] 1351 && this.datum_params[2] == dest.datum_params[2] 1352 && this.datum_params[3] == dest.datum_params[3] 1353 && this.datum_params[4] == dest.datum_params[4] 1354 && this.datum_params[5] == dest.datum_params[5] 1355 && this.datum_params[6] == dest.datum_params[6]); 1356 } else if ( this.datum_type == Proj4js.common.PJD_GRIDSHIFT || 1357 dest.datum_type == Proj4js.common.PJD_GRIDSHIFT ) { 1358 alert("ERROR: Grid shift transformations are not implemented."); 1359 return false 1360 } else { 1361 return true; // datums are equal 1362 } 1363 }, // cs_compare_datums() 1364 1365 /* 1366 * The function Convert_Geodetic_To_Geocentric converts geodetic coordinates 1367 * (latitude, longitude, and height) to geocentric coordinates (X, Y, Z), 1368 * according to the current ellipsoid parameters. 1369 * 1370 * Latitude : Geodetic latitude in radians (input) 1371 * Longitude : Geodetic longitude in radians (input) 1372 * Height : Geodetic height, in meters (input) 1373 * X : Calculated Geocentric X coordinate, in meters (output) 1374 * Y : Calculated Geocentric Y coordinate, in meters (output) 1375 * Z : Calculated Geocentric Z coordinate, in meters (output) 1376 * 1377 */ 1378 geodetic_to_geocentric : function(p) { 1379 var Longitude = p.x; 1380 var Latitude = p.y; 1381 var Height = p.z ? p.z : 0; //Z value not always supplied 1382 var X; // output 1383 var Y; 1384 var Z; 1385 1386 var Error_Code=0; // GEOCENT_NO_ERROR; 1387 var Rn; /* Earth radius at location */ 1388 var Sin_Lat; /* Math.sin(Latitude) */ 1389 var Sin2_Lat; /* Square of Math.sin(Latitude) */ 1390 var Cos_Lat; /* Math.cos(Latitude) */ 1391 1392 /* 1393 ** Don't blow up if Latitude is just a little out of the value 1394 ** range as it may just be a rounding issue. Also removed longitude 1395 ** test, it should be wrapped by Math.cos() and Math.sin(). NFW for PROJ.4, Sep/2001. 1396 */ 1397 if( Latitude < -Proj4js.common.HALF_PI && Latitude > -1.001 * Proj4js.common.HALF_PI ) { 1398 Latitude = -Proj4js.common.HALF_PI; 1399 } else if( Latitude > Proj4js.common.HALF_PI && Latitude < 1.001 * Proj4js.common.HALF_PI ) { 1400 Latitude = Proj4js.common.HALF_PI; 1401 } else if ((Latitude < -Proj4js.common.HALF_PI) || (Latitude > Proj4js.common.HALF_PI)) { 1402 /* Latitude out of range */ 1403 Proj4js.reportError('geocent:lat out of range:'+Latitude); 1404 return null; 1405 } 1406 1407 if (Longitude > Proj4js.common.PI) Longitude -= (2*Proj4js.common.PI); 1408 Sin_Lat = Math.sin(Latitude); 1409 Cos_Lat = Math.cos(Latitude); 1410 Sin2_Lat = Sin_Lat * Sin_Lat; 1411 Rn = this.a / (Math.sqrt(1.0e0 - this.es * Sin2_Lat)); 1412 X = (Rn + Height) * Cos_Lat * Math.cos(Longitude); 1413 Y = (Rn + Height) * Cos_Lat * Math.sin(Longitude); 1414 Z = ((Rn * (1 - this.es)) + Height) * Sin_Lat; 1415 1416 p.x = X; 1417 p.y = Y; 1418 p.z = Z; 1419 return Error_Code; 1420 }, // cs_geodetic_to_geocentric() 1421 1422 1423 geocentric_to_geodetic : function (p) { 1424/* local defintions and variables */ 1425/* end-criterium of loop, accuracy of sin(Latitude) */ 1426var genau = 1.E-12; 1427var genau2 = (genau*genau); 1428var maxiter = 30; 1429 1430 var P; /* distance between semi-minor axis and location */ 1431 var RR; /* distance between center and location */ 1432 var CT; /* sin of geocentric latitude */ 1433 var ST; /* cos of geocentric latitude */ 1434 var RX; 1435 var RK; 1436 var RN; /* Earth radius at location */ 1437 var CPHI0; /* cos of start or old geodetic latitude in iterations */ 1438 var SPHI0; /* sin of start or old geodetic latitude in iterations */ 1439 var CPHI; /* cos of searched geodetic latitude */ 1440 var SPHI; /* sin of searched geodetic latitude */ 1441 var SDPHI; /* end-criterium: addition-theorem of sin(Latitude(iter)-Latitude(iter-1)) */ 1442 var At_Pole; /* indicates location is in polar region */ 1443 var iter; /* # of continous iteration, max. 30 is always enough (s.a.) */ 1444 1445 var X = p.x; 1446 var Y = p.y; 1447 var Z = p.z ? p.z : 0.0; //Z value not always supplied 1448 var Longitude; 1449 var Latitude; 1450 var Height; 1451 1452 At_Pole = false;
1453 P = Math.sqrt(X*X+Y*Y); 1454 RR = Math.sqrt(X*X+Y*Y+Z*Z); 1455 1456/* special cases for latitude and longitude */ 1457 if (P/this.a < genau) { 1458 1459/* special case, if P=0. (X=0., Y=0.) */ 1460 At_Pole = true; 1461 Longitude = 0.0; 1462 1463/* if (X,Y,Z)=(0.,0.,0.) then Height becomes semi-minor axis 1464 * of ellipsoid (=center of mass), Latitude becomes PI/2 */ 1465 if (RR/this.a < genau) { 1466 Latitude = Proj4js.common.HALF_PI; 1467 Height = -this.b; 1468 return; 1469 } 1470 } else { 1471/* ellipsoidal (geodetic) longitude 1472 * interval: -PI < Longitude <= +PI */ 1473 Longitude=Math.atan2(Y,X); 1474 } 1475 1476/* -------------------------------------------------------------- 1477 * Following iterative algorithm was developped by 1478 * "Institut f�r Erdmessung", University of Hannover, July 1988. 1479 * Internet: www.ife.uni-hannover.de 1480 * Iterative computation of CPHI,SPHI and Height. 1481 * Iteration of CPHI and SPHI to 10**-12 radian resp. 1482 * 2*10**-7 arcsec. 1483 * -------------------------------------------------------------- 1484 */ 1485 CT = Z/RR; 1486 ST = P/RR; 1487 RX = 1.0/Math.sqrt(1.0-this.es*(2.0-this.es)*ST*ST); 1488 CPHI0 = ST*(1.0-this.es)*RX; 1489 SPHI0 = CT*RX; 1490 iter = 0; 1491 1492/* loop to find sin(Latitude) resp. Latitude 1493 * until |sin(Latitude(iter)-Latitude(iter-1))| < genau */ 1494 do 1495 { 1496 iter++; 1497 RN = this.a/Math.sqrt(1.0-this.es*SPHI0*SPHI0); 1498 1499/* ellipsoidal (geodetic) height */ 1500 Height = P*CPHI0+Z*SPHI0-RN*(1.0-this.es*SPHI0*SPHI0); 1501 1502 RK = this.es*RN/(RN+Height); 1503 RX = 1.0/Math.sqrt(1.0-RK*(2.0-RK)*ST*ST); 1504 CPHI = ST*(1.0-RK)*RX; 1505 SPHI = CT*RX; 1506 SDPHI = SPHI*CPHI0-CPHI*SPHI0; 1507 CPHI0 = CPHI; 1508 SPHI0 = SPHI; 1509 } 1510 while (SDPHI*SDPHI > genau2 && iter < maxiter); 1511 1512/* ellipsoidal (geodetic) latitude */ 1513 Latitude=Math.atan(SPHI/Math.abs(CPHI)); 1514 1515 p.x = Longitude; 1516 p.y = Latitude; 1517 p.z = Height; 1518 return p; 1519 }, // cs_geocentric_to_geodetic() 1520 1521 /** Convert_Geocentric_To_Geodetic 1522 * The method used here is derived from 'An Improved Algorithm for 1523 * Geocentric to Geodetic Coordinate Conversion', by Ralph Toms, Feb 1996 1524 */ 1525 geocentric_to_geodetic_noniter : function (p) { 1526 var X = p.x; 1527 var Y = p.y; 1528 var Z = p.z ? p.z : 0; //Z value not always supplied 1529 var Longitude; 1530 var Latitude; 1531 var Height; 1532 1533 var W; /* distance from Z axis */ 1534 var W2; /* square of distance from Z axis */ 1535 var T0; /* initial estimate of vertical component */ 1536 var T1; /* corrected estimate of vertical component */ 1537 var S0; /* initial estimate of horizontal component */ 1538 var S1; /* corrected estimate of horizontal component */ 1539 var Sin_B0; /* Math.sin(B0), B0 is estimate of Bowring aux variable */ 1540 var Sin3_B0; /* cube of Math.sin(B0) */ 1541 var Cos_B0; /* Math.cos(B0) */ 1542 var Sin_p1; /* Math.sin(phi1), phi1 is estimated latitude */ 1543 var Cos_p1; /* Math.cos(phi1) */ 1544 var Rn; /* Earth radius at location */ 1545 var Sum; /* numerator of Math.cos(phi1) */ 1546 var At_Pole; /* indicates location is in polar region */ 1547 1548 X = parseFloat(X); // cast from string to float 1549 Y = parseFloat(Y); 1550 Z = parseFloat(Z); 1551 1552 At_Pole = false; 1553 if (X != 0.0) 1554 { 1555 Longitude = Math.atan2(Y,X); 1556 } 1557 else 1558 { 1559 if (Y > 0) 1560 { 1561 Longitude = Proj4js.common.HALF_PI; 1562 } 1563 else if (Y < 0) 1564 { 1565 Longitude = -Proj4js.common.HALF_PI; 1566 } 1567 else 1568 { 1569 At_Pole = true; 1570 Longitude = 0.0; 1571 if (Z > 0.0) 1572 { /* north pole */ 1573 Latitude = Proj4js.common.HALF_PI; 1574 } 1575 else if (Z < 0.0) 1576 { /* south pole */ 1577 Latitude = -Proj4js.common.HALF_PI; 1578 } 1579 else 1580 { /* center of earth */ 1581 Latitude = Proj4js.common.HALF_PI; 1582 Height = -this.b; 1583 return; 1584 } 1585 } 1586 } 1587 W2 = X*X + Y*Y; 1588 W = Math.sqrt(W2); 1589 T0 = Z * Proj4js.common.AD_C; 1590 S0 = Math.sqrt(T0 * T0 + W2); 1591 Sin_B0 = T0 / S0; 1592 Cos_B0 = W / S0; 1593 Sin3_B0 = Sin_B0 * Sin_B0 * Sin_B0; 1594 T1 = Z + this.b * this.ep2 * Sin3_B0; 1595 Sum = W - this.a * this.es * Cos_B0 * Cos_B0 * Cos_B0; 1596 S1 = Math.sqrt(T1*T1 + Sum * Sum); 1597 Sin_p1 = T1 / S1; 1598 Cos_p1 = Sum / S1; 1599 Rn = this.a / Math.sqrt(1.0 - this.es * Sin_p1 * Sin_p1); 1600 if (Cos_p1 >= Proj4js.common.COS_67P5) 1601 { 1602 Height = W / Cos_p1 - Rn; 1603 } 1604 else if (Cos_p1 <= -Proj4js.common.COS_67P5) 1605 { 1606 Height = W / -Cos_p1 - Rn; 1607 } 1608 else 1609 { 1610 Height = Z / Sin_p1 + Rn * (this.es - 1.0); 1611 } 1612 if (At_Pole == false) 1613 { 1614 Latitude = Math.atan(Sin_p1 / Cos_p1); 1615 } 1616 1617 p.x = Longitude; 1618 p.y = Latitude; 1619 p.z = Height; 1620 return p; 1621 }, // geocentric_to_geodetic_noniter() 1622 1623 /****************************************************************/ 1624 // pj_geocentic_to_wgs84( p ) 1625 // p = point to transform in geocentric coordinates (x,y,z) 1626 geocentric_to_wgs84 : function ( p ) { 1627 1628 if( this.datum_type == Proj4js.common.PJD_3PARAM ) 1629 { 1630 // if( x[io] == HUGE_VAL ) 1631 // continue; 1632 p.x += this.datum_params[0]; 1633 p.y += this.datum_params[1]; 1634 p.z += this.datum_params[2]; 1635 1636 } 1637 else if (this.datum_type == Proj4js.common.PJD_7PARAM) 1638 { 1639 var Dx_BF =this.datum_params[0]; 1640 var Dy_BF =this.datum_params[1]; 1641 var Dz_BF =this.datum_params[2]; 1642 var Rx_BF =this.datum_params[3]; 1643 var Ry_BF =this.datum_params[4]; 1644 var Rz_BF =this.datum_params[5]; 1645 var M_BF =this.datum_params[6]; 1646 // if( x[io] == HUGE_VAL ) 1647 // continue; 1648 var x_out = M_BF*( p.x - Rz_BF*p.y + Ry_BF*p.z) + Dx_BF; 1649 var y_out = M_BF*( Rz_BF*p.x + p.y - Rx_BF*p.z) + Dy_BF; 1650 var z_out = M_BF*(-Ry_BF*p.x + Rx_BF*p.y + p.z) + Dz_B
1650F; 1651 p.x = x_out; 1652 p.y = y_out; 1653 p.z = z_out; 1654 } 1655 }, // cs_geocentric_to_wgs84 1656 1657 /****************************************************************/ 1658 // pj_geocentic_from_wgs84() 1659 // coordinate system definition, 1660 // point to transform in geocentric coordinates (x,y,z) 1661 geocentric_from_wgs84 : function( p ) { 1662 1663 if( this.datum_type == Proj4js.common.PJD_3PARAM ) 1664 { 1665 //if( x[io] == HUGE_VAL ) 1666 // continue; 1667 p.x -= this.datum_params[0]; 1668 p.y -= this.datum_params[1]; 1669 p.z -= this.datum_params[2]; 1670 1671 } 1672 else if (this.datum_type == Proj4js.common.PJD_7PARAM) 1673 { 1674 var Dx_BF =this.datum_params[0]; 1675 var Dy_BF =this.datum_params[1]; 1676 var Dz_BF =this.datum_params[2]; 1677 var Rx_BF =this.datum_params[3]; 1678 var Ry_BF =this.datum_params[4]; 1679 var Rz_BF =this.datum_params[5]; 1680 var M_BF =this.datum_params[6]; 1681 var x_tmp = (p.x - Dx_BF) / M_BF; 1682 var y_tmp = (p.y - Dy_BF) / M_BF; 1683 var z_tmp = (p.z - Dz_BF) / M_BF; 1684 //if( x[io] == HUGE_VAL ) 1685 // continue; 1686 1687 p.x = x_tmp + Rz_BF*y_tmp - Ry_BF*z_tmp; 1688 p.y = -Rz_BF*x_tmp + y_tmp + Rx_BF*z_tmp; 1689 p.z = Ry_BF*x_tmp - Rx_BF*y_tmp + z_tmp; 1690 } //cs_geocentric_from_wgs84() 1691 } 1692}); 1693 1694/** point object, nothing fancy, just allows values to be 1695 passed back and forth by reference rather than by value. 1696 Other point classes may be used as long as they have 1697 x and y properties, which will get modified in the transform method. 1698*/ 1699Proj4js.Point = Proj4js.Class({ 1700 1701 /** 1702 * Constructor: Proj4js.Point 1703 * 1704 * Parameters: 1705 * - x {float} or {Array} either the first coordinates component or 1706 * the full coordinates 1707 * - y {float} the second component 1708 * - z {float} the third component, optional. 1709 */ 1710 initialize : function(x,y,z) { 1711 if (typeof x == 'object') { 1712 this.x = x[0]; 1713 this.y = x[1]; 1714 this.z = x[2] || 0.0; 1715 } else if (typeof x == 'string' && typeof y == 'undefined') { 1716 var coords = x.split(','); 1717 this.x = parseFloat(coords[0]); 1718 this.y = parseFloat(coords[1]); 1719 this.z = parseFloat(coords[2]) || 0.0; 1720 } else { 1721 this.x = x; 1722 this.y = y; 1723 this.z = z || 0.0; 1724 } 1725 }, 1726 1727 /** 1728 * APIMethod: clone 1729 * Build a copy of a Proj4js.Point object. 1730 * 1731 * Return: 1732 * {Proj4js}.Point the cloned point. 1733 */ 1734 clone : function() { 1735 return new Proj4js.Point(this.x, this.y, this.z); 1736 }, 1737 1738 /** 1739 * APIMethod: toString 1740 * Return a readable string version of the point 1741 * 1742 * Return: 1743 * {String} String representation of Proj4js.Point object. 1744 * (ex. <i>"x=5,y=42"</i>) 1745 */ 1746 toString : function() { 1747 return ("x=" + this.x + ",y=" + this.y); 1748 }, 1749 1750 /** 1751 * APIMethod: toShortString 1752 * Return a short string version of the point. 1753 * 1754 * Return: 1755 * {String} Shortened String representation of Proj4js.Point object. 1756 * (ex. <i>"5, 42"</i>) 1757 */ 1758 toShortString : function() { 1759 return (this.x + ", " + this.y); 1760 } 1761}); 1762 1763Proj4js.PrimeMeridian = { 1764 "greenwich": 0.0, //"0dE", 1765 "lisbon": -9.131906111111, //"9d07'54.862\"W", 1766 "paris": 2.337229166667, //"2d20'14.025\"E", 1767 "bogota": -74.080916666667, //"74d04'51.3\"W", 1768 "madrid": -3.687938888889, //"3d41'16.58\"W", 1769 "rome": 12.452333333333, //"12d27'8.4\"E", 1770 "bern": 7.439583333333, //"7d26'22.5\"E", 1771 "jakarta": 106.807719444444, //"106d48'27.79\"E", 1772 "ferro": -17.666666666667, //"17d40'W", 1773 "brussels": 4.367975, //"4d22'4.71\"E", 1774 "stockholm": 18.058277777778, //"18d3'29.8\"E", 1775 "athens": 23.7163375, //"23d42'58.815\"E", 1776 "oslo": 10.722916666667 //"10d43'22.5\"E" 1777}; 1778 1779Proj4js.Ellipsoid = { 1780 "MERIT": {a:6378137.0, rf:298.257, ellipseName:"MERIT 1983"}, 1781 "SGS85": {a:6378136.0, rf:298.257, ellipseName:"Soviet Geodetic System 85"}, 1782 "GRS80": {a:6378137.0, rf:298.257222101, ellipseName:"GRS 1980(IUGG, 1980)"}, 1783 "IAU76": {a:6378140.0, rf:298.257, ellipseName:"IAU 1976"}, 1784 "airy": {a:6377563.396, b:6356256.910, ellipseName:"Airy 1830"}, 1785 "APL4.": {a:6378137, rf:298.25, ellipseName:"Appl. Physics. 1965"}, 1786 "NWL9D": {a:6378145.0, rf:298.25, ellipseName:"Naval Weapons Lab., 1965"}, 1787 "mod_airy": {a:6377340.189, b:6356034.446, ellipseName:"Modified Airy"}, 1788 "andrae": {a:6377104.43, rf:300.0, ellipseName:"Andrae 1876 (Den., Iclnd.)"}, 1789 "aust_SA": {a:6378160.0, rf:298.25, ellipseName:"Australian Natl & S. Amer. 1969"}, 1790 "GRS67": {a:6378160.0, rf:298.2471674270, ellipseName:"GRS 67(IUGG 1967)"}, 1791 "bessel": {a:6377397.155, rf:299.1528128, ellipseName:"Bessel 1841"},
1792 "bess_nam": {a:6377483.865, rf:299.1528128, ellipseName:"Bessel 1841 (Namibia)"}, 1793 "clrk66": {a:6378206.4, b:6356583.8, ellipseName:"Clarke 1866"}, 1794 "clrk80": {a:6378249.145, rf:293.4663, ellipseName:"Clarke 1880 mod."}, 1795 "CPM": {a:6375738.7, rf:334.29, ellipseName:"Comm. des Poids et Mesures 1799"}, 1796 "delmbr": {a:6376428.0, rf:311.5, ellipseName:"Delambre 1810 (Belgium)"}, 1797 "engelis": {a:6378136.05, rf:298.2566, ellipseName:"Engelis 1985"}, 1798 "evrst30": {a:6377276.345, rf:300.8017, ellipseName:"Everest 1830"}, 1799 "evrst48": {a:6377304.063, rf:300.8017, ellipseName:"Everest 1948"}, 1800 "evrst56": {a:6377301.243, rf:300.8017, ellipseName:"Everest 1956"}, 1801 "evrst69": {a:6377295.664, rf:300.8017, ellipseName:"Everest 1969"}, 1802 "evrstSS": {a:6377298.556, rf:300.8017, ellipseName:"Everest (Sabah & Sarawak)"}, 1803 "fschr60": {a:6378166.0, rf:298.3, ellipseName:"Fischer (Mercury Datum) 1960"}, 1804 "fschr60m": {a:6378155.0, rf:298.3, ellipseName:"Fischer 1960"}, 1805 "fschr68": {a:6378150.0, rf:298.3, ellipseName:"Fischer 1968"}, 1806 "helmert": {a:6378200.0, rf:298.3, ellipseName:"Helmert 1906"}, 1807 "hough": {a:6378270.0, rf:297.0, ellipseName:"Hough"}, 1808 "intl": {a:6378388.0, rf:297.0, ellipseName:"International 1909 (Hayford)"}, 1809 "kaula": {a:6378163.0, rf:298.24, ellipseName:"Kaula 1961"}, 1810 "lerch": {a:6378139.0, rf:298.257, ellipseName:"Lerch 1979"}, 1811 "mprts": {a:6397300.0, rf:191.0, ellipseName:"Maupertius 1738"}, 1812 "new_intl": {a:6378157.5, b:6356772.2, ellipseName:"New International 1967"}, 1813 "plessis": {a:6376523.0, rf:6355863.0, ellipseName:"Plessis 1817 (France)"}, 1814 "krass": {a:6378245.0, rf:298.3, ellipseName:"Krassovsky, 1942"}, 1815 "SEasia": {a:6378155.0, b:6356773.3205, ellipseName:"Southeast Asia"}, 1816 "walbeck": {a:6376896.0, b:6355834.8467, ellipseName:"Walbeck"}, 1817 "WGS60": {a:6378165.0, rf:298.3, ellipseName:"WGS 60"}, 1818 "WGS66": {a:6378145.0, rf:298.25, ellipseName:"WGS 66"}, 1819 "WGS72": {a:6378135.0, rf:298.26, ellipseName:"WGS 72"}, 1820 "WGS84": {a:6378137.0, rf:298.257223563, ellipseName:"WGS 84"}, 1821 "sphere": {a:6370997.0, b:6370997.0, ellipseName:"Normal Sphere (r=6370997)"} 1822}; 1823 1824Proj4js.Datum = { 1825 "WGS84": {towgs84: "0,0,0", ellipse: "WGS84", datumName: "WGS84"}, 1826 "GGRS87": {towgs84: "-199.87,74.79,246.62", ellipse: "GRS80", datumName: "Greek_Geodetic_Reference_System_1987"}, 1827 "NAD83": {towgs84: "0,0,0", ellipse: "GRS80", datumName: "North_American_Datum_1983"}, 1828 "NAD27": {nadgrids: "@conus,@alaska,@ntv2_0.gsb,@ntv1_can.dat", ellipse: "clrk66", datumName: "North_American_Datum_1927"}, 1829 "potsdam": {towgs84: "606.0,23.0,413.0", ellipse: "bessel", datumName: "Potsdam Rauenberg 1950 DHDN"}, 1830 "carthage": {towgs84: "-263.0,6.0,431.0", ellipse: "clark80", datumName: "Carthage 1934 Tunisia"}, 1831 "hermannskogel": {towgs84: "653.0,-212.0,449.0", ellipse: "bessel", datumName: "Hermannskogel"}, 1832 "ire65": {towgs84: "482.530,-130.596,564.557,-1.042,-0.214,-0.631,8.15", ellipse: "mod_airy", datumName: "Ireland 1965"}, 1833 "nzgd49": {towgs84: "59.47,-5.04,187.44,0.47,-0.1,1.024,-4.5993", ellipse: "intl", datumName: "New Zealand Geodetic Datum 1949"}, 1834 "OSGB36": {towgs84: "446.448,-125.157,542.060,0.1502,0.2470,0.8421,-20.4894", ellipse: "airy", datumName: "Airy 1830"} 1835}; 1836 1837Proj4js.WGS84 = new Proj4js.Proj('WGS84'); 1838Proj4js.Datum['OSB36'] = Proj4js.Datum['OSGB36']; //as returned from spatialreference.org 1839 1840//lookup table to go from the projection name in WKT to the Proj4js projection name 1841//build this out as required 1842Proj4js.wktProjections = { 1843 "Lambert Tangential Conformal Conic Projection": "lcc", 1844 "Mercator": "merc", 1845 "Popular Visualisation Pseudo Mercator": "merc", 1846 "Mercator_1SP": "merc", 1847 "Transverse_Mercator": "tmerc", 1848 "Transverse Mercator": "tmerc", 1849 "Lambert Azimuthal Equal Area": "laea", 1850 "Universal Transverse Mercator System": "utm" 1851}; 1852 1853 1854/* ====================================================================== 1855 projCode/aea.js 1856 ====================================================================== */ 1857 1858/******************************************************************************* 1859NAME ALBERS CONICAL EQUAL AREA 1860 1861PURPOSE: Transforms input longitude and latitude to Easting and Northing 1862 for the Albers Conical Equal Area projection. The longitude 1863 and latitude must be in radians. The Easting and Northing 1864 values will be returned in meters. 1865 1866PROGRAMMER DATE 1867---------- ----
1868T. Mittan, Feb, 1992 1869 1870ALGORITHM REFERENCES 1871 18721. Snyder, John P., "Map Projections--A Working Manual", U.S. Geological 1873 Survey Professional Paper 1395 (Supersedes USGS Bulletin 1532), United 1874 State Government Printing Office, Washington D.C., 1987. 1875 18762. Snyder, John P. and Voxland, Philip M., "An Album of Map Projections", 1877 U.S. Geological Survey Professional Paper 1453 , United State Government 1878 Printing Office, Washington D.C., 1989. 1879*******************************************************************************/ 1880 1881 1882Proj4js.Proj.aea = { 1883 init : function() { 1884 1885 if (Math.abs(this.lat1 + this.lat2) < Proj4js.common.EPSLN) { 1886 Proj4js.reportError("aeaInitEqualLatitudes"); 1887 return; 1888 } 1889 this.temp = this.b / this.a; 1890 this.es = 1.0 - Math.pow(this.temp,2); 1891 this.e3 = Math.sqrt(this.es); 1892 1893 this.sin_po=Math.sin(this.lat1); 1894 this.cos_po=Math.cos(this.lat1); 1895 this.t1=this.sin_po; 1896 this.con = this.sin_po; 1897 this.ms1 = Proj4js.common.msfnz(this.e3,this.sin_po,this.cos_po); 1898 this.qs1 = Proj4js.common.qsfnz(this.e3,this.sin_po,this.cos_po); 1899 1900 this.sin_po=Math.sin(this.lat2); 1901 this.cos_po=Math.cos(this.lat2); 1902 this.t2=this.sin_po; 1903 this.ms2 = Proj4js.common.msfnz(this.e3,this.sin_po,this.cos_po); 1904 this.qs2 = Proj4js.common.qsfnz(this.e3,this.sin_po,this.cos_po); 1905 1906 this.sin_po=Math.sin(this.lat0); 1907 this.cos_po=Math.cos(this.lat0); 1908 this.t3=this.sin_po; 1909 this.qs0 = Proj4js.common.qsfnz(this.e3,this.sin_po,this.cos_po); 1910 1911 if (Math.abs(this.lat1 - this.lat2) > Proj4js.common.EPSLN) { 1912 this.ns0 = (this.ms1 * this.ms1 - this.ms2 *this.ms2)/ (this.qs2 - this.qs1); 1913 } else { 1914 this.ns0 = this.con; 1915 } 1916 this.c = this.ms1 * this.ms1 + this.ns0 * this.qs1; 1917 this.rh = this.a * Math.sqrt(this.c - this.ns0 * this.qs0)/this.ns0; 1918 }, 1919 1920/* Albers Conical Equal Area forward equations--mapping lat,long to x,y 1921 -------------------------------------------------------------------*/ 1922 forward: function(p){ 1923 1924 var lon=p.x; 1925 var lat=p.y; 1926 1927 this.sin_phi=Math.sin(lat); 1928 this.cos_phi=Math.cos(lat); 1929 1930 var qs = Proj4js.common.qsfnz(this.e3,this.sin_phi,this.cos_phi); 1931 var rh1 =this.a * Math.sqrt(this.c - this.ns0 * qs)/this.ns0; 1932 var theta = this.ns0 * Proj4js.common.adjust_lon(lon - this.long0); 1933 var x = rh1 * Math.sin(theta) + this.x0; 1934 var y = this.rh - rh1 * Math.cos(theta) + this.y0; 1935 1936 p.x = x; 1937 p.y = y; 1938 return p; 1939 }, 1940 1941 1942 inverse: function(p) { 1943 var rh1,qs,con,theta,lon,lat; 1944 1945 p.x -= this.x0; 1946 p.y = this.rh - p.y + this.y0; 1947 if (this.ns0 >= 0) { 1948 rh1 = Math.sqrt(p.x *p.x + p.y * p.y); 1949 con = 1.0; 1950 } else { 1951 rh1 = -Math.sqrt(p.x * p.x + p.y *p.y); 1952 con = -1.0; 1953 } 1954 theta = 0.0; 1955 if (rh1 != 0.0) { 1956 theta = Math.atan2(con * p.x, con * p.y); 1957 } 1958 con = rh1 * this.ns0 / this.a; 1959 qs = (this.c - con * con) / this.ns0; 1960 if (this.e3 >= 1e-10) { 1961 con = 1 - .5 * (1.0 -this.es) * Math.log((1.0 - this.e3) / (1.0 + this.e3))/this.e3; 1962 if (Math.abs(Math.abs(con) - Math.abs(qs)) > .0000000001 ) { 1963 lat = this.phi1z(this.e3,qs); 1964 } else { 1965 if (qs >= 0) { 1966 lat = .5 * Proj4js.common.PI; 1967 } else { 1968 lat = -.5 * Proj4js.common.PI; 1969 } 1970 } 1971 } else { 1972 lat = this.phi1z(this.e3,qs); 1973 } 1974 1975 lon = Proj4js.common.adjust_lon(theta/this.ns0 + this.long0); 1976 p.x = lon; 1977 p.y = lat; 1978 return p; 1979 }, 1980 1981/* Function to compute phi1, the latitude for the inverse of the 1982 Albers Conical Equal-Area projection. 1983-------------------------------------------*/ 1984 phi1z: function (eccent,qs) { 1985 var sinphi, cosphi, con, com, dphi; 1986 var phi = Proj4js.common.asinz(.5 * qs); 1987 if (eccent < Proj4js.common.EPSLN) return phi; 1988 1989 var eccnts = eccent * eccent; 1990 for (var i = 1; i <= 25; i++) { 1991 sinphi = Math.sin(phi); 1992 cosphi = Math.cos(phi); 1993 con = eccent * sinphi; 1994 com = 1.0 - con * con; 1995 dphi = .5 * com * com / cosphi * (qs / (1.0 - eccnts) - sinphi / com + .5 / eccent * Math.log((1.0 - con) / (1.0 + con))); 1996 phi = phi + dphi; 1997 if (Math.abs(dphi) <= 1e-7) return phi; 1998 } 1999 Proj4js.reportError("aea:phi1z:Convergence error"); 2000 return null; 2001 } 2002 2003}; 2004 2005 2006 2007/* ====================================================================== 2008 projCode/sterea.js 2009 ====================================================================== */ 2010 2011 2012Proj4js.Proj.sterea = { 2013 dependsOn : 'gauss', 2014 2015 init : function() { 2016 Proj4js.Proj['gauss'].init.apply(this); 2017 if (!this.rc) { 2018 Proj4js.reportError("sterea:init:E_ERROR_0"); 2019 return; 2020 } 2021 this.sinc0 = Math.sin(this.phic0); 2022 this.cosc0 = Math.cos(this.phic0); 2023 this.R2 = 2.0 * this.rc; 2024 if (!this.title) this.title = "Oblique Stereographic Alternative"; 2025 }, 2026 2027 forward : function(p) { 2028 var sinc, cosc, cosl, k; 2029 p.x = Proj4js.common.adjust_lon(p.x-this.long0); /* adjust del longitude */ 2030 Proj4js.Proj['gauss'].forward.apply(this, [p]); 2031 sinc = Math.sin(p.y); 2032 cosc = Math.cos(p.y);
2033 cosl = Math.cos(p.x); 2034 k = this.k0 * this.R2 / (1.0 + this.sinc0 * sinc + this.cosc0 * cosc * cosl); 2035 p.x = k * cosc * Math.sin(p.x); 2036 p.y = k * (this.cosc0 * sinc - this.sinc0 * cosc * cosl); 2037 p.x = this.a * p.x + this.x0; 2038 p.y = this.a * p.y + this.y0; 2039 return p; 2040 }, 2041 2042 inverse : function(p) { 2043 var sinc, cosc, lon, lat, rho; 2044 p.x = (p.x - this.x0) / this.a; /* descale and de-offset */ 2045 p.y = (p.y - this.y0) / this.a; 2046 2047 p.x /= this.k0; 2048 p.y /= this.k0; 2049 if ( (rho = Math.sqrt(p.x*p.x + p.y*p.y)) ) { 2050 var c = 2.0 * Math.atan2(rho, this.R2); 2051 sinc = Math.sin(c); 2052 cosc = Math.cos(c); 2053 lat = Math.asin(cosc * this.sinc0 + p.y * sinc * this.cosc0 / rho); 2054 lon = Math.atan2(p.x * sinc, rho * this.cosc0 * cosc - p.y * this.sinc0 * sinc); 2055 } else { 2056 lat = this.phic0; 2057 lon = 0.; 2058 } 2059 2060 p.x = lon; 2061 p.y = lat; 2062 Proj4js.Proj['gauss'].inverse.apply(this,[p]); 2063 p.x = Proj4js.common.adjust_lon(p.x + this.long0); /* adjust longitude to CM */ 2064 return p; 2065 } 2066}; 2067 2068/* ====================================================================== 2069 projCode/poly.js 2070 ====================================================================== */ 2071 2072/* Function to compute, phi4, the latitude for the inverse of the 2073 Polyconic projection. 2074------------------------------------------------------------*/ 2075function phi4z (eccent,e0,e1,e2,e3,a,b,c,phi) { 2076 var sinphi, sin2ph, tanphi, ml, mlp, con1, con2, con3, dphi, i; 2077 2078 phi = a; 2079 for (i = 1; i <= 15; i++) { 2080 sinphi = Math.sin(phi); 2081 tanphi = Math.tan(phi); 2082 c = tanphi * Math.sqrt (1.0 - eccent * sinphi * sinphi); 2083 sin2ph = Math.sin (2.0 * phi); 2084 /* 2085 ml = e0 * *phi - e1 * sin2ph + e2 * sin (4.0 * *phi); 2086 mlp = e0 - 2.0 * e1 * cos (2.0 * *phi) + 4.0 * e2 * cos (4.0 * *phi); 2087 */ 2088 ml = e0 * phi - e1 * sin2ph + e2 * Math.sin (4.0 * phi) - e3 * Math.sin (6.0 * phi); 2089 mlp = e0 - 2.0 * e1 * Math.cos (2.0 * phi) + 4.0 * e2 * Math.cos (4.0 * phi) - 6.0 * e3 * Math.cos (6.0 * phi); 2090 con1 = 2.0 * ml + c * (ml * ml + b) - 2.0 * a * (c * ml + 1.0); 2091 con2 = eccent * sin2ph * (ml * ml + b - 2.0 * a * ml) / (2.0 *c); 2092 con3 = 2.0 * (a - ml) * (c * mlp - 2.0 / sin2ph) - 2.0 * mlp; 2093 dphi = con1 / (con2 + con3); 2094 phi += dphi; 2095 if (Math.abs(dphi) <= .0000000001 ) return(phi); 2096 } 2097 Proj4js.reportError("phi4z: No convergence"); 2098 return null; 2099} 2100 2101 2102/* Function to compute the constant e4 from the input of the eccentricity 2103 of the spheroid, x. This constant is used in the Polar Stereographic 2104 projection. 2105--------------------------------------------------------------------*/ 2106function e4fn(x) { 2107 var con, com; 2108 con = 1.0 + x; 2109 com = 1.0 - x; 2110 return (Math.sqrt((Math.pow(con,con))*(Math.pow(com,com)))); 2111} 2112 2113 2114 2115 2116 2117/******************************************************************************* 2118NAME POLYCONIC 2119 2120PURPOSE: Transforms input longitude and latitude to Easting and 2121 Northing for the Polyconic projection. The 2122 longitude and latitude must be in radians. The Easting 2123 and Northing values will be returned in meters. 2124 2125PROGRAMMER DATE 2126---------- ---- 2127T. Mittan Mar, 1993 2128 2129ALGORITHM REFERENCES 2130 21311. Snyder, John P., "Map Projections--A Working Manual", U.S. Geological 2132 Survey Professional Paper 1395 (Supersedes USGS Bulletin 1532), United 2133 State Government Printing Office, Washington D.C., 1987. 2134 21352. Snyder, John P. and Voxland, Philip M., "An Album of Map Projections", 2136 U.S. Geological Survey Professional Paper 1453 , United State Government 2137 Printing Office, Washington D.C., 1989. 2138*******************************************************************************/ 2139 2140Proj4js.Proj.poly = { 2141 2142 /* Initialize the POLYCONIC projection 2143 ----------------------------------*/ 2144 init: function() { 2145 var temp; /* temporary variable */ 2146 if (this.lat0 == 0) this.lat0 = 90;//this.lat0 ca 2147 2148 /* Place parameters in static storage for common use 2149 -------------------------------------------------*/ 2150 this.temp = this.b / this.a; 2151 this.es = 1.0 - Math.pow(this.temp,2);// devait etre dans tmerc.js mais n y est pas donc je commente sinon retour de valeurs nulles 2152 this.e = Math.sqrt(this.es); 2153 this.e0 = Proj4js.common.e0fn(this.es); 2154 this.e1 = Proj4js.common.e1fn(this.es); 2155 this.e2 = Proj4js.common.e2fn(this.es); 2156 this.e3 = Proj4js.common.e3fn(this.es); 2157 this.ml0 = Proj4js.common.mlfn(this.e0, this.e1,this.e2, this.e3, this.lat0);//si que des zeros le calcul ne se fait pas 2158 //if (!this.ml0) {this.ml0=0;} 2159 }, 2160 2161 2162 /* Polyconic forward equations--mapping lat,long to x,y 2163 ---------------------------------------------------*/ 2164 forward: function(p) { 2165 var sinphi, cosphi; /* sin and cos value */ 2166 var al; /* temporary values */ 2167 var c; /* temporary values */ 2168 var con, ml; /* cone constant, small m */ 2169 var ms; /* small m */ 2170 var x,y; 2171 2172 var lon=p.x; 2173 var lat=p.y; 2174 2175 con = Proj4js.common.adjust_lon(lon - this.long0); 2176 if (Math.abs(lat) <= .0000001) { 2177 x = this.x0 + this.a * con; 2178 y = this.y0 - this.a * this.ml0; 2179 } else { 2180 sinphi = Math.sin(lat); 2181 cosphi = Math.cos(lat); 2182 2183 ml = Proj4js.common.mlfn(this.e0, this.e1, this.e2, this.e3, lat); 2184 ms = Proj4js.common.msfnz(this.e,sinphi,cosphi); 2185 con = sinphi; 2186 x = this.x0 + this.a * ms * Math.sin(con)/sinphi; 2187 y = this.y0 + this.a * (ml - this.ml0 + ms * (1.0 - Math.cos(con))/sinphi); 2188 } 2189 2190 p.x=x; 2191 p.y=y; 2192 return p; 2193 }, 2194 2195 2196 /* Inverse equations 2197 -----------------*/ 2198 inverse: function(p) { 2199 var sin_phi, cos_phi; /* sin and cos value */ 2200 var al; /* temporary values */ 2201 var b; /* temporary values */ 2202 var c; /* temporary values */ 2203 var con, ml; /* cone constant, small m */ 2204 var iflg; /* error flag */ 2205 var lon,lat; 2206 p.x -= this.x0; 2207 p.y -= this.y0; 2208 al = this.ml0 + p.y/this.a; 2209 iflg = 0; 2210 2211 if (Math.abs(al) <= .0000001) { 2212 lon = p.x/this.a + this.long0; 2213 lat = 0.0; 2214 } else { 2215 b = al * al + (p.x/this.a) * (p.x/this.a); 2216 iflg = phi4z(this.es,this.e0,this.e1,this.e2,this.e3,this.al,b,c,lat); 2217 if (iflg != 1) return(iflg); 2218 lon = Proj4js.common.adjust_lon((Proj4js.common.asinz(p.x * c / this.a) / Math.sin(lat)) + this.long0); 2219 } 2220 2221 p.x=lon; 2222 p.y=lat; 2223 return p; 2224 } 2225}; 2226 2227 2228 2229/* ====================================================================== 2230 projCode/equi.js 2231 ====================================================================== */ 2232
2233/******************************************************************************* 2234NAME EQUIRECTANGULAR 2235 2236PURPOSE: Transforms input longitude and latitude to Easting and 2237 Northing for the Equirectangular projection. The 2238 longitude and latitude must be in radians. The Easting 2239 and Northing values will be returned in meters. 2240 2241PROGRAMMER DATE 2242---------- ---- 2243T. Mittan Mar, 1993 2244 2245ALGORITHM REFERENCES 2246 22471. Snyder, John P., "Map Projections--A Working Manual", U.S. Geological 2248 Survey Professional Paper 1395 (Supersedes USGS Bulletin 1532), United 2249 State Government Printing Office, Washington D.C., 1987. 2250 22512. Snyder, John P. and Voxland, Philip M., "An Album of Map Projections", 2252 U.S. Geological Survey Professional Paper 1453 , United State Government 2253 Printing Office, Washington D.C., 1989. 2254*******************************************************************************/ 2255Proj4js.Proj.equi = { 2256 2257 init: function() { 2258 if(!this.x0) this.x0=0; 2259 if(!this.y0) this.y0=0; 2260 if(!this.lat0) this.lat0=0; 2261 if(!this.long0) this.long0=0; 2262 ///this.t2; 2263 }, 2264 2265 2266 2267/* Equirectangular forward equations--mapping lat,long to x,y 2268 ---------------------------------------------------------*/ 2269 forward: function(p) { 2270 2271 var lon=p.x; 2272 var lat=p.y; 2273 2274 var dlon = Proj4js.common.adjust_lon(lon - this.long0); 2275 var x = this.x0 +this. a * dlon *Math.cos(this.lat0); 2276 var y = this.y0 + this.a * lat; 2277 2278 this.t1=x; 2279 this.t2=Math.cos(this.lat0); 2280 p.x=x; 2281 p.y=y; 2282 return p; 2283 }, //equiFwd() 2284 2285 2286 2287/* Equirectangular inverse equations--mapping x,y to lat/long 2288 ---------------------------------------------------------*/ 2289 inverse: function(p) { 2290 2291 p.x -= this.x0; 2292 p.y -= this.y0; 2293 var lat = p.y /this. a; 2294 2295 if ( Math.abs(lat) > Proj4js.common.HALF_PI) { 2296 Proj4js.reportError("equi:Inv:DataError"); 2297 } 2298 var lon = Proj4js.common.adjust_lon(this.long0 + p.x / (this.a * Math.cos(this.lat0))); 2299 p.x=lon; 2300 p.y=lat; 2301 }//equiInv() 2302}; 2303 2304 2305/* ====================================================================== 2306 projCode/merc.js 2307 ====================================================================== */ 2308 2309/******************************************************************************* 2310NAME MERCATOR 2311 2312PURPOSE: Transforms input longitude and latitude to Easting and 2313 Northing for the Mercator projection. The 2314 longitude and latitude must be in radians. The Easting 2315 and Northing values will be returned in meters. 2316 2317PROGRAMMER DATE 2318---------- ---- 2319D. Steinwand, EROS Nov, 1991 2320T. Mittan Mar, 1993 2321 2322ALGORITHM REFERENCES 2323 23241. Snyder, John P., "Map Projections--A Working Manual", U.S. Geological 2325 Survey Professional Paper 1395 (Supersedes USGS Bulletin 1532), United 2326 State Government Printing Office, Washington D.C., 1987. 2327 23282. Snyder, John P. and Voxland, Philip M., "An Album of Map Projections", 2329 U.S. Geological Survey Professional Paper 1453 , United State Government 2330 Printing Office, Washington D.C., 1989. 2331*******************************************************************************/ 2332 2333//static double r_major = a; /* major axis */ 2334//static double r_minor = b; /* minor axis */ 2335//static double lon_center = long0; /* Center longitude (projection center) */ 2336//static double lat_origin = lat0; /* center latitude */ 2337//static double e,es; /* eccentricity constants */ 2338//static double m1; /* small value m */ 2339//static double false_northing = y0; /* y offset in meters */ 2340//static double false_easting = x0; /* x offset in meters */ 2341//scale_fact = k0 2342 2343Proj4js.Proj.merc = { 2344 init : function() { 2345 //?this.temp = this.r_minor / this.r_major; 2346 //this.temp = this.b / this.a; 2347 //this.es = 1.0 - Math.sqrt(this.temp); 2348 //this.e = Math.sqrt( this.es ); 2349 //?this.m1 = Math.cos(this.lat_origin) / (Math.sqrt( 1.0 - this.es * Math.sin(this.lat_origin) * Math.sin(this.lat_origin))); 2350 //this.m1 = Math.cos(0.0) / (Math.sqrt( 1.0 - this.es * Math.sin(0.0) * Math.sin(0.0))); 2351 if (this.lat_ts) { 2352 if (this.sphere) { 2353 this.k0 = Math.cos(this.lat_ts); 2354 } else { 2355 this.k0 = Proj4js.common.msfnz(this.es, Math.sin(this.lat_ts), Math.cos(this.lat_ts)); 2356 } 2357 } 2358 }, 2359 2360/* Mercator forward equations--mapping lat,long to x,y 2361 --------------------------------------------------*/ 2362 2363 forward : function(p) { 2364 //alert("ll2m coords : "+coords); 2365 var lon = p.x; 2366 var lat = p.y; 2367 // convert to radians 2368 if ( lat*Proj4js.common.R2D > 90.0 && 2369 lat*Proj4js.common.R2D < -90.0 && 2370 lon*Proj4js.common.R2D > 180.0 && 2371 lon*Proj4js.common.R2D < -180.0) { 2372 Proj4js.reportError("merc:forward: llInputOutOfRange: "+ lon +" : " + lat); 2373 return null; 2374 } 2375 2376 var x,y; 2377 if(Math.abs( Math.abs(lat) - Proj4js.common.HALF_PI) <= Proj4js.common.EPSLN) { 2378 Proj4js.reportError("merc:forward: ll2mAtPoles"); 2379 return null; 2380 } else { 2381 if (this.sphere) { 2382 x = this.x0 + this.a * this.k0 * Proj4js.common.adjust_lon(lon - this.long0); 2383 y = this.y0 + this.a * this.k0 * Math.log(Math.tan(Proj4js.common.FORTPI + 0.5*lat)); 2384 } else { 2385 var sinphi = Math.sin(lat); 2386 var ts = Proj4js.common.tsfnz(this.e,lat,sinphi); 2387 x = this.x0 + this.a * this.k0 * Proj4js.common.adjust_lon(lon - this.long0); 2388 y = this.y0 - this.a * this.k0 * Math.log(ts); 2389 } 2390 p.x = x; 2391 p.y = y; 2392 return p; 2393 } 2394 }, 2395 2396 2397 /* Mercator inverse equations--mapping x,y to lat/long 2398 --------------------------------------------------*/ 2399 inverse : function(p) { 2400 2401 var x = p.x - this.x0; 2402 var y = p.y - this.y0; 2403 var lon,lat; 2404 2405 if (this.sphere) { 2406 lat = Proj4js.common.HALF_PI - 2.0 * Math.atan(Math.exp(-y / this.a * this.k0)); 2407 } else { 2408 var ts = Math.exp(-y / (this.a * this.k0)); 2409 lat = Proj4js.common.phi2z(this.e,ts); 2410 if(lat == -9999) { 2411 Proj4js.reportError("merc:inverse: lat = -9999"); 2412 return null; 2413 } 2414 } 2415 lon = Proj4js.common.adjust_lon(this.long0+ x / (this.a * this.k0)); 2416 2417 p.x = lon; 2418 p.y = lat; 2419 return p; 2420 } 2421}; 2422 2423 2424/* ====================================================================== 2425 projCode/utm.js 2426 ====================================================================== */ 2427
2428/******************************************************************************* 2429NAME TRANSVERSE MERCATOR 2430 2431PURPOSE: Transforms input longitude and latitude to Easting and 2432 Northing for the Transverse Mercator projection. The 2433 longitude and latitude must be in radians. The Easting 2434 and Northing values will be returned in meters. 2435 2436ALGORITHM REFERENCES 2437 24381. Snyder, John P., "Map Projections--A Working Manual", U.S. Geological 2439 Survey Professional Paper 1395 (Supersedes USGS Bulletin 1532), United 2440 State Government Printing Office, Washington D.C., 1987. 2441 24422. Snyder, John P. and Voxland, Philip M., "An Album of Map Projections", 2443 U.S. Geological Survey Professional Paper 1453 , United State Government 2444 Printing Office, Washington D.C., 1989. 2445*******************************************************************************/ 2446 2447 2448/** 2449 Initialize Transverse Mercator projection 2450*/ 2451 2452Proj4js.Proj.utm = { 2453 dependsOn : 'tmerc', 2454 2455 init : function() { 2456 if (!this.zone) { 2457 Proj4js.reportError("utm:init: zone must be specified for UTM"); 2458 return; 2459 } 2460 this.lat0 = 0.0; 2461 this.long0 = ((6 * Math.abs(this.zone)) - 183) * Proj4js.common.D2R; 2462 this.x0 = 500000.0; 2463 this.y0 = this.utmSouth ? 10000000.0 : 0.0; 2464 this.k0 = 0.9996; 2465 2466 Proj4js.Proj['tmerc'].init.apply(this); 2467 this.forward = Proj4js.Proj['tmerc'].forward; 2468 this.inverse = Proj4js.Proj['tmerc'].inverse; 2469 } 2470}; 2471/* ====================================================================== 2472 projCode/eqdc.js 2473 ====================================================================== */ 2474 2475/******************************************************************************* 2476NAME EQUIDISTANT CONIC 2477 2478PURPOSE: Transforms input longitude and latitude to Easting and Northing 2479 for the Equidistant Conic projection. The longitude and 2480 latitude must be in radians. The Easting and Northing values 2481 will be returned in meters. 2482 2483PROGRAMMER DATE 2484---------- ---- 2485T. Mittan Mar, 1993 2486 2487ALGORITHM REFERENCES 2488 24891. Snyder, John P., "Map Projections--A Working Manual", U.S. Geological 2490 Survey Professional Paper 1395 (Supersedes USGS Bulletin 1532), United 2491 State Government Printing Office, Washington D.C., 1987. 2492 24932. Snyder, John P. and Voxland, Philip M., "An Album of Map Projections", 2494 U.S. Geological Survey Professional Paper 1453 , United State Government 2495 Printing Office, Washington D.C., 1989. 2496*******************************************************************************/ 2497 2498/* Variables common to all subroutines in this code file 2499 -----------------------------------------------------*/ 2500 2501Proj4js.Proj.eqdc = { 2502 2503/* Initialize the Equidistant Conic projection 2504 ------------------------------------------*/ 2505 init: function() { 2506 2507 /* Place parameters in static storage for common use 2508 -------------------------------------------------*/ 2509 2510 if(!this.mode) this.mode=0;//chosen default mode 2511 this.temp = this.b / this.a; 2512 this.es = 1.0 - Math.pow(this.temp,2); 2513 this.e = Math.sqrt(this.es); 2514 this.e0 = Proj4js.common.e0fn(this.es); 2515 this.e1 = Proj4js.common.e1fn(this.es); 2516 this.e2 = Proj4js.common.e2fn(this.es); 2517 this.e3 = Proj4js.common.e3fn(this.es); 2518 2519 this.sinphi=Math.sin(this.lat1); 2520 this.cosphi=Math.cos(this.lat1); 2521 2522 this.ms1 = Proj4js.common.msfnz(this.e,this.sinphi,this.cosphi); 2523 this.ml1 = Proj4js.common.mlfn(this.e0, this.e1, this.e2,this.e3, this.lat1); 2524 2525 /* format B 2526 ---------*/ 2527 if (this.mode != 0) { 2528 if (Math.abs(this.lat1 + this.lat2) < Proj4js.common.EPSLN) { 2529 Proj4js.reportError("eqdc:Init:EqualLatitudes"); 2530 //return(81); 2531 } 2532 this.sinphi=Math.sin(this.lat2); 2533 this.cosphi=Math.cos(this.lat2); 2534 2535 this.ms2 = Proj4js.common.msfnz(this.e,this.sinphi,this.cosphi); 2536 this.ml2 = Proj4js.common.mlfn(this.e0, this.e1, this.e2, this.e3, this.lat2); 2537 if (Math.abs(this.lat1 - this.lat2) >= Proj4js.common.EPSLN) { 2538 this.ns = (this.ms1 - this.ms2) / (this.ml2 - this.ml1); 2539 } else { 2540 this.ns = this.sinphi; 2541 } 2542 } else { 2543 this.ns = this.sinphi; 2544 } 2545 this.g = this.ml1 + this.ms1/this.ns; 2546 this.ml0 = Proj4js.common.mlfn(this.e0, this.e1,this. e2, this.e3, this.lat0); 2547 this.rh = this.a * (this.g - this.ml0); 2548 }, 2549 2550 2551/* Equidistant Conic forward equations--mapping lat,long to x,y 2552 -----------------------------------------------------------*/ 2553 forward: function(p) { 2554 var lon=p.x; 2555 var lat=p.y; 2556 2557 /* Forward equations 2558 -----------------*/ 2559 var ml = Proj4js.common.mlfn(this.e0, this.e1, this.e2, this.e3, lat); 2560 var rh1 = this.a * (this.g - ml); 2561 var theta = this.ns * Proj4js.common.adjust_lon(lon - this.long0); 2562 2563 var x = this.x0 + rh1 * Math.sin(theta); 2564 var y = this.y0 + this.rh - rh1 * Math.cos(theta); 2565 p.x=x; 2566 p.y=y; 2567 return p; 2568 }, 2569 2570/* Inverse equations 2571 -----------------*/ 2572 inverse: function(p) { 2573 p.x -= this.x0; 2574 p.y = this.rh - p.y + this.y0; 2575 var con, rh1; 2576 if (this.ns >= 0) { 2577 rh1 = Math.sqrt(p.x *p.x + p.y * p.y); 2578 con = 1.0; 2579 } else {
2580 rh1 = -Math.sqrt(p.x *p. x +p. y * p.y); 2581 con = -1.0; 2582 } 2583 var theta = 0.0; 2584 if (rh1 != 0.0) theta = Math.atan2(con *p.x, con *p.y); 2585 var ml = this.g - rh1 /this.a; 2586 var lat = this.phi3z(ml,this.e0,this.e1,this.e2,this.e3); 2587 var lon = Proj4js.common.adjust_lon(this.long0 + theta / this.ns); 2588 2589 p.x=lon; 2590 p.y=lat; 2591 return p; 2592 }, 2593 2594/* Function to compute latitude, phi3, for the inverse of the Equidistant 2595 Conic projection. 2596-----------------------------------------------------------------*/ 2597 phi3z: function(ml,e0,e1,e2,e3) { 2598 var phi; 2599 var dphi; 2600 2601 phi = ml; 2602 for (var i = 0; i < 15; i++) { 2603 dphi = (ml + e1 * Math.sin(2.0 * phi) - e2 * Math.sin(4.0 * phi) + e3 * Math.sin(6.0 * phi))/ e0 - phi; 2604 phi += dphi; 2605 if (Math.abs(dphi) <= .0000000001) { 2606 return phi; 2607 } 2608 } 2609 Proj4js.reportError("PHI3Z-CONV:Latitude failed to converge after 15 iterations"); 2610 return null; 2611 } 2612 2613 2614}; 2615/* ====================================================================== 2616 projCode/tmerc.js 2617 ====================================================================== */ 2618 2619/******************************************************************************* 2620NAME TRANSVERSE MERCATOR 2621 2622PURPOSE: Transforms input longitude and latitude to Easting and 2623 Northing for the Transverse Mercator projection. The 2624 longitude and latitude must be in radians. The Easting 2625 and Northing values will be returned in meters. 2626 2627ALGORITHM REFERENCES 2628 26291. Snyder, John P., "Map Projections--A Working Manual", U.S. Geological 2630 Survey Professional Paper 1395 (Supersedes USGS Bulletin 1532), United 2631 State Government Printing Office, Washington D.C., 1987. 2632 26332. Snyder, John P. and Voxland, Philip M., "An Album of Map Projections", 2634 U.S. Geological Survey Professional Paper 1453 , United State Government 2635 Printing Office, Washington D.C., 1989. 2636*******************************************************************************/ 2637 2638 2639/** 2640 Initialize Transverse Mercator projection 2641*/ 2642 2643Proj4js.Proj.tmerc = { 2644 init : function() { 2645 this.e0 = Proj4js.common.e0fn(this.es); 2646 this.e1 = Proj4js.common.e1fn(this.es); 2647 this.e2 = Proj4js.common.e2fn(this.es); 2648 this.e3 = Proj4js.common.e3fn(this.es); 2649 this.ml0 = this.a * Proj4js.common.mlfn(this.e0, this.e1, this.e2, this.e3, this.lat0); 2650 }, 2651 2652 /** 2653 Transverse Mercator Forward - long/lat to x/y 2654 long/lat in radians 2655 */ 2656 forward : function(p) { 2657 var lon = p.x; 2658 var lat = p.y; 2659 2660 var delta_lon = Proj4js.common.adjust_lon(lon - this.long0); // Delta longitude 2661 var con; // cone constant 2662 var x, y; 2663 var sin_phi=Math.sin(lat); 2664 var cos_phi=Math.cos(lat); 2665 2666 if (this.sphere) { /* spherical form */ 2667 var b = cos_phi * Math.sin(delta_lon); 2668 if ((Math.abs(Math.abs(b) - 1.0)) < .0000000001) { 2669 Proj4js.reportError("tmerc:forward: Point projects into infinity"); 2670 return(93); 2671 } else { 2672 x = .5 * this.a * this.k0 * Math.log((1.0 + b)/(1.0 - b)); 2673 con = Math.acos(cos_phi * Math.cos(delta_lon)/Math.sqrt(1.0 - b*b)); 2674 if (lat < 0) con = - con; 2675 y = this.a * this.k0 * (con - this.lat0); 2676 } 2677 } else { 2678 var al = cos_phi * delta_lon; 2679 var als = Math.pow(al,2); 2680 var c = this.ep2 * Math.pow(cos_phi,2); 2681 var tq = Math.tan(lat); 2682 var t = Math.pow(tq,2); 2683 con = 1.0 - this.es * Math.pow(sin_phi,2); 2684 var n = this.a / Math.sqrt(con); 2685 var ml = this.a * Proj4js.common.mlfn(this.e0, this.e1, this.e2, this.e3, lat); 2686 2687 x = this.k0 * 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.ep2))) + this.x0; 2688 y = this.k0 * (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.ep2))))) + this.y0; 2689 2690 } 2691 p.x = x; p.y = y; 2692 return p; 2693 }, // tmercFwd() 2694 2695 /** 2696 Transverse Mercator Inverse - x/y to long/lat 2697 */ 2698 inverse : function(p) { 2699 var con, phi; /* temporary angles */ 2700 var delta_phi; /* difference between longitudes */ 2701 var i; 2702 var max_iter = 6; /* maximun number of iterations */ 2703 var lat, lon; 2704 2705 if (this.sphere) { /* spherical form */ 2706 var f = Math.exp(p.x/(this.a * this.k0)); 2707 var g = .5 * (f - 1/f); 2708 var temp = this.lat0 + p.y/(this.a * this.k0); 2709 var h = Math.cos(temp); 2710 con = Math.sqrt((1.0 - h * h)/(1.0 + g * g)); 2711 lat = Proj4js.common.asinz(con); 2712 if (temp < 0) 2713 lat = -lat; 2714 if ((g == 0) && (h == 0)) { 2715 lon = this.long0; 2716 } else { 2717 lon = Proj4js.common.adjust_lon(Math.atan2(g,h) + this.long0); 2718 } 2719 } else { // ellipsoidal form 2720 var x = p.x - this.x0; 2721 var y = p.y - this.y0; 2722 2723 con = (this.ml0 + y / this.k0) / this.a; 2724 phi = con; 2725 for (i=0;true;i++) { 2726 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; 2727 phi += delta_phi; 2728 if (Math.abs(delta_phi) <= Proj4js.common.EPSLN) break; 2729 if (i >= max_iter) { 2730 Proj4js.reportError("tmerc:inverse: Latitude failed to converge"); 2731 return(95); 2732 } 2733 } // for() 2734 if (Math.abs(phi) < Proj4js.common.HALF_PI) { 2735 // sincos(phi, &sin_phi, &cos_phi); 2736 var sin_phi=Math.sin(phi); 2737 var cos_phi=Math.cos(phi); 2738 var tan_phi = Math.tan(phi); 2739 var c = this.ep2 * Math.pow(cos_phi,2); 2740 var cs = Math.pow(c,2); 2741 var t = Math.pow(tan_phi,2); 2742 var ts = Math.pow(t,2); 2743 con = 1.0 - this.es * Math.pow(sin_phi,2); 2744 var n = this.a / Math.sqrt(con); 2745 var r = n * (1.0 - this.es) / con; 2746 var d = x / (n * this.k0); 2747 var ds = Math.pow(d,2); 2748 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.ep2 - ds / 30.0 * (61.0 + 90.0 * t + 298.0 * c + 45.0 * ts - 252.0 * this.ep2 - 3.0 * cs))); 2749 lon = Proj4js.common.adjust_lon(this.long0 + (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.ep2 + 24.0 * ts))) / cos_phi)); 2750 } else { 2751 lat = Proj4js.common.HALF_PI * Proj4js.common.sign(y); 2752 lon = this.long0; 2753 } 2754 } 2755 p.x = lon; 2756 p.y = lat; 2757 return p; 2758 }
2758 // tmercInv() 2759}; 2760/* ====================================================================== 2761 defs/GOOGLE.js 2762 ====================================================================== */ 2763 2764Proj4js.defs["GOOGLE"]="+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs"; 2765Proj4js.defs["EPSG:900913"]=Proj4js.defs["GOOGLE"]; 2766/* ====================================================================== 2767 projCode/gstmerc.js 2768 ====================================================================== */ 2769 2770Proj4js.Proj.gstmerc = { 2771 init : function() { 2772 2773 // array of: a, b, lon0, lat0, k0, x0, y0 2774 var temp= this.b / this.a; 2775 this.e= Math.sqrt(1.0 - temp*temp); 2776 this.lc= this.long0; 2777 this.rs= Math.sqrt(1.0+this.e*this.e*Math.pow(Math.cos(this.lat0),4.0)/(1.0-this.e*this.e)); 2778 var sinz= Math.sin(this.lat0); 2779 var pc= Math.asin(sinz/this.rs); 2780 var sinzpc= Math.sin(pc); 2781 this.cp= Proj4js.common.latiso(0.0,pc,sinzpc)-this.rs*Proj4js.common.latiso(this.e,this.lat0,sinz); 2782 this.n2= this.k0*this.a*Math.sqrt(1.0-this.e*this.e)/(1.0-this.e*this.e*sinz*sinz); 2783 this.xs= this.x0; 2784 this.ys= this.y0-this.n2*pc; 2785 2786 if (!this.title) this.title = "Gauss Schreiber transverse mercator"; 2787 }, 2788 2789 2790 // forward equations--mapping lat,long to x,y 2791 // ----------------------------------------------------------------- 2792 forward : function(p) { 2793 2794 var lon= p.x; 2795 var lat= p.y; 2796 2797 var L= this.rs*(lon-this.lc); 2798 var Ls= this.cp+(this.rs*Proj4js.common.latiso(this.e,lat,Math.sin(lat))); 2799 var lat1= Math.asin(Math.sin(L)/Proj4js.common.cosh(Ls)); 2800 var Ls1= Proj4js.common.latiso(0.0,lat1,Math.sin(lat1)); 2801 p.x= this.xs+(this.n2*Ls1); 2802 p.y= this.ys+(this.n2*Math.atan(Proj4js.common.sinh(Ls)/Math.cos(L))); 2803 return p; 2804 }, 2805 2806 // inverse equations--mapping x,y to lat/long 2807 // ----------------------------------------------------------------- 2808 inverse : function(p) { 2809 2810 var x= p.x; 2811 var y= p.y; 2812 2813 var L= Math.atan(Proj4js.common.sinh((x-this.xs)/this.n2)/Math.cos((y-this.ys)/this.n2)); 2814 var lat1= Math.asin(Math.sin((y-this.ys)/this.n2)/Proj4js.common.cosh((x-this.xs)/this.n2)); 2815 var LC= Proj4js.common.latiso(0.0,lat1,Math.sin(lat1)); 2816 p.x= this.lc+L/this.rs; 2817 p.y= Proj4js.common.invlatiso(this.e,(LC-this.cp)/this.rs); 2818 return p; 2819 } 2820 2821}; 2822/* ====================================================================== 2823 projCode/ortho.js 2824 ====================================================================== */ 2825 2826/******************************************************************************* 2827NAME ORTHOGRAPHIC 2828 2829PURPOSE: Transforms input longitude and latitude to Easting and 2830 Northing for the Orthographic projection. The 2831 longitude and latitude must be in radians. The Easting 2832 and Northing values will be returned in meters. 2833 2834PROGRAMMER DATE 2835---------- ---- 2836T. Mittan Mar, 1993 2837 2838ALGORITHM REFERENCES 2839 28401. Snyder, John P., "Map Projections--A Working Manual", U.S. Geological 2841 Survey Professional Paper 1395 (Supersedes USGS Bulletin 1532), United 2842 State Government Printing Office, Washington D.C., 1987. 2843 28442. Snyder, John P. and Voxland, Philip M., "An Album of Map Projections", 2845 U.S. Geological Survey Professional Paper 1453 , United State Government 2846 Printing Office, Washington D.C., 1989. 2847*******************************************************************************/ 2848 2849Proj4js.Proj.ortho = { 2850 2851 /* Initialize the Orthographic projection 2852 -------------------------------------*/ 2853 init: function(def) { 2854 //double temp; /* temporary variable */ 2855 2856 /* Place parameters in static storage for common use 2857 -------------------------------------------------*/; 2858 this.sin_p14=Math.sin(this.lat0); 2859 this.cos_p14=Math.cos(this.lat0); 2860 }, 2861 2862 2863 /* Orthographic forward equations--mapping lat,long to x,y 2864 ---------------------------------------------------*/ 2865 forward: function(p) { 2866 var sinphi, cosphi; /* sin and cos value */ 2867 var dlon; /* delta longitude value */ 2868 var coslon; /* cos of longitude */ 2869 var ksp; /* scale factor */ 2870 var g; 2871 var lon=p.x; 2872 var lat=p.y; 2873 /* Forward equations 2874 -----------------*/ 2875 dlon = Proj4js.common.adjust_lon(lon - this.long0); 2876 2877 sinphi=Math.sin(lat); 2878 cosphi=Math.cos(lat); 2879 2880 coslon = Math.cos(dlon); 2881 g = this.sin_p14 * sinphi + this.cos_p14 * cosphi * coslon; 2882 ksp = 1.0; 2883 if ((g > 0) || (Math.abs(g) <= Proj4js.common.EPSLN)) { 2884 var x = this.a * ksp * cosphi * Math.sin(dlon); 2885 var y = this.y0 + this.a * ksp * (this.cos_p14 * sinphi - this.sin_p14 * cosphi * coslon); 2886 } else { 2887 Proj4js.reportError("orthoFwdPointError"); 2888 } 2889 p.x=x; 2890 p.y=y; 2891 return p; 2892 }, 2893 2894 2895 inverse: function(p) { 2896 var rh; /* height above ellipsoid */ 2897 var z; /* angle */ 2898 var sinz,cosz; /* sin of z and cos of z */ 2899 var temp; 2900 var con; 2901 var lon , lat; 2902 /* Inverse equations 2903 -----------------*/ 2904 p.x -= this.x0; 2905 p.y -= this.y0; 2906 rh = Math.sqrt(p.x * p.x + p.y * p.y); 2907 if (rh > this.a + .0000001) { 2908 Proj4js.reportError("orthoInvDataError"); 2909 } 2910 z = Proj4js.common.asinz(rh / this.a); 2911 2912 sinz=Math.sin(z); 2913 cosz=Math.cos(z); 2914 2915 lon = this.long0; 2916 if (Math.abs(rh) <= Proj4js.common.EPSLN) { 2917 lat = this.lat0; 2918 } 2919 lat = Proj4js.common.asinz(cosz * this.sin_p14 + (p.y * sinz * this.cos_p14)/rh); 2920 con = Math.abs(this.lat0) - Proj4js.common.HALF_PI; 2921 if (Math.abs(con) <= Proj4js.common.EPSLN) { 2922 if (this.lat0 >= 0) { 2923 lon = Proj4js.common.adjust_lon(this.long0 + Math.atan2(p.x, -p.y)); 2924 } else { 2925 lon = Proj4js.common.adjust_lon(this.long0 -Math.atan2(-p.x, p.y)); 2926 } 2927 } 2928 con = cosz - this.sin_p14 * Math.sin(lat); 2929 p.x=lon; 2930 p.y=lat; 2931 return p; 2932 } 2933}; 2934 2935 2936/* ====================================================================== 2937 projCode/krovak.js 2938 ====================================================================== */ 2939
vendor: 5,085 bytes, lines 2940-3078
2940/** 2941 NOTES: According to EPSG the full Krovak projection method should have 2942 the following parameters. Within PROJ.4 the azimuth, and pseudo 2943 standard parallel are hardcoded in the algorithm and can't be 2944 altered from outside. The others all have defaults to match the 2945 common usage with Krovak projection. 2946 2947 lat_0 = latitude of centre of the projection 2948 2949 lon_0 = longitude of centre of the projection 2950 2951 ** = azimuth (true) of the centre line passing through the centre of the projection 2952 2953 ** = latitude of pseudo standard parallel 2954 2955 k = scale factor on the pseudo standard parallel 2956 2957 x_0 = False Easting of the centre of the projection at the apex of the cone 2958 2959 y_0 = False Northing of the centre of the projection at the apex of the cone 2960 2961 **/ 2962 2963Proj4js.Proj.krovak = { 2964 2965 init: function() { 2966 /* we want Bessel as fixed ellipsoid */ 2967 this.a = 6377397.155; 2968 this.es = 0.006674372230614; 2969 this.e = Math.sqrt(this.es); 2970 /* if latitude of projection center is not set, use 49d30'N */ 2971 if (!this.lat0) { 2972 this.lat0 = 0.863937979737193; 2973 } 2974 if (!this.long0) { 2975 this.long0 = 0.7417649320975901 - 0.308341501185665; 2976 } 2977 /* if scale not set default to 0.9999 */ 2978 if (!this.k0) { 2979 this.k0 = 0.9999; 2980 } 2981 this.s45 = 0.785398163397448; /* 45° */ 2982 this.s90 = 2 * this.s45; 2983 this.fi0 = this.lat0; /* Latitude of projection centre 49° 30' */ 2984 /* Ellipsoid Bessel 1841 a = 6377397.155m 1/f = 299.1528128, 2985 e2=0.006674372230614; 2986 */ 2987 this.e2 = this.es; /* 0.006674372230614; */ 2988 this.e = Math.sqrt(this.e2); 2989 this.alfa = Math.sqrt(1. + (this.e2 * Math.pow(Math.cos(this.fi0), 4)) / (1. - this.e2)); 2990 this.uq = 1.04216856380474; /* DU(2, 59, 42, 42.69689) */ 2991 this.u0 = Math.asin(Math.sin(this.fi0) / this.alfa); 2992 this.g = Math.pow( (1. + this.e * Math.sin(this.fi0)) / (1. - this.e * Math.sin(this.fi0)) , this.alfa * this.e / 2. ); 2993 this.k = Math.tan( this.u0 / 2. + this.s45) / Math.pow (Math.tan(this.fi0 / 2. + this.s45) , this.alfa) * this.g; 2994 this.k1 = this.k0; 2995 this.n0 = this.a * Math.sqrt(1. - this.e2) / (1. - this.e2 * Math.pow(Math.sin(this.fi0), 2)); 2996 this.s0 = 1.37008346281555; /* Latitude of pseudo standard parallel 78° 30'00" N */ 2997 this.n = Math.sin(this.s0); 2998 this.ro0 = this.k1 * this.n0 / Math.tan(this.s0); 2999 this.ad = this.s90 - this.uq; 3000 }, 3001 3002 /* ellipsoid */ 3003 /* calculate xy from lat/lon */ 3004 /* Constants, identical to inverse transform function */ 3005 forward: function(p) { 3006 var gfi, u, deltav, s, d, eps, ro; 3007 var lon = p.x; 3008 var lat = p.y; 3009 var delta_lon = Proj4js.common.adjust_lon(lon - this.long0); // Delta longitude 3010 /* Transformation */ 3011 gfi = Math.pow ( ((1. + this.e * Math.sin(lat)) / (1. - this.e * Math.sin(lat))) , (this.alfa * this.e / 2.)); 3012 u= 2. * (Math.atan(this.k * Math.pow( Math.tan(lat / 2. + this.s45), this.alfa) / gfi)-this.s45); 3013 deltav = - delta_lon * this.alfa; 3014 s = Math.asin(Math.cos(this.ad) * Math.sin(u) + Math.sin(this.ad) * Math.cos(u) * Math.cos(deltav)); 3015 d = Math.asin(Math.cos(u) * Math.sin(deltav) / Math.cos(s)); 3016 eps = this.n * d; 3017 ro = this.ro0 * Math.pow(Math.tan(this.s0 / 2. + this.s45) , this.n) / Math.pow(Math.tan(s / 2. + this.s45) , this.n); 3018 /* x and y are reverted! */ 3019 //p.y = ro * Math.cos(eps) / a; 3020 //p.x = ro * Math.sin(eps) / a; 3021 p.y = ro * Math.cos(eps) / 1.0; 3022 p.x = ro * Math.sin(eps) / 1.0; 3023 3024 if(this.czech) { 3025 p.y *= -1.0; 3026 p.x *= -1.0; 3027 } 3028 return (p); 3029 }, 3030 3031 /* calculate lat/lon from xy */ 3032 inverse: function(p) { 3033 /* Constants, identisch wie in der Umkehrfunktion */ 3034 var u, deltav, s, d, eps, ro, fi1; 3035 var ok; 3036 3037 /* Transformation */ 3038 /* revert y, x*/ 3039 var tmp = p.x; 3040 p.x=p.y; 3041 p.y=tmp; 3042 if(this.czech) { 3043 p.y *= -1.0; 3044 p.x *= -1.0; 3045 } 3046 ro = Math.sqrt(p.x * p.x + p.y * p.y); 3047 eps = Math.atan2(p.y, p.x); 3048 d = eps / Math.sin(this.s0); 3049 s = 2. * (Math.atan( Math.pow(this.ro0 / ro, 1. / this.n) * Math.tan(this.s0 / 2. + this.s45)) - this.s45); 3050 u = Math.asin(Math.cos(this.ad) * Math.sin(s) - Math.sin(this.ad) * Math.cos(s) * Math.cos(d)); 3051 deltav = Math.asin(Math.cos(s) * Math.sin(d) / Math.cos(u)); 3052 p.x = this.long0 - deltav / this.alfa; 3053 /* ITERATION FOR lat */ 3054 fi1 = u; 3055 ok = 0; 3056 var iter = 0; 3057 do { 3058 p.y = 2. * ( Math.atan( Math.pow( this.k, -1. / this.alfa) * 3059 Math.pow( Math.tan(u / 2. + this.s45) , 1. / this.alfa) * 3060 Math.pow( (1. + this.e * Math.sin(fi1)) / (1. - this.e * Math.sin(fi1)) , this.e / 2.) 3061 ) - this.s45); 3062 if (Math.abs(fi1 - p.y) < 0.0000000001) ok=1; 3063 fi1 = p.y; 3064 iter += 1; 3065 } while (ok==0 && iter < 15); 3066 if (iter >= 15) { 3067 Proj4js.reportError("PHI3Z-CONV:Latitude failed to converge after 15 iterations"); 3068 //console.log('iter:', iter); 3069 return null; 3070 } 3071 3072 return (p); 3073 } 3074}; 3075/* ====================================================================== 3076 projCode/somerc.js 3077 ====================================================================== */ 3078
3079/******************************************************************************* 3080NAME SWISS OBLIQUE MERCATOR 3081 3082PURPOSE: Swiss projection. 3083WARNING: X and Y are inverted (weird) in the swiss coordinate system. Not 3084 here, since we want X to be horizontal and Y vertical. 3085 3086ALGORITHM REFERENCES 30871. "Formules et constantes pour le Calcul pour la 3088 projection cylindrique conforme àaxe oblique et pour la transformation entre 3089 des systèmes de référence". 3090 http://www.swisstopo.admin.ch/internet/swisstopo/fr/home/topics/survey/sys/refsys/switzerland.parsysrelated1.31216.downloadList.77004.DownloadFile.tmp/swissprojectionfr.pdf 3091 3092*******************************************************************************/ 3093 3094Proj4js.Proj.somerc = { 3095 3096 init: function() { 3097 var phy0 = this.lat0; 3098 this.lambda0 = this.long0; 3099 var sinPhy0 = Math.sin(phy0); 3100 var semiMajorAxis = this.a; 3101 var invF = this.rf; 3102 var flattening = 1 / invF; 3103 var e2 = 2 * flattening - Math.pow(flattening, 2); 3104 var e = this.e = Math.sqrt(e2); 3105 this.R = this.k0 * semiMajorAxis * Math.sqrt(1 - e2) / (1 - e2 * Math.pow(sinPhy0, 2.0)); 3106 this.alpha = Math.sqrt(1 + e2 / (1 - e2) * Math.pow(Math.cos(phy0), 4.0)); 3107 this.b0 = Math.asin(sinPhy0 / this.alpha); 3108 this.K = Math.log(Math.tan(Math.PI / 4.0 + this.b0 / 2.0)) 3109 - this.alpha 3110 * Math.log(Math.tan(Math.PI / 4.0 + phy0 / 2.0)) 3111 + this.alpha 3112 * e / 2 3113 * Math.log((1 + e * sinPhy0) 3114 / (1 - e * sinPhy0)); 3115 }, 3116 3117 3118 forward: function(p) { 3119 var Sa1 = Math.log(Math.tan(Math.PI / 4.0 - p.y / 2.0)); 3120 var Sa2 = this.e / 2.0 3121 * Math.log((1 + this.e * Math.sin(p.y)) 3122 / (1 - this.e * Math.sin(p.y))); 3123 var S = -this.alpha * (Sa1 + Sa2) + this.K; 3124 3125 // spheric latitude 3126 var b = 2.0 * (Math.atan(Math.exp(S)) - Math.PI / 4.0); 3127 3128 // spheric longitude 3129 var I = this.alpha * (p.x - this.lambda0); 3130 3131 // psoeudo equatorial rotation 3132 var rotI = Math.atan(Math.sin(I) 3133 / (Math.sin(this.b0) * Math.tan(b) + 3134 Math.cos(this.b0) * Math.cos(I))); 3135 3136 var rotB = Math.asin(Math.cos(this.b0) * Math.sin(b) - 3137 Math.sin(this.b0) * Math.cos(b) * Math.cos(I)); 3138 3139 p.y = this.R / 2.0 3140 * Math.log((1 + Math.sin(rotB)) / (1 - Math.sin(rotB))) 3141 + this.y0; 3142 p.x = this.R * rotI + this.x0; 3143 return p; 3144 }, 3145 3146 inverse: function(p) { 3147 var Y = p.x - this.x0; 3148 var X = p.y - this.y0; 3149 3150 var rotI = Y / this.R; 3151 var rotB = 2 * (Math.atan(Math.exp(X / this.R)) - Math.PI / 4.0); 3152 3153 var b = Math.asin(Math.cos(this.b0) * Math.sin(rotB) 3154 + Math.sin(this.b0) * Math.cos(rotB) * Math.cos(rotI)); 3155 var I = Math.atan(Math.sin(rotI) 3156 / (Math.cos(this.b0) * Math.cos(rotI) - Math.sin(this.b0) 3157 * Math.tan(rotB))); 3158 3159 var lambda = this.lambda0 + I / this.alpha; 3160 3161 var S = 0.0; 3162 var phy = b; 3163 var prevPhy = -1000.0; 3164 var iteration = 0; 3165 while (Math.abs(phy - prevPhy) > 0.0000001) 3166 { 3167 if (++iteration > 20) 3168 { 3169 Proj4js.reportError("omercFwdInfinity"); 3170 return; 3171 } 3172 //S = Math.log(Math.tan(Math.PI / 4.0 + phy / 2.0)); 3173 S = 1.0 3174 / this.alpha 3175 * (Math.log(Math.tan(Math.PI / 4.0 + b / 2.0)) - this.K) 3176 + this.e 3177 * Math.log(Math.tan(Math.PI / 4.0 3178 + Math.asin(this.e * Math.sin(phy)) 3179 / 2.0)); 3180 prevPhy = phy; 3181 phy = 2.0 * Math.atan(Math.exp(S)) - Math.PI / 2.0; 3182 } 3183 3184 p.x = lambda; 3185 p.y = phy; 3186 return p; 3187 } 3188}; 3189/* ====================================================================== 3190 projCode/stere.js 3191 ====================================================================== */ 3192 3193 3194// Initialize the Stereographic projection 3195 3196Proj4js.Proj.stere = { 3197 ssfn_: function(phit, sinphi, eccen) { 3198 sinphi *= eccen; 3199 return (Math.tan (.5 * (Proj4js.common.HALF_PI + phit)) * Math.pow((1. - sinphi) / (1. + sinphi), .5 * eccen)); 3200 }, 3201 TOL: 1.e-8, 3202 NITER: 8, 3203 CONV: 1.e-10, 3204 S_POLE: 0, 3205 N_POLE: 1, 3206 OBLIQ: 2, 3207 EQUIT: 3, 3208 3209 init: function() { 3210 this.phits = this.lat_ts ? this.lat_ts : Proj4js.common.HALF_PI; 3211 var t = Math.abs(this.lat0); 3212 if ((Math.abs(t) - Proj4js.common.HALF_PI) < Proj4js.common.EPSLN) { 3213 this.mode = this.lat0 < 0. ? this.S_POLE : this.N_POLE; 3214 } else { 3215 this.mode = t > Proj4js.common.EPSLN ? this.OBLIQ : this.EQUIT; 3216 } 3217 this.phits = Math.abs(this.phits); 3218 if (this.es) { 3219 var X; 3220 3221 switch (this.mode) { 3222 case this.N_POLE:
3223 case this.S_POLE: 3224 if (Math.abs(this.phits - Proj4js.common.HALF_PI) < Proj4js.common.EPSLN) { 3225 this.akm1 = 2. * this.k0 / Math.sqrt(Math.pow(1+this.e,1+this.e)*Math.pow(1-this.e,1-this.e)); 3226 } else { 3227 t = Math.sin(this.phits); 3228 this.akm1 = Math.cos(this.phits) / Proj4js.common.tsfnz(this.e, this.phits, t); 3229 t *= this.e; 3230 this.akm1 /= Math.sqrt(1. - t * t); 3231 } 3232 break; 3233 case this.EQUIT: 3234 this.akm1 = 2. * this.k0; 3235 break; 3236 case this.OBLIQ: 3237 t = Math.sin(this.lat0); 3238 X = 2. * Math.atan(this.ssfn_(this.lat0, t, this.e)) - Proj4js.common.HALF_PI; 3239 t *= this.e; 3240 this.akm1 = 2. * this.k0 * Math.cos(this.lat0) / Math.sqrt(1. - t * t); 3241 this.sinX1 = Math.sin(X); 3242 this.cosX1 = Math.cos(X); 3243 break; 3244 } 3245 } else { 3246 switch (this.mode) { 3247 case this.OBLIQ: 3248 this.sinph0 = Math.sin(this.lat0); 3249 this.cosph0 = Math.cos(this.lat0); 3250 case this.EQUIT: 3251 this.akm1 = 2. * this.k0; 3252 break; 3253 case this.S_POLE: 3254 case this.N_POLE: 3255 this.akm1 = Math.abs(this.phits - Proj4js.common.HALF_PI) >= Proj4js.common.EPSLN ? 3256 Math.cos(this.phits) / Math.tan(Proj4js.common.FORTPI - .5 * this.phits) : 3257 2. * this.k0 ; 3258 break; 3259 } 3260 } 3261 }, 3262 3263// Stereographic forward equations--mapping lat,long to x,y 3264 forward: function(p) { 3265 var lon = p.x; 3266 lon = Proj4js.common.adjust_lon(lon - this.long0); 3267 var lat = p.y; 3268 var x, y; 3269 3270 if (this.sphere) { 3271 var sinphi, cosphi, coslam, sinlam; 3272 3273 sinphi = Math.sin(lat); 3274 cosphi = Math.cos(lat); 3275 coslam = Math.cos(lon); 3276 sinlam = Math.sin(lon); 3277 switch (this.mode) { 3278 case this.EQUIT: 3279 y = 1. + cosphi * coslam; 3280 if (y <= Proj4js.common.EPSLN) { 3281 Proj4js.reportError("stere:forward:Equit"); 3282 } 3283 y = this.akm1 / y; 3284 x = y * cosphi * sinlam; 3285 y *= sinphi; 3286 break; 3287 case this.OBLIQ: 3288 y = 1. + this.sinph0 * sinphi + this.cosph0 * cosphi * coslam; 3289 if (y <= Proj4js.common.EPSLN) { 3290 Proj4js.reportError("stere:forward:Obliq"); 3291 } 3292 y = this.akm1 / y; 3293 x = y * cosphi * sinlam; 3294 y *= this.cosph0 * sinphi - this.sinph0 * cosphi * coslam; 3295 break; 3296 case this.N_POLE: 3297 coslam = -coslam; 3298 lat = -lat; 3299 //Note no break here so it conitnues through S_POLE 3300 case this.S_POLE: 3301 if (Math.abs(lat - Proj4js.common.HALF_PI) < this.TOL) { 3302 Proj4js.reportError("stere:forward:S_POLE"); 3303 } 3304 y = this.akm1 * Math.tan(Proj4js.common.FORTPI + .5 * lat); 3305 x = sinlam * y; 3306 y *= coslam; 3307 break; 3308 } 3309 } else { 3310 coslam = Math.cos(lon); 3311 sinlam = Math.sin(lon); 3312 sinphi = Math.sin(lat); 3313 var sinX, cosX; 3314 if (this.mode == this.OBLIQ || this.mode == this.EQUIT) { 3315 var Xt = 2. * Math.atan(this.ssfn_(lat, sinphi, this.e)); 3316 sinX = Math.sin(Xt - Proj4js.common.HALF_PI); 3317 cosX = Math.cos(Xt); 3318 } 3319 switch (this.mode) { 3320 case this.OBLIQ: 3321 var A = this.akm1 / (this.cosX1 * (1. + this.sinX1 * sinX + this.cosX1 * cosX * coslam)); 3322 y = A * (this.cosX1 * sinX - this.sinX1 * cosX * coslam); 3323 x = A * cosX; 3324 break; 3325 case this.EQUIT: 3326 var A = 2. * this.akm1 / (1. + cosX * coslam); 3327 y = A * sinX; 3328 x = A * cosX; 3329 break; 3330 case this.S_POLE: 3331 lat = -lat; 3332 coslam = - coslam; 3333 sinphi = -sinphi; 3334 case this.N_POLE: 3335 x = this.akm1 * Proj4js.common.tsfnz(this.e, lat, sinphi); 3336 y = - x * coslam; 3337 break; 3338 } 3339 x = x * sinlam; 3340 } 3341 p.x = x*this.a + this.x0; 3342 p.y = y*this.a + this.y0; 3343 return p; 3344 }, 3345 3346 3347//* Stereographic inverse equations--mapping x,y to lat/long 3348 inverse: function(p) { 3349 var x = (p.x - this.x0)/this.a; /* descale and de-offset */ 3350 var y = (p.y - this.y0)/this.a; 3351 var lon, lat; 3352 3353 var cosphi, sinphi, tp=0.0, phi_l=0.0, rho, halfe=0.0, pi2=0.0; 3354 var i; 3355 3356 if (this.sphere) { 3357 var c, rh, sinc, cosc; 3358 3359 rh = Math.sqrt(x*x + y*y); 3360 c = 2. * Math.atan(rh / this.akm1); 3361 sinc = Math.sin(c); 3362 cosc = Math.cos(c); 3363 lon = 0.; 3364 switch (this.mode) { 3365 case this.EQUIT: 3366 if (Math.abs(rh) <= Proj4js.common.EPSLN) { 3367 lat = 0.; 3368 } else { 3369 lat = Math.asin(y * sinc / rh); 3370 } 3371 if (cosc != 0. || x != 0.) lon = Math.atan2(x * sinc, cosc * rh); 3372 break; 3373 case this.OBLIQ: 3374 if (Math.abs(rh) <= Proj4js.common.EPSLN) { 3375 lat = this.phi0; 3376 } else { 3377 lat = Math.asin(cosc * this.sinph0 + y * sinc * this.cosph0 / rh); 3378 } 3379 c = cosc - this.sinph0 * Math.sin(lat); 3380 if (c != 0. || x != 0.) { 3381 lon = Math.atan2(x * sinc * this.cosph0, c * rh); 3382 } 3383 break; 3384 case this.N_POLE: 3385 y = -y; 3386 case this.S_POLE: 3387 if (Math.abs(rh) <= Proj4js.common.EPSLN) { 3388 lat = this.phi0; 3389 } else { 3390 lat = Math.asin(this.mode == this.S_POLE ? -cosc : cosc); 3391 } 3392 lon = (x == 0. && y == 0.) ? 0. : Math.atan2(x, y); 3393 break; 3394 } 3395 p.x = Proj4js.common.adjust_lon(lon + this.long0); 3396 p.y = lat; 3397 } else { 3398 rho = Math.sqrt(x*x + y*y); 3399 switch (this.mode) { 3400 case this.OBLIQ: 3401 case this.EQUIT: 3402 tp = 2. * Math.atan2(rho * this.cosX1 , this.akm1); 3403 cosphi = Math.cos(tp); 3404 sinphi = Math.sin(tp); 3405 if( rho == 0.0 ) { 3406 phi_l = Math.asin(cosphi * this.sinX1); 3407 } else { 3408 phi_l = Math.asin(cosphi * this.sinX1 + (y * sinphi * this.cosX1 / rho)); 3409 } 3410 3411 tp = Math.tan(.5 * (Proj4js.common.HALF_PI + phi_l)); 3412 x *= sinphi; 3413 y = rho * this.cosX1 * cosphi - y * this.sinX1* sinphi; 3414 pi2 = Proj4js.common.HALF_PI; 3415 halfe = .5 * this.e; 3416 break; 3417 case this.N_POLE: 3418 y = -y; 3419 case this.S_POLE: 3420 tp = - rho / this.akm1;
3421 phi_l = Proj4js.common.HALF_PI - 2. * Math.atan(tp); 3422 pi2 = -Proj4js.common.HALF_PI; 3423 halfe = -.5 * this.e; 3424 break; 3425 } 3426 for (i = this.NITER; i--; phi_l = lat) { //check this 3427 sinphi = this.e * Math.sin(phi_l); 3428 lat = 2. * Math.atan(tp * Math.pow((1.+sinphi)/(1.-sinphi), halfe)) - pi2; 3429 if (Math.abs(phi_l - lat) < this.CONV) { 3430 if (this.mode == this.S_POLE) lat = -lat; 3431 lon = (x == 0. && y == 0.) ? 0. : Math.atan2(x, y); 3432 p.x = Proj4js.common.adjust_lon(lon + this.long0); 3433 p.y = lat; 3434 return p; 3435 } 3436 } 3437 } 3438 } 3439}; 3440/* ====================================================================== 3441 projCode/nzmg.js 3442 ====================================================================== */ 3443 3444/******************************************************************************* 3445NAME NEW ZEALAND MAP GRID 3446 3447PURPOSE: Transforms input longitude and latitude to Easting and 3448 Northing for the New Zealand Map Grid projection. The 3449 longitude and latitude must be in radians. The Easting 3450 and Northing values will be returned in meters. 3451 3452 3453ALGORITHM REFERENCES 3454 34551. Department of Land and Survey Technical Circular 1973/32 3456 http://www.linz.govt.nz/docs/miscellaneous/nz-map-definition.pdf 3457 34582. OSG Technical Report 4.1 3459 http://www.linz.govt.nz/docs/miscellaneous/nzmg.pdf 3460 3461 3462IMPLEMENTATION NOTES 3463 3464The two references use different symbols for the calculated values. This 3465implementation uses the variable names similar to the symbols in reference [1]. 3466 3467The alogrithm uses different units for delta latitude and delta longitude. 3468The delta latitude is assumed to be in units of seconds of arc x 10^-5. 3469The delta longitude is the usual radians. Look out for these conversions. 3470 3471The algorithm is described using complex arithmetic. There were three 3472options: 3473 * find and use a Javascript library for complex arithmetic 3474 * write my own complex library 3475 * expand the complex arithmetic by hand to simple arithmetic 3476 3477This implementation has expanded the complex multiplication operations 3478into parallel simple arithmetic operations for the real and imaginary parts. 3479The imaginary part is way over to the right of the display; this probably 3480violates every coding standard in the world, but, to me, it makes it much 3481more obvious what is going on. 3482 3483The following complex operations are used: 3484 - addition 3485 - multiplication 3486 - division 3487 - complex number raised to integer power 3488 - summation 3489 3490A summary of complex arithmetic operations: 3491 (from http://en.wikipedia.org/wiki/Complex_arithmetic) 3492 addition: (a + bi) + (c + di) = (a + c) + (b + d)i 3493 subtraction: (a + bi) - (c + di) = (a - c) + (b - d)i 3494 multiplication: (a + bi) x (c + di) = (ac - bd) + (bc + ad)i 3495 division: (a + bi) / (c + di) = [(ac + bd)/(cc + dd)] + [(bc - ad)/(cc + dd)]i 3496 3497The algorithm needs to calculate summations of simple and complex numbers. This is 3498implemented using a for-loop, pre-loading the summed value to zero. 3499 3500The algorithm needs to calculate theta^2, theta^3, etc while doing a summation. 3501There are three possible implementations: 3502 - use Math.pow in the summation loop - except for complex numbers 3503 - precalculate the values before running the loop 3504 - calculate theta^n = theta^(n-1) * theta during the loop 3505This implementation uses the third option for both real and complex arithmetic. 3506 3507For example 3508 psi_n = 1; 3509 sum = 0; 3510 for (n = 1; n <=6; n++) { 3511 psi_n1 = psi_n * psi; // calculate psi^(n+1) 3512 psi_n = psi_n1; 3513 sum = sum + A[n] * psi_n; 3514 } 3515 3516 3517TEST VECTORS 3518 3519NZMG E, N: 2487100.638 6751049.719 metres 3520NZGD49 long, lat: 172.739194 -34.444066 degrees 3521 3522NZMG E, N: 2486533.395 6077263.661 metres 3523NZGD49 long, lat: 172.723106 -40.512409 degrees 3524 3525NZMG E, N: 2216746.425 5388508.765 metres 3526NZGD49 long, lat: 169.172062 -46.651295 degrees 3527 3528Note that these test vectors convert from NZMG metres to lat/long referenced 3529to NZGD49, not the more usual WGS84. The difference is about 70m N/S and about 353010m E/W. 3531 3532These test vectors are provided in reference [1]. Many more test 3533vectors are available in 3534 http://www.linz.govt.nz/docs/topography/topographicdata/placenamesdatabase/nznamesmar08.zip 3535which is a catalog of names on the 260-series maps. 3536 3537 3538EPSG CODES 3539 3540NZMG EPSG:27200 3541NZGD49 EPSG:4272 3542 3543http://spatialreference.org/ defines these as 3544 Proj4js.defs["EPSG:4272"] = "+proj=longlat +ellps=intl +datum=nzgd49 +no_defs "; 3545 Proj4js.defs["EPSG:27200"] = "+proj=nzmg +lat_0=-41 +lon_0=173 +x_0=2510000 +y_0=6023150 +ellps=intl +datum=nzgd49 +units=m +no_defs "; 3546 3547 3548LICENSE 3549 Copyright: Stephen Irons 2008 3550 Released under terms of the LGPL as per: http://www.gnu.org/copyleft/lesser.html 3551 3552*******************************************************************************/ 3553 3554 3555/** 3556 Initialize New Zealand Map Grip projection 3557*/ 3558 3559Proj4js.Proj.nzmg = { 3560 3561 /** 3562 * iterations: Number of iterations to refine inverse transform. 3563 * 0 -> km accuracy 3564 * 1 -> m accuracy -- suitable for most mapping applications 3565 * 2 -> mm accuracy 3566 */ 3567 iterations: 1, 3568 3569 init : function() { 3570 this.A = new Array(); 3571 this.A[1] = +0.6399175073; 3572 this.A[2] = -0.1358797613; 3573 this.A[3] = +0.063294409; 3574 this.A[4] = -0.02526853; 3575 this.A[5] = +0.0117879; 3576 this.A[6] = -0.0055161; 3577 this.A[7] = +0.0026906; 3578 this.A[8] = -0.001333; 3579 this.A[9] = +0.00067; 3580 this.A[10] = -0.00034; 3581 3582 this.B_re = new Array(); this.B_im = new Array(); 3583 this.B_re[1] = +0.7557853228; this.B_im[1] = 0.0; 3584 this.B_re[2] = +0.249204646; this.B_im[2] = +0.003371507; 3585 this.B_re[3] = -0.001541739; this.B_im[3] = +0.041058560; 3586 this.B_re[4] = -0.10162907; this.B_im[4] = +0.01727609; 3587 this.B_re[5] = -0.26623489; this.B_im[5] = -0.36249218; 3588 this.B_re[6] = -0.6870983; this.B_im[6] = -1.1651967; 3589 3590 this.C_re = new Array(); this.C_im = new Array(); 3591 this.C_re[1] = +1.3231270439; this.C_im[1] = 0.0; 3592 this.C_re[2] = -0.577245789; this.C_im[2] = -0.007809598; 3593 this.C_re[3] = +0.508307513; this.C_im[3] = -0.112208952; 3594 this.C_re[4] = -0.15094762; this.C_im[4] = +0.18200602; 3595 this.C_re[5] = +1.01418179; this.C_im[5] = +1.64497696; 3596 this.C_re[6] = +1.9660549; this.C_im[6] = +2.5127645; 3597 3598 this.D = new Array(); 3599 this.D[1] = +1.5627014243; 3600 this.D[2] = +0.5185406398; 3601 this.D[3] = -0.03333098; 3602 this.D[4] = -0.1052906; 3603 this.D[5] = -0.0368594;
3604 this.D[6] = +0.007317; 3605 this.D[7] = +0.01220; 3606 this.D[8] = +0.00394; 3607 this.D[9] = -0.0013; 3608 }, 3609 3610 /** 3611 New Zealand Map Grid Forward - long/lat to x/y 3612 long/lat in radians 3613 */ 3614 forward : function(p) { 3615 var lon = p.x; 3616 var lat = p.y; 3617 3618 var delta_lat = lat - this.lat0; 3619 var delta_lon = lon - this.long0; 3620 3621 // 1. Calculate d_phi and d_psi ... // and d_lambda 3622 // For this algorithm, delta_latitude is in seconds of arc x 10-5, so we need to scale to those units. Longitude is radians. 3623 var d_phi = delta_lat / Proj4js.common.SEC_TO_RAD * 1E-5; var d_lambda = delta_lon; 3624 var d_phi_n = 1; // d_phi^0 3625 3626 var d_psi = 0; 3627 for (var n = 1; n <= 10; n++) { 3628 d_phi_n = d_phi_n * d_phi; 3629 d_psi = d_psi + this.A[n] * d_phi_n; 3630 } 3631 3632 // 2. Calculate theta 3633 var th_re = d_psi; var th_im = d_lambda; 3634 3635 // 3. Calculate z 3636 var th_n_re = 1; var th_n_im = 0; // theta^0 3637 var th_n_re1; var th_n_im1; 3638 3639 var z_re = 0; var z_im = 0; 3640 for (var n = 1; n <= 6; n++) { 3641 th_n_re1 = th_n_re*th_re - th_n_im*th_im; th_n_im1 = th_n_im*th_re + th_n_re*th_im; 3642 th_n_re = th_n_re1; th_n_im = th_n_im1; 3643 z_re = z_re + this.B_re[n]*th_n_re - this.B_im[n]*th_n_im; z_im = z_im + this.B_im[n]*th_n_re + this.B_re[n]*th_n_im; 3644 } 3645 3646 // 4. Calculate easting and northing 3647 p.x = (z_im * this.a) + this.x0; 3648 p.y = (z_re * this.a) + this.y0; 3649 3650 return p; 3651 }, 3652 3653 3654 /** 3655 New Zealand Map Grid Inverse - x/y to long/lat 3656 */ 3657 inverse : function(p) { 3658 3659 var x = p.x; 3660 var y = p.y; 3661 3662 var delta_x = x - this.x0; 3663 var delta_y = y - this.y0; 3664 3665 // 1. Calculate z 3666 var z_re = delta_y / this.a; var z_im = delta_x / this.a; 3667 3668 // 2a. Calculate theta - first approximation gives km accuracy 3669 var z_n_re = 1; var z_n_im = 0; // z^0 3670 var z_n_re1; var z_n_im1; 3671 3672 var th_re = 0; var th_im = 0; 3673 for (var n = 1; n <= 6; n++) { 3674 z_n_re1 = z_n_re*z_re - z_n_im*z_im; z_n_im1 = z_n_im*z_re + z_n_re*z_im; 3675 z_n_re = z_n_re1; z_n_im = z_n_im1; 3676 th_re = th_re + this.C_re[n]*z_n_re - this.C_im[n]*z_n_im; th_im = th_im + this.C_im[n]*z_n_re + this.C_re[n]*z_n_im; 3677 } 3678 3679 // 2b. Iterate to refine the accuracy of the calculation 3680 // 0 iterations gives km accuracy 3681 // 1 iteration gives m accuracy -- good enough for most mapping applications 3682 // 2 iterations bives mm accuracy 3683 for (var i = 0; i < this.iterations; i++) { 3684 var th_n_re = th_re; var th_n_im = th_im; 3685 var th_n_re1; var th_n_im1; 3686 3687 var num_re = z_re; var num_im = z_im; 3688 for (var n = 2; n <= 6; n++) { 3689 th_n_re1 = th_n_re*th_re - th_n_im*th_im; th_n_im1 = th_n_im*th_re + th_n_re*th_im; 3690 th_n_re = th_n_re1; th_n_im = th_n_im1; 3691 num_re = num_re + (n-1)*(this.B_re[n]*th_n_re - this.B_im[n]*th_n_im); num_im = num_im + (n-1)*(this.B_im[n]*th_n_re + this.B_re[n]*th_n_im); 3692 } 3693 3694 th_n_re = 1; th_n_im = 0; 3695 var den_re = this.B_re[1]; var den_im = this.B_im[1]; 3696 for (var n = 2; n <= 6; n++) { 3697 th_n_re1 = th_n_re*th_re - th_n_im*th_im; th_n_im1 = th_n_im*th_re + th_n_re*th_im; 3698 th_n_re = th_n_re1; th_n_im = th_n_im1; 3699 den_re = den_re + n * (this.B_re[n]*th_n_re - this.B_im[n]*th_n_im); den_im = den_im + n * (this.B_im[n]*th_n_re + this.B_re[n]*th_n_im); 3700 } 3701 3702 // Complex division 3703 var den2 = den_re*den_re + den_im*den_im; 3704 th_re = (num_re*den_re + num_im*den_im) / den2; th_im = (num_im*den_re - num_re*den_im) / den2; 3705 } 3706 3707 // 3. Calculate d_phi ... // and d_lambda 3708 var d_psi = th_re; var d_lambda = th_im; 3709 var d_psi_n = 1; // d_psi^0 3710 3711 var d_phi = 0; 3712 for (var n = 1; n <= 9; n++) { 3713 d_psi_n = d_psi_n * d_psi; 3714 d_phi = d_phi + this.D[n] * d_psi_n; 3715 } 3716 3717 // 4. Calculate latitude and longitude 3718 // d_phi is calcuated in second of arc * 10^-5, so we need to scale back to radians. d_lambda is in radians. 3719 var lat = this.lat0 + (d_phi * Proj4js.common.SEC_TO_RAD * 1E5); 3720 var lon = this.long0 + d_lambda; 3721 3722 p.x = lon; 3723 p.y = lat; 3724 3725 return p; 3726 } 3727}; 3728/* ====================================================================== 3729 projCode/mill.js 3730 ====================================================================== */ 3731
3732/******************************************************************************* 3733NAME MILLER CYLINDRICAL 3734 3735PURPOSE: Transforms input longitude and latitude to Easting and 3736 Northing for the Miller Cylindrical projection. The 3737 longitude and latitude must be in radians. The Easting 3738 and Northing values will be returned in meters. 3739 3740PROGRAMMER DATE 3741---------- ---- 3742T. Mittan March, 1993 3743 3744This function was adapted from the Lambert Azimuthal Equal Area projection 3745code (FORTRAN) in the General Cartographic Transformation Package software 3746which is available from the U.S. Geological Survey National Mapping Division. 3747 3748ALGORITHM REFERENCES 3749 37501. "New Equal-Area Map Projections for Noncircular Regions", John P. Snyder, 3751 The American Cartographer, Vol 15, No. 4, October 1988, pp. 341-355. 3752 37532. Snyder, John P., "Map Projections--A Working Manual", U.S. Geological 3754 Survey Professional Paper 1395 (Supersedes USGS Bulletin 1532), United 3755 State Government Printing Office, Washington D.C., 1987. 3756 37573. "Software Documentation for GCTP General Cartographic Transformation 3758 Package", U.S. Geological Survey National Mapping Division, May 1982. 3759*******************************************************************************/ 3760 3761Proj4js.Proj.mill = { 3762 3763/* Initialize the Miller Cylindrical projection 3764 -------------------------------------------*/ 3765 init: function() { 3766 //no-op 3767 }, 3768 3769 3770 /* Miller Cylindrical forward equations--mapping lat,long to x,y 3771 ------------------------------------------------------------*/ 3772 forward: function(p) { 3773 var lon=p.x; 3774 var lat=p.y; 3775 /* Forward equations 3776 -----------------*/ 3777 var dlon = Proj4js.common.adjust_lon(lon -this.long0); 3778 var x = this.x0 + this.a * dlon; 3779 var y = this.y0 + this.a * Math.log(Math.tan((Proj4js.common.PI / 4.0) + (lat / 2.5))) * 1.25; 3780 3781 p.x=x; 3782 p.y=y; 3783 return p; 3784 },//millFwd() 3785 3786 /* Miller Cylindrical inverse equations--mapping x,y to lat/long 3787 ------------------------------------------------------------*/ 3788 inverse: function(p) { 3789 p.x -= this.x0; 3790 p.y -= this.y0; 3791 3792 var lon = Proj4js.common.adjust_lon(this.long0 + p.x /this.a); 3793 var lat = 2.5 * (Math.atan(Math.exp(0.8*p.y/this.a)) - Proj4js.common.PI / 4.0); 3794 3795 p.x=lon; 3796 p.y=lat; 3797 return p; 3798 }//millInv() 3799}; 3800/* ====================================================================== 3801 projCode/gnom.js 3802 ====================================================================== */ 3803 3804/***************************************************************************** 3805NAME GNOMONIC 3806 3807PURPOSE: Transforms input longitude and latitude to Easting and 3808 Northing for the Gnomonic Projection. 3809 Implementation based on the existing sterea and ortho 3810 implementations. 3811 3812PROGRAMMER DATE 3813---------- ---- 3814Richard Marsden November 2009 3815 3816ALGORITHM REFERENCES 3817 38181. Snyder, John P., "Flattening the Earth - Two Thousand Years of Map 3819 Projections", University of Chicago Press 1993 3820 38212. Wolfram Mathworld "Gnomonic Projection" 3822 http://mathworld.wolfram.com/GnomonicProjection.html 3823 Accessed: 12th November 2009 3824******************************************************************************/ 3825 3826Proj4js.Proj.gnom = { 3827 3828 /* Initialize the Gnomonic projection 3829 -------------------------------------*/ 3830 init: function(def) { 3831 3832 /* Place parameters in static storage for common use 3833 -------------------------------------------------*/ 3834 this.sin_p14=Math.sin(this.lat0); 3835 this.cos_p14=Math.cos(this.lat0); 3836 // Approximation for projecting points to the horizon (infinity) 3837 this.infinity_dist = 1000 * this.a; 3838 this.rc = 1; 3839 }, 3840 3841 3842 /* Gnomonic forward equations--mapping lat,long to x,y 3843 ---------------------------------------------------*/ 3844 forward: function(p) { 3845 var sinphi, cosphi; /* sin and cos value */ 3846 var dlon; /* delta longitude value */ 3847 var coslon; /* cos of longitude */ 3848 var ksp; /* scale factor */ 3849 var g; 3850 var x, y; 3851 var lon=p.x; 3852 var lat=p.y; 3853 /* Forward equations 3854 -----------------*/ 3855 dlon = Proj4js.common.adjust_lon(lon - this.long0); 3856 3857 sinphi=Math.sin(lat); 3858 cosphi=Math.cos(lat); 3859 3860 coslon = Math.cos(dlon); 3861 g = this.sin_p14 * sinphi + this.cos_p14 * cosphi * coslon; 3862 ksp = 1.0; 3863 if ((g > 0) || (Math.abs(g) <= Proj4js.common.EPSLN)) { 3864 x = this.x0 + this.a * ksp * cosphi * Math.sin(dlon) / g; 3865 y = this.y0 + this.a * ksp * (this.cos_p14 * sinphi - this.sin_p14 * cosphi * coslon) / g; 3866 } else { 3867 Proj4js.reportError("orthoFwdPointError"); 3868 3869 // Point is in the opposing hemisphere and is unprojectable 3870 // We still need to return a reasonable point, so we project 3871 // to infinity, on a bearing 3872 // equivalent to the northern hemisphere equivalent 3873 // This is a reasonable approximation for short shapes and lines that 3874 // straddle the horizon.
3875 3876 x = this.x0 + this.infinity_dist * cosphi * Math.sin(dlon); 3877 y = this.y0 + this.infinity_dist * (this.cos_p14 * sinphi - this.sin_p14 * cosphi * coslon); 3878 3879 } 3880 p.x=x; 3881 p.y=y; 3882 return p; 3883 }, 3884 3885 3886 inverse: function(p) { 3887 var rh; /* Rho */ 3888 var z; /* angle */ 3889 var sinc, cosc; 3890 var c; 3891 var lon , lat; 3892 3893 /* Inverse equations 3894 -----------------*/ 3895 p.x = (p.x - this.x0) / this.a; 3896 p.y = (p.y - this.y0) / this.a; 3897 3898 p.x /= this.k0; 3899 p.y /= this.k0; 3900 3901 if ( (rh = Math.sqrt(p.x * p.x + p.y * p.y)) ) { 3902 c = Math.atan2(rh, this.rc); 3903 sinc = Math.sin(c); 3904 cosc = Math.cos(c); 3905 3906 lat = Proj4js.common.asinz(cosc*this.sin_p14 + (p.y*sinc*this.cos_p14) / rh); 3907 lon = Math.atan2(p.x*sinc, rh*this.cos_p14*cosc - p.y*this.sin_p14*sinc); 3908 lon = Proj4js.common.adjust_lon(this.long0+lon); 3909 } else { 3910 lat = this.phic0; 3911 lon = 0.0; 3912 } 3913 3914 p.x=lon; 3915 p.y=lat; 3916 return p; 3917 } 3918}; 3919 3920 3921/* ====================================================================== 3922 projCode/sinu.js 3923 ====================================================================== */ 3924 3925/******************************************************************************* 3926NAME SINUSOIDAL 3927 3928PURPOSE: Transforms input longitude and latitude to Easting and 3929 Northing for the Sinusoidal projection. The 3930 longitude and latitude must be in radians. The Easting 3931 and Northing values will be returned in meters. 3932 3933PROGRAMMER DATE 3934---------- ---- 3935D. Steinwand, EROS May, 1991 3936 3937This function was adapted from the Sinusoidal projection code (FORTRAN) in the 3938General Cartographic Transformation Package software which is available from 3939the U.S. Geological Survey National Mapping Division. 3940 3941ALGORITHM REFERENCES 3942 39431. Snyder, John P., "Map Projections--A Working Manual", U.S. Geological 3944 Survey Professional Paper 1395 (Supersedes USGS Bulletin 1532), United 3945 State Government Printing Office, Washington D.C., 1987. 3946 39472. "Software Documentation for GCTP General Cartographic Transformation 3948 Package", U.S. Geological Survey National Mapping Division, May 1982. 3949*******************************************************************************/ 3950 3951Proj4js.Proj.sinu = { 3952 3953 /* Initialize the Sinusoidal projection 3954 ------------------------------------*/ 3955 init: function() { 3956 /* Place parameters in static storage for common use 3957 -------------------------------------------------*/ 3958 3959 3960 if (!this.sphere) { 3961 this.en = Proj4js.common.pj_enfn(this.es); 3962 } else { 3963 this.n = 1.; 3964 this.m = 0.; 3965 this.es = 0; 3966 this.C_y = Math.sqrt((this.m + 1.) / this.n); 3967 this.C_x = this.C_y/(this.m + 1.); 3968 } 3969 3970 }, 3971 3972 /* Sinusoidal forward equations--mapping lat,long to x,y 3973 -----------------------------------------------------*/ 3974 forward: function(p) { 3975 var x,y,delta_lon; 3976 var lon=p.x; 3977 var lat=p.y; 3978 /* Forward equations 3979 -----------------*/ 3980 lon = Proj4js.common.adjust_lon(lon - this.long0); 3981 3982 if (this.sphere) { 3983 if (!this.m) { 3984 lat = this.n != 1. ? Math.asin(this.n * Math.sin(lat)): lat; 3985 } else { 3986 var k = this.n * Math.sin(lat); 3987 for (var i = Proj4js.common.MAX_ITER; i ; --i) { 3988 var V = (this.m * lat + Math.sin(lat) - k) / (this.m + Math.cos(lat)); 3989 lat -= V; 3990 if (Math.abs(V) < Proj4js.common.EPSLN) break; 3991 } 3992 } 3993 x = this.a * this.C_x * lon * (this.m + Math.cos(lat)); 3994 y = this.a * this.C_y * lat; 3995 3996 } else { 3997 3998 var s = Math.sin(lat); 3999 var c = Math.cos(lat); 4000 y = this.a * Proj4js.common.pj_mlfn(lat, s, c, this.en); 4001 x = this.a * lon * c / Math.sqrt(1. - this.es * s * s); 4002 } 4003 4004 p.x=x; 4005 p.y=y; 4006 return p; 4007 }, 4008 4009 inverse: function(p) { 4010 var lat,temp,lon; 4011 4012 /* Inverse equations 4013 -----------------*/ 4014 p.x -= this.x0; 4015 p.y -= this.y0; 4016 lat = p.y / this.a; 4017 4018 if (this.sphere) { 4019 4020 p.y /= this.C_y; 4021 lat = this.m ? Math.asin((this.m * p.y + Math.sin(p.y)) / this.n) : 4022 ( this.n != 1. ? Math.asin(Math.sin(p.y) / this.n) : p.y ); 4023 lon = p.x / (this.C_x * (this.m + Math.cos(p.y))); 4024 4025 } else { 4026 lat = Proj4js.common.pj_inv_mlfn(p.y/this.a, this.es, this.en) 4027 var s = Math.abs(lat); 4028 if (s < Proj4js.common.HALF_PI) { 4029 s = Math.sin(lat); 4030 temp = this.long0 + p.x * Math.sqrt(1. - this.es * s * s) /(this.a * Math.cos(lat)); 4031 //temp = this.long0 + p.x / (this.a * Math.cos(lat)); 4032 lon = Proj4js.common.adjust_lon(temp); 4033 } else if ((s - Proj4js.common.EPSLN) < Proj4js.common.HALF_PI) { 4034 lon = this.long0; 4035 } 4036 4037 } 4038 4039 p.x=lon; 4040 p.y=lat; 4041 return p; 4042 } 4043}; 4044 4045 4046/* ====================================================================== 4047 projCode/vandg.js 4048 ====================================================================== */ 4049
4050/******************************************************************************* 4051NAME VAN DER GRINTEN 4052 4053PURPOSE: Transforms input Easting and Northing to longitude and 4054 latitude for the Van der Grinten projection. The 4055 Easting and Northing must be in meters. The longitude 4056 and latitude values will be returned in radians. 4057 4058PROGRAMMER DATE 4059---------- ---- 4060T. Mittan March, 1993 4061 4062This function was adapted from the Van Der Grinten projection code 4063(FORTRAN) in the General Cartographic Transformation Package software 4064which is available from the U.S. Geological Survey National Mapping Division. 4065 4066ALGORITHM REFERENCES 4067 40681. "New Equal-Area Map Projections for Noncircular Regions", John P. Snyder, 4069 The American Cartographer, Vol 15, No. 4, October 1988, pp. 341-355. 4070 40712. Snyder, John P., "Map Projections--A Working Manual", U.S. Geological 4072 Survey Professional Paper 1395 (Supersedes USGS Bulletin 1532), United 4073 State Government Printing Office, Washington D.C., 1987. 4074 40753. "Software Documentation for GCTP General Cartographic Transformation 4076 Package", U.S. Geological Survey National Mapping Division, May 1982. 4077*******************************************************************************/ 4078 4079Proj4js.Proj.vandg = { 4080 4081/* Initialize the Van Der Grinten projection 4082 ----------------------------------------*/ 4083 init: function() { 4084 this.R = 6370997.0; //Radius of earth 4085 }, 4086 4087 forward: function(p) { 4088 4089 var lon=p.x; 4090 var lat=p.y; 4091 4092 /* Forward equations 4093 -----------------*/ 4094 var dlon = Proj4js.common.adjust_lon(lon - this.long0); 4095 var x,y; 4096 4097 if (Math.abs(lat) <= Proj4js.common.EPSLN) { 4098 x = this.x0 + this.R * dlon; 4099 y = this.y0; 4100 } 4101 var theta = Proj4js.common.asinz(2.0 * Math.abs(lat / Proj4js.common.PI)); 4102 if ((Math.abs(dlon) <= Proj4js.common.EPSLN) || (Math.abs(Math.abs(lat) - Proj4js.common.HALF_PI) <= Proj4js.common.EPSLN)) { 4103 x = this.x0; 4104 if (lat >= 0) { 4105 y = this.y0 + Proj4js.common.PI * this.R * Math.tan(.5 * theta); 4106 } else { 4107 y = this.y0 + Proj4js.common.PI * this.R * - Math.tan(.5 * theta); 4108 } 4109 // return(OK); 4110 } 4111 var al = .5 * Math.abs((Proj4js.common.PI / dlon) - (dlon / Proj4js.common.PI)); 4112 var asq = al * al; 4113 var sinth = Math.sin(theta); 4114 var costh = Math.cos(theta); 4115 4116 var g = costh / (sinth + costh - 1.0); 4117 var gsq = g * g; 4118 var m = g * (2.0 / sinth - 1.0); 4119 var msq = m * m; 4120 var con = Proj4js.common.PI * this.R * (al * (g - msq) + Math.sqrt(asq * (g - msq) * (g - msq) - (msq + asq) * (gsq - msq))) / (msq + asq); 4121 if (dlon < 0) { 4122 con = -con; 4123 } 4124 x = this.x0 + con; 4125 con = Math.abs(con / (Proj4js.common.PI * this.R)); 4126 if (lat >= 0) { 4127 y = this.y0 + Proj4js.common.PI * this.R * Math.sqrt(1.0 - con * con - 2.0 * al * con); 4128 } else { 4129 y = this.y0 - Proj4js.common.PI * this.R * Math.sqrt(1.0 - con * con - 2.0 * al * con); 4130 } 4131 p.x = x; 4132 p.y = y; 4133 return p; 4134 }, 4135 4136/* Van Der Grinten inverse equations--mapping x,y to lat/long 4137 ---------------------------------------------------------*/ 4138 inverse: function(p) { 4139 var lon, lat; 4140 var xx,yy,xys,c1,c2,c3; 4141 var al,asq; 4142 var a1; 4143 var m1; 4144 var con; 4145 var th1; 4146 var d; 4147 4148 /* inverse equations 4149 -----------------*/ 4150 p.x -= this.x0; 4151 p.y -= this.y0; 4152 con = Proj4js.common.PI * this.R; 4153 xx = p.x / con; 4154 yy =p.y / con; 4155 xys = xx * xx + yy * yy; 4156 c1 = -Math.abs(yy) * (1.0 + xys); 4157 c2 = c1 - 2.0 * yy * yy + xx * xx; 4158 c3 = -2.0 * c1 + 1.0 + 2.0 * yy * yy + xys * xys; 4159 d = yy * yy / c3 + (2.0 * c2 * c2 * c2 / c3 / c3 / c3 - 9.0 * c1 * c2 / c3 /c3) / 27.0; 4160 a1 = (c1 - c2 * c2 / 3.0 / c3) / c3; 4161 m1 = 2.0 * Math.sqrt( -a1 / 3.0); 4162 con = ((3.0 * d) / a1) / m1; 4163 if (Math.abs(con) > 1.0) { 4164 if (con >= 0.0) { 4165 con = 1.0; 4166 } else { 4167 con = -1.0; 4168 } 4169 } 4170 th1 = Math.acos(con) / 3.0; 4171 if (p.y >= 0) { 4172 lat = (-m1 *Math.cos(th1 + Proj4js.common.PI / 3.0) - c2 / 3.0 / c3) * Proj4js.common.PI; 4173 } else { 4174 lat = -(-m1 * Math.cos(th1 + Proj4js.common.PI / 3.0) - c2 / 3.0 / c3) * Proj4js.common.PI; 4175 } 4176 4177 if (Math.abs(xx) < Proj4js.common.EPSLN) { 4178 lon = this.long0; 4179 } 4180 lon = Proj4js.common.adjust_lon(this.long0 + Proj4js.common.PI * (xys - 1.0 + Math.sqrt(1.0 + 2.0 * (
4180xx * xx - yy * yy) + xys * xys)) / 2.0 / xx); 4181 4182 p.x=lon; 4183 p.y=lat; 4184 return p; 4185 } 4186}; 4187/* ====================================================================== 4188 projCode/cea.js 4189 ====================================================================== */ 4190 4191/******************************************************************************* 4192NAME LAMBERT CYLINDRICAL EQUAL AREA 4193 4194PURPOSE: Transforms input longitude and latitude to Easting and 4195 Northing for the Lambert Cylindrical Equal Area projection. 4196 This class of projection includes the Behrmann and 4197 Gall-Peters Projections. The 4198 longitude and latitude must be in radians. The Easting 4199 and Northing values will be returned in meters. 4200 4201PROGRAMMER DATE 4202---------- ---- 4203R. Marsden August 2009 4204Winwaed Software Tech LLC, http://www.winwaed.com 4205 4206This function was adapted from the Miller Cylindrical Projection in the Proj4JS 4207library. 4208 4209Note: This implementation assumes a Spherical Earth. The (commented) code 4210has been included for the ellipsoidal forward transform, but derivation of 4211the ellispoidal inverse transform is beyond me. Note that most of the 4212Proj4JS implementations do NOT currently support ellipsoidal figures. 4213Therefore this is not seen as a problem - especially this lack of support 4214is explicitly stated here. 4215 4216ALGORITHM REFERENCES 4217 42181. "Cartographic Projection Procedures for the UNIX Environment - 4219 A User's Manual" by Gerald I. Evenden, USGS Open File Report 90-284 4220 and Release 4 Interim Reports (2003) 4221 42222. Snyder, John P., "Flattening the Earth - Two Thousand Years of Map 4223 Projections", Univ. Chicago Press, 1993 4224*******************************************************************************/ 4225 4226Proj4js.Proj.cea = { 4227 4228/* Initialize the Cylindrical Equal Area projection 4229 -------------------------------------------*/ 4230 init: function() { 4231 //no-op 4232 }, 4233 4234 4235 /* Cylindrical Equal Area forward equations--mapping lat,long to x,y 4236 ------------------------------------------------------------*/ 4237 forward: function(p) { 4238 var lon=p.x; 4239 var lat=p.y; 4240 /* Forward equations 4241 -----------------*/ 4242 var dlon = Proj4js.common.adjust_lon(lon -this.long0); 4243 var x = this.x0 + this.a * dlon * Math.cos(this.lat_ts); 4244 var y = this.y0 + this.a * Math.sin(lat) / Math.cos(this.lat_ts); 4245 /* Elliptical Forward Transform 4246 Not implemented due to a lack of a matchign inverse function 4247 { 4248 var Sin_Lat = Math.sin(lat); 4249 var Rn = this.a * (Math.sqrt(1.0e0 - this.es * Sin_Lat * Sin_Lat )); 4250 x = this.x0 + this.a * dlon * Math.cos(this.lat_ts); 4251 y = this.y0 + Rn * Math.sin(lat) / Math.cos(this.lat_ts); 4252 } 4253 */ 4254 4255 4256 p.x=x; 4257 p.y=y; 4258 return p; 4259 },//ceaFwd() 4260 4261 /* Cylindrical Equal Area inverse equations--mapping x,y to lat/long 4262 ------------------------------------------------------------*/ 4263 inverse: function(p) { 4264 p.x -= this.x0; 4265 p.y -= this.y0; 4266 4267 var lon = Proj4js.common.adjust_lon( this.long0 + (p.x / this.a) / Math.cos(this.lat_ts) ); 4268 4269 var lat = Math.asin( (p.y/this.a) * Math.cos(this.lat_ts) ); 4270 4271 p.x=lon; 4272 p.y=lat; 4273 return p; 4274 }//ceaInv() 4275}; 4276/* ====================================================================== 4277 projCode/eqc.js 4278 ====================================================================== */ 4279 4280/* similar to equi.js FIXME proj4 uses eqc */ 4281Proj4js.Proj.eqc = { 4282 init : function() { 4283 4284 if(!this.x0) this.x0=0; 4285 if(!this.y0) this.y0=0; 4286 if(!this.lat0) this.lat0=0; 4287 if(!this.long0) this.long0=0; 4288 if(!this.lat_ts) this.lat_ts=0; 4289 if (!this.title) this.title = "Equidistant Cylindrical (Plate Carre)"; 4290 4291 this.rc= Math.cos(this.lat_ts); 4292 }, 4293 4294 4295 // forward equations--mapping lat,long to x,y 4296 // ----------------------------------------------------------------- 4297 forward : function(p) { 4298 4299 var lon= p.x; 4300 var lat= p.y; 4301 4302 var dlon = Proj4js.common.adjust_lon(lon - this.long0); 4303 var dlat = Proj4js.common.adjust_lat(lat - this.lat0 ); 4304 p.x= this.x0 + (this.a*dlon*this.rc); 4305 p.y= this.y0 + (this.a*dlat ); 4306 return p; 4307 }, 4308 4309 // inverse equations--mapping x,y to lat/long 4310 // ----------------------------------------------------------------- 4311 inverse : function(p) { 4312 4313 var x= p.x; 4314 var y= p.y; 4315 4316 p.x= Proj4js.common.adjust_lon(this.long0 + ((x - this.x0)/(this.a*this.rc))); 4317 p.y= Proj4js.common.adjust_lat(this.lat0 + ((y - this.y0)/(this.a ))); 4318 return p; 4319 } 4320 4321}; 4322/* ====================================================================== 4323 projCode/cass.js 4324 ====================================================================== */ 4325
4326/******************************************************************************* 4327NAME CASSINI 4328 4329PURPOSE: Transforms input longitude and latitude to Easting and 4330 Northing for the Cassini projection. The 4331 longitude and latitude must be in radians. The Easting 4332 and Northing values will be returned in meters. 4333 Ported from PROJ.4. 4334 4335 4336ALGORITHM REFERENCES 4337 43381. Snyder, John P., "Map Projections--A Working Manual", U.S. Geological 4339 Survey Professional Paper 1395 (Supersedes USGS Bulletin 1532), United 4340 State Government Printing Office, Washington D.C., 1987. 4341 43422. Snyder, John P. and Voxland, Philip M., "An Album of Map Projections", 4343 U.S. Geological Survey Professional Paper 1453 , United State Government 4344*******************************************************************************/ 4345 4346 4347//Proj4js.defs["EPSG:28191"] = "+proj=cass +lat_0=31.73409694444445 +lon_0=35.21208055555556 +x_0=170251.555 +y_0=126867.909 +a=6378300.789 +b=6356566.435 +towgs84=-275.722,94.7824,340.894,-8.001,-4.42,-11.821,1 +units=m +no_defs"; 4348 4349// Initialize the Cassini projection 4350// ----------------------------------------------------------------- 4351 4352Proj4js.Proj.cass = { 4353 init : function() { 4354 if (!this.sphere) { 4355 this.en = Proj4js.common.pj_enfn(this.es) 4356 this.m0 = Proj4js.common.pj_mlfn(this.lat0, Math.sin(this.lat0), Math.cos(this.lat0), this.en); 4357 } 4358 }, 4359 4360 C1: .16666666666666666666, 4361 C2: .00833333333333333333, 4362 C3: .04166666666666666666, 4363 C4: .33333333333333333333, 4364 C5: .06666666666666666666, 4365 4366 4367/* Cassini forward equations--mapping lat,long to x,y 4368 -----------------------------------------------------------------------*/ 4369 forward: function(p) { 4370 4371 /* Forward equations 4372 -----------------*/ 4373 var x,y; 4374 var lam=p.x; 4375 var phi=p.y; 4376 lam = Proj4js.common.adjust_lon(lam - this.long0); 4377 4378 if (this.sphere) { 4379 x = Math.asin(Math.cos(phi) * Math.sin(lam)); 4380 y = Math.atan2(Math.tan(phi) , Math.cos(lam)) - this.phi0; 4381 } else { 4382 //ellipsoid 4383 this.n = Math.sin(phi); 4384 this.c = Math.cos(phi); 4385 y = Proj4js.common.pj_mlfn(phi, this.n, this.c, this.en); 4386 this.n = 1./Math.sqrt(1. - this.es * this.n * this.n); 4387 this.tn = Math.tan(phi); 4388 this.t = this.tn * this.tn; 4389 this.a1 = lam * this.c; 4390 this.c *= this.es * this.c / (1 - this.es); 4391 this.a2 = this.a1 * this.a1; 4392 x = this.n * this.a1 * (1. - this.a2 * this.t * (this.C1 - (8. - this.t + 8. * this.c) * this.a2 * this.C2)); 4393 y -= this.m0 - this.n * this.tn * this.a2 * (.5 + (5. - this.t + 6. * this.c) * this.a2 * this.C3); 4394 } 4395 4396 p.x = this.a*x + this.x0; 4397 p.y = this.a*y + this.y0; 4398 return p; 4399 },//cassFwd() 4400 4401/* Inverse equations 4402 -----------------*/ 4403 inverse: function(p) { 4404 p.x -= this.x0; 4405 p.y -= this.y0; 4406 var x = p.x/this.a; 4407 var y = p.y/this.a; 4408 var phi, lam; 4409 4410 if (this.sphere) { 4411 this.dd = y + this.lat0; 4412 phi = Math.asin(Math.sin(this.dd) * Math.cos(x)); 4413 lam = Math.atan2(Math.tan(x), Math.cos(this.dd)); 4414 } else { 4415 /* ellipsoid */ 4416 var ph1 = Proj4js.common.pj_inv_mlfn(this.m0 + y, this.es, this.en); 4417 this.tn = Math.tan(ph1); 4418 this.t = this.tn * this.tn; 4419 this.n = Math.sin(ph1); 4420 this.r = 1. / (1. - this.es * this.n * this.n); 4421 this.n = Math.sqrt(this.r); 4422 this.r *= (1. - this.es) * this.n; 4423 this.dd = x / this.n; 4424 this.d2 = this.dd * this.dd; 4425 phi = ph1 - (this.n * this.tn / this.r) * this.d2 * (.5 - (1. + 3. * this.t) * this.d2 * this.C3); 4426 lam = this.dd * (1. + this.t * this.d2 * (-this.C4 + (1. + 3. * this.t) * this.d2 * this.C5)) / Math.cos(ph1); 4427 } 4428 p.x = Proj4js.common.adjust_lon(this.long0+lam); 4429 p.y = phi; 4430 return p; 4431 }//cassInv() 4432 4433} 4434/* ====================================================================== 4435 projCode/gauss.js 4436 ====================================================================== */ 4437 4438 4439Proj4js.Proj.gauss = { 4440 4441 init : function() { 4442 var sphi = Math.sin(this.lat0); 4443 var cphi = Math.cos(this.lat0); 4444 cphi *= cphi; 4445 this.rc = Math.sqrt(1.0 - this.es) / (1.0 - this.es * sphi * sphi); 4446 this.C = Math.sqrt(1.0 + this.es * cphi * cphi / (1.0 - this.es)); 4447 this.phic0 = Math.asin(sphi / this.C); 4448 this.ratexp = 0.5 * this.C * this.e; 4449 this.K = Math.tan(0.5 * this.phic0 + Proj4js.common.FORTPI) / (Math.pow(Math.tan(0.5*this.lat0 + Proj4js.common.FORTPI), this.C) * Proj4js.common.srat(this.e*sphi, this.ratexp)); 4450 }, 4451 4452 forward : function(p) { 4453 var lon = p.x; 4454 var lat = p.y; 4455 4456 p.y = 2.0 * Math.atan( this.K * Math.pow(Math.tan(0.5 * lat + Proj4js.common.FORTPI), this.C) * Proj4js.common.srat(this.e * Math.sin(lat), this.ratexp) ) - Proj4js.common.HALF_PI; 4457 p.x = this.C * lon; 4458 return p; 4459 }, 4460 4461 inverse : function(p) { 4462 var DEL_TOL = 1e-14; 4463 var lon = p.x / this.C; 4464 var lat = p.y; 4465 var num = Math.pow(Math.tan(0.5 * lat + Proj4js.common.FORTPI)/this.K, 1./this.C); 4466 for (var i = Proj4js.common.MAX_ITER; i>0; --i) { 4467 lat = 2.0 * Math.atan(num * Proj4js.common.srat(this.e * Math.sin(p.y), -0.5 * this.e)) - Proj4js.common.HALF_PI; 4468 if (Math.abs(lat - p.y) < DEL_TOL) break; 4469 p.y = lat; 4470 } 4471 /* convergence failed */ 4472 if (!i) { 4473 Proj4js.reportError("gauss:inverse:convergence failed"); 4474 return null; 4475 } 4476 p.x = lon; 4477 p.y = lat; 4478 return p; 4479 } 4480}; 4481 4482/* ====================================================================== 4483 projCode/omerc.js 4484 ====================================================================== */ 4485
4486/******************************************************************************* 4487NAME OBLIQUE MERCATOR (HOTINE) 4488 4489PURPOSE: Transforms input longitude and latitude to Easting and 4490 Northing for the Oblique Mercator projection. The 4491 longitude and latitude must be in radians. The Easting 4492 and Northing values will be returned in meters. 4493 4494PROGRAMMER DATE 4495---------- ---- 4496T. Mittan Mar, 1993 4497 4498ALGORITHM REFERENCES 4499 45001. Snyder, John P., "Map Projections--A Working Manual", U.S. Geological 4501 Survey Professional Paper 1395 (Supersedes USGS Bulletin 1532), United 4502 State Government Printing Office, Washington D.C., 1987. 4503 45042. Snyder, John P. and Voxland, Philip M., "An Album of Map Projections", 4505 U.S. Geological Survey Professional Paper 1453 , United State Government 4506 Printing Office, Washington D.C., 1989. 4507*******************************************************************************/ 4508 4509Proj4js.Proj.omerc = { 4510 4511 /* Initialize the Oblique Mercator projection 4512 ------------------------------------------*/ 4513 init: function() { 4514 if (!this.mode) this.mode=0; 4515 if (!this.lon1) {this.lon1=0;this.mode=1;} 4516 if (!this.lon2) this.lon2=0; 4517 if (!this.lat2) this.lat2=0; 4518 4519 /* Place parameters in static storage for common use 4520 -------------------------------------------------*/ 4521 var temp = this.b/ this.a; 4522 var es = 1.0 - Math.pow(temp,2); 4523 var e = Math.sqrt(es); 4524 4525 this.sin_p20=Math.sin(this.lat0); 4526 this.cos_p20=Math.cos(this.lat0); 4527 4528 this.con = 1.0 - this.es * this.sin_p20 * this.sin_p20; 4529 this.com = Math.sqrt(1.0 - es); 4530 this.bl = Math.sqrt(1.0 + this.es * Math.pow(this.cos_p20,4.0)/(1.0 - es)); 4531 this.al = this.a * this.bl * this.k0 * this.com / this.con; 4532 if (Math.abs(this.lat0) < Proj4js.common.EPSLN) { 4533 this.ts = 1.0; 4534 this.d = 1.0; 4535 this.el = 1.0; 4536 } else { 4537 this.ts = Proj4js.common.tsfnz(this.e,this.lat0,this.sin_p20); 4538 this.con = Math.sqrt(this.con); 4539 this.d = this.bl * this.com / (this.cos_p20 * this.con); 4540 if ((this.d * this.d - 1.0) > 0.0) { 4541 if (this.lat0 >= 0.0) { 4542 this.f = this.d + Math.sqrt(this.d * this.d - 1.0); 4543 } else { 4544 this.f = this.d - Math.sqrt(this.d * this.d - 1.0); 4545 } 4546 } else { 4547 this.f = this.d; 4548 } 4549 this.el = this.f * Math.pow(this.ts,this.bl); 4550 } 4551 4552 //this.longc=52.60353916666667; 4553 4554 if (this.mode != 0) { 4555 this.g = .5 * (this.f - 1.0/this.f); 4556 this.gama = Proj4js.common.asinz(Math.sin(this.alpha) / this.d); 4557 this.longc= this.longc - Proj4js.common.asinz(this.g * Math.tan(this.gama))/this.bl; 4558 4559 /* Report parameters common to format B 4560 -------------------------------------*/ 4561 //genrpt(azimuth * R2D,"Azimuth of Central Line: "); 4562 //cenlon(lon_origin); 4563 // cenlat(lat_origin); 4564 4565 this.con = Math.abs(this.lat0); 4566 if ((this.con > Proj4js.common.EPSLN) && (Math.abs(this.con - Proj4js.common.HALF_PI) > Proj4js.common.EPSLN)) { 4567 this.singam=Math.sin(this.gama); 4568 this.cosgam=Math.cos(this.gama); 4569 4570 this.sinaz=Math.sin(this.alpha); 4571 this.cosaz=Math.cos(this.alpha); 4572 4573 if (this.lat0>= 0) { 4574 this.u = (this.al / this.bl) * Math.atan(Math.sqrt(this.d*this.d - 1.0)/this.cosaz); 4575 } else { 4576 this.u = -(this.al / this.bl) *Math.atan(Math.sqrt(this.d*this.d - 1.0)/this.cosaz); 4577 } 4578 } else { 4579 Proj4js.reportError("omerc:Init:DataError"); 4580 } 4581 } else { 4582 this.sinphi =Math. sin(this.at1); 4583 this.ts1 = Proj4js.common.tsfnz(this.e,this.lat1,this.sinphi); 4584 this.sinphi = Math.sin(this.lat2); 4585 this.ts2 = Proj4js.common.tsfnz(this.e,this.lat2,this.sinphi); 4586 this.h = Math.pow(this.ts1,this.bl); 4587 this.l = Math.pow(this.ts2,this.bl); 4588 this.f = this.el/this.h; 4589 this.g = .5 * (this.f - 1.0/this.f); 4590 this.j = (this.el * this.el - this.l * this.h)/(this.el * this.el + this.l * this.h); 4591 this.p = (this.l - this.h) / (this.l + this.h); 4592 this.dlon = this.lon1 - this.lon2; 4593 if (this.dlon < -Proj4js.common.PI) this.lon2 = this.lon2 - 2.0 * Proj4js.common.PI; 4594 if (this.dlon > Proj4js.common.PI) this.lon2 = this.lon2 + 2.0 * Proj4js.common.PI; 4595 this.dlon = this.lon1 - this.lon2; 4596 this.longc = .5 * (this.lon1 + this.lon2) -Math.atan(this.j * Math.tan(.5 * this.bl * this.dlon)/this.p)/this.bl; 4597 this.dlon = Proj4js.common.adjust_lon(this.lon1 - this.longc); 4598 this.gama = Math.atan(Math.sin(this.bl * this.dlon)/this.g); 4599 this.alpha = Proj4js.common.asinz(this.d * Math.sin(this.gama)); 4600 4601 /* Report parameters common to format A 4602 -------------------------------------*/ 4603 4604 if (Math.abs(this.lat1 - this.lat2) <= Proj4js.common.EPSLN) { 4605 Proj4js.reportError("omercInitDataError"); 4606 //return(202); 4607 } else { 4608 this.con = Math.abs(this.lat1); 4609 } 4610 if ((this.con <= Proj4js.common.EPSLN) || (Math.abs(this.con - Proj4js.common.HALF_PI) <= Proj4js.common.EPSLN)) { 4611 Proj4js.reportError("omercInitDataError"); 4612 //return(202); 4613 } else { 4614 if (Math.abs(Math.abs(this.lat0) - Proj4js.common.HALF_PI) <= Proj4js.common.EPSLN) { 4615 Proj4js.reportError("omercInitDataError"); 4616 //return(202); 4617 } 4618 } 4619
4620 this.singam=Math.sin(this.gam); 4621 this.cosgam=Math.cos(this.gam); 4622 4623 this.sinaz=Math.sin(this.alpha); 4624 this.cosaz=Math.cos(this.alpha); 4625 4626 4627 if (this.lat0 >= 0) { 4628 this.u = (this.al/this.bl) * Math.atan(Math.sqrt(this.d * this.d - 1.0)/this.cosaz); 4629 } else { 4630 this.u = -(this.al/this.bl) * Math.atan(Math.sqrt(this.d * this.d - 1.0)/this.cosaz); 4631 } 4632 } 4633 }, 4634 4635 4636 /* Oblique Mercator forward equations--mapping lat,long to x,y 4637 ----------------------------------------------------------*/ 4638 forward: function(p) { 4639 var theta; /* angle */ 4640 var sin_phi, cos_phi;/* sin and cos value */ 4641 var b; /* temporary values */ 4642 var c, t, tq; /* temporary values */ 4643 var con, n, ml; /* cone constant, small m */ 4644 var q,us,vl; 4645 var ul,vs; 4646 var s; 4647 var dlon; 4648 var ts1; 4649 4650 var lon=p.x; 4651 var lat=p.y; 4652 /* Forward equations 4653 -----------------*/ 4654 sin_phi = Math.sin(lat); 4655 dlon = Proj4js.common.adjust_lon(lon - this.longc); 4656 vl = Math.sin(this.bl * dlon); 4657 if (Math.abs(Math.abs(lat) - Proj4js.common.HALF_PI) > Proj4js.common.EPSLN) { 4658 ts1 = Proj4js.common.tsfnz(this.e,lat,sin_phi); 4659 q = this.el / (Math.pow(ts1,this.bl)); 4660 s = .5 * (q - 1.0 / q); 4661 t = .5 * (q + 1.0/ q); 4662 ul = (s * this.singam - vl * this.cosgam) / t; 4663 con = Math.cos(this.bl * dlon); 4664 if (Math.abs(con) < .0000001) { 4665 us = this.al * this.bl * dlon; 4666 } else { 4667 us = this.al * Math.atan((s * this.cosgam + vl * this.singam) / con)/this.bl; 4668 if (con < 0) us = us + Proj4js.common.PI * this.al / this.bl; 4669 } 4670 } else { 4671 if (lat >= 0) { 4672 ul = this.singam; 4673 } else { 4674 ul = -this.singam; 4675 } 4676 us = this.al * lat / this.bl; 4677 } 4678 if (Math.abs(Math.abs(ul) - 1.0) <= Proj4js.common.EPSLN) { 4679 //alert("Point projects into infinity","omer-for"); 4680 Proj4js.reportError("omercFwdInfinity"); 4681 //return(205); 4682 } 4683 vs = .5 * this.al * Math.log((1.0 - ul)/(1.0 + ul)) / this.bl; 4684 us = us - this.u; 4685 var x = this.x0 + vs * this.cosaz + us * this.sinaz; 4686 var y = this.y0 + us * this.cosaz - vs * this.sinaz; 4687 4688 p.x=x; 4689 p.y=y; 4690 return p; 4691 }, 4692 4693 inverse: function(p) { 4694 var delta_lon; /* Delta longitude (Given longitude - center */ 4695 var theta; /* angle */ 4696 var delta_theta; /* adjusted longitude */ 4697 var sin_phi, cos_phi;/* sin and cos value */ 4698 var b; /* temporary values */ 4699 var c, t, tq; /* temporary values */ 4700 var con, n, ml; /* cone constant, small m */ 4701 var vs,us,q,s,ts1; 4702 var vl,ul,bs; 4703 var lon, lat; 4704 var flag; 4705 4706 /* Inverse equations 4707 -----------------*/ 4708 p.x -= this.x0; 4709 p.y -= this.y0; 4710 flag = 0; 4711 vs = p.x * this.cosaz - p.y * this.sinaz; 4712 us = p.y * this.cosaz + p.x * this.sinaz; 4713 us = us + this.u; 4714 q = Math.exp(-this.bl * vs / this.al); 4715 s = .5 * (q - 1.0/q); 4716 t = .5 * (q + 1.0/q); 4717 vl = Math.sin(this.bl * us / this.al); 4718 ul = (vl * this.cosgam + s * this.singam)/t; 4719 if (Math.abs(Math.abs(ul) - 1.0) <= Proj4js.common.EPSLN) 4720 { 4721 lon = this.longc; 4722 if (ul >= 0.0) { 4723 lat = Proj4js.common.HALF_PI; 4724 } else { 4725 lat = -Proj4js.common.HALF_PI; 4726 } 4727 } else { 4728 con = 1.0 / this.bl; 4729 ts1 =Math.pow((this.el / Math.sqrt((1.0 + ul) / (1.0 - ul))),con); 4730 lat = Proj4js.common.phi2z(this.e,ts1); 4731 //if (flag != 0) 4732 //return(flag); 4733 //~ con = Math.cos(this.bl * us /al); 4734 theta = this.longc - Math.atan2((s * this.cosgam - vl * this.singam) , con)/this.bl; 4735 lon = Proj4js.common.adjust_lon(theta); 4736 } 4737 p.x=lon; 4738 p.y=lat; 4739 return p; 4740 } 4741}; 4742/* ====================================================================== 4743 projCode/lcc.js 4744 ====================================================================== */ 4745 4746/******************************************************************************* 4747NAME LAMBERT CONFORMAL CONIC 4748 4749PURPOSE: Transforms input longitude and latitude to Easting and 4750 Northing for the Lambert Conformal Conic projection. The 4751 longitude and latitude must be in radians. The Easting 4752 and Northing values will be returned in meters. 4753 4754 4755ALGORITHM REFERENCES 4756 47571. Snyder, John P., "Map Projections--A Working Manual", U.S. Geological 4758 Survey Professional Paper 1395 (Supersedes USGS Bulletin 1532), United 4759 State Government Printing Office, Washington D.C., 1987. 4760 47612. Snyder, John P. and Voxland, Philip M., "An Album of Map Projections", 4762 U.S. Geological Survey Professional Paper 1453 , United State Government 4763*******************************************************************************/ 4764 4765 4766//<2104> +proj=lcc +lat_1=10.16666666666667 +lat_0=10.16666666666667 +lon_0=-71.60561777777777 +k_0=1 +x0=-17044 +x0=-23139.97 +ellps=intl +units=m +no_defs no_defs 4767 4768// Initialize the Lambert Conformal conic projection 4769// ----------------------------------------------------------------- 4770 4771//Proj4js.Proj.lcc = Class.create(); 4772Proj4js.Proj.lcc = { 4773 init : function() { 4774 4775 // array of: r_maj,r_min,lat1,lat2,c_lon,c_lat,false_east,fal
4775se_north 4776 //double c_lat; /* center latitude */ 4777 //double c_lon; /* center longitude */ 4778 //double lat1; /* first standard parallel */ 4779 //double lat2; /* second standard parallel */ 4780 //double r_maj; /* major axis */ 4781 //double r_min; /* minor axis */ 4782 //double false_east; /* x offset in meters */ 4783 //double false_north; /* y offset in meters */ 4784 4785 if (!this.lat2){this.lat2=this.lat0;}//if lat2 is not defined 4786 if (!this.k0) this.k0 = 1.0; 4787 4788 // Standard Parallels cannot be equal and on opposite sides of the equator 4789 if (Math.abs(this.lat1+this.lat2) < Proj4js.common.EPSLN) { 4790 Proj4js.reportError("lcc:init: Equal Latitudes"); 4791 return; 4792 } 4793 4794 var temp = this.b / this.a; 4795 this.e = Math.sqrt(1.0 - temp*temp); 4796 4797 var sin1 = Math.sin(this.lat1); 4798 var cos1 = Math.cos(this.lat1); 4799 var ms1 = Proj4js.common.msfnz(this.e, sin1, cos1); 4800 var ts1 = Proj4js.common.tsfnz(this.e, this.lat1, sin1); 4801 4802 var sin2 = Math.sin(this.lat2); 4803 var cos2 = Math.cos(this.lat2); 4804 var ms2 = Proj4js.common.msfnz(this.e, sin2, cos2); 4805 var ts2 = Proj4js.common.tsfnz(this.e, this.lat2, sin2); 4806 4807 var ts0 = Proj4js.common.tsfnz(this.e, this.lat0, Math.sin(this.lat0)); 4808 4809 if (Math.abs(this.lat1 - this.lat2) > Proj4js.common.EPSLN) { 4810 this.ns = Math.log(ms1/ms2)/Math.log(ts1/ts2); 4811 } else { 4812 this.ns = sin1; 4813 } 4814 this.f0 = ms1 / (this.ns * Math.pow(ts1, this.ns)); 4815 this.rh = this.a * this.f0 * Math.pow(ts0, this.ns); 4816 if (!this.title) this.title = "Lambert Conformal Conic"; 4817 }, 4818 4819 4820 // Lambert Conformal conic forward equations--mapping lat,long to x,y 4821 // ----------------------------------------------------------------- 4822 forward : function(p) { 4823 4824 var lon = p.x; 4825 var lat = p.y; 4826 4827 // convert to radians 4828 if ( lat <= 90.0 && lat >= -90.0 && lon <= 180.0 && lon >= -180.0) { 4829 //lon = lon * Proj4js.common.D2R; 4830 //lat = lat * Proj4js.common.D2R; 4831 } else { 4832 Proj4js.reportError("lcc:forward: llInputOutOfRange: "+ lon +" : " + lat); 4833 return null; 4834 } 4835 4836 var con = Math.abs( Math.abs(lat) - Proj4js.common.HALF_PI); 4837 var ts, rh1; 4838 if (con > Proj4js.common.EPSLN) { 4839 ts = Proj4js.common.tsfnz(this.e, lat, Math.sin(lat) ); 4840 rh1 = this.a * this.f0 * Math.pow(ts, this.ns); 4841 } else { 4842 con = lat * this.ns; 4843 if (con <= 0) { 4844 Proj4js.reportError("lcc:forward: No Projection"); 4845 return null; 4846 } 4847 rh1 = 0; 4848 } 4849 var theta = this.ns * Proj4js.common.adjust_lon(lon - this.long0); 4850 p.x = this.k0 * (rh1 * Math.sin(theta)) + this.x0; 4851 p.y = this.k0 * (this.rh - rh1 * Math.cos(theta)) + this.y0; 4852 4853 return p; 4854 }, 4855 4856 // Lambert Conformal Conic inverse equations--mapping x,y to lat/long 4857 // ----------------------------------------------------------------- 4858 inverse : function(p) { 4859 4860 var rh1, con, ts; 4861 var lat, lon; 4862 var x = (p.x - this.x0)/this.k0; 4863 var y = (this.rh - (p.y - this.y0)/this.k0); 4864 if (this.ns > 0) { 4865 rh1 = Math.sqrt (x * x + y * y); 4866 con = 1.0; 4867 } else { 4868 rh1 = -Math.sqrt (x * x + y * y); 4869 con = -1.0; 4870 } 4871 var theta = 0.0; 4872 if (rh1 != 0) { 4873 theta = Math.atan2((con * x),(con * y)); 4874 } 4875 if ((rh1 != 0) || (this.ns > 0.0)) { 4876 con = 1.0/this.ns; 4877 ts = Math.pow((rh1/(this.a * this.f0)), con); 4878 lat = Proj4js.common.phi2z(this.e, ts); 4879 if (lat == -9999) return null; 4880 } else { 4881 lat = -Proj4js.common.HALF_PI; 4882 } 4883 lon = Proj4js.common.adjust_lon(theta/this.ns + this.long0); 4884 4885 p.x = lon; 4886 p.y = lat; 4887 return p; 4888 } 4889}; 4890 4891 4892 4893 4894/* ====================================================================== 4895 projCode/laea.js 4896 ====================================================================== */ 4897
4898/******************************************************************************* 4899NAME LAMBERT AZIMUTHAL EQUAL-AREA 4900 4901PURPOSE: Transforms input longitude and latitude to Easting and 4902 Northing for the Lambert Azimuthal Equal-Area projection. The 4903 longitude and latitude must be in radians. The Easting 4904 and Northing values will be returned in meters. 4905 4906PROGRAMMER DATE 4907---------- ---- 4908D. Steinwand, EROS March, 1991 4909 4910This function was adapted from the Lambert Azimuthal Equal Area projection 4911code (FORTRAN) in the General Cartographic Transformation Package software 4912which is available from the U.S. Geological Survey National Mapping Division. 4913 4914ALGORITHM REFERENCES 4915 49161. "New Equal-Area Map Projections for Noncircular Regions", John P. Snyder, 4917 The American Cartographer, Vol 15, No. 4, October 1988, pp. 341-355. 4918 49192. Snyder, John P., "Map Projections--A Working Manual", U.S. Geological 4920 Survey Professional Paper 1395 (Supersedes USGS Bulletin 1532), United 4921 State Government Printing Office, Washington D.C., 1987. 4922 49233. "Software Documentation for GCTP General Cartographic Transformation 4924 Package", U.S. Geological Survey National Mapping Division, May 1982. 4925*******************************************************************************/ 4926 4927Proj4js.Proj.laea = { 4928 S_POLE: 1, 4929 N_POLE: 2, 4930 EQUIT: 3, 4931 OBLIQ: 4, 4932 4933 4934/* Initialize the Lambert Azimuthal Equal Area projection 4935 ------------------------------------------------------*/ 4936 init: function() { 4937 var t = Math.abs(this.lat0); 4938 if (Math.abs(t - Proj4js.common.HALF_PI) < Proj4js.common.EPSLN) { 4939 this.mode = this.lat0 < 0. ? this.S_POLE : this.N_POLE; 4940 } else if (Math.abs(t) < Proj4js.common.EPSLN) { 4941 this.mode = this.EQUIT; 4942 } else { 4943 this.mode = this.OBLIQ; 4944 } 4945 if (this.es > 0) { 4946 var sinphi; 4947 4948 this.qp = Proj4js.common.qsfnz(this.e, 1.0); 4949 this.mmf = .5 / (1. - this.es); 4950 this.apa = this.authset(this.es); 4951 switch (this.mode) { 4952 case this.N_POLE: 4953 case this.S_POLE: 4954 this.dd = 1.; 4955 break; 4956 case this.EQUIT: 4957 this.rq = Math.sqrt(.5 * this.qp); 4958 this.dd = 1. / this.rq; 4959 this.xmf = 1.; 4960 this.ymf = .5 * this.qp; 4961 break; 4962 case this.OBLIQ: 4963 this.rq = Math.sqrt(.5 * this.qp); 4964 sinphi = Math.sin(this.lat0); 4965 this.sinb1 = Proj4js.common.qsfnz(this.e, sinphi) / this.qp; 4966 this.cosb1 = Math.sqrt(1. - this.sinb1 * this.sinb1); 4967 this.dd = Math.cos(this.lat0) / (Math.sqrt(1. - this.es * sinphi * sinphi) * this.rq * this.cosb1); 4968 this.ymf = (this.xmf = this.rq) / this.dd; 4969 this.xmf *= this.dd; 4970 break; 4971 } 4972 } else { 4973 if (this.mode == this.OBLIQ) { 4974 this.sinph0 = Math.sin(this.lat0); 4975 this.cosph0 = Math.cos(this.lat0); 4976 } 4977 } 4978 }, 4979 4980/* Lambert Azimuthal Equal Area forward equations--mapping lat,long to x,y 4981 -----------------------------------------------------------------------*/ 4982 forward: function(p) { 4983 4984 /* Forward equations 4985 -----------------*/ 4986 var x,y; 4987 var lam=p.x; 4988 var phi=p.y; 4989 lam = Proj4js.common.adjust_lon(lam - this.long0); 4990 4991 if (this.sphere) { 4992 var coslam, cosphi, sinphi; 4993 4994 sinphi = Math.sin(phi); 4995 cosphi = Math.cos(phi); 4996 coslam = Math.cos(lam); 4997 switch (this.mode) { 4998 case this.OBLIQ: 4999 case this.EQUIT: 5000 y = (this.mode == this.EQUIT) ? 1. + cosphi * coslam : 1. + this.sinph0 * sinphi + this.cosph0 * cosphi * coslam; 5001 if (y <= Proj4js.common.EPSLN) { 5002 Proj4js.reportError("laea:fwd:y less than eps"); 5003 return null; 5004 } 5005 y = Math.sqrt(2. / y); 5006 x = y * cosphi * Math.sin(lam); 5007 y *= (this.mode == this.EQUIT) ? sinphi : this.cosph0 * sinphi - this.sinph0 * cosphi * coslam; 5008 break; 5009 case this.N_POLE: 5010 coslam = -coslam; 5011 case this.S_POLE: 5012 if (Math.abs(phi + this.phi0) < Proj4js.common.EPSLN) { 5013 Proj4js.reportError("laea:fwd:phi < eps"); 5014 return null; 5015 } 5016 y = Proj4js.common.FORTPI - phi * .5; 5017 y = 2. * ((this.mode == this.S_POLE) ? Math.cos(y) : Math.sin(y)); 5018 x = y * Math.sin(lam); 5019 y *= coslam; 5020 break; 5021 } 5022 } else { 5023 var coslam, sinlam, sinphi, q, sinb=0.0, cosb=0.0, b=0.0; 5024
5025 coslam = Math.cos(lam); 5026 sinlam = Math.sin(lam); 5027 sinphi = Math.sin(phi); 5028 q = Proj4js.common.qsfnz(this.e, sinphi); 5029 if (this.mode == this.OBLIQ || this.mode == this.EQUIT) { 5030 sinb = q / this.qp; 5031 cosb = Math.sqrt(1. - sinb * sinb); 5032 } 5033 switch (this.mode) { 5034 case this.OBLIQ: 5035 b = 1. + this.sinb1 * sinb + this.cosb1 * cosb * coslam; 5036 break; 5037 case this.EQUIT: 5038 b = 1. + cosb * coslam; 5039 break; 5040 case this.N_POLE: 5041 b = Proj4js.common.HALF_PI + phi; 5042 q = this.qp - q; 5043 break; 5044 case this.S_POLE: 5045 b = phi - Proj4js.common.HALF_PI; 5046 q = this.qp + q; 5047 break; 5048 } 5049 if (Math.abs(b) < Proj4js.common.EPSLN) { 5050 Proj4js.reportError("laea:fwd:b < eps"); 5051 return null; 5052 } 5053 switch (this.mode) { 5054 case this.OBLIQ: 5055 case this.EQUIT: 5056 b = Math.sqrt(2. / b); 5057 if (this.mode == this.OBLIQ) { 5058 y = this.ymf * b * (this.cosb1 * sinb - this.sinb1 * cosb * coslam); 5059 } else { 5060 y = (b = Math.sqrt(2. / (1. + cosb * coslam))) * sinb * this.ymf; 5061 } 5062 x = this.xmf * b * cosb * sinlam; 5063 break; 5064 case this.N_POLE: 5065 case this.S_POLE: 5066 if (q >= 0.) { 5067 x = (b = Math.sqrt(q)) * sinlam; 5068 y = coslam * ((this.mode == this.S_POLE) ? b : -b); 5069 } else { 5070 x = y = 0.; 5071 } 5072 break; 5073 } 5074 } 5075 5076 //v 1.0 5077 /* 5078 var sin_lat=Math.sin(lat); 5079 var cos_lat=Math.cos(lat); 5080 5081 var sin_delta_lon=Math.sin(delta_lon); 5082 var cos_delta_lon=Math.cos(delta_lon); 5083 5084 var g =this.sin_lat_o * sin_lat +this.cos_lat_o * cos_lat * cos_delta_lon; 5085 if (g == -1.0) { 5086 Proj4js.reportError("laea:fwd:Point projects to a circle of radius "+ 2.0 * R); 5087 return null; 5088 } 5089 var ksp = this.a * Math.sqrt(2.0 / (1.0 + g)); 5090 var x = ksp * cos_lat * sin_delta_lon + this.x0; 5091 var y = ksp * (this.cos_lat_o * sin_lat - this.sin_lat_o * cos_lat * cos_delta_lon) + this.y0; 5092 */ 5093 p.x = this.a*x + this.x0; 5094 p.y = this.a*y + this.y0; 5095 return p; 5096 },//lamazFwd() 5097 5098/* Inverse equations 5099 -----------------*/ 5100 inverse: function(p) { 5101 p.x -= this.x0; 5102 p.y -= this.y0; 5103 var x = p.x/this.a; 5104 var y = p.y/this.a; 5105 var lam, phi; 5106 5107 if (this.sphere) { 5108 var cosz=0.0, rh, sinz=0.0; 5109 5110 rh = Math.sqrt(x*x + y*y); 5111 phi = rh * .5; 5112 if (phi > 1.) { 5113 Proj4js.reportError("laea:Inv:DataError"); 5114 return null; 5115 } 5116 phi = 2. * Math.asin(phi); 5117 if (this.mode == this.OBLIQ || this.mode == this.EQUIT) { 5118 sinz = Math.sin(phi); 5119 cosz = Math.cos(phi); 5120 } 5121 switch (this.mode) { 5122 case this.EQUIT: 5123 phi = (Math.abs(rh) <= Proj4js.common.EPSLN) ? 0. : Math.asin(y * sinz / rh); 5124 x *= sinz; 5125 y = cosz * rh; 5126 break; 5127 case this.OBLIQ: 5128 phi = (Math.abs(rh) <= Proj4js.common.EPSLN) ? this.phi0 : Math.asin(cosz * this.sinph0 + y * sinz * this.cosph0 / rh); 5129 x *= sinz * this.cosph0; 5130 y = (cosz - Math.sin(phi) * this.sinph0) * rh; 5131 break; 5132 case this.N_POLE: 5133 y = -y; 5134 phi = Proj4js.common.HALF_PI - phi; 5135 break; 5136 case this.S_POLE: 5137 phi -= Proj4js.common.HALF_PI; 5138 break; 5139 } 5140 lam = (y == 0. && (this.mode == this.EQUIT || this.mode == this.OBLIQ)) ? 0. : Math.atan2(x, y); 5141 } else { 5142 var cCe, sCe, q, rho, ab=0.0; 5143 5144 switch (this.mode) { 5145 case this.EQUIT: 5146 case this.OBLIQ: 5147 x /= this.dd; 5148 y *= this.dd; 5149 rho = Math.sqrt(x*x + y*y); 5150 if (rho < Proj4js.common.EPSLN) { 5151 p.x = 0.; 5152 p.y = this.phi0; 5153 return p; 5154 } 5155 sCe = 2. * Math.asin(.5 * rho / this.rq);
5156 cCe = Math.cos(sCe); 5157 x *= (sCe = Math.sin(sCe)); 5158 if (this.mode == this.OBLIQ) { 5159 ab = cCe * this.sinb1 + y * sCe * this.cosb1 / rho 5160 q = this.qp * ab; 5161 y = rho * this.cosb1 * cCe - y * this.sinb1 * sCe; 5162 } else { 5163 ab = y * sCe / rho; 5164 q = this.qp * ab; 5165 y = rho * cCe; 5166 } 5167 break; 5168 case this.N_POLE: 5169 y = -y; 5170 case this.S_POLE: 5171 q = (x * x + y * y); 5172 if (!q ) { 5173 p.x = 0.; 5174 p.y = this.phi0; 5175 return p; 5176 } 5177 /* 5178 q = this.qp - q; 5179 */ 5180 ab = 1. - q / this.qp; 5181 if (this.mode == this.S_POLE) { 5182 ab = - ab; 5183 } 5184 break; 5185 } 5186 lam = Math.atan2(x, y); 5187 phi = this.authlat(Math.asin(ab), this.apa); 5188 } 5189 5190 /* 5191 var Rh = Math.Math.sqrt(p.x *p.x +p.y * p.y); 5192 var temp = Rh / (2.0 * this.a); 5193 5194 if (temp > 1) { 5195 Proj4js.reportError("laea:Inv:DataError"); 5196 return null; 5197 } 5198 5199 var z = 2.0 * Proj4js.common.asinz(temp); 5200 var sin_z=Math.sin(z); 5201 var cos_z=Math.cos(z); 5202 5203 var lon =this.long0; 5204 if (Math.abs(Rh) > Proj4js.common.EPSLN) { 5205 var lat = Proj4js.common.asinz(this.sin_lat_o * cos_z +this. cos_lat_o * sin_z *p.y / Rh); 5206 var temp =Math.abs(this.lat0) - Proj4js.common.HALF_PI; 5207 if (Math.abs(temp) > Proj4js.common.EPSLN) { 5208 temp = cos_z -this.sin_lat_o * Math.sin(lat); 5209 if(temp!=0.0) lon=Proj4js.common.adjust_lon(this.long0+Math.atan2(p.x*sin_z*this.cos_lat_o,temp*Rh)); 5210 } else if (this.lat0 < 0.0) { 5211 lon = Proj4js.common.adjust_lon(this.long0 - Math.atan2(-p.x,p.y)); 5212 } else { 5213 lon = Proj4js.common.adjust_lon(this.long0 + Math.atan2(p.x, -p.y)); 5214 } 5215 } else { 5216 lat = this.lat0; 5217 } 5218 */ 5219 //return(OK); 5220 p.x = Proj4js.common.adjust_lon(this.long0+lam); 5221 p.y = phi; 5222 return p; 5223 },//lamazInv() 5224 5225/* determine latitude from authalic latitude */ 5226 P00: .33333333333333333333, 5227 P01: .17222222222222222222, 5228 P02: .10257936507936507936, 5229 P10: .06388888888888888888, 5230 P11: .06640211640211640211, 5231 P20: .01641501294219154443, 5232 5233 authset: function(es) { 5234 var t; 5235 var APA = new Array(); 5236 APA[0] = es * this.P00; 5237 t = es * es; 5238 APA[0] += t * this.P01; 5239 APA[1] = t * this.P10; 5240 t *= es; 5241 APA[0] += t * this.P02; 5242 APA[1] += t * this.P11; 5243 APA[2] = t * this.P20; 5244 return APA; 5245 }, 5246 5247 authlat: function(beta, APA) { 5248 var t = beta+beta; 5249 return(beta + APA[0] * Math.sin(t) + APA[1] * Math.sin(t+t) + APA[2] * Math.sin(t+t+t)); 5250 } 5251 5252}; 5253 5254 5255 5256/* ====================================================================== 5257 projCode/aeqd.js 5258 ====================================================================== */ 5259 5260Proj4js.Proj.aeqd = { 5261 5262 init : function() { 5263 this.sin_p12=Math.sin(this.lat0); 5264 this.cos_p12=Math.cos(this.lat0); 5265 }, 5266 5267 forward: function(p) { 5268 var lon=p.x; 5269 var lat=p.y; 5270 var ksp; 5271 5272 var sinphi=Math.sin(p.y); 5273 var cosphi=Math.cos(p.y); 5274 var dlon = Proj4js.common.adjust_lon(lon - this.long0); 5275 var coslon = Math.cos(dlon); 5276 var g = this.sin_p12 * sinphi + this.cos_p12 * cosphi * coslon; 5277 if (Math.abs(Math.abs(g) - 1.0) < Proj4js.common.EPSLN) { 5278 ksp = 1.0; 5279 if (g < 0.0) { 5280 Proj4js.reportError("aeqd:Fwd:PointError"); 5281 return; 5282 } 5283 } else { 5284 var z = Math.acos(g); 5285 ksp = z/Math.sin(z); 5286 } 5287 p.x = this.x0 + this.a * ksp * cosphi * Math.sin(dlon); 5288 p.y = this.y0 + this.a * ksp * (this.cos_p12 * sinphi - this.sin_p12 * cosphi * coslon); 5289 return p; 5290 }, 5291 5292 inverse: function(p){ 5293 p.x -= this.x0; 5294 p.y -= this.y0; 5295 5296 var rh = Math.sqrt(p.x * p.x + p.y *p.y); 5297 if (rh > (2.0 * Proj4js.common.HALF_PI * this.a)) { 5298 Proj4js.reportError("aeqdInvDataError"); 5299 return; 5300 } 5301 var z = rh / this.a; 5302 5303 var sinz=Math.sin(z); 5304 var cosz=Math.cos(z); 5305 5306 var lon = this.long0; 5307 var lat; 5308 if (Math.abs(rh) <= Proj4js.common.EPSLN) { 5309 lat = this.lat0; 5310 } else { 5311 lat = Proj4js.common.asinz(cosz * this.sin_p12 + (p.y * sinz * this.cos_p12) / rh); 5312 var con = Math.abs(this.lat0) - Proj4js.common.HALF_PI; 5313 if (Math.abs(con) <= Proj4js.common.EPSLN) { 5314 if (this.lat0 >= 0.0) { 5315 lon = Proj4js.common.adjust_lon(this.long0 + Math.atan2(p.x , -p.y)); 5316 } else { 5317 lon = Proj4js.common.adjust_lon(this.long0 - Math.atan2(-p.x , p.y)); 5318 } 5319 } else { 5320 con = cosz - this.sin_p12 * Math.sin(lat); 5321 if ((Math.abs(con) < Proj4js.common.EPSLN) && (Math.abs(p.x) < Proj4js.common.EPSLN)) { 5322 //no-op, just keep the lon value as is 5323 } else { 5324 var temp = Math.atan2((p.x * sinz * this.cos_p12), (con * rh)); 5325 lon = Proj4js.common.adjust_lon(this.long0 + Math.atan2((p.x * sinz * this.cos_p12), (con * rh))); 5326 } 5327 } 5328 } 5329 5330 p.x = lon; 5331 p.y = lat; 5332 return p; 5333 } 5334}; 5335/* ====================================================================== 5336 projCode/moll.js 5337 ====================================================================== */ 5338
5339/******************************************************************************* 5340NAME MOLLWEIDE 5341 5342PURPOSE: Transforms input longitude and latitude to Easting and 5343 Northing for the MOllweide projection. The 5344 longitude and latitude must be in radians. The Easting 5345 and Northing values will be returned in meters. 5346 5347PROGRAMMER DATE 5348---------- ---- 5349D. Steinwand, EROS May, 1991; Updated Sept, 1992; Updated Feb, 1993 5350S. Nelson, EDC Jun, 2993; Made corrections in precision and 5351 number of iterations. 5352 5353ALGORITHM REFERENCES 5354 53551. Snyder, John P. and Voxland, Philip M., "An Album of Map Projections", 5356 U.S. Geological Survey Professional Paper 1453 , United State Government 5357 Printing Office, Washington D.C., 1989. 5358 53592. Snyder, John P., "Map Projections--A Working Manual", U.S. Geological 5360 Survey Professional Paper 1395 (Supersedes USGS Bulletin 1532), United 5361 State Government Printing Office, Washington D.C., 1987. 5362*******************************************************************************/ 5363 5364Proj4js.Proj.moll = { 5365 5366 /* Initialize the Mollweide projection 5367 ------------------------------------*/ 5368 init: function(){ 5369 //no-op 5370 }, 5371 5372 /* Mollweide forward equations--mapping lat,long to x,y 5373 ----------------------------------------------------*/ 5374 forward: function(p) { 5375 5376 /* Forward equations 5377 -----------------*/ 5378 var lon=p.x; 5379 var lat=p.y; 5380 5381 var delta_lon = Proj4js.common.adjust_lon(lon - this.long0); 5382 var theta = lat; 5383 var con = Proj4js.common.PI * Math.sin(lat); 5384 5385 /* Iterate using the Newton-Raphson method to find theta 5386 -----------------------------------------------------*/ 5387 for (var i=0;true;i++) { 5388 var delta_theta = -(theta + Math.sin(theta) - con)/ (1.0 + Math.cos(theta)); 5389 theta += delta_theta; 5390 if (Math.abs(delta_theta) < Proj4js.common.EPSLN) break; 5391 if (i >= 50) { 5392 Proj4js.reportError("moll:Fwd:IterationError"); 5393 //return(241); 5394 } 5395 } 5396 theta /= 2.0; 5397 5398 /* If the latitude is 90 deg, force the x coordinate to be "0 + false easting" 5399 this is done here because of precision problems with "cos(theta)" 5400 --------------------------------------------------------------------------*/ 5401 if (Proj4js.common.PI/2 - Math.abs(lat) < Proj4js.common.EPSLN) delta_lon =0; 5402 var x = 0.900316316158 * this.a * delta_lon * Math.cos(theta) + this.x0; 5403 var y = 1.4142135623731 * this.a * Math.sin(theta) + this.y0; 5404 5405 p.x=x; 5406 p.y=y; 5407 return p; 5408 }, 5409 5410 inverse: function(p){ 5411 var theta; 5412 var arg; 5413 5414 /* Inverse equations 5415 -----------------*/ 5416 p.x-= this.x0; 5417 //~ p.y -= this.y0; 5418 var arg = p.y / (1.4142135623731 * this.a); 5419 5420 /* Because of division by zero problems, 'arg' can not be 1.0. Therefore 5421 a number very close to one is used instead. 5422 -------------------------------------------------------------------*/ 5423 if(Math.abs(arg) > 0.999999999999) arg=0.999999999999; 5424 var theta =Math.asin(arg); 5425 var lon = Proj4js.common.adjust_lon(this.long0 + (p.x / (0.900316316158 * this.a * Math.cos(theta)))); 5426 if(lon < (-Proj4js.common.PI)) lon= -Proj4js.common.PI; 5427 if(lon > Proj4js.common.PI) lon= Proj4js.common.PI; 5428 arg = (2.0 * theta + Math.sin(2.0 * theta)) / Proj4js.common.PI; 5429 if(Math.abs(arg) > 1.0)arg=1.0; 5430 var lat = Math.asin(arg); 5431 //return(OK); 5432 5433 p.x=lon; 5434 p.y=lat; 5435 return p; 5436 } 5437}; 5438
Line numbers count LF bytes from the start of the resource, as the search results do. Vendor segments are library code the classifier recognised; they are stored but not indexed. Bytes are shown as Latin1 characters, one per byte.