1# 2# This file is part of: 3# 4# gpsman --- GPS Manager: a manager for GPS receiver data 5# 6# Copyright (c) 1998-2013 Miguel Filgueiras migfilg@t-online.de 7# 8# This program is free software; you can redistribute it and/or modify 9# it under the terms of the GNU General Public License as published by 10# the Free Software Foundation; either version 3 of the License, or 11# (at your option) any later version. 12# 13# This program is distributed in the hope that it will be useful, 14# but WITHOUT ANY WARRANTY; without even the implied warranty of 15# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 16# GNU General Public License for more details. 17# 18# You should have received a copy of the GNU General Public License 19# along with this program. 20# 21# File: projs_main.tcl 22# Last change: 6 October 2013 23# 24 25# Includes contributions by 26# - Sandor Laky (laky.sandor_at_freemail.hu) marked "SL contribution" 27# - Miguel Filgueiras (to code by others) marked "MF contribution/change" 28 29### most of this code used to be at the end of file projections.tcl 30 31### general initialization procedure for all main projections 32# see also the comments before the procedures for UTM in projections.tcl 33 34proc ProjInit {proj data datum ps} { 35 # changes here must be reflected on proc BadProjSetup below 36 # initialize projection parameters 37 # $ps is a list of latd,longd,datum to be projected 38 # should call proc ProjParams if its parameters may be changed by 39 # the user and if $ASKPROJPARAMS is true 40 # should call proc to compute auxiliary parameters if any 41 # must return projection of first position in $ps 42 global $data TMPARAM LNTFPARAMS ASKPROJPARAMS 43 44 foreach "latd longd posdatum" [lindex $ps 0] { break } 45 if { $posdatum != $datum } { 46 foreach "latd longd" \ 47 [ToDatum $latd $longd $posdatum $datum] { break } 48 } 49 set ${data}(datum) $datum 50 set askauxcomp 1 51 switch $proj { 52 UTM { 53 # only the first position is used 54 # no parameters can be changed by the user 55 set askauxcomp 0 56 set p [DegreesToUTM $latd $longd $datum] 57 scan [lindex $p 0] %0d ze 58 set zn [lindex $p 1] 59 set ${data}(UTMzone) [list $ze $zn] 60 set ${data}(m_0) [expr -183+6.0*$ze] 61 set res [list [lindex $p 2] [lindex $p 3]] 62 } 63 TM { 64 # take averages of latitudes and of longitudes as central 65 # coordinates, and use k0=0.9996 66 foreach "${data}(lat0) ${data}(long0)" \ 67 [AverageLatLong $latd $longd $datum $ps] {} 68 set ${data}(k0) 0.9996 69 if { $ASKPROJPARAMS } { 70 ProjParams change TM $data 71 } 72 set res [ConvToTM $latd $longd [set ${data}(lat0)] \ 73 [set ${data}(long0)] [set ${data}(k0)] $datum] 74 set askauxcomp 0 75 } 76 BMN - CTR - GKK - KKJP - TWG { 77 # special cases of TM 78 # take averages of longitudes to find long0, and use lat0=0 79 set long0 [AverageLong $longd $datum $ps] 80 foreach "pr lmin lmax delta k0" $TMPARAM($proj) {} 81 if { $long0 < $lmin } { 82 set long0 $lmin 83 } elseif { $long0 > $lmax } { 84 set long0 $lmax 85 } else { 86 set long0 [expr int(1.0*($long0-$lmin+0.5*$delta)/ \ 87 $delta)*$delta+$lmin] 88 } 89 set ${data}(${pr}long0) $long0 90 set ${data}(m_0) $k0 91 if { $ASKPROJPARAMS } { 92 ProjParams change $proj $data 93 } 94 set res [ConvToTM $latd $longd 0 [set ${data}(${pr}long0)] \ 95 $k0 $datum] 96 set askauxcomp 0 97 } 98 LCC1 - Stereogr - SOM { 99 # Lambert Conformal Conic single parallel, 100 # Stereographic and Swiss Oblique Mercator projections 101 # take averages of latitudes and of longitudes as central 102 # coordinates, and use k0=1 103 # in Stereographic central latitude determines the case: polar, 104 # oblique, equatorial 105 foreach "${data}(lat0) ${data}(long0)" \ 106 [AverageLatLong $latd $longd $datum $ps] {} 107 set ${data}(k0) 1.0 108 } 109 LCC2 { 110 # Lambert Conformal Conic projection, two parallels 111 # take extremes of latitudes as latitudes of standard parallels, 112 # averages of latitudes and longitudes as false origin coordinates 113 foreach "${data}(latF) ${data}(longF) \ 114 ${data}(lat1) ${data}(lat2)" \ 115 [AverageExtrLatLong $latd $longd $datum $ps] {} 116 } 117 LambNTF { 118 # as LCC2 but with zones with specific parameters 119 # only the first position is used 120 # no parameters can be changed by the user 121 set askauxcomp 0 122 set p [DegreesToLambNTF $latd $longd $datum] 123 if { [set z [lindex $p 0]] == "--" } { 124 set z II 125 } 126 global LNTFz$z 127 array set $data [array get LNTFz$z] 128 set res [ProjLambNTFPoint $data $latd $longd $datum] 129 } 130 AEA { 131 # Albers Equal Area Conic projection 132 # take extremes of latitudes as latitudes of standard parallels, 133 # average of longitudes as central longitude, 0 for central 134 # latitude 135 foreach "x ${data}(long0) lamax lamin" \ 136 [AverageExtrLatLong $latd $longd $datum $ps] {} 137 if { abs($lamax) < abs($lamin) } { 138 set x $lamax ; set lamax $lamin ; set lamin $x 139 } 140 set ${data}(lat1) $lamax ; set ${data}(lat2) $lamin 141 set ${data}(lat0) 0 142 } 143 LEAC { 144 # Lambert Equal Area Conic projection 145 # take extreme of latitudes as latitude of standard parallel, 146 # average of longitudes as central longitude, North polar 147 # aspect if extreme latitude is nonnegative, 0 for central 148 # latitude 149 foreach "x ${data}(long0) lamax lamin" \ 150 [AverageExtrLatLong $latd $longd $datum $ps] {} 151 if { abs($lamax) < abs($lamin) } { 152 set lamax $lamin 153 } 154 set ${data}(lat1) $lamax 155 if { $lamax < 0 } { 156 set ${data}(polasp) south 157 } else { set ${data}(polasp) north } 158 set ${data}(lat0) 0 159 } 160 Merc1 { 161 # Mercator projection, single parallel, ellipsoidal case 162 # take average of longitudes as central longitude, and use k0=1 163 set ${data}(long0) [AverageLong $longd $datum $ps] 164 set ${data}(k0) 1.0 165 } 166 Merc2 { 167 # Mercator projection, two parallels 168 # take average of longitudes as central longitude and extreme 169 # latitude (in absolute value) as first standard parallel 170 foreach "x ${data}(long0) lamin lamax" \ 171 [AverageExtrLatLong $latd $longd $datum $ps] {} 172 if { abs($lamin) > $lamax } { set lamax [expr abs($lamin)] } 173 set ${data}(lat1) $lamax 174 } 175 SphMerc - CS { 176 # Mercator projection, single parallel, spherical case 177 # and Cassini-Soldner projection 178 # take average of longitudes and longitudes as central coordinates 179 foreach "${data}(lat0) ${data}(long0)" \ 180 [AverageLatLong $latd $longd $datum $ps] {} 181 } 182 APOLY - EqCyl { 183 # American Polyconic and Equidistant Cylindrical projections 184 # take averages of latitudes as standard latitude 185 set ${data}(lat0) [AverageLat $latd $datum $ps] 186 } 187 Schreiber { 188 # Schreiber projection 189 # all parameters are fixed, except the datum that must use the 190 # "Bessel 1841" ellipsoid 191 if { [EllipsdOf $datum] != "Bessel 1841" } { 192 GMMessage "Datum: Rijks Driehoeksmeting" 193 set datum "Rijks Driehoeksmeting" 194 } 195 ProjSchreiberComputeAux $data $datum 196 set res [ProjSchreiberPoint $data $latd $longd $datum] 197 set askauxcomp 0 198 } 199 EOV { 200 # Hungarian grid projection: no parameters needed 201 set res [ProjEOVPoint {} $latd $longd $datum] 202 } 203 default { 204 BUG "no initialization proc for projection" 205 set askauxcomp 0 206 } 207 } 208 if { $askauxcomp } { 209 if { $ASKPROJPARAMS } { 210 ProjParams change $proj $data 211 } 212 Proj${proj}ComputeAux $data $datum 213 set res [Proj${proj}Point $data $latd $longd $datum] 214 } 215 return $res 216} 217 218proc BadProjSetup {data proj latd longd datum} { 219 # this must be kept in line with proc ProjInit (above) 220 # initialize projection for use in command mode 221 # $data is array for the parameters with external parameters and 222 # datum defined, except in the case of "UTM" and "EOV" that only 223 # have datum 224 # return 0 on success, 1 on error 225 global ASKPROJPARAMS $data 226 227 set ASKPROJPARAMS 0 228 set aux 0 229 switch $proj { 230 UTM { 231 if { [catch {set ${data}(UTMzone)}] } { 232 set p [DegreesToUTM $latd $longd $datum] 233 scan [lindex $p 0] %0d ze 234 set zn [lindex $p 1] 235 set ${data}(UTMzone) [list $ze $zn] 236 set ${data}(m_0) [expr -183+6.0*$ze] 237 } 238 } 239 LambNTF { 240 # as LCC2 but with zones with specific parameters 241 # no parameters can be changed by the user 242 set p [DegreesToLambNTF $latd $longd $datum] 243 if { [set z [lindex $p 0]] == "--" } { 244 set z II 245 } 246 global LNTFz$z 247 catch {unset $data} 248 array set $data [array get LNTFz$z] 249 } 250 Schreiber { 251 # Schreiber projection 252 # all parameters are fixed, except the datum that must use the 253 # "Bessel 1841" ellipsoid 254 if { [EllipsdOf [set ${data}(datum)]] != "Bessel 1841" } { 255 set ${data}(datum) "Rijks Driehoeksmeting" 256 } 257 } 258 TM - BMN - CTR - EOV - GKK - KKJP - TWG { 259 } 260 LCC1 - Stereogr - SOM - LCC2 - AEA - 261 LEAC - Merc1 - Merc2 - SphMerc - CS - APOLY - EqCyl { 262 set aux 1 263 } 264 default { 265 return 1 266 } 267 } 268 if { $aux } { 269 Proj${proj}ComputeAux $data [set ${data}(datum)] 270 } 271 return 0 272} 273 274proc AverageExtrLatLong {latd longd datum ps} { 275 # compute the averages of latitudes and longitudes and the minimum and 276 # maximum latitude for $latd,$longd and all positions in $ps except 277 # the first 278 # return list with the two averages followed by the maximum and the minimum 279 280 set n 1 281 set sla [set lamx [set lamn $latd]] ; set slo $longd 282 foreach p [lreplace $ps 0 0] { 283 foreach "la lo posdatum" $p { break } 284 if { $posdatum != $datum } { 285 foreach "la lo" \ 286 [ToDatum $la $lo $posdatum $datum] { break } 287 } 288 set sla [expr $sla+$la] ; set slo [expr $slo+$lo] 289 if { $la > $lamx } { set lamx $la } 290 if { $la < $lamn } { set lamn $la } 291 incr n 292 } 293 return [list [expr 1.0*$sla/$n] [expr 1.0*$slo/$n] $lamx $lamn] 294} 295 296proc AverageLatLong {latd longd datum ps} { 297 # compute the averages of latitudes and longitudes for $latd,$longd and 298 # all positions in $ps except the first 299 # return list with the two averages 300 301 set n 1 302 foreach p [lreplace $ps 0 0] { 303 foreach "la lo posdatum" $p { break } 304 if { $posdatum != $datum } { 305 foreach "la lo" \ 306 [ToDatum $la $lo $posdatum $datum] { break } 307 } 308 set latd [expr $latd+$la] ; set longd [expr $longd+$lo] 309 incr n 310 } 311 return [list [expr 1.0*$latd/$n] [expr 1.0*$longd/$n]] 312} 313 314proc AverageLat {latd datum ps} { 315 # compute the average of latitudes for $latd and 316 # all positions in $ps except the first 317 318 set n 1 319 foreach p [lreplace $ps 0 0] { 320 set la [lindex $p 0] ; set pdatum [lindex $p 2] 321 if { $pdatum != $datum } { 322 set la [lindex [ToDatum $la [lindex $p 1] $pdatum $datum] 0] 323 } 324 set latd [expr $la+$latd] 325 incr n 326 } 327 return [expr 1.0*$latd/$n] 328} 329 330proc AverageLong {longd datum ps} { 331 # compute the average of longitudes for $longd and 332 # all positions in $ps except the first 333 334 set n 1 335 foreach p [lreplace $ps 0 0] { 336 set lo [lindex $p 1] ; set pdatum [lindex $p 2] 337 if { $pdatum != $datum } { 338 set lo [lindex [ToDatum [lindex $p 0] $lo $pdatum $datum] 1] 339 } 340 set longd [expr $longd+$lo] 341 incr n 342 } 343 return [expr 1.0*$longd/$n] 344} 345 346## converting longitudes to -180, 180 (inclusive) or -360, 360 (inclusive) 347 348proc NormalLong {long} { 349 # convert longitude in decimal degrees to the -180,180 (inclusive) range 350 351 if { abs([set x [expr abs($long)*0.00555555555555555556]])-1 > 1e-14 } { 352 set x [expr 0.5*($x+1)] 353 set x [expr (($x-floor($x))-0.5)*360] 354 } else { set x [expr abs($long)] } 355 if { $long < 0 && abs($x-180) > 1e-14 } { return [expr -$x] } 356 return $x 357} 358 359proc NormalLongCentred {long longc} { 360 # convert longitude in decimal degrees to the closest value to $longc in 361 # the -360,360 (inclusive) range 362 363 set long [NormalLong $long] 364 if { $long <= $longc-180 } { 365 set long [expr $long+360] 366 } elseif { $long > $longc+180 } { 367 set long [expr $long-360] 368 } 369 return $long 370} 371 372## projection procedures for UTM projection 373 374# there may be auxiliary parameters to be computed from the main ones 375 376proc ProjUTMComputeAux {data args} { 377 # compute auxiliary parameters from main parameters 378 # $data is name of global array for the parameters 379 # $args is not used but is needed because of uniform way of calling 380 # these procedures 381 global $data 382 383 set z [set ${data}(UTMzone)] 384 # set the zone to a list of number and letter if it is a string 385 if { [regexp {^([0-9]+)([A-Z])$} $z x ze zn] } { 386 set ${data}(UTMzone) [list $ze $zn] 387 } else { set ze [lindex $z 0] } 388 set ${data}(m_0) [expr -183+6.0*$ze] 389 return 390} 391 392# projection proc is a function of latd,longd,datum to planar Cartesian 393# coordinates x,y in the terrain 394 395proc ProjUTMPoint {data latd longd datum} { 396 # compute planar Cartesian coordinates for given position 397 global $data 398 399 set prdatum [set ${data}(datum)] 400 if { $datum != $prdatum } { 401 foreach "latd longd" \ 402 [ToDatum $latd $longd $datum $prdatum] { break } 403 } 404 foreach "pze pzn px py" [DegreesToUTM $latd $longd $prdatum] {} 405 foreach "mze mzn" [set ${data}(UTMzone)] {} 406 if { $mze!=$pze || $mzn != $pzn } { 407 foreach "px py" [CompUTMOnZone $latd $longd [set ${data}(m_0)]] { 408 break 409 } 410 } 411 return [list $px $py] 412} 413 414# inverse projection proc is a function of planar Cartesian coordinates in the 415# terrain to list with latitude and longitude in the datum of the projection 416 417proc ProjUTMInvert {data x y} { 418 # return list with latitude and longitude for point at $x,$y in the terrain 419 global $data 420 421 foreach "ze zn" [set ${data}(UTMzone)] {} 422 set datum [set ${data}(datum)] 423 set p [UTMToDegrees $ze $zn $x $y $datum] 424 return $p 425} 426 427## Transverse Mercator projection 428 429proc ProjTMPoint {data latd longd datum} { 430 # compute planar Cartesian coordinates for given position 431 global $data 432 433 set prdatum [set ${data}(datum)] 434 if { $datum != $prdatum } { 435 foreach "latd longd" \ 436 [ToDatum $latd $longd $datum $prdatum] { break } 437 } 438 return [ConvToTM $latd $longd [set ${data}(lat0)] \ 439 [set ${data}(long0)] [set ${data}(k0)] $prdatum] 440} 441 442proc ProjTMInvert {data x y} { 443 # return list with latitude and longitude for point at $x,$y in the terrain 444 global $data 445 446 set datum [set ${data}(datum)] 447 set p [ConvFromTM $x $y [set ${data}(lat0)] [set ${data}(long0)] \ 448 [set ${data}(k0)] $datum] 449 return $p 450} 451 452## German Grid (Gauss-Krueger-Koordinatensystem) projection 453# Information provided by Andreas Lange (Andreas.C.Lange_at_GMX.de) 454# This is a Transverse Mercator projection with zone in [0-6], lat0=0, 455# lon0=zone*3 (0, 3, 6, 9, 12, 15E), and scale factor at central meridian 456# of k0=1.0 457 458# Basic Finnish Grid (KKJP) projection 459# Similar with zone in [1-4], lon0=zone*3+18 (21, 24, 27, 30E) 460 461# Taiwan Grid projection 462# Information provided by Dan Jacobson (jidanni_at_yahoo.com.tw) 463# Similar with zone in [1-6], lon0=zone*2+113 (115, 117, ..., 125E), 464# and k0=0.9999 465 466# Carta Tecnica Regionale (CTR), Italian projection 467# "Le carte topografiche CTR ed il loro uso GPS" 468# (http://www.gpscomefare.com/guide/tutorialgps/mapdatum.htm) 469# May 2003 (information kindly sent by Alessandro Palmas) 470# Similar with zone in [1-2], lon0=zone*6+3 (9, 15), k0=0.9996 471 472# Austrian BMN grid projection 473# information kindly sent by Alessandro Palmas, July 2003 474# Similar with zones M28, M31, M34, lon0=zone_index*3+10.333333, k0=1 475 476proc ProjGKKComputeAux {data args} { 477 # compute auxiliary parameters from main parameters 478 # $data is name of global array for the parameters 479 # $args is not used but is needed because of uniform way of calling 480 # these procedures 481 global $data 482 483 set ${data}(m_0) 1.0 484 return 485} 486 487proc ProjKKJPComputeAux {data args} { 488 # compute auxiliary parameters from main parameters 489 # $data is name of global array for the parameters 490 # $args is not used but is needed because of uniform way of calling 491 # these procedures 492 global $data 493 494 set ${data}(m_0) 1.0 495 return 496} 497 498proc ProjTWGComputeAux {data args} { 499 # compute auxiliary parameters from main parameters 500 # $data is name of global array for the parameters 501 # $args is not used but is needed because of uniform way of calling 502 # these procedures 503 global $data 504 505 set ${data}(m_0) 0.9999 506 return 507} 508 509proc ProjCTRComputeAux {data args} { 510 # compute auxiliary parameters from main parameters 511 # $data is name of global array for the parameters 512 # $args is not used but is needed because of uniform way of calling 513 # these procedures 514 global $data 515 516 set ${data}(m_0) 0.9996 517 return 518} 519 520proc ProjBMNComputeAux {data args} { 521 # compute auxiliary parameters from main parameters 522 # $data is name of global array for the parameters 523 # $args is not used but is needed because of uniform way of calling 524 # these procedures 525 global $data 526 527 set ${data}(m_0) 1.0 528 return 529} 530 531proc ProjGKKPoint {data latd longd datum} { 532 # compute planar Cartesian coordinates for given position 533 534 return [Proj_TMZ_Point gkk $data $latd $longd $datum] 535} 536 537proc ProjKKJPPoint {data latd longd datum} { 538 # compute planar Cartesian coordinates for given position 539 540 return [Proj_TMZ_Point kkjp $data $latd $longd $datum] 541} 542 543proc ProjTWGPoint {data latd longd datum} { 544 # compute planar Cartesian coordinates for given position 545 546 return [Proj_TMZ_Point twg $data $latd $longd $datum] 547} 548 549proc ProjCTRPoint {data latd longd datum} { 550 # compute planar Cartesian coordinates for given position 551 552 return [Proj_TMZ_Point ctr $data $latd $longd $datum] 553} 554 555proc ProjBMNPoint {data latd longd datum} { 556 # compute planar Cartesian coordinates for given position 557 558 return [Proj_TMZ_Point bmn $data $latd $longd $datum] 559} 560 561proc Proj_TMZ_Point {pr data latd longd datum} { 562 # compute planar Cartesian coordinates for given position 563 global $data 564 565 set prdatum [set ${data}(datum)] 566 if { $datum != $prdatum } { 567 foreach "latd longd" \ 568 [ToDatum $latd $longd $datum $prdatum] { break } 569 } 570 return [ConvToTM $latd $longd 0 [set ${data}(${pr}long0)] \ 571 [set ${data}(m_0)] $prdatum] 572} 573 574proc ProjGKKInvert {data x y} { 575 # return list with latitude and longitude for point at $x,$y in the terrain 576 577 return [Proj_TMZ_Invert gkk $data $x $y] 578} 579 580proc ProjKKJPInvert {data x y} { 581 # return list with latitude and longitude for point at $x,$y in the terrain 582 583 return [Proj_TMZ_Invert kkjp $data $x $y] 584} 585 586proc ProjTWGInvert {data x y} { 587 # return list with latitude and longitude for point at $x,$y in the terrain 588 589 return [Proj_TMZ_Invert twg $data $x $y] 590} 591 592proc ProjCTRInvert {data x y} { 593 # return list with latitude and longitude for point at $x,$y in the terrain 594 595 return [Proj_TMZ_Invert ctr $data $x $y] 596} 597 598proc ProjBMNInvert {data x y} { 599 # return list with latitude and longitude for point at $x,$y in the terrain 600 601 return [Proj_TMZ_Invert bmn $data $x $y] 602} 603 604proc Proj_TMZ_Invert {pr data x y} { 605 # return list with latitude and longitude for point at $x,$y in the terrain 606 global $data 607 608 set datum [set ${data}(datum)] 609 set p [ConvFromTM $x $y 0 [set ${data}(${pr}long0)] \ 610 [set ${data}(m_0)] $datum] 611 return $p 612} 613 614## Lambert Conic Conformal projections 615 616# single standard parallel 617 618proc ProjLCC1ComputeAux {data datum} { 619 # compute auxiliary parameters from main parameters 620 # $data is name of global array for the parameters 621 global $data 622 623 set d [EllipsdData $datum] 624 set a [lindex $d 0] ; set f [lindex $d 1] 625 set es [expr $f*(2-$f)] 626 set e [set ${data}(m_e) [expr sqrt($es)]] 627 set phi0 [expr [set ${data}(lat0)]*0.01745329251994329576] 628 if { [set n [set ${data}(lcc_n) [expr sin($phi0)]]] < 0 } { 629 set ${data}(lcc_sn) -1 630 } else { set ${data}(lcc_sn) 1 } 631 set m0 [expr cos($phi0)/sqrt(1-$es*$n*$n)] 632 set t0n [expr pow([ExpRedLat $phi0 $n $e],$n)] 633 set F [expr $m0/($n*$t0n)] 634 set aFk0 [set ${data}(m_a) [expr $a*$F*[set ${data}(k0)]]] 635 set ${data}(lcc_rho0F) [expr $aFk0*$t0n] 636 return 637} 638 639proc ProjLCC1Point {data latd longd datum} { 640 # compute planar Cartesian coordinates for given position 641 global $data 642 643 set prdatum [set ${data}(datum)] 644 if { $datum != $prdatum } { 645 foreach "latd longd" \ 646 [ToDatum $latd $longd $datum $prdatum] { break } 647 } 648 set long0 [set ${data}(long0)] 649 set longd [NormalLongCentred $longd $long0] 650 set phi [expr $latd*0.01745329251994329576] 651 set theta [expr [set ${data}(lcc_n)]* \ 652 ($longd-$long0)*0.01745329251994329576] 653 set t [ExpRedLat $phi [expr sin($phi)] [set ${data}(m_e)]] 654 set rho [expr [set ${data}(m_a)]* \ 655 pow($t,[set ${data}(lcc_n)])] 656 set x [expr $rho*sin($theta)] 657 set y [expr [set ${data}(lcc_rho0F)]-$rho*cos($theta)] 658 return [list $x $y] 659} 660 661proc ProjLCC1Invert {data x y} { 662 # return list with latitude and longitude for point at $x,$y in the terrain 663 global $data 664 665 set z [expr [set ${data}(lcc_rho0F)]-$y] 666 set thetap [expr atan(1.0*$x/$z)] 667 set rhop [expr [set ${data}(lcc_sn)]*[Hypot $x $z]] 668 set tp [expr pow($rhop/([set ${data}(m_a)]), \ 669 1.0/[set ${data}(lcc_n)])] 670 set longd [expr $thetap/[set ${data}(lcc_n)]*57.29577951308232087684+ \ 671 [set ${data}(long0)]] 672 set longd [NormalLong $longd] 673 set latd [expr [LatFromRedLat $tp [set ${data}(m_e)]]* \ 674 57.29577951308232087684] 675 return [list $latd $longd] 676} 677 678# two standard parallels 679 680proc ProjLCC2ComputeAux {data datum} { 681 # compute auxiliary parameters from main parameters 682 # $data is name of global array for the parameters 683 global $data MESS 684 685 set d [EllipsdData $datum] 686 set a [lindex $d 0] ; set f [lindex $d 1] 687 set es [expr $f*(2-$f)] 688 set e [set ${data}(m_e) [expr sqrt($es)]] 689 set phi1 [expr [set ${data}(lat1)]*0.01745329251994329576] 690 set phi2 [expr [set ${data}(lat2)]*0.01745329251994329576] 691 set phiF [expr [set ${data}(latF)]*0.01745329251994329576] 692 set s1 [expr sin($phi1)] ; set s2 [expr sin($phi2)] 693 set t1 [ExpRedLat $phi1 $s1 $e] 694 set t2 [ExpRedLat $phi2 $s2 $e] 695 if { [catch {set m1 [expr cos($phi1)/sqrt(1-$es*$s1*$s1)]}] || \ 696 [catch {set m2 [expr cos($phi2)/sqrt(1-$es*$s2*$s2)]}] || \ 697 [catch {set n [expr (log($m1)-log($m2))/(log($t1)-log($t2))]}] } { 698 GMMessage $MESS(badProjargs) 699 ProjParams change LCC2 $data 700 ProjLCC2ComputeAux $data $datum 701 return 702 } 703 set ${data}(lcc_n) $n 704 if { $n < 0 } { 705 set ${data}(lcc_sn) -1 706 } else { set ${data}(lcc_sn) 1 } 707 set F [expr $m1/($n*pow($t1,$n))] 708 set aF [set ${data}(m_a) [expr $a*$F]] 709 set tF [ExpRedLat $phiF [expr sin($phiF)] $e] 710 set ${data}(lcc_rho0F) [expr $aF*pow($tF,$n)] 711 return 712} 713 714proc ProjLCC2Point {data latd longd datum} { 715 # compute planar Cartesian coordinates for given position 716 global $data 717 718 set prdatum [set ${data}(datum)] 719 if { $datum != $prdatum } { 720 foreach "latd longd" \ 721 [ToDatum $latd $longd $datum $prdatum] { break } 722 } 723 set longF [set ${data}(longF)] 724 set longd [NormalLongCentred $longd $longF] 725 set phi [expr $latd*0.01745329251994329576] 726 set theta [expr [set ${data}(lcc_n)]* \ 727 ($longd-$longF)*0.01745329251994329576] 728 set t [ExpRedLat $phi [expr sin($phi)] [set ${data}(m_e)]] 729 set rho [expr [set ${data}(m_a)]* \ 730 pow($t,[set ${data}(lcc_n)])] 731 set x [expr $rho*sin($theta)] 732 set y [expr [set ${data}(lcc_rho0F)]-$rho*cos($theta)] 733 return [list $x $y] 734} 735 736proc ProjLCC2Invert {data x y} { 737 # return list with latitude and longitude for point at $x,$y in the terrain 738 global $data 739 740 set z [expr [set ${data}(lcc_rho0F)]-$y] 741 set thetap [expr atan(1.0*$x/$z)] 742 set rhop [expr [set ${data}(lcc_sn)]*[Hypot $x $z]] 743 set tp [expr pow($rhop/([set ${data}(m_a)]), \ 744 1.0/[set ${data}(lcc_n)])] 745 set longd [expr $thetap/[set ${data}(lcc_n)]*57.29577951308232087684+ \ 746 [set ${data}(longF)]] 747 set longd [NormalLong $longd] 748 set latd [expr [LatFromRedLat $tp \ 749 [set ${data}(m_e)]]*57.29577951308232087684] 750 return [list $latd $longd] 751} 752 753# French NTF projection 754 755proc ProjLambNTFComputeAux {data datum} { 756 global $data 757 758 set z [set $data(NTFzone)] 759 global LNTFz$z 760 array set $data [array get LNTFz$z] 761 return 762} 763 764proc ProjLambNTFPoint {data latd longd datum} { 765 766 return [ProjLCC2Point $data $latd $longd $datum] 767} 768 769proc ProjLambNTFInvert {data x y} { 770 771 return [ProjLCC2Invert $data $x $y] 772} 773 774# ancillary procs for LCC projections 775 776proc ExpRedLat {phi sinphi e} { 777 # compute exp of reduced latitude for $phi 778 # $e is first eccentricity of ellipsoid 779 780 set se [expr $sinphi*$e] 781 return [expr tan((1.5707963267948966-$phi)/2.0) / \ 782 pow((1.0-$se)/(1.0+$se), $e/2.0)] 783} 784 785proc LatFromRedLat {ts e} { 786 # compute latitude from reduced latitude 787 788 set eccnth [expr $e/2.0] 789 set phi [expr 1.5707963267948966-2*atan($ts)] 790 set i 16 ; set d 1 791 while { abs($d) > 1e-10 && [incr i -1] } { 792 set con [expr $e*sin($phi)] 793 set d [expr 1.5707963267948966- \ 794 2*atan($ts*pow((1-$con)/(1+$con),$eccnth))-$phi] 795 set phi [expr $phi+$d] 796 } 797 return $phi 798} 799 800## Albers Equal Area Conic and Lambert Equal Area Conic projection 801# adapted from PROJ4.3.3 802 803proc ProjAEAComputeAux {data datum} { 804 # compute auxiliary parameters from main parameters 805 # $data is name of global array for the parameters 806 global $data 807 808 set phi1 [expr [set ${data}(lat1)]*0.01745329251994329576] 809 set phi2 [expr [set ${data}(lat2)]*0.01745329251994329576] 810 ComputeAux_AEA_LEAC AEA $data $datum $phi1 $phi2 811 return 812} 813 814proc ProjAEAPoint {data latd longd datum} { 815 # compute planar Cartesian coordinates for given position 816 global $data MESS 817 818 set prdatum [set ${data}(datum)] 819 if { $datum != $prdatum } { 820 foreach "latd longd" \ 821 [ToDatum $latd $longd $datum $prdatum] { break } 822 } 823 set long0 [set ${data}(long0)] 824 set longd [NormalLongCentred $longd $long0] 825 set a [set ${data}(m_a)] ; set n [set ${data}(m_1)] 826 set phi [expr $latd*0.01745329251994329576] 827 set lam [expr $n*($longd-$long0)*0.01745329251994329576] 828 set rho [expr [set ${data}(m_3)]-$n* \ 829 [SmallQ [expr sin($phi)] [set ${data}(m_e)] \ 830 [set ${data}(m_0)]]] 831 if { $rho < 0 } { 832 GMMessage $MESS(outrngproj) 833 return [list 0 0] 834 } 835 set rho [expr [set ${data}(m_4)]*sqrt($rho)] 836 set x [expr $a*$rho*sin($lam)] 837 set y [expr $a*([set ${data}(m_5)]-$rho*cos($lam))] 838 return [list $x $y] 839} 840 841proc ProjAEAInvert {data x y} { 842 # return list with latitude and longitude for point at $x,$y in the terrain 843 global $data MESS TXT 844 845 set a [set ${data}(m_a)] ; set n [set ${data}(m_1)] 846 set datum [set ${data}(datum)] 847 set x [expr $x/$a] 848 set y [expr [set ${data}(m_5)]-$y/$a] 849 if { [set rho [Hypot $x $y]] == 0 } { 850 set lam 0 851 if { $n > 0 } { 852 set phi 1.5707963267948966 853 } else { set phi -1.5707963267948966 } 854 } else { 855 if { $n < 0 } { 856 set rho [expr -$rho] ; set x [expr -$x] ; set y [expr -$y] 857 } 858 set phi [expr 1.0*$rho/[set ${data}(m_4)]] 859 set phi [expr 1.0*([set ${data}(m_3)]-$phi*$phi)/$n] 860 if { abs([set ${data}(m_2)]-abs($phi)) > 1e-7 } { 861 set phi [LatAngle $phi [set ${data}(m_e)] \ 862 [set ${data}(m_0)]] 863 if { $phi > 1e19 } { 864 ## return 0 180 on error... 865 GMMessage [format $MESS(badinvproj) $TXT(PRJAEA)/$TXT(PRJLEAC)] 866 return [list 0 180] 867 } 868 } elseif { $phi < 0 } { 869 set phi -1.5707963267948966 870 } else { set phi 1.5707963267948966 } 871 set lam [expr 1.0*atan2($x,$y)/$n] 872 } 873 set latd [expr $phi*57.29577951308232087684] 874 set longd [expr ($lam*57.29577951308232087684)+[set ${data}(long0)]] 875 return [list $latd [NormalLong $longd]] 876} 877 878proc ProjLEACComputeAux {data datum} { 879 # compute auxiliary parameters from main parameters 880 # $data is name of global array for the parameters 881 global $data 882 883 set phi2 [expr [set ${data}(lat1)]*0.01745329251994329576] 884 if { [set ${data}(polasp)] == "south" } { 885 set phi1 -1.5707963267948966 886 } else { set phi1 1.5707963267948966 } 887 ComputeAux_AEA_LEAC LEAC $data $datum $phi1 $phi2 888 return 889} 890 891proc ProjLEACPoint {data latd longd datum} { 892 893 return [ProjAEAPoint $data $latd $longd $datum] 894} 895 896proc ProjLEACInvert {data x y} { 897 898 return [ProjAEAInvert $data $x $y] 899} 900 901# ancillary procedures for AEA and LEAC projections 902 903proc ComputeAux_AEA_LEAC {proj data datum phi1 phi2} { 904 global $data MESS 905 906 if { abs($phi1+$phi2) < 1e-10 } { 907 GMMessage $MESS(badProjargs) 908 ProjParams change $proj $data 909 ProjAEAComputeAux $data $datum 910 return 911 } 912 set d [EllipsdData $datum] 913 set ${data}(m_a) [lindex $d 0] ; set f [lindex $d 1] 914 set es [expr $f*(2-$f)] 915 set e [set ${data}(m_e) [expr sqrt($es)]] 916 set oes [set ${data}(m_0) [expr 1-$es]] 917 set n [set ${data}(m_1) [expr sin($phi1)]] 918 set cosp [expr cos($phi1)] 919 set m1 [expr $cosp/sqrt(1.0-$es*$n*$n)] 920 set ml1 [SmallQ $n $e $oes] 921 if { abs($phi1-$phi2) > 1e-10 } { 922 # secant 923 set s2 [expr sin($phi2)] 924 set m2 [expr cos($phi2)/sqrt(1.0-$es*$s2*$s2)] 925 set ml2 [SmallQ $s2 $e $oes] 926 set n [set ${data}(m_1) [expr ($m1*$m1-$m2*$m2)/($ml2-$ml1)]] 927 } 928 set ${data}(m_2) [expr 1.0-0.5*$oes*log((1.0-$e)/(1.0+$e))/$e] 929 set c [set ${data}(m_3) [expr $m1*$m1+$n*$ml1]] 930 set dd [set ${data}(m_4) [expr 1.0/$n]] 931 set phi0 [expr [set ${data}(lat0)]*0.01745329251994329576] 932 set ${data}(m_5) [expr $dd*sqrt($c-$n*[SmallQ sin($phi0) $e $oes])] 933 return 934} 935 936proc LatAngle {qs te toes} { 937 # compute latitude angle phi-1 938 939 set phi [Aasin [expr 0.5*$qs]] 940 if { $te < 1e-7 } { return $phi } 941 for { set i 15 } { $i > 0 } { incr i -1 } { 942 set sinpi [expr sin($phi)] ; set cospi [expr cos($phi)] 943 set con [expr $te*$sinpi] 944 set com [expr 1.0-$con*$con] 945 set dphi [expr 0.5*$com*$com/$cospi* \ 946 ($qs/$toes-$sinpi/$com+0.5/$te* \ 947 log((1.0-$con)/(1.0+$con)))] 948 set phi [expr $phi+$dphi] 949 if { abs($dphi) < 1e-10 } { return $phi } 950 } 951 return 1e20 952} 953 954proc SmallQ {sinphi e oes} { 955 # compute small q 956 957 if { $e < 1e-7 } { return [expr $sinphi+$sinphi] } 958 set con [expr $e*$sinphi] 959 return [expr $oes*($sinphi/(1.0-$con*$con)- \ 960 (0.5/$e)*log((1.0-$con)/(1.0+$con)))] 961} 962 963## Mercator projection: a special case of Lambert Conic Conformal 964# with the equator as standard parallel 965 966# single standard parallel 967 968proc ProjMerc1ComputeAux {data datum} { 969 # compute auxiliary parameters from main parameters 970 # $data is name of global array for the parameters 971 global $data 972 973 set d [EllipsdData $datum] 974 set a [lindex $d 0] ; set f [lindex $d 1] 975 set ${data}(m_a) [expr 1.0*$a*[set ${data}(k0)]] 976 set es [expr $f*(2-$f)] 977 set ${data}(m_e) [expr sqrt($es)] 978 set e4 [expr $es*$es] ; set e6 [expr $e4*$es] 979 set ${data}(m_1) [expr $es* \ 980 (0.5+0.20833333333333333333*$es+0.08333333333333333333*$e4+ \ 981 0.03611111111111111111*$e6)] 982 set ${data}(m_2) [expr $es* \ 983 (0.14583333333333333333*$es+0.12083333333333333333*$e4+ \ 984 0.07039930555555555555*$e6)] 985 set ${data}(m_3) [expr $es* \ 986 (0.05833333333333333333*$e4+0.07232142857142857142*$e6)] 987 set ${data}(m_4) [expr 0.02653149801587301587*$es*$e6] 988 return 989} 990 991proc ProjMerc1Point {data latd longd datum} { 992 # compute planar Cartesian coordinates for given position 993 global $data 994 995 set prdatum [set ${data}(datum)] 996 if { $datum != $prdatum } { 997 foreach "latd longd" \ 998 [ToDatum $latd $longd $datum $prdatum] { break } 999 } 1000 set long0 [set ${data}(long0)] 1001 set longd [NormalLongCentred $longd $long0] 1002 set phi [expr $latd*0.01745329251994329576] 1003 set se [expr sin($phi)*[set ${data}(m_e)]] 1004 set mm [expr tan((1.5707963267948966+$phi)/2.0)* \ 1005 pow((1.0-$se)/(1.0+$se), [set ${data}(m_e)]/2.0)] 1006 set x [expr [set ${data}(m_a)]*($longd-$long0)*0.01745329251994329576] 1007 set y [expr [set ${data}(m_a)]*log($mm)] 1008 return [list $x $y] 1009} 1010 1011proc ProjMerc1Invert {data x y} { 1012 # return list with latitude and longitude for point at $x,$y in the terrain 1013 global $data 1014 1015 set ksi [expr 1.5707963267948966- \ 1016 2*atan(exp(-$y/[set ${data}(m_a)]))] 1017 set latd [expr ($ksi+[set ${data}(m_1)]*sin($ksi+$ksi)+ \ 1018 [set ${data}(m_2)]*sin(4*$ksi)+ \ 1019 [set ${data}(m_3)]*sin(6*$ksi)+ \ 1020 [set ${data}(m_4)]*sin(8*$ksi))* \ 1021 57.29577951308232087684] 1022 set longd [expr $x/[set ${data}(m_a)]*57.29577951308232087684+ \ 1023 [set ${data}(long0)]] 1024 return [list $latd [NormalLong $longd]] 1025} 1026 1027# two standard parallels (symmetrical with respect to the equator) 1028 1029proc ProjMerc2ComputeAux {data datum} { 1030 # compute auxiliary parameters from main parameters 1031 # $data is name of global array for the parameters 1032 global $data 1033 1034 set d [EllipsdData $datum] 1035 set a [lindex $d 0] ; set f [lindex $d 1] 1036 set es [expr $f*(2-$f)] 1037 set ${data}(m_e) [expr sqrt($es)] 1038 set phi1 [expr abs([set ${data}(lat1)])*0.01745329251994329576] 1039 set s [expr sin($phi1)] 1040 set k0 [expr cos($phi1)/sqrt(1-$es*$s*$s)] 1041 set ${data}(m_a) [expr 1.0*$a*$k0] 1042 set e4 [expr $es*$es] ; set e6 [expr $e4*$es] 1043 set ${data}(m_1) [expr $es* \ 1044 (0.5+0.20833333333333333333*$es+0.08333333333333333333*$e4+ \ 1045 0.03611111111111111111*$e6)] 1046 set ${data}(m_2) [expr $es* \ 1047 (0.14583333333333333333*$es+0.12083333333333333333*$e4+ \ 1048 0.07039930555555555555*$e6)] 1049 set ${data}(m_3) [expr $es* \ 1050 (0.05833333333333333333*$e4+0.07232142857142857142*$e6)] 1051 set ${data}(m_4) [expr 0.02653149801587301587*$es*$e6] 1052 return 1053} 1054 1055proc ProjMerc2Point {data latd longd datum} { 1056 # compute planar Cartesian coordinates for given position 1057 global $data 1058 1059 return [ProjMerc1Point $data $latd $longd $datum] 1060} 1061 1062proc ProjMerc2Invert {data x y} { 1063 # return list with latitude and longitude for point at $x,$y in the terrain 1064 global $data 1065 1066 return [ProjMerc1Invert $data $x $y] 1067} 1068 1069## Mercator projection: a special case of Lambert Conic Conformal 1070# with the equator as standard parallel 1071 1072# single standard parallel, spherical case 1073# adapted from libproj 4.6.1 1074 1075proc ProjSphMercComputeAux {data datum} { 1076 # compute auxiliary parameters from main parameters 1077 # $data is name of global array for the parameters 1078 global $data 1079 1080 set a [lindex [EllipsdData $datum] 0] 1081 set lat0 [set ${data}(lat0)] 1082 if { $lat0 > 1e-15 } { 1083 set k0 [expr cos($lat0*0.01745329251994329576)] 1084 } else { set k0 1 } 1085 set ${data}(m_a) [expr 1.0*$a*$k0] 1086 return 1087} 1088 1089proc ProjSphMercPoint {data latd longd datum} { 1090 # compute planar Cartesian coordinates for given position 1091 global $data 1092 1093 set prdatum [set ${data}(datum)] 1094 if { $datum != $prdatum } { 1095 foreach "latd longd" \ 1096 [ToDatum $latd $longd $datum $prdatum] { break } 1097 } 1098 set long0 [set ${data}(long0)] 1099 set longd [NormalLongCentred $longd $long0] 1100 set x [expr [set ${data}(m_a)]*($longd-$long0)*0.01745329251994329576] 1101 set y [expr [set ${data}(m_a)]* \ 1102 log(tan(0.78539816339744830962+$latd*0.00872664625997164788))] 1103 return [list $x $y] 1104} 1105 1106proc ProjSphMercInvert {data x y} { 1107 # return list with latitude and longitude for point at $x,$y in the terrain 1108 global $data 1109 1110 set latd [expr (1.5707963267948966-2*atan(exp(-$y/[set ${data}(m_a)]))) * \ 1111 57.29577951308232087684] 1112 set longd [expr $x/[set ${data}(m_a)]*57.29577951308232087684+ \ 1113 [set ${data}(long0)]] 1114 return [list $latd [NormalLong $longd]] 1115} 1116 1117## Cassini-Soldner projection 1118 1119proc ProjCSComputeAux {data datum} { 1120 # compute auxiliary parameters from main parameters 1121 # $data is name of global array for the parameters 1122 global $data 1123 1124 set d [EllipsdData $datum] 1125 set a [set ${data}(m_a) [lindex $d 0]] ; set f [lindex $d 1] 1126 set es [set ${data}(m_e) [expr $f*(2-$f)]] 1127 set phi0 [expr [set ${data}(lat0)]*0.01745329251994329576] 1128 set e4 [expr $es*$es] ; set e6 [expr $es*$e4] 1129 set m1 [set ${data}(m_1) [expr 1-0.25*$es-0.046875*$e4-0.01953125*$e6]] 1130 set m2 [set ${data}(m_2) [expr 0.375*$es+0.09375*$e4+0.0439453125*$e6]] 1131 set m3 [set ${data}(m_3) [expr 0.05859375*$e4+0.0439453125*$e6]] 1132 set m4 [set ${data}(m_4) [expr 0.01139322916666666667*$e6]] 1133 set ${data}(m_0) [expr $a*($m1*$phi0-$m2*sin($phi0+$phi0)+ \ 1134 $m3*sin(4*$phi0)-$m4*sin(6*$phi0))] 1135 set se [expr sqrt(1-$es)] 1136 set e1 [expr (1.0-$se)/(1+$se)] ; set e12 [expr $e1*$e1] 1137 set ${data}(m_5) [expr $e1*(1.5-0.84375*$e12)] 1138 set ${data}(m_6) [expr $e12*(1.3125-1.71875*$e12)] 1139 set ${data}(m_7) [expr 1.57291666666666666667*$e12*$e1] 1140 set ${data}(m_8) [expr 2.142578125*$e12*$e12] 1141 return 1142} 1143 1144proc ProjCSPoint {data latd longd datum} { 1145 # compute planar Cartesian coordinates for given position 1146 global $data 1147 1148 set prdatum [set ${data}(datum)] 1149 if { $datum != $prdatum } { 1150 foreach "latd longd" \ 1151 [ToDatum $latd $longd $datum $prdatum] { break } 1152 } 1153 set long0 [set ${data}(long0)] 1154 set longd [NormalLongCentred $longd $long0] 1155 set a [set ${data}(m_a)] ; set es [set ${data}(m_e)] 1156 set phi [expr $latd*0.01745329251994329576] 1157 set cphi [expr cos($phi)] ; set sphi [expr sin($phi)] 1158 set A [expr ($longd-$long0)*0.01745329251994329576*$cphi] 1159 set A2 [expr $A*$A] ; set A4 [expr $A2*$A2] 1160 set t [expr tan($phi)] ; set T [expr $t*$t] 1161 set C [expr $es*$cphi*$cphi/(1.0-$es)] 1162 set v [expr $a/sqrt(1-$es*$sphi*$sphi)] 1163 set M [expr $a*([set ${data}(m_1)]*$phi- \ 1164 [set ${data}(m_2)]*sin($phi+$phi)+ \ 1165 [set ${data}(m_3)]*sin(4*$phi)- \ 1166 [set ${data}(m_4)]*sin(6*$phi))] 1167 set x [expr $v*$A*(1-$T*$A2/6.0-(8-$T+8*$C)*$T*$A4/120.0)] 1168 set y [expr $M-[set ${data}(m_0)]+$v*$t*(0.5*$A2+(5-$T+6*$C)*$A4/24.0)] 1169 return [list $x $y] 1170} 1171 1172proc ProjCSInvert {data x y} { 1173 # return list with latitude and longitude for point at $x,$y in the terrain 1174 global $data 1175 1176 set a [set ${data}(m_a)] ; set es [set ${data}(m_e)] 1177 set miu1 [expr ([set ${data}(m_0)]+$y)/($a*[set ${data}(m_1)])] 1178 set phi1 [expr $miu1+[set ${data}(m_5)]*sin($miu1+$miu1)+ \ 1179 [set ${data}(m_6)]*sin(4*$miu1)+ \ 1180 [set ${data}(m_7)]*sin(6*$miu1)+ \ 1181 [set ${data}(m_8)]*sin(8*$miu1)] 1182 set sphi1 [expr sin($phi1)] 1183 set t [expr 1-$es*$sphi1*$sphi1] 1184 set v1 [expr $a/sqrt($t)] ; set rho1 [expr $v1*(1-$es)/$t] 1185 set t [expr tan($phi1)] ; set T1 [expr $t*$t] 1186 set D [expr $x/$v1] ; set D2 [expr $D*$D] 1187 set tt [expr (1+3*$T1)*$D2] 1188 set lat [expr ($phi1-$v1*$t/$rho1*$D2*(0.5-$tt/24.0))* \ 1189 57.29577951308232087684] 1190 set long [expr $D*(1-$T1*$D2/3.0+$tt*$D2/15.0)/cos($phi1)* \ 1191 57.29577951308232087684+[set ${data}(long0)]] 1192 return [list $lat [NormalLong $long]] 1193} 1194 1195## American polyconic projection 1196# adapted from PROJ4.0 1197 1198proc ProjAPOLYComputeAux {data datum} { 1199 # compute auxiliary parameters from main parameters 1200 # $data is name of global array for the parameters 1201 global $data 1202 1203 set d [EllipsdData $datum] 1204 set ${data}(m_a) [lindex $d 0] 1205 set f [lindex $d 1] 1206 set es [set ${data}(m_5) [expr $f*(2-$f)]] 1207 MeridDistParams $data $es 1208 set phi0 [expr [set ${data}(lat0)]*0.01745329251994329576] 1209 set ${data}(m_6) \ 1210 [MeridionalDist $data $phi0 [expr sin($phi0)] [expr cos($phi0)]] 1211 return 1212} 1213 1214proc ProjAPOLYPoint {data latd longd datum} { 1215 # compute planar Cartesian coordinates for given position 1216 global $data 1217 1218 set prdatum [set ${data}(datum)] 1219 if { $datum != $prdatum } { 1220 foreach "latd longd" \ 1221 [ToDatum $latd $longd $datum $prdatum] { break } 1222 } 1223 set a [set ${data}(m_a)] 1224 set phi [expr $latd*0.01745329251994329576] 1225 set lam [expr $longd*0.01745329251994329576] 1226 if { abs($phi) <= 1e-10 } { 1227 set x [expr $lam*$a] ; set y [expr -$a*[set ${data}(m_6)]] 1228 } else { 1229 set sp [expr sin($phi)] 1230 if { abs([set cp [expr cos($phi)]]) > 1e-10 } { 1231 set ms [expr $cp/sqrt(1-[set ${data}(m_5)]*$sp*$sp)/$sp] 1232 } else { set ms 0 } 1233 set t [expr $lam*$sp] 1234 set x [expr $a*$ms*sin($t)] 1235 set y [expr $a*(([MeridionalDist $data $phi $sp $cp]- \ 1236 [set ${data}(m_6)])+$ms*(1-cos($t)))] 1237 } 1238 return [list $x $y] 1239} 1240 1241proc ProjAPOLYInvert {data x y} { 1242 # return list with latitude and longitude for point at $x,$y in the terrain 1243 global $data MESS TXT 1244 1245 set a [set ${data}(m_a)] 1246 set x [expr 1.0*$x/$a] ; set y [expr 1.0*$y/$a] 1247 set y [expr $y+[set ${data}(m_6)]] 1248 if { abs($y) <= 1e-10 } { 1249 set phi 0 ; set lam $x 1250 } else { 1251 set r [expr $x*$x+$y*$y] ; set phi $y 1252 set es [set ${data}(m_5)] 1253 for { set i 0 } { $i < 20 } { incr i } { 1254 set sp [expr sin($phi)] ; set cp [expr cos($phi)] 1255 set s2ph [expr $sp*$cp] 1256 if { abs($cp) < 1e-12 } { 1257 GMMessage [format $MESS(badinvproj) $TXT(PRJAPOLY)] 1258 set i 0 1259 break 1260 } 1261 set mlp [expr sqrt(1-$es*$sp*$sp)] 1262 set c [expr $sp*$mlp/$cp] 1263 set ml [MeridionalDist $data $phi $sp $cp] 1264 set mlb [expr $ml*$ml+$r] 1265 set mlp [expr (1-$es)/($mlp*$mlp*$mlp)] 1266 set dphi [expr ($ml+$ml+$c*$mlb-2.0*$y*($c*$ml+1))/ \ 1267 ($es*$s2ph*($mlb-2*$y*$ml)/$c+ \ 1268 2*($y-$ml)*($c*$mlp-1.0/$s2ph)-$mlp-$mlp)] 1269 set phi [expr $phi+$dphi] 1270 if { abs($dphi) <= 1e-12 } { break } 1271 } 1272 if { $i == 0 } { 1273 GMMessage [format $MESS(badinvproj) $TXT(PRJAPOLY)] 1274 } 1275 set c [expr sin($phi)] 1276 set lam [expr [Aasin [expr $x*tan($phi)*sqrt(1-$es*$c*$c)]]/sin($phi)] 1277 } 1278 set longd [expr $lam*57.29577951308232087684] 1279 set latd [expr $phi*57.29577951308232087684] 1280 return [list $latd $longd] 1281} 1282 1283## Equidistant Cylindrical projection 1284# adapted from libproj 4.6.1 1285 1286proc ProjEqCylComputeAux {data datum} { 1287 # compute auxiliary parameters from main parameters 1288 # $data is name of global array for the parameters 1289 global $data 1290 1291 set a [set ${data}(m_a) [lindex [EllipsdData $datum] 0]] 1292 set ${data}(m_0) [expr cos([set ${data}(lat0)]*0.01745329251994329576)*$a] 1293 return 1294} 1295 1296proc ProjEqCylPoint {data latd longd datum} { 1297 # compute planar Cartesian coordinates for given position 1298 global $data 1299 1300 set prdatum [set ${data}(datum)] 1301 if { $datum != $prdatum } { 1302 foreach "latd longd" \ 1303 [ToDatum $latd $longd $datum $prdatum] { break } 1304 } 1305 set x [expr [set ${data}(m_0)]*$longd*0.01745329251994329576] 1306 set y [expr [set ${data}(m_a)]*($latd-[set ${data}(lat0)])* \ 1307 0.01745329251994329576] 1308 return [list $x $y] 1309} 1310 1311proc ProjEqCylInvert {data x y} { 1312 # return list with latitude and longitude for point at $x,$y in the terrain 1313 global $data 1314 1315 set latd [expr $y/[set ${data}(m_a)]*57.29577951308232087684+ \ 1316 [set ${data}(lat0)]] 1317 set longd [expr $x/[set ${data}(m_0)]*57.29577951308232087684] 1318 return [list $latd [NormalLong $longd]] 1319} 1320 1321## auxiliary procs for polyconic projection 1322# meridional distance for ellipsoid and inverse 1323 1324proc MeridDistParams {data es} { 1325 # compute parameters for meridional distance and inverse functions 1326 # saving them at indices m_0, ..., m_4 of array $data 1327 global $data 1328 1329 set t [expr $es*(0.046875+$es*(0.01953125+$es*0.01068115234375))] 1330 set ${data}(m_0) [expr 1-$es*(0.25+$t)] 1331 set ${data}(m_1) [expr $es*(0.75-$t)] 1332 set t [expr $es*$es] 1333 set ${data}(m_2) [expr $t*(0.46875-$es* \ 1334 (0.01302083333333333333+$es*0.00712076822916666666))] 1335 set t [expr $es*$t] 1336 set ${data}(m_3) \ 1337 [expr $t*(0.36458333333333333333-$es*0.00569661458333333333)] 1338 set ${data}(m_4) [expr $t*$es*0.3076171875] 1339 return 1340} 1341 1342proc MeridionalDist {data phi sphi cphi} { 1343 # compute meridional distance for given latitude 1344 # assume parameters for transformation are in array $data at indices 1345 # m_0, ..., m_4 1346 global $data 1347 1348 set cphi [expr $cphi*$sphi] 1349 set sphi [expr $sphi*$sphi] 1350 return [expr [set ${data}(m_0)]*$phi-$cphi* \ 1351 ([set ${data}(m_1)]+$sphi* \ 1352 ([set ${data}(m_2)]+$sphi* \ 1353 ([set ${data}(m_3)]+$sphi*[set ${data}(m_4)])))] 1354} 1355 1356proc LatFromMeridDist {data d} { 1357 # compute latitude from given meridional distance 1358 # assume parameters for transformation and e^2 are in array $data at 1359 # indices m_0, ..., m_4 and m_5 1360 global $data MESS 1361 1362 set es [set ${data}(m_5)] 1363 set k [expr 1.0/(1-$es)] 1364 set phi $d 1365 for { set i 0 } { $i < 10 } { incr i } { 1366 set s [expr sin($phi)] 1367 set t [expr 1.0-$es*$s*$s] 1368 set dphi [expr [MeridionalDist $data $phi $s [expr cos($phi)]]* \ 1369 ($t*sqrt($t))*$k] 1370 set phi [expr $phi-$dphi] 1371 if { abs($dphi) < 1e-11 } { return $phi } 1372 } 1373 GMMessage $MESS(badinvmdist) 1374 return $phi 1375} 1376 1377## Stereographic projection 1378# 1379 1380proc ProjStereogrComputeAux {data datum} { 1381 # compute auxiliary parameters from main parameters 1382 # $data is name of global array for the parameters 1383 global $data 1384 1385 set d [EllipsdData $datum] 1386 set ${data}(m_a) [lindex $d 0] ; set f [lindex $d 1] 1387 set e [set ${data}(m_e) [expr sqrt($f*(2-$f))]] 1388 set phi0 [expr [set ${data}(lat0)]*0.01745329251994329576] 1389 set k0 [set ${data}(k0)] 1390 set t [expr abs($phi0)] 1391 if { abs($t-1.5707963267948966) < 1e-10 } { 1392 if { $phi0 < 0 } { 1393 set mode s_pole 1394 } else { set mode n_pole } 1395 set ${data}(m_0) [expr 2*$k0/sqrt(pow(1+$e,1+$e)*pow(1-$e,1-$e))] 1396 } elseif { $t > 1e-10 } { 1397 set mode obliq 1398 set t [expr sin($phi0)] 1399 set X [expr 2*atan([SSFunc_ $phi0 $t $e])-1.5707963267948966] 1400 set t [expr $t*$e] 1401 set ${data}(m_0) [expr 2*$k0*cos($phi0)/sqrt(1.0-$t*$t)] 1402 set ${data}(m_1) [expr sin($X)] 1403 set ${data}(m_2) [expr cos($X)] 1404 } else { 1405 set mode equat 1406 set ${data}(m_0) [expr 2*$k0] 1407 set ${data}(m_1) 0 1408 set ${data}(m_2) 1.0 1409 } 1410 set ${data}(m_3) $mode 1411 return 1412} 1413 1414proc ProjStereogrPoint {data latd longd datum} { 1415 # compute planar Cartesian coordinates for given position 1416 global $data 1417 1418 set prdatum [set ${data}(datum)] 1419 if { $datum != $prdatum } { 1420 foreach "latd longd" \ 1421 [ToDatum $latd $longd $datum $prdatum] { break } 1422 } 1423 set long0 [set ${data}(long0)] 1424 set longd [NormalLongCentred $longd $long0] 1425 set e [set ${data}(m_e)] 1426 set akm1 [set ${data}(m_0)] 1427 set phi [expr $latd*0.01745329251994329576] 1428 set lam [expr ($longd-$long0)*0.01745329251994329576] 1429 set coslam [expr cos($lam)] ; set sinlam [expr sin($lam)] 1430 set sinphi [expr sin($phi)] 1431 switch [set ${data}(m_3)] { 1432 obliq - 1433 equat { 1434 set X [expr 2*atan([SSFunc_ $phi $sinphi $e])-1.5707963267948966] 1435 set sinX [expr sin($X)] ; set cosX [expr cos($X)] 1436 set sinX1 [set ${data}(m_1)] 1437 set cosX1 [set ${data}(m_2)] 1438 set A [expr $akm1/($cosX1*(1+$sinX1*$sinX+$cosX1*$cosX*$coslam))] 1439 set y [expr $A*($cosX1*$sinX-$sinX1*$cosX*$coslam)] 1440 set x [expr $A*$cosX] 1441 } 1442 n_pole { 1443 set x [expr $akm1*[ExpRedLat $phi $sinphi $e]] 1444 set y [expr -$x*$coslam] 1445 } 1446 s_pole { 1447 set x [expr $akm1*[ExpRedLat [expr -$phi] [expr -$sinphi] $e]] 1448 set y [expr $x*$coslam] 1449 } 1450 } 1451 set a [set ${data}(m_a)] 1452 set x [expr $x*$sinlam*$a] ; set y [expr $y*$a] 1453 return [list $x $y] 1454} 1455 1456proc ProjStereogrInvert {data x y} { 1457 # return list with latitude and longitude for point at $x,$y in the terrain 1458 global $data MESS TXT 1459 1460 set e [set ${data}(m_e)] ; set a [set ${data}(m_a)] 1461 set akm1 [set ${data}(m_0)] 1462 set long0 [set ${data}(long0)] ; set datum [set ${data}(datum)] 1463 set x [expr $x/$a] ; set y [expr $y/$a] 1464 if { abs($x)<1e-6 && abs($y)<1e-6 } { 1465 return [list [set ${data}(lat0)] $long0] 1466 } 1467 set rho [Hypot $x $y] 1468 set spole 0 1469 switch [set m [set ${data}(m_3)]] { 1470 obliq - 1471 equat { 1472 set sinX1 [set ${data}(m_1)] 1473 set cosX1 [set ${data}(m_2)] 1474 set tp [expr 2*atan2($rho*$cosX1,$akm1)] 1475 set cosphi [expr cos($tp)] ; set sinphi [expr sin($tp)] 1476 if { $rho == 0.0 } { 1477 set phi_l [Aasin [expr $cosphi*$sinX1]] 1478 } else { 1479 set phi_l [Aasin [expr $cosphi*$sinX1+($y*$sinphi*$cosX1/$rho)]] 1480 } 1481 set tp [expr tan(0.5*(1.5707963267948966+$phi_l))] 1482 set x [expr $x*$sinphi] 1483 set y [expr $rho*$cosX1*$cosphi-$y*$sinX1*$sinphi] 1484 set halfpi 1.5707963267948966 1485 set halfe [expr 0.5*$e] 1486 } 1487 n_pole { 1488 set y [expr -$y] 1489 set tp [expr -1.0*$rho/$akm1] 1490 set phi_l [expr 1.5707963267948966-2*atan($tp)] 1491 set halfpi -1.5707963267948966 1492 set halfe [expr -0.5*$e] 1493 } 1494 s_pole { 1495 set spole 1 1496 set tp [expr -1.0*$rho/$akm1] 1497 set phi_l [expr 1.5707963267948966-2*atan($tp)] 1498 set halfpi -1.5707963267948966 1499 set halfe [expr -0.5*$e] 1500 } 1501 } 1502 for { set i 8 } { [incr i -1] } { set phi_l $phi } { 1503 set sinphi [expr $e*sin($phi_l)] 1504 set phi [expr 2*atan($tp*pow((1+$sinphi)/(1-$sinphi),$halfe))-$halfpi] 1505 if { abs($phi_l-$phi) < 1e-10 } { 1506 if { $spole } { set phi [expr -$phi] } 1507 if { $x == 0 && $y == 0 } { 1508 set longd 0 1509 } else { set longd [expr atan2($x,$y)*57.29577951308232087684] } 1510 set latd [expr $phi*57.29577951308232087684] 1511 set longd [expr $longd+$long0] 1512 return [list $latd $longd] 1513 } 1514 } 1515 ## return 0 180 on error... 1516 GMMessage [format $MESS(badinvproj) $TXT(PRJStereogr)] 1517 return [list 0 180] 1518} 1519 1520# Schreiber double projection: used for the RDG (Netherlands) grid 1521# this is a special case of the oblique Stereographic projection 1522# 1523 1524proc ProjSchreiberComputeAux {data datum} { 1525 # compute auxiliary parameters from main parameters 1526 # $data is name of global array for the parameters 1527 global $data 1528 1529 set phi0 [set ${data}(m_0) 0.91029672689323934677] 1530 set lam0 [set ${data}(m_1) 0.09403203751960005358] 1531 set e [set ${data}(m_2) 0.08169683122252750299] 1532 set k0 0.9999079 1533 # 2*k*R 1534 set ${data}(m_3) [expr 12765289.142*$k0] 1535 set n 1.00047585668 ; set m 0.003773953832 1536 set tau0 [expr log([SSFunc_ $phi0 [expr sin($phi0)] $e])] 1537 set B0 [expr 2*atan(pow(2.7182818284590452354,$n*$tau0+$m))- \ 1538 1.5707963267948966] 1539 set ${data}(m_4) [expr sin($B0)] 1540 set ${data}(m_5) [expr cos($B0)] 1541 return 1542} 1543 1544proc ProjSchreiberPoint {data latd longd datum} { 1545 # compute planar Cartesian coordinates for given position 1546 global $data 1547 1548 set prdatum [set ${data}(datum)] 1549 if { $datum != $prdatum } { 1550 foreach "latd longd" \ 1551 [ToDatum $latd $longd $datum $prdatum] { break } 1552 } 1553 # just for the sake of those not knowing that this projection 1554 # was designed for use in the Netherlands... 1555 set long0 5.38763888888888888564 1556 set longd [NormalLongCentred $longd $long0] 1557 set n 1.00047585668 ; set m 0.003773953832 1558 set e [set ${data}(m_2)] 1559 set phi [expr $latd*0.01745329251994329576] 1560 set lam [expr $n*($longd*0.01745329251994329576-[set ${data}(m_1)])] 1561 set coslam [expr cos($lam)] ; set sinlam [expr sin($lam)] 1562 set sinphi [expr sin($phi)] 1563 set tau [expr log([SSFunc_ $phi $sinphi $e])] 1564 set B [expr 2*atan(pow(2.7182818284590452354,$n*$tau+$m))- \ 1565 1.5707963267948966] 1566 set sinB [expr sin($B)] ; set cosB [expr cos($B)] 1567 set sinB0 [set ${data}(m_4)] 1568 set cosB0 [set ${data}(m_5)] 1569 set tkR [set ${data}(m_3)] 1570 set A [expr $tkR/(1+$sinB*$sinB0+$cosB*$cosB0*$coslam)] 1571 set x [expr $A*$sinlam*$cosB] 1572 set y [expr $A*($sinB*$cosB0-$cosB*$sinB0*$coslam)] 1573 return [list $x $y] 1574} 1575 1576proc ProjSchreiberInvert {data x y} { 1577 # return list with latitude and longitude for point at $x,$y in the terrain 1578 global $data MESS TXT 1579 1580 set datum [set ${data}(datum)] 1581 set long0 [expr [set ${data}(m_1)]*57.29577951308232087684] 1582 if { abs($x)<1e-6 && abs($y)<1e-6 } { 1583 set lat0 [expr [set ${data}(m_0)]*57.29577951308232087684] 1584 return [list $lat0 $long0] 1585 } 1586 set n 1.00047585668 ; set m 0.003773953832 1587 set sinB0 [set ${data}(m_4)] 1588 set cosB0 [set ${data}(m_5)] 1589 set tkR [set ${data}(m_3)] 1590 set rho [Hypot $x $y] 1591 set Psi [expr 2*atan2($rho,$tkR)] 1592 set sinPsi [expr sin($Psi)] 1593 set B [Aasin [expr cos($Psi)*$sinB0+$y*$sinPsi*$cosB0/$rho]] 1594 set lam [Aasin [expr $x*$sinPsi/$rho/cos($B)]] 1595 set longd [expr $lam/$n*57.29577951308232087684+$long0] 1596 set tau [expr (log(tan(0.5*(1.5707963267948966+$B)))-$m)/$n] 1597 set etau [expr pow(2.7182818284590452354,$tau)] 1598 set phi [expr 2*atan($etau)-1.5707963267948966] 1599 set e [set ${data}(m_2)] ; set he [expr 0.5*$e] 1600 for { set i 0 } { $i < 50 } { incr i } { 1601 set esp [expr $e*sin($phi)] 1602 set phi_ [expr 2*atan($etau*pow((1+$esp)/(1-$esp),$he))- \ 1603 1.5707963267948966] 1604 if { abs($phi_-$phi) < 1e-10 } { 1605 set latd [expr $phi*57.29577951308232087684] 1606 return [list $latd $longd] 1607 } 1608 set phi $phi_ 1609 } 1610 ## return 0 180 on error... 1611 GMMessage [format $MESS(badinvproj) $TXT(PRJSchreiber)] 1612 return [list 0 180] 1613} 1614 1615## ancillary procs for Stereographic and Schreiber projections 1616 1617proc SSFunc_ {phit sinphi eccen} { 1618 1619 set sinphi [expr $sinphi*$eccen] 1620 return [expr tan(0.5*(1.5707963267948966+$phit))* \ 1621 pow((1-$sinphi)/(1+$sinphi),0.5*$eccen)] 1622} 1623 1624# Swiss Oblique Mercator projection 1625 1626proc ProjSOMComputeAux {data datum} { 1627 # compute auxiliary parameters from main parameters 1628 # $data is name of global array for the parameters 1629 global $data 1630 1631 set d [EllipsdData $datum] 1632 set a [lindex $d 0] ; set f [lindex $d 1] 1633 set es [expr $f*(2-$f)] 1634 set e [set ${data}(m_e) [expr sqrt($es)]] 1635 set he [set ${data}(m_0) [expr 0.5*$e]] 1636 set oes [expr 1-$es] 1637 set roes [set ${data}(m_6) [expr 1.0/$oes]] 1638 set phi0 [expr [set ${data}(lat0)]*0.01745329251994329576] 1639 set cp [expr cos($phi0)] ; set cp [expr $cp*$cp] 1640 set c [set ${data}(m_1) [expr sqrt(1+$es*$cp*$cp*$roes)]] 1641 set sp [expr sin($phi0)] 1642 set sp0 [set ${data}(m_3) [expr $sp/$c]] 1643 set p0 [Aasin $sp0] 1644 set ${data}(m_2) [expr cos($p0)] 1645 set sp [expr $sp*$e] 1646 set ${data}(m_4) [expr log(tan(0.78539816339744833+0.5*$p0))- \ 1647 $c*(log(tan(0.78539816339744833+0.5*$phi0))- \ 1648 $he*log((1.0+$sp)/(1.0-$sp)))] 1649 set ${data}(m_5) [expr $a*[set ${data}(k0)]*sqrt($oes)/(1.0-$sp*$sp)] 1650 return 1651} 1652 1653proc ProjSOMPoint {data latd longd datum} { 1654 # compute planar Cartesian coordinates for given position 1655 global $data 1656 1657 set prdatum [set ${data}(datum)] 1658 if { $datum != $prdatum } { 1659 foreach "latd longd" \ 1660 [ToDatum $latd $longd $datum $prdatum] { break } 1661 } 1662 # just for the sake of those not knowing that this projection 1663 # was designed for use in Switzerland... 1664 set long0 [set ${data}(long0)] 1665 set longd [NormalLongCentred $longd $long0] 1666 set phi [expr $latd*0.01745329251994329576] 1667 set sp [expr [set ${data}(m_e)]*sin($phi)] 1668 set c [set ${data}(m_1)] 1669 set phip [expr 2*atan(exp($c*(log(tan(0.78539816339744833+0.5*$phi))- \ 1670 [set ${data}(m_0)]*log((1.0+$sp)/(1.0-$sp)))+ \ 1671 [set ${data}(m_4)]))-1.5707963267948966] 1672 set lam [expr $c*($longd-$long0)*0.01745329251994329576] 1673 set cp [expr cos($phip)] 1674 set phipp [Aasin [expr [set ${data}(m_2)]*sin($phip)- \ 1675 [set ${data}(m_3)]*$cp*cos($lam)]] 1676 set akR [set ${data}(m_5)] 1677 set x [expr $akR*[Aasin [expr $cp*sin($lam)/cos($phipp)]]] 1678 set y [expr $akR*log(tan(0.78539816339744833+0.5*$phipp))] 1679 return [list $x $y] 1680} 1681 1682proc ProjSOMInvert {data x y} { 1683 # return list with latitude and longitude for point at $x,$y in the terrain 1684 global $data MESS TXT 1685 1686 set datum [set ${data}(datum)] 1687 set akR [set ${data}(m_5)] 1688 set phipp [expr 2*(atan(exp(1.0*$y/$akR))-0.78539816339744833)] 1689 set lampp [expr 1.0*$x/$akR] 1690 set cp [expr cos($phipp)] 1691 set phip [Aasin [expr [set ${data}(m_2)]*sin($phipp)+ \ 1692 [set ${data}(m_3)]*$cp*cos($lampp)]] 1693 set lamp [Aasin [expr $cp*sin($lampp)/cos($phip)]] 1694 set e [set ${data}(m_e)] 1695 set he [set ${data}(m_0)] ; set c [set ${data}(m_1)] 1696 set roes [set ${data}(m_6)] 1697 set con [expr ([set ${data}(m_4)]- \ 1698 log(tan(0.78539816339744833+0.5*$phip)))/$c] 1699 set i 10 1700 while { [incr i -1] } { 1701 set esp [expr $e*sin($phip)] 1702 set delp [expr ($con+log(tan(0.78539816339744833+0.5*$phip))- \ 1703 $he*log((1.0+$esp)/(1.0-$esp)))*(1.0-$esp*$esp)*cos($phip)* \ 1704 $roes] 1705 set phip [expr $phip-$delp] 1706 if { abs($delp) < 1e-10 } { break } 1707 } 1708 if { $i } { 1709 set latd [expr $phip*57.29577951308232087684] 1710 set longd [expr 57.29577951308232087684*$lamp/$c+ \ 1711 [set ${data}(long0)]] 1712 return [list $latd $longd] 1713 } 1714 ## return 0 180 on error... 1715 GMMessage [format $MESS(badinvproj) $TXT(PRJSOM)] 1716 return [list 0 180] 1717} 1718 1719# ancillary procedures for Swiss Oblique Mercator and others 1720 1721proc Aasin {v} { 1722 # arcsin function 1723 global MESS 1724 1725 set av [expr abs($v)] 1726 if { $av >= 1 } { 1727 if { $av > 1.00000000000001 } { 1728 GMMessage [format $MESS(badargtofunc) asin] 1729 } 1730 if { $v < 0 } { 1731 return -1.5707963267948966 1732 } 1733 return 1.5707963267948966 1734 } 1735 return [expr asin($v)] 1736} 1737 1738# ancillary procedures for LCC1, LCC2, Stereographic and Schreiber projections 1739 1740proc Hypot {x y} { 1741 # computation of sqrt($x*$x+$y*$y) avoiding overflows 1742 1743 if { $x < 0 } { 1744 set x [expr -$x] 1745 } elseif { $x == 0 } { 1746 if { $y < 0 } { return [expr -$y] } 1747 return $y 1748 } 1749 if { $y < 0 } { 1750 set y [expr -$y] 1751 } elseif { $y == 0 } { return $x } 1752 if { $x < $y } { 1753 set x [expr 1.0*$x/$y] 1754 return [expr $y*sqrt(1+$x*$x)] 1755 } 1756 set y [expr 1.0*$y/$x] 1757 return [expr $x*sqrt(1+$y*$y)] 1758} 1759 1760# 1761# SL contribution 1762# 1763# Implementation of the EOV (Hungarian National Projection) 1764# 1765# 2008-07-15 1766# 1767# Based on publications by Jozsef Varga Phd. 1768# http://www.agt.bme.hu/staff_h/varga/Osszes/Dok3uj.htm 1769# http://www.agt.bme.hu/staff_h/varga/gps/kezdoknek.html 1770# 1771 1772proc ProjEOVPoint {data latd longd datum} { 1773 # compute planar Cartesian coordinates for given position 1774 # although this projection computes x for northing and 1775 # y for easting, for compatibility with the mapping procedures 1776 # this procedure returns the usual (easting, northing) pair, 1777 # and the grid coordinates are computed by inversing it 1778 # $data is not used 1779 1780 # MF change: use the 3 parameter conversion to the HD72 if needed 1781 if { $datum != "Hungarian Datum 1972" } { 1782 foreach "latd longd" \ 1783 [ToDatum $latd $longd $datum "Hungarian Datum 1972"] { break } 1784 } 1785 # MF changes after this point: for efficiency, original code commented out 1786 # HD72 Phi, Lambda to radians 1787 set Phi [expr $latd*0.01745329251994329576] 1788 set Lambda [expr $longd*0.01745329251994329576] 1789 #--- 1790 1791 # HD72 Phi, Lambda to spherical phi, lambda 1792 ## GRS67 semi-major / semi-minor axis 1793 # set a 6378160.0 1794 # set b 6356774.516 1795 ## EOV parameters 1796 set k 1.003110007693 1797 set n 1.000719704936 1798 # set Lambda_0 [expr (19.0 + 02.0/60.0 + 54.8584/3600.0) / 360.0 * 2.0*$pi] 1799 # set Lambda_0 0.33246029532469185632 1800 ## GRS67 eccentricity 1801 # set epsilon [expr sqrt(($a*$a-$b*$b)/($a*$a))] 1802 set epsilon 0.0818205680555 1803 # set A [expr pow(tan($pi/4.0+$Phi/2.0), $n)] 1804 set A [expr pow(tan(0.785398163397448+0.5*$Phi), $n)] 1805 # set B [expr pow((1.0-$epsilon*sin($Phi))/(1.0+$epsilon*sin($Phi)), 1806 # $n*$epsilon/2.0)] 1807 set sinPhi [expr sin($Phi)] 1808 set B [expr pow((1.0-$epsilon*$sinPhi)/(1.0+$epsilon*$sinPhi), \ 1809 0.5*$n*$epsilon)] 1810 # set phi [expr 2.0 * atan($k*$A*$B) - $pi/2.0] 1811 set phi [expr 2.0 * atan($k*$A*$B) - 1.570796326794897] 1812 # set lambda [expr $n*($Lambda-$Lambda_0)] 1813 set lambda [expr $n*($Lambda-0.33246029532469185632)] 1814 1815 # Spherical phi, lambda to auxiliary phi, lambda 1816 ## EOV starting point 1817 # set phi_0 [expr (47.0 + 06.0/60.0 + 00.0000/3600.0) / 360.0 * 2.0*$pi] 1818 set phi_0 0.8220500776893292307101144 1819 # set phi_a [expr asin(sin($phi)*cos($phi_0) 1820 # - cos($phi)*sin($phi_0)*cos($lambda))] 1821 set cosphi [expr cos($phi)] 1822 set phi_a [expr asin(sin($phi)*cos($phi_0)- \ 1823 $cosphi*sin($phi_0)*cos($lambda))] 1824 # set lambda_a [expr asin((cos($phi)*sin($lambda))/(cos($phi_a)))] 1825 set lambda_a [expr asin(($cosphi*sin($lambda))/cos($phi_a))] 1826 1827 # Auxiliary phi, lambda to cartesian x, y (x == northing && y == easting) 1828 ## Radius of the "new" Gauss-sphere 1829 # set R 6379743.001 1830 ## Reduction (scale factor) 1831 # set m_0 0.99993 1832 set Rm_0 6379296.41898993 1833 # set x [expr $R * $m_0 * log(tan($pi/4+$phi_a/2))] 1834 # set y [expr $R * $m_0 * $lambda_a] 1835 set x [expr $Rm_0 * log(tan(0.78539816339744830961+0.5*$phi_a))] 1836 set y [expr $Rm_0 * $lambda_a] 1837 1838 # MF change: the false easting and false northing are added by the 1839 # grid procedures; return the "raw" easting and northing 1840 return [list $y $x] 1841} 1842 1843proc ProjEOVInvert {data x y} { 1844 # return list with latitude and longitude for point at $x,$y in the terrain 1845 # where $x is the easting and $y the northing (against the EOV convention) 1846 # with no false easting/northing added 1847 # $data is not used 1848 # return 0, 180 if computation of latitude diverges 1849 global MESS TXT 1850 1851 # Return to the EOV convention ($y == easting && $x == northing) 1852 set z $x ; set x $y ; set y $z 1853 1854 # MF changes after this point: for efficiency, original code commented out 1855 # Cartesian y, x to auxiliary phi, lambda 1856 # set R 6379743.001 1857 ## Reduction (scale factor) 1858 # set m_0 0.99993 1859 set Rm_0 6379296.41898993 1860 ## Pi :-P 1861 # set pi 3.1415926535897932384626434 1862 # set phi_a [expr 2*atan(exp($x/($R*$m_0)))-$pi/2] 1863 set phi_a [expr 2*atan(exp($x/$Rm_0))-1.57079632679489661923] 1864 # set lambda_a [expr $y/($R*$m_0)] 1865 set lambda_a [expr $y/($Rm_0)] 1866 1867 # Auxiliary phi, lambda to spherical phi, lambda 1868 ## EOV starting point 1869 # set phi_0 [expr (47.0 + 06.0/60.0 + 00.0000/3600.0) / 360.0 * 2.0*$pi] 1870 # set phi_0 0.82205007768932923029 1871 set cosphi_a [expr cos($phi_a)] 1872 # set phi [expr asin( sin($phi_a)*cos($phi_0) + 1873 # cos($phi_a)*sin($phi_0)*cos($lambda_a) )] 1874 set phi [expr asin( sin($phi_a)*0.680720868959 + \ 1875 $cosphi_a*0.732542898787*cos($lambda_a) )] 1876 # set lambda [expr asin( (cos($phi_a)*sin($lambda_a)) / cos($phi) )] 1877 set lambda [expr asin( ($cosphi_a*sin($lambda_a)) / cos($phi) )] 1878 1879 # Spherical phi, lambda to HD72 Phi, Lambda 1880 ## GRS67 semi-major / semi-minor axis 1881 # set a 6378160.0 1882 # set b 6356774.516 1883 ## EOV parameters 1884 set k 1.003110007693 1885 set n 1.000719704936 1886 # set Lambda_0 [expr (19.0 + 02.0/60.0 + 54.8584/3600.0) / 360.0 * 2.0*$pi] 1887 # set Lambda_0 0.33246029532469185632 1888 ## GRS67 eccentricity 1889 # set epsilon [expr sqrt(($a*$a-$b*$b)/($a*$a))] 1890 set epsilon 0.0818205680555 1891 ## Iteration to calculate Phi 1892 set Phi $phi 1893 set i 0 1894 # set A [expr tan($pi/4+$phi/2)] 1895 set A [expr tan(0.785398163397448+0.5*$phi)] 1896 set invn [expr 1/$n] 1897 while 1 { 1898 set B [expr $k *pow( (1-$epsilon*sin($Phi)) / (1+$epsilon*sin($Phi)), \ 1899 $n*$epsilon/2)] 1900 # set Phi_new [expr 2 * atan(pow($A/$B, 1/$n)) - $pi/2] 1901 set Phi_new [expr 2 * atan(pow(1.0*$A/$B, $invn)) - 1.570796326794897] 1902 if { abs($Phi-$Phi_new) <= 1.5E-9 } { 1903 break 1904 } 1905 if { [incr i] > 20} { 1906 GMMessage [format $MESS(badinvproj) $TXT(PRJEOV)] 1907 return [list 0 180] 1908 } 1909 set Phi $Phi_new 1910 } 1911 ## Calculate Lambda 1912 # set Lambda [expr $Lambda_0 + $lambda/$n] 1913 set Lambda [expr 0.33246029532469185632 + $lambda/$n] 1914 1915 # MF change: results will be for the HD72 1916 # Convert the result from radian to DD.DDDDDDDD... 1917 set la [expr $Phi*57.29577951308232087684] 1918 set lo [expr $Lambda*57.29577951308232087684] 1919 #--- 1920 return [list $la $lo] 1921} 1922 1923 1924 1925 1926