1# mapproj.tcl -- 2# 3# Package for map projections. 4# 5# Copyright (c) 2007 by Kevin B. Kenny. 6# 7# See the file "license.terms" for information on usage and redistribution 8# of this file, and for a DISCLAIMER OF ALL WARRANTIES. 9# 10# RCS: @(#) $Id: mapproj.tcl,v 1.1 2007/08/24 22:36:35 kennykb Exp $ 11#------------------------------------------------------------------------------ 12 13package require Tcl 8.4 14package require math::interpolate 1.0 15package require math::special 0.2.1 16 17package provide mapproj 1.0 18 19# ::mapproj -- 20# 21# Namespace holding the procedures and values. 22 23namespace eval ::mapproj { 24 25 # Cylindrical projections 26 27 namespace export toPlateCarree fromPlateCarree 28 namespace export toCassini fromCassini 29 namespace export toCylindricalEqualArea fromCylindricalEqualArea 30 namespace export toMercator fromMercator 31 32 # Named cylindric equal-area projections with specific 33 # standard parallels 34 35 namespace export toLambertCylindricalEqualArea \ 36 fromLambertCylindricalEqualArea 37 namespace export toBehrmann fromBehrmann 38 namespace export toTrystanEdwards fromTrystanEdwards 39 namespace export toHoboDyer fromHoboDyer 40 namespace export toGallPeters fromGallPeters 41 namespace export toBalthasart fromBalthasart 42 43 # Pseudocylindrical projections - equal area 44 45 namespace export toSinusoidal fromSinusoidal 46 namespace export toMillerCylindrical fromMillerCylindrical 47 namespace export toMollweide fromMollweide 48 namespace export toEckertIV fromEckertIV toEckertVI fromEckertVI 49 50 # Pseudocylindrical projections - compromise 51 52 namespace export toRobinson fromRobinson 53 54 # Azimuthal projections 55 56 namespace export toAzimuthalEquidistant fromAzimuthalEquidistant 57 namespace export toLambertAzimuthalEqualArea fromLambertAzimuthalEqualArea 58 namespace export toStereographic fromStereographic 59 namespace export toOrthographic fromOrthographic 60 namespace export toGnomonic fromGnomonic 61 62 # Pseudo-azimuthal projections 63 64 namespace export toHammer fromHammer 65 66 # Conic projections 67 68 namespace export toConicEquidistant fromConicEquidistant 69 namespace export toLambertConformalConic fromLambertConformalConic 70 namespace export toAlbersEqualAreaConic fromAlbersEqualAreaConic 71 72 # Miscellaneous projections 73 74 namespace export toPeirceQuincuncial fromPeirceQuincuncial 75 76 # Fundamental constants 77 78 variable pi [expr {acos(-1.0)}] 79 variable twopi [expr {2.0 * $pi}] 80 variable halfpi [expr {0.5 * $pi}] 81 variable quarterpi [expr {0.25 * $pi}] 82 variable threequarterpi [expr {0.75 * $pi}] 83 variable mquarterpi [expr {-0.25 * $pi}] 84 variable mthreequarterpi [expr {-0.75 * $pi}] 85 variable radian [expr {180. / $pi}] 86 variable degree [expr {$pi / 180.}] 87 variable sqrt2 [expr {sqrt(2.)}] 88 variable sqrt8 [expr {2. * $sqrt2}] 89 variable halfSqrt3 [expr {sqrt(3.) / 2.}] 90 variable halfSqrt2 [expr {sqrt(2.) / 2.}] 91 variable EckertIVK1 [expr {2.0 / sqrt($pi * (4.0 + $pi))}] 92 variable EckertIVK2 [expr {2.0 * sqrt($pi / (4.0 + $pi))}] 93 variable EckertVIK1 [expr {sqrt(2.0 + $pi)}] 94 variable PeirceQuincuncialScale 3.7081493546027438 ;# 2*K(1/2) 95 variable PeirceQuincuncialLimit 1.8540746773013719 ;# K(1/2) 96 97 # Table of parallel length and distance from equator for the 98 # Robinson projection 99 100 variable RobinsonLatitude { 101 -90.0 -85.0 -80.0 -75.0 -70.0 -65.0 102 -60.0 -55.0 -50.0 -45.0 -40.0 -35.0 103 -30.0 -25.0 -20.0 -15.0 -10.0 -5.0 104 0.0 5.0 10.0 15.0 20.0 25.0 30.0 105 35.0 40.0 45.0 50.0 55.0 60.0 106 65.0 70.0 75.0 80.0 85.0 90.0 107 } 108 variable RobinsonPLEN { 109 0.5322 0.5722 0.6213 0.6732 0.7186 0.7597 110 0.7986 0.8350 0.8679 0.8962 0.9216 0.9427 111 0.9600 0.9730 0.9822 0.9900 0.9954 0.9986 112 1.0000 0.9986 0.9954 0.9900 0.9822 0.9730 0.9600 113 0.9427 0.9216 0.8962 0.8679 0.8350 0.7986 114 0.7597 0.7186 0.6732 0.6213 0.5722 0.5322 115 } 116 variable RobinsonPDFE { 117 -1.0000 -0.9761 -0.9394 -0.8936 -0.8435 -0.7903 118 -0.7346 -0.6769 -0.6176 -0.5571 -0.4958 -0.4340 119 -0.3720 -0.3100 -0.2480 -0.1860 -0.1240 -0.0620 120 0.0000 0.0620 0.1240 0.1860 0.2480 0.3100 0.3720 121 0.4340 0.4958 0.5571 0.6176 0.6769 0.7346 122 0.7903 0.8435 0.8936 0.9394 0.9761 1.0000 123 } 124 125 # Interpolation tables for Robinson 126 127 variable RobinsonSplinePLEN \ 128 [math::interpolate::prepare-cubic-splines \ 129 $RobinsonLatitude $RobinsonPLEN] 130 variable RobinsonSplinePDFE \ 131 [math::interpolate::prepare-cubic-splines \ 132 $RobinsonLatitude $RobinsonPDFE] 133 variable RobinsonM [expr {0.5072 * $pi}] 134 135 namespace import ::math::special::cn 136 137} 138 139# ::mapproj::ellF - 140# 141# Computes the Legendre incomplete elliptic integral of the 142# first kind: 143# 144# F(phi, k) = \integral_0^phi dtheta/sqrt(1 - k**2 sin**2 theta) 145# 146# 147# Parameters: 148# phi -- Limit of integration; angle around the ellipse 149# k -- Eccentricity 150# 151# Results: 152# Returns F(phi, k) 153# 154# Notes: 155# We compute this integral in terms of the Carlson elliptic integral 156# ellRF(x, y, z). 157 158proc ::mapproj::ellF {phi k} { 159 return [ellFaux [expr {cos($phi)}] [expr {sin($phi)}] $k] 160} 161 162# ::mapproj::ellFaux - 163# 164# Computes the Legendre incomplete elliptic integral of the 165# first kind when circular functions of the 'phi' argument 166# are already available. 167# 168# Parameters: 169# cos_phi - Cosine of the argument 170# sin_phi - Sine of the argument 171# k - Parameter 172# 173# Results: 174# Returns F(atan(sin_phi/cos_phi), k) 175 176proc ::mapproj::ellFaux {cos_phi sin_phi k} { 177 set rf [ellRF [expr {$cos_phi * $cos_phi}] \ 178 [expr {1.0 - $k * $k * $sin_phi * $sin_phi}] \ 179 1.0] 180 return [expr {$sin_phi * $rf}] 181} 182 183# ::mapproj::ellRF -- 184# 185# Computes the Carlson incomplete elliptic integral of the 186# first kind: 187# 188# RF(x, y, z) = 1/2 * integral_0^inf dt/sqrt((t+x)*(t+y)*(t+z)) 189# 190# Parameters: 191# x, y, z -- Interchangeable parameters of the integral 192# 193# Results: 194# Returns the value of RF 195 196proc ::mapproj::ellRF {x y z} { 197 if {$x < 0.0 || $y < 0.0 || $z < 0.0} { 198 return -code error "Negative argument to Carlson's ellRF" \ 199 -errorCode "ellRF negArgument" 200 } 201 set delx 1.0; set dely 1.0; set delz 1.0 202 while {abs($delx) > 0.0025 || abs($dely) > 0.0025 || abs($delz) > 0.0025} { 203 set sx [expr {sqrt($x)}] 204 set sy [expr {sqrt($y)}] 205 set sz [expr {sqrt($z)}] 206 set len [expr {$sx * ($sy + $sz) + $sy * $sz}] 207 set x [expr {0.25 * ($x + $len)}] 208 set y [expr {0.25 * ($y + $len)}] 209 set z [expr {0.25 * ($z + $len)}] 210 set mean [expr {($x + $y + $z) / 3.0}] 211 set delx [expr {($mean - $x) / $mean}] 212 set dely [expr {($mean - $y) / $mean}] 213 set delz [expr {($mean - $z) / $mean}] 214 } 215 set e2 [expr {$delx * $dely - $delz * $delz}] 216 set e3 [expr {$delx * $dely * $delz}] 217 return [expr {(1.0 + ($e2 / 24.0 - 0.1 - 3.0 * $e3 / 44.0) * $e2 218 + $e3 / 14.) / sqrt($mean)}] 219} 220 221# ::mapproj::toPlateCarree -- 222# 223# Project a latitude and longitude onto the plate carr�e. 224# 225# Parameters: 226# phi_0 -- Latitude of the center of the sheet in degrees 227# lambda_0 -- Longitude of the center of sheet in degrees 228# lambda -- Longitude of the point to be projected in degrees 229# phi -- Latitude of the point to be projected in degrees 230# 231# Results: 232# Returns x and y co-ordinates. Units of x and y are Earth radii; 233# scale is true at the Equator. 234 235proc ::mapproj::toPlateCarree {lambda_0 phi_0 lambda phi} { 236 variable degree 237 set lambda [expr {$lambda - $lambda_0 + 180.}] 238 if {$lambda < 0.0 || $lambda > 360.0} { 239 set lambda [expr {$lambda - 360. * floor($lambda / 360.)}] 240 } 241 set lambda [expr {$lambda - 180.}] 242 set x [expr {$lambda * $degree}] 243 set y [expr {($phi - $phi_0) * $degree}] 244 return [list $x $y] 245} 246 247# ::mapproj::fromPlateCarree -- 248# 249# Solve a plate carr�e projection for the 250# latitude and longitude represented by a point on the map. 251# 252# Parameters: 253# phi_0 -- Latitude of the center of the sheet in degrees 254# lambda_0 -- Longitude of the center of projection 255# x,y -- normalized x and y co-ordinates of a point on the map 256# 257# Results: 258# Returns a list consisting of the longitude and latitude in degrees. 259 260proc ::mapproj::fromPlateCarree {phi_0 lambda_0 x y} { 261 variable radian 262 set lambda [expr {$lambda_0 + $x * $radian + 180.}] 263 if {$lambda < 0.0 || $lambda > 360.0} { 264 set lambda [expr {$lambda - 360. * floor($lambda / 360.)}] 265 } 266 set lambda [expr {$lambda - 180.}] 267 set phi [expr {$y * $radian + $phi_0}] 268 return [list $lambda $phi] 269} 270 271# mapproj::toCylindricalEqualArea -- 272# 273# Project a latitude and longitude into cylindrical equal-area 274# co-ordinates. 275# 276# Parameters: 277# phi_1 -- Standard latitude in degrees 278# phi_0 -- Latitude of the center of the sheet in degrees 279# lambda_0 -- Longitude of the center of projection in degrees 280# lambda -- Longitude of the point to be projected in degrees 281# phi -- Latitude of the point to be projected in degrees 282# 283# Results: 284# Returns x and y co-ordinates. Units of x and y are Earth radii; 285# scale is true at the reference latitude 286 287proc ::mapproj::toCylindricalEqualArea {phi_1 lambda_0 phi_0 lambda phi} { 288 variable degree 289 set cos_phi_s [expr {cos($phi_1 * $degree)}] 290 set lambda [expr {$lambda - $lambda_0 + 180.}] 291 if {$lambda < 0.0 || $lambda > 360.0} { 292 set lambda [expr {$lambda - 360. * floor($lambda / 360.)}] 293 } 294 set lambda [expr {$lambda - 180.}] 295 set x [expr {$lambda * $degree * $cos_phi_s}] 296 set y0 [expr {sin($phi_0 * $degree) / $cos_phi_s}] 297 set y [expr {sin($phi * $degree) / $cos_phi_s}] 298 return [list $x [expr {$y - $y0}]] 299} 300 301# ::mapproj::fromCylindricalEqualArea -- 302# 303# Solve a cylindrical equal area map projection for the 304# latitude and longitude represented by a point on the map. 305# 306# Parameters: 307# phi_1 -- Standard latitude in degrees 308# phi_0 -- Latitude of the center of the sheet in degrees 309# lambda_0 -- Longitude of the center of sheet in degrees 310# x,y -- normalized x and y co-ordinates of a point on the map 311# 312# Results: 313# Returns a list consisting of the longitude and latitude in degrees. 314 315proc ::mapproj::fromCylindricalEqualArea {phi_1 lambda_0 phi_0 x y} { 316 variable degree 317 variable radian 318 set cos_phi_s [expr {cos($phi_1 * $degree)}] 319 set lambda [expr {$lambda_0 + $x / $cos_phi_s * $radian + 180.}] 320 if {$lambda < 0.0 || $lambda > 360.0} { 321 set lambda [expr {$lambda - 360. * floor($lambda / 360.)}] 322 } 323 set lambda [expr {$lambda - 180.}] 324 set y0 [expr {sin($phi_0 * $degree) / $cos_phi_s}] 325 set phi [expr {asin(($y + $y0) * $cos_phi_s) * $radian}] 326 return [list $lambda $phi] 327} 328 329# ::mapproj::toMercator -- 330# 331# Project a latitude and longitude into the Mercator projection 332# co-ordinates. 333# 334# Parameters: 335# lambda_0 -- Longitude of the center of projection in degrees 336# phi_0 -- Latitude of the center of sheet in degrees 337# lambda -- Longitude of the point to be projected in degrees 338# phi -- Latitude of the point to be projected in degrees 339# 340# Results: 341# Returns x and y co-ordinates in Earth radii. Scale is true 342# at the Equator and increases without bounds toward the Poles. 343 344proc ::mapproj::toMercator {lambda_0 phi_0 lambda phi} { 345 variable trace; if {[info exists trace]} { puts [info level 0] } 346 variable degree 347 variable quarterpi 348 set lambda [expr {$lambda - $lambda_0 + 180.}] 349 if {$lambda < 0.0 || $lambda > 360.0} { 350 set lambda [expr {$lambda - 360. * floor($lambda / 360.)}] 351 } 352 set lambda [expr {$lambda - 180.}] 353 set x [expr {$lambda * $degree}] 354 set y [expr {log(tan($quarterpi + 0.5 * $phi * $degree))}] 355 set y0 [expr {log(tan($quarterpi + 0.5 * $phi_0 * $degree))}] 356 if {[info exists trace]} { puts "[info level 0] -> $x [expr {$y - $y0}]" } 357 return [list $x [expr {$y - $y0}]] 358} 359 360# ::mapproj::fromMercator -- 361# 362# Converts Mercator map co-ordinates to latitude and longitude. 363# 364# Parameters: 365# lambda_0 -- Longitude of the center of projection 366# phi_0 -- Latitude of the center of sheet in degrees 367# x,y -- normalized x and y co-ordinates of a point on the map 368# 369# Results: 370# Returns a list consisting of the longitude and latitude in degrees. 371 372proc ::mapproj::fromMercator {lambda_0 phi_0 x y} { 373 variable quarterpi 374 variable degree 375 variable radian 376 set y0 [expr {log(tan($quarterpi + 0.5 * $phi_0 * $degree))}] 377 set lambda [expr {$lambda_0 + $x * $radian + 180.}] 378 if {$lambda < 0.0 || $lambda > 360.0} { 379 set lambda [expr {$lambda - 360. * floor($lambda / 360.)}] 380 } 381 set lambda [expr {$lambda - 180.}] 382 set phi [expr {$radian * 2.0 * atan(exp($y + $y0)) - 90.0}] 383 return [list $lambda $phi] 384} 385 386# ::mapproj::toMillerCylindrical -- 387# 388# Project a latitude and longitude into the Miller Cylindrical projection 389# co-ordinates. 390# 391# Parameters: 392# lambda_0 -- Longitude of the center of projection in degrees 393# lambda -- Longitude of the point to be projected in degrees 394# phi -- Latitude of the point to be projected in degrees 395# 396# Results: 397# Returns x and y co-ordinates. The x co-ordinate ranges from 398# -$pi to pi. 399 400proc ::mapproj::toMillerCylindrical {lambda_0 lambda phi} { 401 foreach {x y} [toMercator $lambda_0 0.0 \ 402 $lambda [expr {0.8 * $phi}]] break 403 set y [expr {1.25 * $y}] 404 return [list $x $y] 405} 406 407# ::mapproj::fromMillerCylindrical -- 408# 409# Converts Miller Cylindrical projected map co-ordinates 410# to latitude and longitude. 411# 412# Parameters: 413# lambda_0 -- Longitude of the center of projection 414# x,y -- normalized x and y co-ordinates of a point on the map 415# 416# Results: 417# Returns a list consisting of the longitude and latitude in degrees. 418 419proc ::mapproj::fromMillerCylindrical {lambda_0 x y} { 420 foreach {lambda phi} [fromMercator $lambda_0 0.0 \ 421 $x [expr {0.8 * $y}]] break 422 return [list $lambda [expr {1.25 * $phi}]] 423} 424 425# ::mapproj::toSinusoidal -- 426# 427# Project a latitude and longitude into the sinusoidal 428# (Sanson-Flamsteed) projection. 429# 430# Parameters: 431# lambda_0 -- Longitude of the center of projection in degrees 432# phi_0 -- Latitude of the center of the sheet, in degrees 433# lambda -- Longitude of the point to be projected in degrees 434# phi -- Latitude of the point to be projected in degrees 435# 436# Results: 437# Returns x and y co-ordinates, in Earth radii. 438# Scale is true along the Equator and central meridian. 439 440proc ::mapproj::toSinusoidal {lambda_0 phi_0 lambda phi} { 441 variable degree 442 variable quarterpi 443 set lambda [expr {$lambda - $lambda_0 + 180.}] 444 if {$lambda < 0.0 || $lambda > 360.0} { 445 set lambda [expr {$lambda - 360. * floor($lambda / 360.)}] 446 } 447 set lambda [expr {($lambda - 180.) * $degree}] 448 set phi [expr {$phi * $degree}] 449 set x [expr {$lambda * cos($phi)}] 450 set phi [expr {$phi - $phi_0 * $degree}] 451 return [list $x $phi] 452} 453 454# ::mapproj::fromSinusoidal -- 455# 456# Converts sinusoidal (Sanson-Flamsteed) map co-ordinates 457# to latitude and longitude. 458# 459# Parameters: 460# lambda_0 -- Longitude of the center of projection 461# phi_0 -- Latitude of the center of the sheet, in degrees 462# x,y -- normalized x and y co-ordinates of a point on the map 463# 464# Results: 465# Returns a list consisting of the longitude and latitude in degrees. 466 467proc ::mapproj::fromSinusoidal {lambda_0 phi_0 x y} { 468 variable degree 469 variable radian 470 set y [expr {$y + $phi_0 * $degree}] 471 set phi [expr {$y * $radian}] 472 set lambda [expr {180. + $lambda_0 + $radian * $x / cos($y)}] 473 if {$lambda < 0.0 || $lambda > 360.0} { 474 set lambda [expr {$lambda - 360. * floor($lambda / 360.)}] 475 } 476 set lambda [expr {$lambda - 180.}] 477 return [list $lambda $phi] 478} 479 480# ::mapproj::toMollweide -- 481# 482# Project a latitude and longitude into the Mollweide projection 483# co-ordinates. 484# 485# Parameters: 486# lambda_0 -- Longitude of the center of projection in degrees 487# lambda -- Longitude of the point to be projected in degrees 488# phi -- Latitude of the point to be projected in degrees 489# 490# Results: 491# Returns x and y co-ordinates, in Earth radii. 492# Scale is true along the 40 deg 44 min parallels 493 494proc ::mapproj::toMollweide {lambda_0 lambda phi} { 495 variable degree 496 variable pi 497 variable halfpi 498 variable sqrt2 499 variable sqrt8 500 set lambda [expr {$lambda - $lambda_0 + 180.}] 501 if {$lambda < 0.0 || $lambda > 360.0} { 502 set lambda [expr {$lambda - 360. * floor($lambda / 360.)}] 503 } 504 set lambda [expr {($lambda - 180.) * $degree}] 505 set phi [expr {$phi * $degree}] 506 set theta [expr {2.0 * asin(2.0 * $phi / $pi)}] 507 set diff 1.0 508 set pisinphi [expr {$pi * sin($phi)}] 509 while {abs($diff) >= 1.0e-4} { 510 set diff [expr {($theta + sin($theta) - $pisinphi) 511 / (1.0 + cos($theta))}] 512 set theta [expr {$theta - $diff}] 513 } 514 set theta [expr {0.5 * $theta}] 515 set x [expr {$sqrt8 * $lambda * cos($theta) / $pi}] 516 set y [expr {$sqrt2 * sin($theta)}] 517 return [list $x $y] 518} 519 520# ::mapproj::fromMollweide -- 521# 522# Converts Mollweide projected map co-ordinates to latitude 523# and longitude. 524# 525# Parameters: 526# lambda_0 -- Longitude of the center of projection 527# x,y -- normalized x and y co-ordinates of a point on the map 528# 529# Results: 530# Returns a list consisting of the longitude and latitude in degrees. 531 532proc ::mapproj::fromMollweide {lambda_0 x y} { 533 variable pi 534 variable radian 535 variable degree 536 variable halfpi 537 variable sqrt2 538 variable sqrt8 539 set theta [expr {asin($y / $sqrt2)}] 540 set lambda [expr {$lambda_0 + $radian * $pi * $x / 541 ($sqrt8 * cos($theta)) + 180.}] 542 set phi [expr {asin((2.0 * $theta + sin(2.0 * $theta)) / $pi)}] 543 if {$lambda < 0.0 || $lambda > 360.0} { 544 set lambda [expr {$lambda - 360. * floor($lambda / 360.)}] 545 } 546 set lambda [expr {$lambda - 180.}] 547 return [list $lambda [expr {$phi * $radian}]] 548} 549 550# ::mapproj::toEckertIV -- 551# 552# Project a latitude and longitude into the Eckert IV projection 553# co-ordinates. 554# 555# Parameters: 556# lambda_0 -- Longitude of the center of projection in degrees 557# lambda -- Longitude of the point to be projected in degrees 558# phi -- Latitude of the point to be projected in degrees 559# 560# Results: 561# Returns x and y co-ordinates. Scale is true along the 40 deg 30 min 562# parallels. 563 564proc ::mapproj::toEckertIV {lambda_0 lambda phi} { 565 variable degree 566 variable pi 567 variable halfpi 568 variable EckertIVK1 569 variable EckertIVK2 570 set lambda [expr {$lambda - $lambda_0 + 180.}] 571 if {$lambda < 0.0 || $lambda > 360.0} { 572 set lambda [expr {$lambda - 360. * floor($lambda / 360.)}] 573 } 574 set lambda [expr {($lambda - 180.) * $degree}] 575 set phi [expr {$phi * $degree}] 576 set theta [expr {$phi / 2}] 577 set diff 1.0 578 set A [expr {(2.0 + $halfpi) * sin($phi)}] 579 while {abs($diff) >= 1.0e-4} { 580 set costheta [expr {cos($theta)}] 581 set sintheta [expr {sin($theta)}] 582 set diff \ 583 [expr {($theta + $sintheta * $costheta + 2.0 * sin($theta) - $A) 584 / (2.0 * $costheta * (1.0 + $costheta))}] 585 set theta [expr {$theta - $diff}] 586 } 587 set x [expr {$EckertIVK1 * $lambda * (1.0 + cos($theta))}] 588 set y [expr {$EckertIVK2 * sin($theta)}] 589 return [list $x $y] 590} 591 592# ::mapproj::fromEckertIV -- 593# 594# Converts Eckert IV projected map co-ordinates to latitude 595# and longitude. 596# 597# Parameters: 598# lambda_0 -- Longitude of the center of projection 599# x,y -- normalized x and y co-ordinates of a point on the map 600# 601# Results: 602# Returns a list consisting of the longitude and latitude in degrees. 603 604proc ::mapproj::fromEckertIV {lambda_0 x y} { 605 variable pi 606 variable radian 607 variable degree 608 variable halfpi 609 variable sqrt2 610 variable EckertIVK1 611 variable EckertIVK2 612 set sintheta [expr {$y / $EckertIVK2}] 613 set costheta [expr {sqrt(1.0 - $sintheta * $sintheta)}] 614 set theta [expr {atan2($sintheta, $costheta)}] 615 set phi [expr {asin(($theta + $sintheta*$costheta + 2.*$sintheta) 616 / (2. + $halfpi))}] 617 set lambda [expr {180.0 + $lambda_0 618 + $radian / $EckertIVK1 * $x / (1.0 + $costheta)}] 619 if {$lambda < 0.0 || $lambda > 360.0} { 620 set lambda [expr {$lambda - 360. * floor($lambda / 360.)}] 621 } 622 set lambda [expr {$lambda - 180.}] 623 return [list $lambda [expr {$phi * $radian}]] 624} 625 626# ::mapproj::toEckertVI -- 627# 628# Project a latitude and longitude into the Eckert IV projection 629# co-ordinates. 630# 631# Parameters: 632# lambda_0 -- Longitude of the center of projection in degrees 633# lambda -- Longitude of the point to be projected in degrees 634# phi -- Latitude of the point to be projected in degrees 635# 636# Results: 637# Returns x and y co-ordinates. Scale is true along the 40 deg 30 min 638# parallels. 639 640proc ::mapproj::toEckertVI {lambda_0 lambda phi} { 641 variable degree 642 variable pi 643 variable halfpi 644 variable EckertVIK1 645 set lambda [expr {$lambda - $lambda_0 + 180.}] 646 if {$lambda < 0.0 || $lambda > 360.0} { 647 set lambda [expr {$lambda - 360. * floor($lambda / 360.)}] 648 } 649 set lambda [expr {($lambda - 180.) * $degree}] 650 set phi [expr {$phi * $degree}] 651 set theta [expr {$phi / 2}] 652 set diff 1.0 653 set A [expr {(1.0 + $halfpi) * sin($phi)}] 654 while {abs($diff) >= 1.0e-4} { 655 set costheta [expr {cos($theta)}] 656 set sintheta [expr {sin($theta)}] 657 set diff \ 658 [expr {($theta + $sintheta - $A) 659 / (1.0 + $costheta)}] 660 set theta [expr {$theta - $diff}] 661 } 662 set x [expr {$lambda * (1.0 + cos($theta)) / $EckertVIK1}] 663 set y [expr {2.0 * $theta / $EckertVIK1}] 664 return [list $x $y] 665} 666 667# ::mapproj::fromEckertVI -- 668# 669# Converts Eckert IV projected map co-ordinates to latitude 670# and longitude. 671# 672# Parameters: 673# lambda_0 -- Longitude of the center of projection 674# x,y -- normalized x and y co-ordinates of a point on the map 675# 676# Results: 677# Returns a list consisting of the longitude and latitude in degrees. 678 679proc ::mapproj::fromEckertVI {lambda_0 x y} { 680 variable pi 681 variable radian 682 variable degree 683 variable halfpi 684 variable sqrt2 685 variable EckertVIK1 686 puts [info level 0] 687 set theta [expr {0.5 * $EckertVIK1 * $y}] 688 puts [list theta = $theta] 689 set phi [expr {asin(($theta + sin($theta)) / (1.0 + $halfpi))}] 690 puts [list phi = $phi] 691 set lambda [expr {180.0 + $lambda_0 + $radian * $EckertVIK1 * $x 692 / (1 + cos($theta))}] 693 if {$lambda < 0.0 || $lambda > 360.0} { 694 set lambda [expr {$lambda - 360. * floor($lambda / 360.)}] 695 } 696 set lambda [expr {$lambda - 180.}] 697 return [list $lambda [expr {$phi * $radian}]] 698} 699 700# ::mapproj::toRobinson -- 701# 702# Project a latitude and longitude into the Robinson projection. 703# 704# Parameters: 705# lambda_0 -- Longitude of the center of projection in degrees 706# lambda -- Longitude of the point to be projected in degrees 707# phi -- Latitude of the point to be projected in degrees 708# 709# Results: 710# Returns x and y co-ordinates, in Earth radii. 711# Scale is true along the Equator. 712 713proc ::mapproj::toRobinson {lambda_0 lambda phi} { 714 variable RobinsonLatitude 715 variable RobinsonSplinePLEN 716 variable RobinsonSplinePDFE 717 variable RobinsonM 718 variable pi 719 variable degree 720 set lambda [expr {$lambda - $lambda_0 + 180.}] 721 if {$lambda < 0.0 || $lambda > 360.0} { 722 set lambda [expr {$lambda - 360. * floor($lambda / 360.)}] 723 } 724 set lambda [expr {$lambda - 180.}] 725 set y [math::interpolate::interp-cubic-splines $RobinsonSplinePDFE $phi] 726 set y [expr {$RobinsonM * $y}] 727 set s [math::interpolate::interp-cubic-splines $RobinsonSplinePLEN $phi] 728 set x [expr {$degree * $s * $lambda}] 729 return [list $x $y] 730} 731 732# ::mapproj::fromRobinson -- 733# 734# Solve the Robinson projection for the 735# latitude and longitude represented by a point on the map. 736# 737# Parameters: 738# lambda_0 -- Longitude of the center of projection 739# x,y -- normalized x and y co-ordinates of a point on the map 740# 741# Results: 742# Returns a list consisting of the longitude and latitude in degrees. 743 744proc ::mapproj::fromRobinson {lambda_0 x y} { 745 variable RobinsonLatitude 746 variable RobinsonPDFE 747 variable RobinsonSplinePLEN 748 variable RobinsonSplinePDFE 749 variable RobinsonM 750 variable radian 751 752 # We know that Robinson latitudes are equally spaced from [-90..90] 753 # at 5-degree intervals. Find the values for RobinsonPDFE that 754 # bracket the y co-ordinate. 755 756 set y [expr {$y / $RobinsonM}] 757 set l 0 758 set u [expr {[llength $RobinsonPDFE] - 1}] 759 while {$l < $u} { 760 set m [expr {($l + $u + 1) / 2}] 761 if {$y >= [lindex $RobinsonPDFE $m]} { 762 set l $m 763 } else { 764 set u [expr {$m - 1}] 765 } 766 } 767 set u [lindex $RobinsonLatitude [expr {$l+1}]] 768 set l [lindex $RobinsonLatitude $l] 769 for {set i 0} {$i < 12} {incr i} { 770 set m [expr {0.5 * ($u + $l)}] 771 set ystar [math::interpolate::interp-cubic-splines \ 772 $RobinsonSplinePDFE $m] 773 if {$ystar < $y} { 774 set l $m 775 } else { 776 set u $m 777 } 778 } 779 puts "latitude $m" 780 set s [math::interpolate::interp-cubic-splines $RobinsonSplinePLEN $m] 781 puts "parallel length $s" 782 set lambda [expr {180.0 + $lambda_0 + $radian * $x / $s}] 783 if {$lambda < 0.0 || $lambda > 360.0} { 784 set lambda [expr {$lambda - 360. * floor($lambda / 360.)}] 785 } 786 set lambda [expr {$lambda - 180.}] 787 return [list $lambda $m] 788} 789 790# ::mapproj::toCassini -- 791# 792# Project a latitude and longitude into the Cassini projection. 793# 794# Parameters: 795# lambda_0 -- Longitude of the center of projection in degrees 796# phi_0 -- Latitude of the center of the sheet 797# lambda -- Longitude of the point to be projected in degrees 798# phi -- Latitude of the point to be projected in degrees 799# 800# Results: 801# Returns x and y co-ordinates, in Earth radii. 802# Scale is true along the central meridian. 803 804proc ::mapproj::toCassini {lambda_0 phi_0 lambda phi} { 805 variable degree 806 variable pi 807 variable twopi 808 variable quarterpi 809 set lambda [expr {$lambda - $lambda_0 + 180.}] 810 if {$lambda < 0.0 || $lambda > 360.0} { 811 set lambda [expr {$lambda - 360. * floor($lambda / 360.)}] 812 } 813 set lambda [expr {($lambda - 180.) * $degree}] 814 set phi [expr {$phi * $degree}] 815 set x [expr {asin(cos($phi) * sin($lambda))}] 816 set y [expr {atan2(tan($phi), cos($lambda)) - $degree * $phi_0}] 817 if {$y < -$pi} { 818 set y [expr {$y + $twopi}] 819 } elseif {$y > $pi} { 820 set y [expr {$y - $twopi}] 821 } 822 return [list $x $y] 823} 824 825# ::mapproj::fromCassini -- 826# 827# Converts Cassini map co-ordinates to latitude and longitude. 828# 829# Parameters: 830# lambda_0 -- Longitude of the center of projection 831# phi_0 -- Latitude of the center of the sheet 832# x,y -- normalized x and y co-ordinates of a point on the map 833# 834# Results: 835# Returns a list consisting of the longitude and latitude in degrees. 836 837proc ::mapproj::fromCassini {lambda_0 phi_0 x y} { 838 variable degree 839 variable radian 840 set y [expr {$y + $degree * $phi_0}] 841 set phi [expr {$radian * asin(cos($x) * sin($y))}] 842 set lambda [expr {180. + $lambda_0 + $radian * atan2(tan($x), cos($y))}] 843 if {$lambda < 0.0 || $lambda > 360.0} { 844 set lambda [expr {$lambda - 360. * floor($lambda / 360.)}] 845 } 846 set lambda [expr {$lambda - 180.}] 847 return [list $lambda $phi] 848} 849 850# ::mapproj::toPeirceQuincuncial 851# 852# Converts geodetic co-ordinates to the Peirce Quincuncial projection. 853# 854# Parameters: 855# lambda_0 - Longitude of the central meridian. (Conventionally, 20.0). 856# lambda - Longitude of the point to be projected in degrees 857# phi - Latitude of the point to be projected in degrees. 858# 859# Results: 860# Returns a list of the x and y co-ordinates. 861 862proc ::mapproj::toPeirceQuincuncial {lambda_0 lambda phi} { 863 variable degree 864 variable halfSqrt2 865 variable pi 866 variable quarterpi 867 variable mquarterpi 868 variable threequarterpi 869 variable mthreequarterpi 870 variable PeirceQuincuncialScale 871 872 # Convert latitude and longitude to radians relative to the 873 # central meridian 874 875 set lambda [expr {$lambda - $lambda_0 + 180.}] 876 if {$lambda < 0.0 || $lambda > 360.0} { 877 set lambda [expr {$lambda - 360. * floor($lambda / 360.)}] 878 } 879 set lambda [expr {($lambda - 180.) * $degree}] 880 set phi [expr {$phi * $degree}] 881 882 # Compute the auxiliary quantities 'm' and 'n'. Set 'm' to match 883 # the sign of 'lambda' and 'n' to be positive if |lambda| > pi/2 884 885 set cos_phiosqrt2 [expr {$halfSqrt2 * cos($phi)}] 886 set cos_lambda [expr {cos($lambda)}] 887 set sin_lambda [expr {sin($lambda)}] 888 set cos_a [expr {$cos_phiosqrt2 * ($sin_lambda + $cos_lambda)}] 889 set cos_b [expr {$cos_phiosqrt2 * ($sin_lambda - $cos_lambda)}] 890 set sin_a [expr {sqrt(1.0 - $cos_a * $cos_a)}] 891 set sin_b [expr {sqrt(1.0 - $cos_b * $cos_b)}] 892 set cos_a_cos_b [expr {$cos_a * $cos_b}] 893 set sin_a_sin_b [expr {$sin_a * $sin_b}] 894 set sin2_m [expr {1.0 + $cos_a_cos_b - $sin_a_sin_b}] 895 set sin2_n [expr {1.0 - $cos_a_cos_b - $sin_a_sin_b}] 896 if {$sin2_m < 0.0} {set sin2_m 0.0} 897 set sin_m [expr {sqrt($sin2_m)}] 898 if {$sin2_m > 1.0} { set sin2_m 1.0 } 899 set cos_m [expr {sqrt(1.0 - $sin2_m)}] 900 if {$sin_lambda < 0.0} { 901 set sin_m [expr {-$sin_m}] 902 } 903 if {$sin2_n < 0.0} { set sin2_n 0.0 } 904 set sin_n [expr {sqrt($sin2_n)}] 905 if {$sin2_n > 1.0} { set sin2_n 1.0 } 906 set cos_n [expr {sqrt(1.0 - $sin2_n)}] 907 if {$cos_lambda > 0.0} { 908 set sin_n [expr {-$sin_n}] 909 } 910 911 # Compute elliptic integrals to map the disc to the square 912 913 set x [ellFaux $cos_m $sin_m $halfSqrt2] 914 set y [ellFaux $cos_n $sin_n $halfSqrt2] 915 916 # Reflect the Southern Hemisphere outward 917 918 if {$phi < 0} { 919 if {$lambda < $mthreequarterpi} { 920 set y [expr {$PeirceQuincuncialScale - $y}] 921 } elseif {$lambda < $mquarterpi} { 922 set x [expr {-$PeirceQuincuncialScale - $x}] 923 } elseif {$lambda < $quarterpi} { 924 set y [expr {-$PeirceQuincuncialScale - $y}] 925 } elseif {$lambda < $threequarterpi} { 926 set x [expr {$PeirceQuincuncialScale - $x}] 927 } else { 928 set y [expr {$PeirceQuincuncialScale - $y}] 929 } 930 } 931 932 # Rotate the square by 45 degrees to fit the screen better 933 934 set X [expr {($x - $y) * $halfSqrt2}] 935 set Y [expr {($x + $y) * $halfSqrt2}] 936 937 return [list $X $Y] 938} 939 940# ::mapproj::fromPeirceQuincuncial -- 941# 942# Converts Peirce Quincuncial map co-ordinates to latitude and longitude. 943# 944# Parameters: 945# lambda_0 -- Longitude of the center of projection 946# x,y -- normalized x and y co-ordinates of a point on the map 947# 948# Results: 949# Returns a list consisting of the longitude and latitude in degrees. 950 951proc ::mapproj::fromPeirceQuincuncial {lambda_0 x y} { 952 variable halfSqrt2 953 variable radian 954 variable pi 955 variable halfpi 956 variable quarterpi 957 variable PeirceQuincuncialScale 958 variable PeirceQuincuncialLimit 959 960 # Rotate x and y 45 degrees 961 962 set X [expr {($x + $y) * $halfSqrt2}] 963 set Y [expr {($y - $x) * $halfSqrt2}] 964 965 # Reflect Southern Hemisphere into the Northern 966 967 set southern 0 968 if {$X < -$PeirceQuincuncialLimit} { 969 set X [expr {-$PeirceQuincuncialScale - $X}] 970 set southern 1 971 } elseif {$X > $PeirceQuincuncialLimit} { 972 set X [expr {$PeirceQuincuncialScale - $X}] 973 set southern 1 974 } elseif {$Y < -$PeirceQuincuncialLimit} { 975 set Y [expr {-$PeirceQuincuncialScale - $Y}] 976 set southern 1 977 } elseif {$Y > $PeirceQuincuncialLimit} { 978 set Y [expr {$PeirceQuincuncialScale - $Y}] 979 set southern 1 980 } 981 982 # Now we know that latitude will be positive. If X is negative, then 983 # longitude will be negative; reflect the Western Hemisphere into the 984 # Eastern. 985 986 set western 0 987 if {$X < 0.0} { 988 set western 1 989 set X [expr {-$X}] 990 } 991 992 # If Y is positive, the point is in the back hemisphere. Reflect 993 # it to the front. 994 995 set back 0 996 if {$Y > 0.0} { 997 set back 1 998 set Y [expr {-$Y}] 999 } 1000 1001 # Finally, constrain longitude to be less than pi/4, by reflecting across 1002 # the 45 degree meridian. 1003 1004 set complement 0 1005 if {$X > -$Y} { 1006 set complement 1 1007 set t [expr {-$X}] 1008 set X [expr {-$Y}] 1009 set Y $t 1010 } 1011 1012 # Compute the elliptic functions to map the plane onto the sphere 1013 1014 set cnx [cn $X $halfSqrt2] 1015 set cny [cn $Y $halfSqrt2] 1016 1017 # Undo the mapping to latitude and longitude 1018 1019 set a1 [expr {acos(-$cnx * $cnx)}] 1020 set a2 [expr {acos($cny * $cny)}] 1021 set b [expr {0.5 * ($a1 + $a2)}] 1022 set a [expr {0.5 * ($a1 - $a2)}] 1023 set cos_a [expr {cos($a)}] 1024 set cos_b [expr {-cos($b)}] 1025 set lambda [expr {$quarterpi - atan2($cos_b, $cos_a)}] 1026 set phi [expr {acos(hypot($cos_b, $cos_a))}] 1027 1028 # Undo the reflections that were done above, to get correct latitude 1029 # and longitude 1030 1031 if {$complement} { 1032 set lambda [expr {$halfpi - $lambda}] 1033 } 1034 if {$back} { 1035 set lambda [expr {$pi - $lambda}] 1036 } 1037 if {$western} { 1038 set lambda [expr {-$lambda}] 1039 } 1040 if {$southern} { 1041 set phi [expr {-$phi}] 1042 } 1043 1044 # Convert latitude and longitude to degrees 1045 1046 set lambda [expr {$lambda * $radian + 180. + $lambda_0}] 1047 if {$lambda < 0.0 || $lambda > 360.0} { 1048 set lambda [expr {$lambda - 360. * floor($lambda / 360.)}] 1049 } 1050 set lambda [expr {$lambda - 180.}] 1051 set phi [expr {$phi * $radian}] 1052 1053 return [list $lambda $phi] 1054} 1055 1056# ::mapproj::RotateCartesianY 1057# 1058# Rotates Cartesian co-ordinates about the y axis 1059# 1060# Parameters: 1061# phi - Angle (in degrees) about which to rotate. 1062# x,y,z - Cartesian co-ordinates 1063# 1064# Results: 1065# Returns a three-element list giving the rotated co-ordinates. 1066 1067proc ::mapproj::RotateCartesianY {phi x y z} { 1068 variable degree 1069 set phi [expr {$degree * $phi}] 1070 set cos_phi [expr {cos($phi)}] 1071 set sin_phi [expr {sin($phi)}] 1072 return [list [expr {$x * $cos_phi - $z * $sin_phi}] \ 1073 $y \ 1074 [expr {$z * $cos_phi + $x * $sin_phi}]] 1075} 1076 1077# ::mapproj::ToCartesian -- 1078# 1079# Converts geodetic co-ordinates to Cartesian 1080# 1081# Parameters: 1082# lambda - Longitude of the point to be projected, in degrees 1083# phi - Latitude of the point to be projected, in degrees 1084# 1085# Results: 1086# Returns a three-element list, x, y, z where x is the component 1087# in the direction of longitude 0, latitude 0, y is the component 1088# in the direction of longitude 90 East, latitude 0, and 1089# z is the component in the direction of the North Pole 1090# 1091# Auxiliary procedure used in several projections to convert 1092# geodetic coordinates to Cartesian range and bearing. 1093 1094proc ::mapproj::ToCartesian {lambda phi} { 1095 variable degree 1096 set lambda [expr {$degree * $lambda}] 1097 set phi [expr {$degree * $phi}] 1098 set cos_phi [expr cos($phi)] 1099 return [list [expr {$cos_phi * cos($lambda)}] \ 1100 [expr {$cos_phi * sin($lambda)}] \ 1101 [expr {sin($phi)}]] 1102} 1103 1104# ::mapproj::CartesianToRangeAndBearing 1105# 1106# Transforms view-relative Cartesian co-ordinates to range and 1107# bearing. 1108# 1109# Parameters: 1110# x,y,z - Cartesian co-ordinates relative to center of Earth; 1111# +x points to the viewer and +z to the "view-up" direction. 1112# 1113# Results: 1114# Returns a three-element list containing, in order, 1115# the cosine (easting) of the bearing, the sine (northing) of the 1116# bearing, and the range. 1117 1118proc ::mapproj::CartesianToRangeAndBearing {x y z} { 1119 set c [expr {hypot($z, $y)}] 1120 if {$c == 0} { 1121 set cos_b 1.0 1122 set sin_b 0.0 1123 } else { 1124 set cos_b [expr {$y / $c}] 1125 set sin_b [expr {$z / $c}] 1126 } 1127 set range [expr {atan2($c, $x)}] 1128 return [list $cos_b $sin_b $range] 1129} 1130 1131# ::mapproj::RangeAndBearingToCartesian -- 1132# 1133# Converts range and bearing to Cartesian co-ordinates. 1134# 1135# Parameters: 1136# cos_b, sin_b -- Cosine (easting) and sine (northing) of the bearing 1137# range - Range, in Earth radii 1138# 1139# Results: 1140# Returns Cartesian co-ordinates relative to center of Earth. 1141# x is toward the station, and z is "view up" 1142 1143proc ::mapproj::RangeAndBearingToCartesian {cos_b sin_b range} { 1144 set c [expr {sin($range)}] 1145 set x [expr {cos($range)}] 1146 set y [expr {$cos_b * $c}] 1147 set z [expr {$sin_b * $c}] 1148 return [list $x $y $z] 1149} 1150 1151 1152# ::mapproj::CartesianToSpherical -- 1153# 1154# Transforms Cartesian x, y, z to spherical co-ordinates 1155# 1156# Parameters: 1157# x, y, z -- Coordinates of a point on the surface of the Earth, 1158# in Earth radii 1159# 1160# Results: 1161# Returns a two-element list comprising longitude and latitude 1162# in radians 1163 1164proc ::mapproj::CartesianToSpherical {x y z} { 1165 return [list [expr {atan2($y, $x)}] [expr {atan2($z, hypot($y, $x))}]] 1166} 1167 1168# ::mapproj::toOrthographic -- 1169# 1170# Transforms latitude and longitude to x and y co-ordinates 1171# on an orthographic projection. Scale is true only at the 1172# point of projection. 1173# 1174# Parameters: 1175# lambda_0, phi_0 -- Longitude and latitude of the center of projection 1176# in degrees 1177# lambda, phi -- Longitude and latitude of the point to be projected 1178# in degrees 1179# 1180# Results: 1181# Returns map x and y co-ordinates, in Earth radii. 1182 1183proc ::mapproj::toOrthographic {lambda_0 phi_0 lambda phi} { 1184 foreach {x y z} [ToCartesian [expr {$lambda-$lambda_0}] $phi] \ 1185 break 1186 foreach {x y z} [RotateCartesianY [expr {-$phi_0}] $x $y $z] \ 1187 break 1188 if {$x < 0} { 1189 return {} 1190 } else { 1191 return [list $y $z] 1192 } 1193} 1194 1195# ::mapproj::fromOrthographic -- 1196# 1197# Transforms x and y on an orthographic projection to latitude 1198# and longitude. 1199# 1200# Parameters: 1201# lambda_0, phi_0 -- Longitude and latitude of the center of projection 1202# in degrees 1203# x, y -- Co-ordinates of the projected point, in Earth radii 1204# 1205# Results: 1206# Returns a two element list containing longitude and latitude 1207# in degrees. 1208 1209proc ::mapproj::fromOrthographic {lambda_0 phi_0 x y} { 1210 variable radian 1211 set r [expr {hypot($x, $y)}] 1212 set alpha [expr {asin($r)}] 1213 set z [expr {sqrt(1.0 - $r*$r)}] 1214 foreach {x y z} [RotateCartesianY $phi_0 $z $x $y] break 1215 foreach {lambda phi} [CartesianToSpherical $x $y $z] break 1216 1217 # Convert latitude and longitude to degrees 1218 1219 set lambda [expr {$lambda * $radian + 180. + $lambda_0}] 1220 if {$lambda < 0.0 || $lambda > 360.0} { 1221 set lambda [expr {$lambda - 360. * floor($lambda / 360.)}] 1222 } 1223 set lambda [expr {$lambda - 180.}] 1224 set phi [expr {$phi * $radian}] 1225 1226 return [list $lambda $phi] 1227} 1228 1229# ::mapproj::toStereographic -- 1230# 1231# Transforms latitude and longitude to x and y co-ordinates 1232# on an orthographic projection. Scale is true only at the 1233# point of projection. 1234# 1235# Parameters: 1236# lambda_0, phi_0 -- Longitude and latitude of the center of projection 1237# in degrees 1238# lambda, phi -- Longitude and latitude of the point to be projected 1239# in degrees 1240# 1241# Results: 1242# Returns map x and y co-ordinates, in Earth radii. 1243 1244proc ::mapproj::toStereographic {lambda_0 phi_0 lambda phi} { 1245 foreach {x y z} [ToCartesian [expr {$lambda-$lambda_0}] $phi] \ 1246 break 1247 foreach {x y z} [RotateCartesianY [expr {-$phi_0}] $x $y $z] \ 1248 break 1249 if {$x < -0.5} { 1250 return {} 1251 } else { 1252 set y [expr {2. * $y / (1. + $x)}] 1253 set z [expr {2. * $z / (1. + $x)}] 1254 return [list $y $z] 1255 } 1256} 1257 1258# ::mapproj::fromStereographic -- 1259# 1260# Transforms x and y on an orthographic projection to latitude 1261# and longitude. 1262# 1263# Parameters: 1264# lambda_0, phi_0 -- Longitude and latitude of the center of projection 1265# in degrees 1266# x, y -- Co-ordinates of the projected point, in Earth radii 1267# 1268# Results: 1269# Returns a two element list containing longitude and latitude 1270# in degrees. 1271 1272proc ::mapproj::fromStereographic {lambda_0 phi_0 x y} { 1273 variable radian 1274 variable halfpi 1275 set denom [expr {4.0 + $x*$x + $y*$y}] 1276 foreach {x y z} [list \ 1277 [expr {(4.0 - $x*$x - $y*$y) / $denom}] \ 1278 [expr {4. * $x / $denom}] \ 1279 [expr {4. * $y / $denom}]] break 1280 1281 foreach {x y z} [RotateCartesianY $phi_0 $x $y $z] break 1282 foreach {lambda phi} [CartesianToSpherical $x $y $z] break 1283 1284 # Convert latitude and longitude to degrees 1285 1286 set lambda [expr {$lambda * $radian + 180. + $lambda_0}] 1287 if {$lambda < 0.0 || $lambda > 360.0} { 1288 set lambda [expr {$lambda - 360. * floor($lambda / 360.)}] 1289 } 1290 set lambda [expr {$lambda - 180.}] 1291 set phi [expr {$phi * $radian}] 1292 1293 return [list $lambda $phi] 1294} 1295 1296# ::mapproj::toGnomonic -- 1297# 1298# Transforms latitude and longitude to x and y co-ordinates 1299# on an orthographic projection. Scale is true only at the 1300# point of projection. 1301# 1302# Parameters: 1303# lambda_0, phi_0 -- Longitude and latitude of the center of projection 1304# in degrees 1305# lambda, phi -- Longitude and latitude of the point to be projected 1306# in degrees 1307# 1308# Results: 1309# Returns map x and y co-ordinates, in Earth radii. 1310 1311proc ::mapproj::toGnomonic {lambda_0 phi_0 lambda phi} { 1312 foreach {x y z} [ToCartesian [expr {$lambda-$lambda_0}] $phi] \ 1313 break 1314 foreach {x y z} [RotateCartesianY [expr {-$phi_0}] $x $y $z] \ 1315 break 1316 if {$x < 0.01} { 1317 return {} 1318 } else { 1319 set y [expr {$y / $x}] 1320 set z [expr {$z / $x}] 1321 return [list $y $z] 1322 } 1323} 1324 1325# ::mapproj::fromGnomonic -- 1326# 1327# Transforms x and y on an orthographic projection to latitude 1328# and longitude. 1329# 1330# Parameters: 1331# lambda_0, phi_0 -- Longitude and latitude of the center of projection 1332# in degrees 1333# x, y -- Co-ordinates of the projected point, in Earth radii 1334# 1335# Results: 1336# Returns a two element list containing longitude and latitude 1337# in degrees. 1338 1339proc ::mapproj::fromGnomonic {lambda_0 phi_0 x y} { 1340 variable radian 1341 variable halfpi 1342 set denom [expr {hypot(1.0, hypot($x, $y))}] 1343 foreach {x y z} [list \ 1344 [expr {1.0 / $denom}] \ 1345 [expr {$x / $denom}] \ 1346 [expr {$y / $denom}]] break 1347 1348 foreach {x y z} [RotateCartesianY $phi_0 $x $y $z] break 1349 foreach {lambda phi} [CartesianToSpherical $x $y $z] break 1350 1351 # Convert latitude and longitude to degrees 1352 1353 set lambda [expr {$lambda * $radian + 180. + $lambda_0}] 1354 if {$lambda < 0.0 || $lambda > 360.0} { 1355 set lambda [expr {$lambda - 360. * floor($lambda / 360.)}] 1356 } 1357 set lambda [expr {$lambda - 180.}] 1358 set phi [expr {$phi * $radian}] 1359 1360 return [list $lambda $phi] 1361} 1362 1363# ::mapproj::toAzimuthalEquidistant -- 1364# 1365# Transforms latitude and longitude to x and y co-ordinates 1366# on an orthographic projection. Scale is true only at the 1367# point of projection. 1368# 1369# Parameters: 1370# lambda_0, phi_0 -- Longitude and latitude of the center of projection 1371# in degrees 1372# lambda, phi -- Longitude and latitude of the point to be projected 1373# in degrees 1374# 1375# Results: 1376# Returns map x and y co-ordinates, in Earth radii. 1377 1378proc ::mapproj::toAzimuthalEquidistant {lambda_0 phi_0 lambda phi} { 1379 foreach {x y z} [ToCartesian [expr {$lambda-$lambda_0}] $phi] \ 1380 break 1381 foreach {x y z} [RotateCartesianY [expr {-$phi_0}] $x $y $z] \ 1382 break 1383 foreach {cs sn range} [CartesianToRangeAndBearing $x $y $z] break 1384 return [list [expr {$cs * $range}] [expr {$sn * $range}]] 1385} 1386 1387# ::mapproj::fromAzimuthalEquidistant -- 1388# 1389# Transforms x and y on an orthographic projection to latitude 1390# and longitude. 1391# 1392# Parameters: 1393# lambda_0, phi_0 -- Longitude and latitude of the center of projection 1394# in degrees 1395# x, y -- Co-ordinates of the projected point, in Earth radii 1396# 1397# Results: 1398# Returns a two element list containing longitude and latitude 1399# in degrees. 1400 1401proc ::mapproj::fromAzimuthalEquidistant {lambda_0 phi_0 x y} { 1402 variable radian 1403 variable halfpi 1404 1405 set range [expr {hypot($y, $x)}] 1406 set cos_b [expr {$x / $range}] 1407 set sin_b [expr {$y / $range}] 1408 foreach {x y z} [RangeAndBearingToCartesian $cos_b $sin_b $range] break 1409 foreach {x y z} [RotateCartesianY $phi_0 $x $y $z] break 1410 foreach {lambda phi} [CartesianToSpherical $x $y $z] break 1411 1412 # Convert latitude and longitude to degrees 1413 1414 set lambda [expr {$lambda * $radian + 180. + $lambda_0}] 1415 if {$lambda < 0.0 || $lambda > 360.0} { 1416 set lambda [expr {$lambda - 360. * floor($lambda / 360.)}] 1417 } 1418 set lambda [expr {$lambda - 180.}] 1419 set phi [expr {$phi * $radian}] 1420 1421 return [list $lambda $phi] 1422} 1423 1424# ::mapproj::toLambertAzimuthalEqualArea -- 1425# 1426# Transforms latitude and longitude to x and y co-ordinates 1427# on an orthographic projection. Scale is true only at the 1428# point of projection. 1429# 1430# Parameters: 1431# lambda_0, phi_0 -- Longitude and latitude of the center of projection 1432# in degrees 1433# lambda, phi -- Longitude and latitude of the point to be projected 1434# in degrees 1435# 1436# Results: 1437# Returns map x and y co-ordinates, in Earth radii. 1438 1439proc ::mapproj::toLambertAzimuthalEqualArea {lambda_0 phi_0 lambda phi} { 1440 foreach {x y z} [ToCartesian [expr {$lambda-$lambda_0}] $phi] \ 1441 break 1442 foreach {x y z} [RotateCartesianY [expr {-$phi_0}] $x $y $z] \ 1443 break 1444 foreach {cs sn range} [CartesianToRangeAndBearing $x $y $z] break 1445 set range [expr {2.0 * sin(0.5 * $range)}] 1446 return [list [expr {$cs * $range}] [expr {$sn * $range}]] 1447} 1448 1449# ::mapproj::fromLambertAzimuthalEqualArea -- 1450# 1451# Transforms x and y on an orthographic projection to latitude 1452# and longitude. 1453# 1454# Parameters: 1455# lambda_0, phi_0 -- Longitude and latitude of the center of projection 1456# in degrees 1457# x, y -- Co-ordinates of the projected point, in Earth radii 1458# 1459# Results: 1460# Returns a two element list containing longitude and latitude 1461# in degrees. 1462 1463proc ::mapproj::fromLambertAzimuthalEqualArea {lambda_0 phi_0 x y} { 1464 variable radian 1465 variable halfpi 1466 1467 set range [expr {hypot($y, $x)}] 1468 set cos_b [expr {$x / $range}] 1469 set sin_b [expr {$y / $range}] 1470 set range [expr {2.0 * asin(0.5 * $range)}] 1471 foreach {x y z} [RangeAndBearingToCartesian $cos_b $sin_b $range] break 1472 foreach {x y z} [RotateCartesianY $phi_0 $x $y $z] break 1473 foreach {lambda phi} [CartesianToSpherical $x $y $z] break 1474 1475 # Convert latitude and longitude to degrees 1476 1477 set lambda [expr {$lambda * $radian + 180. + $lambda_0}] 1478 if {$lambda < 0.0 || $lambda > 360.0} { 1479 set lambda [expr {$lambda - 360. * floor($lambda / 360.)}] 1480 } 1481 set lambda [expr {$lambda - 180.}] 1482 set phi [expr {$phi * $radian}] 1483 1484 return [list $lambda $phi] 1485} 1486 1487# ::mapproj::toHammer -- 1488# 1489# Transforms latitude and longitude to x and y co-ordinates 1490# on an orthographic projection. Scale is true only at the 1491# point of projection. 1492# 1493# Parameters: 1494# lambda_0-- Longitude of the center of projection 1495# in degrees 1496# lambda, phi -- Longitude and latitude of the point to be projected 1497# in degrees 1498# 1499# Results: 1500# Returns map x and y co-ordinates, in Earth radii. 1501 1502proc ::mapproj::toHammer {lambda_0 lambda phi} { 1503 set lambda [expr {$lambda - $lambda_0 + 180.}] 1504 if {$lambda < 0.0 || $lambda > 360.0} { 1505 set lambda [expr {$lambda - 360. * floor($lambda / 360.)}] 1506 } 1507 set lambda [expr {$lambda - 180.0}] 1508 foreach {x y z} [ToCartesian [expr {$lambda/2.}] $phi] \ 1509 break 1510 foreach {cs sn range} [CartesianToRangeAndBearing $x $y $z] break 1511 set range [expr {2.0 * sin(0.5 * $range)}] 1512 return [list [expr {2.0 * $cs * $range}] [expr {$sn * $range}]] 1513} 1514 1515# ::mapproj::fromHammer -- 1516# 1517# Transforms x and y on an orthographic projection to latitude 1518# and longitude. 1519# 1520# Parameters: 1521# lambda_0, phi_0 -- Longitude and latitude of the center of projection 1522# in degrees 1523# x, y -- Co-ordinates of the projected point, in Earth radii 1524# 1525# Results: 1526# Returns a two element list containing longitude and latitude 1527# in degrees. 1528 1529proc ::mapproj::fromHammer {lambda_0 x y} { 1530 variable radian 1531 variable halfpi 1532 1533 set x [expr {0.5 * $x}] 1534 set range [expr {hypot($y, $x)}] 1535 set cos_b [expr {$x / $range}] 1536 set sin_b [expr {$y / $range}] 1537 set range [expr {2.0 * asin(0.5 * $range)}] 1538 foreach {x y z} [RangeAndBearingToCartesian $cos_b $sin_b $range] break 1539 foreach {lambda phi} [CartesianToSpherical $x $y $z] break 1540 1541 # Convert latitude and longitude to degrees 1542 1543 set lambda [expr {2.0 * $lambda * $radian + 180. + $lambda_0}] 1544 if {$lambda < 0.0 || $lambda > 360.0} { 1545 set lambda [expr {$lambda - 360. * floor($lambda / 360.)}] 1546 } 1547 set lambda [expr {$lambda - 180.}] 1548 set phi [expr {$phi * $radian}] 1549 1550 return [list $lambda $phi] 1551} 1552 1553# ::mapproj::toConicEquidistant 1554# 1555# Converts latitude and longitude to map co-ordinates on a 1556# conic equidistant projection. 1557# 1558# Parameters: 1559# lambda_0, phi_0 -- Longitude and latitude of the center of the sheet. 1560# phi_1, phi_2 -- Latitudes of the two standard parallels at which scale 1561# is true 1562# lambda, phi -- Longitude and latitude of the point to be projected 1563# 1564# Results: 1565# Returns a list of map x and y measured in Earth radii. 1566 1567proc ::mapproj::toConicEquidistant {lambda_0 phi_0 phi_1 phi_2 lambda phi} { 1568 variable degree 1569 set phi_0 [expr {$phi_0 * $degree}] 1570 set phi_1 [expr {$phi_1 * $degree}] 1571 set phi_2 [expr {$phi_2 * $degree}] 1572 set lambda [expr {$lambda - $lambda_0 + 180.}] 1573 if {$lambda < 0.0 || $lambda > 360.0} { 1574 set lambda [expr {$lambda - 360. * floor($lambda / 360.)}] 1575 } 1576 set lambda [expr {($lambda - 180.0) * $degree}] 1577 set phi [expr {$phi * $degree}] 1578 set cos_phi_1 [expr {cos($phi_1)}] 1579 set n [expr {($cos_phi_1 - cos($phi_2)) / ($phi_2 - $phi_1)}] 1580 set G [expr {$cos_phi_1 / $n + $phi_1}] 1581 set rho_0 [expr {$G - $phi_0}] 1582 set theta [expr {$n * $lambda}] 1583 set rho [expr {$G - $phi}] 1584 set x [expr {$rho * sin($theta)}] 1585 set y [expr {$rho_0 - $rho * cos($theta)}] 1586 return [list $x $y] 1587} 1588 1589# ::mapproj::fromConicEquidistant -- 1590# 1591# Unprojects map x and y in a conic equidistant projection to 1592# latitude and longitude 1593# 1594# Parameters: 1595# lambda_0, phi_0 -- Longitude and latitude of the center of the sheet. 1596# phi_1, phi_2 -- Latitudes of the two standard parallels at which scale 1597# is true 1598# x, y -- Map co-ordinates in Earth radii 1599# 1600# Results: 1601# Returns a list of longitude and latitude in degrees. 1602 1603proc ::mapproj::fromConicEquidistant {lambda_0 phi_0 phi_1 phi_2 x y} { 1604 variable degree 1605 variable radian 1606 set phi_0 [expr {$phi_0 * $degree}] 1607 set phi_1 [expr {$phi_1 * $degree}] 1608 set phi_2 [expr {$phi_2 * $degree}] 1609 set cos_phi_1 [expr {cos($phi_1)}] 1610 set n [expr {($cos_phi_1 - cos($phi_2)) / ($phi_2 - $phi_1)}] 1611 set G [expr {$cos_phi_1 / $n + $phi_1}] 1612 set rho_0 [expr {$G - $phi_0}] 1613 set rho_0my [expr {$rho_0 - $y}] 1614 set theta [expr {atan2($x, $rho_0my)}] 1615 set rho [expr {sqrt($x*$x + $rho_0my * $rho_0my)}] 1616 if {$n < 0.0} {set rho [expr {-$rho}]} 1617 set phi [expr {($G - $rho) * $radian}] 1618 set lambda [expr {($theta / $n * $radian) + $lambda_0 + 180.}] 1619 if {$lambda < 0.0 || $lambda > 360.0} { 1620 set lambda [expr {$lambda - 360. * floor($lambda / 360.)}] 1621 } 1622 set lambda [expr {$lambda - 180.}] 1623 return [list $lambda $phi] 1624} 1625 1626# ::mapproj::toAlbersEqualAreaConic 1627# 1628# Converts latitude and longitude to map co-ordinates on a 1629# conic equal-area projection. 1630# 1631# Parameters: 1632# lambda_0, phi_0 -- Longitude and latitude of the center of the sheet. 1633# phi_1, phi_2 -- Latitudes of the two standard parallels at which scale 1634# is true 1635# lambda, phi -- Longitude and latitude of the point to be projected 1636# 1637# Results: 1638# Returns a list of map x and y measured in Earth radii. 1639 1640proc ::mapproj::toAlbersEqualAreaConic {lambda_0 phi_0 phi_1 phi_2 lambda phi} { 1641 variable degree 1642 set phi_0 [expr {$phi_0 * $degree}] 1643 set phi_1 [expr {$phi_1 * $degree}] 1644 set phi_2 [expr {$phi_2 * $degree}] 1645 set lambda [expr {$lambda - $lambda_0 + 180.}] 1646 if {$lambda < 0.0 || $lambda > 360.0} { 1647 set lambda [expr {$lambda - 360. * floor($lambda / 360.)}] 1648 } 1649 set lambda [expr {($lambda - 180.0) * $degree}] 1650 set phi [expr {$phi * $degree}] 1651 set cos_phi_1 [expr {cos($phi_1)}] 1652 set sin_phi_1 [expr {sin($phi_1)}] 1653 set n [expr {0.5 * ($sin_phi_1 + sin($phi_2))}] 1654 set theta [expr {$n * $lambda}] 1655 set C [expr {$cos_phi_1 * $cos_phi_1 + 2.0 * $n * $sin_phi_1}] 1656 set rho [expr {sqrt($C - 2.0 * $n * sin($phi)) / $n}] 1657 set rho_0 [expr {sqrt($C - 2.0 * $n * sin($phi_0)) / $n}] 1658 set x [expr {$rho * sin($theta)}] 1659 set y [expr {$rho_0 - $rho * cos($theta)}] 1660 return [list $x $y] 1661} 1662 1663# ::mapproj::fromAlbersEqualAreaConic -- 1664# 1665# Unprojects map x and y in a conic equal-area projection to 1666# latitude and longitude 1667# 1668# Parameters: 1669# lambda_0, phi_0 -- Longitude and latitude of the center of the sheet. 1670# phi_1, phi_2 -- Latitudes of the two standard parallels at which scale 1671# is true 1672# x, y -- Map co-ordinates in Earth radii 1673# 1674# Results: 1675# Returns a list of longitude and latitude in degrees. 1676 1677proc ::mapproj::fromAlbersEqualAreaConic {lambda_0 phi_0 phi_1 phi_2 x y} { 1678 variable degree 1679 variable radian 1680 variable twopi 1681 set phi_0 [expr {$phi_0 * $degree}] 1682 set phi_1 [expr {$phi_1 * $degree}] 1683 set phi_2 [expr {$phi_2 * $degree}] 1684 set cos_phi_1 [expr {cos($phi_1)}] 1685 set sin_phi_1 [expr {sin($phi_1)}] 1686 set n [expr {0.5 * ($sin_phi_1 + sin($phi_2))}] 1687 set C [expr {$cos_phi_1 * $cos_phi_1 + 2.0 * $n * $sin_phi_1}] 1688 set rho_0 [expr {sqrt($C - 2.0 * $n * sin($phi_0)) / $n}] 1689 set theta [expr {atan2($x, $rho_0 - $y)}] 1690 set rho [expr {hypot($x, $rho_0 - $y)}] 1691 set phi [expr {$radian * asin(($C - $rho*$rho*$n*$n) / (2.0 * $n))}] 1692 set lambda [expr {($theta / $n * $radian) + $lambda_0 + 180.}] 1693 if {$lambda < 0.0 || $lambda > 360.0} { 1694 set lambda [expr {$lambda - 360. * floor($lambda / 360.)}] 1695 } 1696 set lambda [expr {$lambda - 180.}] 1697 return [list $lambda $phi] 1698} 1699 1700# ::mapproj::toLambertConformalConic 1701# 1702# Converts latitude and longitude to map co-ordinates on a 1703# conformal conic projection. 1704# 1705# Parameters: 1706# lambda_0, phi_0 -- Longitude and latitude of the center of the sheet. 1707# phi_1, phi_2 -- Latitudes of the two standard parallels at which scale 1708# is true 1709# lambda, phi -- Longitude and latitude of the point to be projected 1710# 1711# Results: 1712# Returns a list of map x and y measured in Earth radii. 1713 1714proc ::mapproj::toLambertConformalConic {lambda_0 phi_0 phi_1 phi_2 lambda phi} { 1715 variable degree 1716 variable quarterpi 1717 set phi_0 [expr {$phi_0 * $degree}] 1718 set phi_1 [expr {$phi_1 * $degree}] 1719 set phi_2 [expr {$phi_2 * $degree}] 1720 set lambda [expr {$lambda - $lambda_0 + 180.}] 1721 if {$lambda < 0.0 || $lambda > 360.0} { 1722 set lambda [expr {$lambda - 360. * floor($lambda / 360.)}] 1723 } 1724 set lambda [expr {($lambda - 180.0) * $degree}] 1725 set phi [expr {$phi * $degree}] 1726 set cos_phi_1 [expr {cos($phi_1)}] 1727 set sin_phi_1 [expr {sin($phi_1)}] 1728 set tan1 [expr {tan($quarterpi + 0.5 * $phi_1)}] 1729 set n [expr {log($cos_phi_1 / cos($phi_2)) 1730 / log(tan($quarterpi + 0.5 * $phi_2) / $tan1)}] 1731 set F [expr {$cos_phi_1 * pow($tan1, $n) / $n}] 1732 set rho [expr {$F * pow(tan($quarterpi + 0.5 * $phi), -$n)}] 1733 set rho_0 [expr {$F * pow(tan($quarterpi + 0.5 * $phi_0), -$n)}] 1734 set x [expr {$rho * sin($n * $lambda)}] 1735 set y [expr {$rho_0 - $rho * cos($n * $lambda)}] 1736 return [list $x $y] 1737} 1738 1739# ::mapproj::fromLambertConformalConic -- 1740# 1741# Unprojects map x and y in a conformal conic projection to 1742# latitude and longitude 1743# 1744# Parameters: 1745# lambda_0, phi_0 -- Longitude and latitude of the center of the sheet. 1746# phi_1, phi_2 -- Latitudes of the two standard parallels at which scale 1747# is true 1748# x, y -- Map co-ordinates in Earth radii 1749# 1750# Results: 1751# Returns a list of longitude and latitude in degrees. 1752 1753proc ::mapproj::fromLambertConformalConic {lambda_0 phi_0 phi_1 phi_2 x y} { 1754 variable degree 1755 variable radian 1756 variable quarterpi 1757 set phi_0 [expr {$phi_0 * $degree}] 1758 set phi_1 [expr {$phi_1 * $degree}] 1759 set phi_2 [expr {$phi_2 * $degree}] 1760 set cos_phi_1 [expr {cos($phi_1)}] 1761 set sin_phi_1 [expr {sin($phi_1)}] 1762 set tan1 [expr {tan($quarterpi + 0.5 * $phi_1)}] 1763 set n [expr {log($cos_phi_1 / cos($phi_2)) 1764 / log(tan($quarterpi + 0.5 * $phi_2) / $tan1)}] 1765 set F [expr {$cos_phi_1 * pow($tan1, $n) / $n}] 1766 set rho_0 [expr {$F * pow(tan($quarterpi + 0.5 * $phi_0), -$n)}] 1767 set y [expr {$rho_0 - $y}] 1768 set rho [expr {sqrt($x*$x + $y*$y)}] 1769 if {$n < 0} { set rho [expr {-$rho}] } 1770 set theta [expr {atan2($x, $y)}] 1771 set phi [expr {$radian * 2 * atan(pow($F / $rho, 1.0 / $n)) - 90.}] 1772 set lambda [expr {($theta / $n * $radian) + $lambda_0 + 180.}] 1773 if {$lambda < 0.0 || $lambda > 360.0} { 1774 set lambda [expr {$lambda - 360. * floor($lambda / 360.)}] 1775 } 1776 set lambda [expr {$lambda - 180.}] 1777 return [list $lambda $phi] 1778} 1779 1780# Define commonly used cylindrical equal-area projections 1781 1782proc ::mapproj::toLambertCylindricalEqualArea {lambda_0 phi_0 lambda phi} { 1783 toCylindricalEqualArea 0.0 $lambda_0 $phi_0 $lambda $phi 1784} 1785proc ::mapproj::fromLambertCylindricalEqualArea {lambda_0 phi_0 x y} { 1786 fromCylindricalEqualArea 0.0 $lambda_0 $phi_0 $x $y 1787} 1788proc ::mapproj::toBehrmann {lambda_0 phi_0 lambda phi} { 1789 toCylindricalEqualArea 30.0 $lambda_0 $phi_0 $lambda $phi 1790} 1791proc ::mapproj::fromBehrmann {lambda_0 phi_0 x y} { 1792 fromCylindricalEqualArea 30.0 $lambda_0 $phi_0 $x $y 1793} 1794proc ::mapproj::toTrystanEdwards {lambda_0 phi_0 lambda phi} { 1795 toCylindricalEqualArea 37.4 $lambda_0 $phi_0 $lambda $phi 1796} 1797proc ::mapproj::fromTrystanEdwards {lambda_0 phi_0 x y} { 1798 fromCylindricalEqualArea 37.4 $lambda_0 $phi_0 $x $y 1799} 1800proc ::mapproj::toHoboDyer {lambda_0 phi_0 lambda phi} { 1801 toCylindricalEqualArea 37.5 $lambda_0 $phi_0 $lambda $phi 1802} 1803proc ::mapproj::fromHoboDyer {lambda_0 phi_0 x y} { 1804 fromCylindricalEqualArea 37.5 $lambda_0 $phi_0 $x $y 1805} 1806proc ::mapproj::toGallPeters {lambda_0 phi_0 lambda phi} { 1807 toCylindricalEqualArea 45.0 $lambda_0 $phi_0 $lambda $phi 1808} 1809proc ::mapproj::fromGallPeters {lambda_0 phi_0 x y} { 1810 fromCylindricalEqualArea 45.0 $lambda_0 $phi_0 $x $y 1811} 1812proc ::mapproj::toBalthasart {lambda_0 phi_0 lambda phi} { 1813 toCylindricalEqualArea 50.0 $lambda_0 $phi_0 $lambda $phi 1814} 1815proc ::mapproj::fromBalthasart {lambda_0 phi_0 x y} { 1816 fromCylindricalEqualArea 50.0 $lambda_0 $phi_0 $x $y 1817} 1818