diff options
author | William Joye <wjoye@cfa.harvard.edu> | 2016-10-27 19:39:39 (GMT) |
---|---|---|
committer | William Joye <wjoye@cfa.harvard.edu> | 2016-10-27 19:39:39 (GMT) |
commit | ea28451286d3ea4a772fa174483f9a7a66bb1ab3 (patch) | |
tree | 6ee9d8a7848333a7ceeee3b13d492e40225f8b86 /tcllib/modules/mapproj | |
parent | b5ca09bae0d6a1edce939eea03594dd56383f2c8 (diff) | |
parent | 7c621da28f07e449ad90c387344f07a453927569 (diff) | |
download | blt-ea28451286d3ea4a772fa174483f9a7a66bb1ab3.zip blt-ea28451286d3ea4a772fa174483f9a7a66bb1ab3.tar.gz blt-ea28451286d3ea4a772fa174483f9a7a66bb1ab3.tar.bz2 |
Merge commit '7c621da28f07e449ad90c387344f07a453927569' as 'tcllib'
Diffstat (limited to 'tcllib/modules/mapproj')
-rwxr-xr-x | tcllib/modules/mapproj/ChangeLog | 47 | ||||
-rwxr-xr-x | tcllib/modules/mapproj/mapproj.man | 308 | ||||
-rwxr-xr-x | tcllib/modules/mapproj/mapproj.tcl | 1817 | ||||
-rwxr-xr-x | tcllib/modules/mapproj/pkgIndex.tcl | 2 |
4 files changed, 2174 insertions, 0 deletions
diff --git a/tcllib/modules/mapproj/ChangeLog b/tcllib/modules/mapproj/ChangeLog new file mode 100755 index 0000000..ce691da --- /dev/null +++ b/tcllib/modules/mapproj/ChangeLog @@ -0,0 +1,47 @@ +2013-02-01 Andreas Kupries <andreas_kupries@users.sourceforge.net> + + * + * Released and tagged Tcllib 1.15 ======================== + * + +2011-12-13 Andreas Kupries <andreas_kupries@users.sourceforge.net> + + * + * Released and tagged Tcllib 1.14 ======================== + * + +2011-01-24 Andreas Kupries <andreas_kupries@users.sourceforge.net> + + * + * Released and tagged Tcllib 1.13 ======================== + * + +2009-12-07 Andreas Kupries <andreas_kupries@users.sourceforge.net> + + * + * Released and tagged Tcllib 1.12 ======================== + * + +2008-12-12 Andreas Kupries <andreas_kupries@users.sourceforge.net> + + * + * Released and tagged Tcllib 1.11.1 ======================== + * + +2008-10-16 Andreas Kupries <andreas_kupries@users.sourceforge.net> + + * + * Released and tagged Tcllib 1.11 ======================== + * + +2007-09-12 Andreas Kupries <andreas_kupries@users.sourceforge.net> + + * + * Released and tagged Tcllib 1.10 ======================== + * + +2007-08-24 Kevin B. Kenny <kennykb@acm.org> + + * mapproj.man: + * mapproj.tcl: + * pkgIndex.tcl: Initial commit. diff --git a/tcllib/modules/mapproj/mapproj.man b/tcllib/modules/mapproj/mapproj.man new file mode 100755 index 0000000..4e642ce --- /dev/null +++ b/tcllib/modules/mapproj/mapproj.man @@ -0,0 +1,308 @@ +[comment {-*- tcl -*- doctools manpage}] +[manpage_begin mapproj n 0.1] +[keywords geodesy] +[keywords map] +[keywords projection] +[copyright {2007 Kevin B. Kenny <kennykb@acm.org>}] +[moddesc {Tcl Library}] +[titledesc {Map projection routines}] +[require Tcl [opt 8.4]] +[require math::interpolate [opt 1.0]] +[require math::special [opt 0.2.1]] +[require mapproj [opt 1.0]] +[description] +The [package mapproj] package provides a set of procedures for +converting between world co-ordinates (latitude and longitude) and map +co-ordinates on a number of different map projections. + +[section Commands] + +The following commands convert between world co-ordinates and +map co-ordinates: + +[list_begin definitions] + +[call [cmd ::mapproj::toPlateCarree] [arg lambda_0] [arg phi_0] [arg lambda] [arg phi]] +Converts to the [emph "plate carr\u00e9e"] (cylindrical equidistant) +projection. +[call [cmd ::mapproj::fromPlateCarree] [arg lambda_0] [arg phi_0] [arg x] [arg y]] +Converts from the [emph "plate carr\u00e9e"] (cylindrical equidistant) +projection. +[call [cmd ::mapproj::toCylindricalEqualArea] [arg lambda_0] [arg phi_0] [arg lambda] [arg phi]] +Converts to the cylindrical equal-area projection. +[call [cmd ::mapproj::fromCylindricalEqualArea] [arg lambda_0] [arg phi_0] [arg x] [arg y]] +Converts from the cylindrical equal-area projection. +[call [cmd ::mapproj::toMercator] [arg lambda_0] [arg phi_0] [arg lambda] [arg phi]] +Converts to the Mercator (cylindrical conformal) projection. +[call [cmd ::mapproj::fromMercator] [arg lambda_0] [arg phi_0] [arg x] [arg y]] +Converts from the Mercator (cylindrical conformal) projection. +[call [cmd ::mapproj::toMillerCylindrical] [arg lambda_0] [arg lambda] [arg phi]] +Converts to the Miller Cylindrical projection. +[call [cmd ::mapproj::fromMillerCylindrical] [arg lambda_0] [arg x] [arg y]] +Converts from the Miller Cylindrical projection. +[call [cmd ::mapproj::toSinusoidal] [arg lambda_0] [arg phi_0] [arg lambda] [arg phi]] +Converts to the sinusoidal (Sanson-Flamsteed) projection. +projection. +[call [cmd ::mapproj::fromSinusoidal] [arg lambda_0] [arg phi_0] [arg x] [arg y]] +Converts from the sinusoidal (Sanson-Flamsteed) projection. +projection. +[call [cmd ::mapproj::toMollweide] [arg lambda_0] [arg lambda] [arg phi]] +Converts to the Mollweide projection. +[call [cmd ::mapproj::fromMollweide] [arg lambda_0] [arg x] [arg y]] +Converts from the Mollweide projection. +[call [cmd ::mapproj::toEckertIV] [arg lambda_0] [arg lambda] [arg phi]] +Converts to the Eckert IV projection. +[call [cmd ::mapproj::fromEckertIV] [arg lambda_0] [arg x] [arg y]] +Converts from the Eckert IV projection. +[call [cmd ::mapproj::toEckertVI] [arg lambda_0] [arg lambda] [arg phi]] +Converts to the Eckert VI projection. +[call [cmd ::mapproj::fromEckertVI] [arg lambda_0] [arg x] [arg y]] +Converts from the Eckert VI projection. +[call [cmd ::mapproj::toRobinson] [arg lambda_0] [arg lambda] [arg phi]] +Converts to the Robinson projection. +[call [cmd ::mapproj::fromRobinson] [arg lambda_0] [arg x] [arg y]] +Converts from the Robinson projection. +[call [cmd ::mapproj::toCassini] [arg lambda_0] [arg phi_0] [arg lambda] [arg phi]] +Converts to the Cassini (transverse cylindrical equidistant) +projection. +[call [cmd ::mapproj::fromCassini] [arg lambda_0] [arg phi_0] [arg x] [arg y]] +Converts from the Cassini (transverse cylindrical equidistant) +projection. +[call [cmd ::mapproj::toPeirceQuincuncial] [arg lambda_0] [arg lambda] [arg phi]] +Converts to the Peirce Quincuncial Projection. +[call [cmd ::mapproj::fromPeirceQuincuncial] [arg lambda_0] [arg x] [arg y]] +Converts from the Peirce Quincuncial Projection. +[call [cmd ::mapproj::toOrthographic] [arg lambda_0] [arg phi_0] [arg lambda] [arg phi]] +Converts to the orthographic projection. +[call [cmd ::mapproj::fromOrthographic] [arg lambda_0] [arg phi_0] [arg x] [arg y]] +Converts from the orthographic projection. +[call [cmd ::mapproj::toStereographic] [arg lambda_0] [arg phi_0] [arg lambda] [arg phi]] +Converts to the stereographic (azimuthal conformal) projection. +[call [cmd ::mapproj::fromStereographic] [arg lambda_0] [arg phi_0] [arg x] [arg y]] +Converts from the stereographic (azimuthal conformal) projection. +[call [cmd ::mapproj::toGnomonic] [arg lambda_0] [arg phi_0] [arg lambda] [arg phi]] +Converts to the gnomonic projection. +[call [cmd ::mapproj::fromGnomonic] [arg lambda_0] [arg phi_0] [arg x] [arg y]] +Converts from the gnomonic projection. +[call [cmd ::mapproj::toAzimuthalEquidistant] [arg lambda_0] [arg phi_0] [arg lambda] [arg phi]] +Converts to the azimuthal equidistant projection. +[call [cmd ::mapproj::fromAzimuthalEquidistant] [arg lambda_0] [arg phi_0] [arg x] [arg y]] +Converts from the azimuthal equidistant projection. +[call [cmd ::mapproj::toLambertAzimuthalEqualArea] [arg lambda_0] [arg phi_0] [arg lambda] [arg phi]] +Converts to the Lambert azimuthal equal-area projection. +[call [cmd ::mapproj::fromLambertAzimuthalEqualArea] [arg lambda_0] [arg phi_0] [arg x] [arg y]] +Converts from the Lambert azimuthal equal-area projection. +[call [cmd ::mapproj::toHammer] [arg lambda_0] [arg lambda] [arg phi]] +Converts to the Hammer projection. +[call [cmd ::mapproj::fromHammer] [arg lambda_0] [arg x] [arg y]] +Converts from the Hammer projection. +[call [cmd ::mapproj::toConicEquidistant] [arg lambda_0] [arg phi_0] [arg phi_1] [arg phi_2] [arg lambda] [arg phi]] +Converts to the conic equidistant projection. +[call [cmd ::mapproj::fromConicEquidistant] [arg lambda_0] [arg phi_0] [arg phi_1] [arg phi_2] [arg x] [arg y]] +Converts from the conic equidistant projection. +[call [cmd ::mapproj::toAlbersEqualAreaConic] [arg lambda_0] [arg phi_0] [arg phi_1] [arg phi_2] [arg lambda] [arg phi]] +Converts to the Albers equal-area conic projection. +[call [cmd ::mapproj::fromAlbersEqualAreaConic] [arg lambda_0] [arg phi_0] [arg phi_1] [arg phi_2] [arg x] [arg y]] +Converts from the Albers equal-area conic projection. +[call [cmd ::mapproj::toLambertConformalConic] [arg lambda_0] [arg phi_0] [arg phi_1] [arg phi_2] [arg lambda] [arg phi]] +Converts to the Lambert conformal conic projection. +[call [cmd ::mapproj::fromLambertConformalConic] [arg lambda_0] [arg phi_0] [arg phi_1] [arg phi_2] [arg x] [arg y]] +Converts from the Lambert conformal conic projection. + +[list_end] + +Among the cylindrical equal-area projections, there are a number of +choices of standard parallels that have names: + +[list_begin definitions] +[call [cmd ::mapproj::toLambertCylindricalEqualArea] [arg lambda_0] [arg phi_0] [arg lambda] [arg phi]] +Converts to the Lambert cylindrical equal area projection. (standard parallel +is the Equator.) +[call [cmd ::mapproj::fromLambertCylindricalEqualArea] [arg lambda_0] [arg phi_0] [arg x] [arg y]] +Converts from the Lambert cylindrical equal area projection. (standard parallel +is the Equator.) +[call [cmd ::mapproj::toBehrmann] [arg lambda_0] [arg phi_0] [arg lambda] [arg phi]] +Converts to the Behrmann cylindrical equal area projection. (standard parallels +are 30 degrees North and South) +[call [cmd ::mapproj::fromBehrmann] [arg lambda_0] [arg phi_0] [arg x] [arg y]] +Converts from the Behrmann cylindrical equal area projection. (standard parallels +are 30 degrees North and South.) +[call [cmd ::mapproj::toTrystanEdwards] [arg lambda_0] [arg phi_0] [arg lambda] [arg phi]] +Converts to the Trystan Edwards cylindrical equal area projection. (standard parallels +are 37.4 degrees North and South) +[call [cmd ::mapproj::fromTrystanEdwards] [arg lambda_0] [arg phi_0] [arg x] [arg y]] +Converts from the Trystan Edwards cylindrical equal area projection. (standard parallels +are 37.4 degrees North and South.) +[call [cmd ::mapproj::toHoboDyer] [arg lambda_0] [arg phi_0] [arg lambda] [arg phi]] +Converts to the Hobo-Dyer cylindrical equal area projection. (standard parallels +are 37.5 degrees North and South) +[call [cmd ::mapproj::fromHoboDyer] [arg lambda_0] [arg phi_0] [arg x] [arg y]] +Converts from the Hobo-Dyer cylindrical equal area projection. (standard parallels +are 37.5 degrees North and South.) +[call [cmd ::mapproj::toGallPeters] [arg lambda_0] [arg phi_0] [arg lambda] [arg phi]] +Converts to the Gall-Peters cylindrical equal area projection. (standard parallels +are 45 degrees North and South) +[call [cmd ::mapproj::fromGallPeters] [arg lambda_0] [arg phi_0] [arg x] [arg y]] +Converts from the Gall-Peters cylindrical equal area projection. (standard parallels +are 45 degrees North and South.) +[call [cmd ::mapproj::toBalthasart] [arg lambda_0] [arg phi_0] [arg lambda] [arg phi]] +Converts to the Balthasart cylindrical equal area projection. (standard parallels +are 50 degrees North and South) +[call [cmd ::mapproj::fromBalthasart] [arg lambda_0] [arg phi_0] [arg x] [arg y]] +Converts from the Balthasart cylindrical equal area projection. (standard parallels +are 50 degrees North and South.) + +[list_end] + +[section Arguments] + +The following arguments are accepted by the projection commands: + +[list_begin definitions] + +[def [arg lambda]] + +Longitude of the point to be projected, in degrees. + +[def [arg phi]] + +Latitude of the point to be projected, in degrees. + +[def [arg lambda_0]] + +Longitude of the center of the sheet, in degrees. For many projections, +this figure is also the reference meridian of the projection. + +[def [arg phi_0]] + +Latitude of the center of the sheet, in degrees. For the azimuthal +projections, this figure is also the latitude of the center of the projection. + +[def [arg phi_1]] + +Latitude of the first reference parallel, for projections that use reference +parallels. + +[def [arg phi_2]] + +Latitude of the second reference parallel, for projections that use reference +parallels. + +[def [arg x]] + +X co-ordinate of a point on the map, in units of Earth radii. + +[def [arg y]] + +Y co-ordinate of a point on the map, in units of Earth radii. + +[list_end] + +[section Results] + +For all of the procedures whose names begin with 'to', the return value +is a list comprising an [emph x] co-ordinate and a [emph y] co-ordinate. +The co-ordinates are relative to the center of the map sheet to be drawn, +measured in Earth radii at the reference location on the map. + +For all of the functions whose names begin with 'from', the return value +is a list comprising the longitude and latitude, in degrees. + +[section {Choosing a projection}] + +This package offers a great many projections, because no single projection +is appropriate to all maps. This section attempts to provide guidance +on how to choose a projection. +[para] +First, consider the type of data that you intend to display on the map. +If the data are [emph directional] ([emph e.g.,] winds, ocean currents, or +magnetic fields) then you need to use a projection that preserves +angles; these are known as [emph conformal] projections. Conformal +projections include the Mercator, the Albers azimuthal equal-area, +the stereographic, and the Peirce Quincuncial projection. If the +data are [emph thematic], describing properties of land or water, such +as temperature, population density, land use, or demographics; then +you need a projection that will show these data with the areas on the map +proportional to the areas in real life. These so-called [emph {equal area}] +projections include the various cylindrical equal area projections, +the sinusoidal projection, the Lambert azimuthal equal-area projection, +the Albers equal-area conic projection, and several of the world-map +projections (Miller Cylindrical, Mollweide, Eckert IV, Eckert VI, Robinson, +and Hammer). If the significant factor in your data is distance from a +central point or line (such as air routes), then you will do best with +an [emph equidistant] projection such as [emph "plate carr\u00e9e"], +Cassini, azimuthal equidistant, or conic equidistant. If direction from +a central point is a critical factor in your data (for instance, +air routes, radio antenna pointing), then you will almost surely want to +use one of the azimuthal projections. Appropriate choices are azimuthal +equidistant, azimuthal equal-area, stereographic, and perhaps orthographic. +[para] +Next, consider how much of the Earth your map will cover, and the general +shape of the area of interest. For maps of the entire Earth, +the cylindrical equal area, Eckert IV and VI, Mollweide, Robinson, and Hammer +projections are good overall choices. The Mercator projection is traditional, +but the extreme distortions of area at high latitudes make it +a poor choice unless a conformal projection is required. The Peirce +projection is a better choice of conformal projection, having less distortion +of landforms. The Miller Cylindrical is a compromise designed to give +shapes similar to the traditional Mercator, but with less polar stretching. +The Peirce Quincuncial projection shows all the continents with acceptable +distortion if a reference meridian close to +20 degrees is chosen. +The Robinson projection yields attractive maps for things like political +divisions, but should be avoided in presenting scientific data, since other +projections have moe useful geometric properties. +[para] +If the map will cover a hemisphere, then choose stereographic, +azimuthal-equidistant, Hammer, or Mollweide projections; these all project +the hemisphere into a circle. +[para] +If the map will cover a large area (at least a few hundred km on a side), +but less than +a hemisphere, then you have several choices. Azimuthal projections +are usually good (choose stereographic, azimuthal equidistant, or +Lambert azimuthal equal-area according to whether shapes, distances from +a central point, or areas are important). Azimuthal projections (and possibly +the Cassini projection) are the only +really good choices for mapping the polar regions. +[para] +If the large area is in one of the temperate zones and is round or has +a primarily east-west extent, then the conic projections are good choices. +Choose the Lambert conformal conic, the conic equidistant, or the Albers +equal-area conic according to whether shape, distance, or area are the +most important parameters. For any of these, the reference parallels +should be chosen at approximately 1/6 and 5/6 of the range of latitudes +to be displayed. For instance, maps of the 48 coterminous United States +are attractive with reference parallels of 28.5 and 45.5 degrees. +[para] +If the large area is equatorial and is round or has a primarily east-west +extent, then the Mercator projection is a good choice for a conformal +projection; Lambert cylindrical equal-area and sinusoidal projections are +good equal-area projections; and the [emph "plate carr\u00e9e"] is a +good equidistant projection. +[para] +Large areas having a primarily North-South aspect, particularly those +spanning the Equator, need some other choices. The Cassini projection +is a good choice for an equidistant projection (for instance, a Cassini +projection with a central meridian of 80 degrees West produces an +attractive map of the Americas). The cylindrical equal-area, Albers +equal-area conic, sinusoidal, Mollweide and Hammer +projections are possible choices for equal-area projections. +A good conformal projection in this situation is the Transverse +Mercator, which alas, is not yet implemented. +[para] +Small areas begin to get into a realm where the ellipticity of the +Earth affects the map scale. This package does not attempt to +handle accurate mapping for large-scale topographic maps. If +slight scale errors are acceptable in your application, then any +of the projections appropriate to large areas should work for +small ones as well. +[para] +There are a few projections that are included for their special +properties. The orthographic projection produces views of the +Earth as seen from space. The gnomonic projection produces a +map on which all great circles (the shortest distance between +two points on the Earth's surface) are rendered as straight lines. +While this projection is useful for navigational planning, it +has extreme distortions of shape and area, and can display +only a limited area of the Earth (substantially less than a hemisphere). +[manpage_end] diff --git a/tcllib/modules/mapproj/mapproj.tcl b/tcllib/modules/mapproj/mapproj.tcl new file mode 100755 index 0000000..1c545f9 --- /dev/null +++ b/tcllib/modules/mapproj/mapproj.tcl @@ -0,0 +1,1817 @@ +# mapproj.tcl -- +# +# Package for map projections. +# +# Copyright (c) 2007 by Kevin B. Kenny. +# +# See the file "license.terms" for information on usage and redistribution +# of this file, and for a DISCLAIMER OF ALL WARRANTIES. +# +# RCS: @(#) $Id: mapproj.tcl,v 1.1 2007/08/24 22:36:35 kennykb Exp $ +#------------------------------------------------------------------------------ + +package require Tcl 8.4 +package require math::interpolate 1.0 +package require math::special 0.2.1 + +package provide mapproj 1.0 + +# ::mapproj -- +# +# Namespace holding the procedures and values. + +namespace eval ::mapproj { + + # Cylindrical projections + + namespace export toPlateCarree fromPlateCarree + namespace export toCassini fromCassini + namespace export toCylindricalEqualArea fromCylindricalEqualArea + namespace export toMercator fromMercator + + # Named cylindric equal-area projections with specific + # standard parallels + + namespace export toLambertCylindricalEqualArea \ + fromLambertCylindricalEqualArea + namespace export toBehrmann fromBehrmann + namespace export toTrystanEdwards fromTrystanEdwards + namespace export toHoboDyer fromHoboDyer + namespace export toGallPeters fromGallPeters + namespace export toBalthasart fromBalthasart + + # Pseudocylindrical projections - equal area + + namespace export toSinusoidal fromSinusoidal + namespace export toMillerCylindrical fromMillerCylindrical + namespace export toMollweide fromMollweide + namespace export toEckertIV fromEckertIV toEckertVI fromEckertVI + + # Pseudocylindrical projections - compromise + + namespace export toRobinson fromRobinson + + # Azimuthal projections + + namespace export toAzimuthalEquidistant fromAzimuthalEquidistant + namespace export toLambertAzimuthalEqualArea fromLambertAzimuthalEqualArea + namespace export toStereographic fromStereographic + namespace export toOrthographic fromOrthographic + namespace export toGnomonic fromGnomonic + + # Pseudo-azimuthal projections + + namespace export toHammer fromHammer + + # Conic projections + + namespace export toConicEquidistant fromConicEquidistant + namespace export toLambertConformalConic fromLambertConformalConic + namespace export toAlbersEqualAreaConic fromAlbersEqualAreaConic + + # Miscellaneous projections + + namespace export toPeirceQuincuncial fromPeirceQuincuncial + + # Fundamental constants + + variable pi [expr {acos(-1.0)}] + variable twopi [expr {2.0 * $pi}] + variable halfpi [expr {0.5 * $pi}] + variable quarterpi [expr {0.25 * $pi}] + variable threequarterpi [expr {0.75 * $pi}] + variable mquarterpi [expr {-0.25 * $pi}] + variable mthreequarterpi [expr {-0.75 * $pi}] + variable radian [expr {180. / $pi}] + variable degree [expr {$pi / 180.}] + variable sqrt2 [expr {sqrt(2.)}] + variable sqrt8 [expr {2. * $sqrt2}] + variable halfSqrt3 [expr {sqrt(3.) / 2.}] + variable halfSqrt2 [expr {sqrt(2.) / 2.}] + variable EckertIVK1 [expr {2.0 / sqrt($pi * (4.0 + $pi))}] + variable EckertIVK2 [expr {2.0 * sqrt($pi / (4.0 + $pi))}] + variable EckertVIK1 [expr {sqrt(2.0 + $pi)}] + variable PeirceQuincuncialScale 3.7081493546027438 ;# 2*K(1/2) + variable PeirceQuincuncialLimit 1.8540746773013719 ;# K(1/2) + + # Table of parallel length and distance from equator for the + # Robinson projection + + variable RobinsonLatitude { + -90.0 -85.0 -80.0 -75.0 -70.0 -65.0 + -60.0 -55.0 -50.0 -45.0 -40.0 -35.0 + -30.0 -25.0 -20.0 -15.0 -10.0 -5.0 + 0.0 5.0 10.0 15.0 20.0 25.0 30.0 + 35.0 40.0 45.0 50.0 55.0 60.0 + 65.0 70.0 75.0 80.0 85.0 90.0 + } + variable RobinsonPLEN { + 0.5322 0.5722 0.6213 0.6732 0.7186 0.7597 + 0.7986 0.8350 0.8679 0.8962 0.9216 0.9427 + 0.9600 0.9730 0.9822 0.9900 0.9954 0.9986 + 1.0000 0.9986 0.9954 0.9900 0.9822 0.9730 0.9600 + 0.9427 0.9216 0.8962 0.8679 0.8350 0.7986 + 0.7597 0.7186 0.6732 0.6213 0.5722 0.5322 + } + variable RobinsonPDFE { + -1.0000 -0.9761 -0.9394 -0.8936 -0.8435 -0.7903 + -0.7346 -0.6769 -0.6176 -0.5571 -0.4958 -0.4340 + -0.3720 -0.3100 -0.2480 -0.1860 -0.1240 -0.0620 + 0.0000 0.0620 0.1240 0.1860 0.2480 0.3100 0.3720 + 0.4340 0.4958 0.5571 0.6176 0.6769 0.7346 + 0.7903 0.8435 0.8936 0.9394 0.9761 1.0000 + } + + # Interpolation tables for Robinson + + variable RobinsonSplinePLEN \ + [math::interpolate::prepare-cubic-splines \ + $RobinsonLatitude $RobinsonPLEN] + variable RobinsonSplinePDFE \ + [math::interpolate::prepare-cubic-splines \ + $RobinsonLatitude $RobinsonPDFE] + variable RobinsonM [expr {0.5072 * $pi}] + + namespace import ::math::special::cn + +} + +# ::mapproj::ellF - +# +# Computes the Legendre incomplete elliptic integral of the +# first kind: +# +# F(phi, k) = \integral_0^phi dtheta/sqrt(1 - k**2 sin**2 theta) +# +# +# Parameters: +# phi -- Limit of integration; angle around the ellipse +# k -- Eccentricity +# +# Results: +# Returns F(phi, k) +# +# Notes: +# We compute this integral in terms of the Carlson elliptic integral +# ellRF(x, y, z). + +proc ::mapproj::ellF {phi k} { + return [ellFaux [expr {cos($phi)}] [expr {sin($phi)}] $k] +} + +# ::mapproj::ellFaux - +# +# Computes the Legendre incomplete elliptic integral of the +# first kind when circular functions of the 'phi' argument +# are already available. +# +# Parameters: +# cos_phi - Cosine of the argument +# sin_phi - Sine of the argument +# k - Parameter +# +# Results: +# Returns F(atan(sin_phi/cos_phi), k) + +proc ::mapproj::ellFaux {cos_phi sin_phi k} { + set rf [ellRF [expr {$cos_phi * $cos_phi}] \ + [expr {1.0 - $k * $k * $sin_phi * $sin_phi}] \ + 1.0] + return [expr {$sin_phi * $rf}] +} + +# ::mapproj::ellRF -- +# +# Computes the Carlson incomplete elliptic integral of the +# first kind: +# +# RF(x, y, z) = 1/2 * integral_0^inf dt/sqrt((t+x)*(t+y)*(t+z)) +# +# Parameters: +# x, y, z -- Interchangeable parameters of the integral +# +# Results: +# Returns the value of RF + +proc ::mapproj::ellRF {x y z} { + if {$x < 0.0 || $y < 0.0 || $z < 0.0} { + return -code error "Negative argument to Carlson's ellRF" \ + -errorCode "ellRF negArgument" + } + set delx 1.0; set dely 1.0; set delz 1.0 + while {abs($delx) > 0.0025 || abs($dely) > 0.0025 || abs($delz) > 0.0025} { + set sx [expr {sqrt($x)}] + set sy [expr {sqrt($y)}] + set sz [expr {sqrt($z)}] + set len [expr {$sx * ($sy + $sz) + $sy * $sz}] + set x [expr {0.25 * ($x + $len)}] + set y [expr {0.25 * ($y + $len)}] + set z [expr {0.25 * ($z + $len)}] + set mean [expr {($x + $y + $z) / 3.0}] + set delx [expr {($mean - $x) / $mean}] + set dely [expr {($mean - $y) / $mean}] + set delz [expr {($mean - $z) / $mean}] + } + set e2 [expr {$delx * $dely - $delz * $delz}] + set e3 [expr {$delx * $dely * $delz}] + return [expr {(1.0 + ($e2 / 24.0 - 0.1 - 3.0 * $e3 / 44.0) * $e2 + + $e3 / 14.) / sqrt($mean)}] +} + +# ::mapproj::toPlateCarree -- +# +# Project a latitude and longitude onto the plate carrée. +# +# Parameters: +# phi_0 -- Latitude of the center of the sheet in degrees +# lambda_0 -- Longitude of the center of sheet in degrees +# lambda -- Longitude of the point to be projected in degrees +# phi -- Latitude of the point to be projected in degrees +# +# Results: +# Returns x and y co-ordinates. Units of x and y are Earth radii; +# scale is true at the Equator. + +proc ::mapproj::toPlateCarree {lambda_0 phi_0 lambda phi} { + variable degree + set lambda [expr {$lambda - $lambda_0 + 180.}] + if {$lambda < 0.0 || $lambda > 360.0} { + set lambda [expr {$lambda - 360. * floor($lambda / 360.)}] + } + set lambda [expr {$lambda - 180.}] + set x [expr {$lambda * $degree}] + set y [expr {($phi - $phi_0) * $degree}] + return [list $x $y] +} + +# ::mapproj::fromPlateCarree -- +# +# Solve a plate carrée projection for the +# latitude and longitude represented by a point on the map. +# +# Parameters: +# phi_0 -- Latitude of the center of the sheet in degrees +# lambda_0 -- Longitude of the center of projection +# x,y -- normalized x and y co-ordinates of a point on the map +# +# Results: +# Returns a list consisting of the longitude and latitude in degrees. + +proc ::mapproj::fromPlateCarree {phi_0 lambda_0 x y} { + variable radian + set lambda [expr {$lambda_0 + $x * $radian + 180.}] + if {$lambda < 0.0 || $lambda > 360.0} { + set lambda [expr {$lambda - 360. * floor($lambda / 360.)}] + } + set lambda [expr {$lambda - 180.}] + set phi [expr {$y * $radian + $phi_0}] + return [list $lambda $phi] +} + +# mapproj::toCylindricalEqualArea -- +# +# Project a latitude and longitude into cylindrical equal-area +# co-ordinates. +# +# Parameters: +# phi_1 -- Standard latitude in degrees +# phi_0 -- Latitude of the center of the sheet in degrees +# lambda_0 -- Longitude of the center of projection in degrees +# lambda -- Longitude of the point to be projected in degrees +# phi -- Latitude of the point to be projected in degrees +# +# Results: +# Returns x and y co-ordinates. Units of x and y are Earth radii; +# scale is true at the reference latitude + +proc ::mapproj::toCylindricalEqualArea {phi_1 lambda_0 phi_0 lambda phi} { + variable degree + set cos_phi_s [expr {cos($phi_1 * $degree)}] + set lambda [expr {$lambda - $lambda_0 + 180.}] + if {$lambda < 0.0 || $lambda > 360.0} { + set lambda [expr {$lambda - 360. * floor($lambda / 360.)}] + } + set lambda [expr {$lambda - 180.}] + set x [expr {$lambda * $degree * $cos_phi_s}] + set y0 [expr {sin($phi_0 * $degree) / $cos_phi_s}] + set y [expr {sin($phi * $degree) / $cos_phi_s}] + return [list $x [expr {$y - $y0}]] +} + +# ::mapproj::fromCylindricalEqualArea -- +# +# Solve a cylindrical equal area map projection for the +# latitude and longitude represented by a point on the map. +# +# Parameters: +# phi_1 -- Standard latitude in degrees +# phi_0 -- Latitude of the center of the sheet in degrees +# lambda_0 -- Longitude of the center of sheet in degrees +# x,y -- normalized x and y co-ordinates of a point on the map +# +# Results: +# Returns a list consisting of the longitude and latitude in degrees. + +proc ::mapproj::fromCylindricalEqualArea {phi_1 lambda_0 phi_0 x y} { + variable degree + variable radian + set cos_phi_s [expr {cos($phi_1 * $degree)}] + set lambda [expr {$lambda_0 + $x / $cos_phi_s * $radian + 180.}] + if {$lambda < 0.0 || $lambda > 360.0} { + set lambda [expr {$lambda - 360. * floor($lambda / 360.)}] + } + set lambda [expr {$lambda - 180.}] + set y0 [expr {sin($phi_0 * $degree) / $cos_phi_s}] + set phi [expr {asin(($y + $y0) * $cos_phi_s) * $radian}] + return [list $lambda $phi] +} + +# ::mapproj::toMercator -- +# +# Project a latitude and longitude into the Mercator projection +# co-ordinates. +# +# Parameters: +# lambda_0 -- Longitude of the center of projection in degrees +# phi_0 -- Latitude of the center of sheet in degrees +# lambda -- Longitude of the point to be projected in degrees +# phi -- Latitude of the point to be projected in degrees +# +# Results: +# Returns x and y co-ordinates in Earth radii. Scale is true +# at the Equator and increases without bounds toward the Poles. + +proc ::mapproj::toMercator {lambda_0 phi_0 lambda phi} { + variable trace; if {[info exists trace]} { puts [info level 0] } + variable degree + variable quarterpi + set lambda [expr {$lambda - $lambda_0 + 180.}] + if {$lambda < 0.0 || $lambda > 360.0} { + set lambda [expr {$lambda - 360. * floor($lambda / 360.)}] + } + set lambda [expr {$lambda - 180.}] + set x [expr {$lambda * $degree}] + set y [expr {log(tan($quarterpi + 0.5 * $phi * $degree))}] + set y0 [expr {log(tan($quarterpi + 0.5 * $phi_0 * $degree))}] + if {[info exists trace]} { puts "[info level 0] -> $x [expr {$y - $y0}]" } + return [list $x [expr {$y - $y0}]] +} + +# ::mapproj::fromMercator -- +# +# Converts Mercator map co-ordinates to latitude and longitude. +# +# Parameters: +# lambda_0 -- Longitude of the center of projection +# phi_0 -- Latitude of the center of sheet in degrees +# x,y -- normalized x and y co-ordinates of a point on the map +# +# Results: +# Returns a list consisting of the longitude and latitude in degrees. + +proc ::mapproj::fromMercator {lambda_0 phi_0 x y} { + variable quarterpi + variable degree + variable radian + set y0 [expr {log(tan($quarterpi + 0.5 * $phi_0 * $degree))}] + set lambda [expr {$lambda_0 + $x * $radian + 180.}] + if {$lambda < 0.0 || $lambda > 360.0} { + set lambda [expr {$lambda - 360. * floor($lambda / 360.)}] + } + set lambda [expr {$lambda - 180.}] + set phi [expr {$radian * 2.0 * atan(exp($y + $y0)) - 90.0}] + return [list $lambda $phi] +} + +# ::mapproj::toMillerCylindrical -- +# +# Project a latitude and longitude into the Miller Cylindrical projection +# co-ordinates. +# +# Parameters: +# lambda_0 -- Longitude of the center of projection in degrees +# lambda -- Longitude of the point to be projected in degrees +# phi -- Latitude of the point to be projected in degrees +# +# Results: +# Returns x and y co-ordinates. The x co-ordinate ranges from +# -$pi to pi. + +proc ::mapproj::toMillerCylindrical {lambda_0 lambda phi} { + foreach {x y} [toMercator $lambda_0 0.0 \ + $lambda [expr {0.8 * $phi}]] break + set y [expr {1.25 * $y}] + return [list $x $y] +} + +# ::mapproj::fromMillerCylindrical -- +# +# Converts Miller Cylindrical projected map co-ordinates +# to latitude and longitude. +# +# Parameters: +# lambda_0 -- Longitude of the center of projection +# x,y -- normalized x and y co-ordinates of a point on the map +# +# Results: +# Returns a list consisting of the longitude and latitude in degrees. + +proc ::mapproj::fromMillerCylindrical {lambda_0 x y} { + foreach {lambda phi} [fromMercator $lambda_0 0.0 \ + $x [expr {0.8 * $y}]] break + return [list $lambda [expr {1.25 * $phi}]] +} + +# ::mapproj::toSinusoidal -- +# +# Project a latitude and longitude into the sinusoidal +# (Sanson-Flamsteed) projection. +# +# Parameters: +# lambda_0 -- Longitude of the center of projection in degrees +# phi_0 -- Latitude of the center of the sheet, in degrees +# lambda -- Longitude of the point to be projected in degrees +# phi -- Latitude of the point to be projected in degrees +# +# Results: +# Returns x and y co-ordinates, in Earth radii. +# Scale is true along the Equator and central meridian. + +proc ::mapproj::toSinusoidal {lambda_0 phi_0 lambda phi} { + variable degree + variable quarterpi + set lambda [expr {$lambda - $lambda_0 + 180.}] + if {$lambda < 0.0 || $lambda > 360.0} { + set lambda [expr {$lambda - 360. * floor($lambda / 360.)}] + } + set lambda [expr {($lambda - 180.) * $degree}] + set phi [expr {$phi * $degree}] + set x [expr {$lambda * cos($phi)}] + set phi [expr {$phi - $phi_0 * $degree}] + return [list $x $phi] +} + +# ::mapproj::fromSinusoidal -- +# +# Converts sinusoidal (Sanson-Flamsteed) map co-ordinates +# to latitude and longitude. +# +# Parameters: +# lambda_0 -- Longitude of the center of projection +# phi_0 -- Latitude of the center of the sheet, in degrees +# x,y -- normalized x and y co-ordinates of a point on the map +# +# Results: +# Returns a list consisting of the longitude and latitude in degrees. + +proc ::mapproj::fromSinusoidal {lambda_0 phi_0 x y} { + variable degree + variable radian + set y [expr {$y + $phi_0 * $degree}] + set phi [expr {$y * $radian}] + set lambda [expr {180. + $lambda_0 + $radian * $x / cos($y)}] + if {$lambda < 0.0 || $lambda > 360.0} { + set lambda [expr {$lambda - 360. * floor($lambda / 360.)}] + } + set lambda [expr {$lambda - 180.}] + return [list $lambda $phi] +} + +# ::mapproj::toMollweide -- +# +# Project a latitude and longitude into the Mollweide projection +# co-ordinates. +# +# Parameters: +# lambda_0 -- Longitude of the center of projection in degrees +# lambda -- Longitude of the point to be projected in degrees +# phi -- Latitude of the point to be projected in degrees +# +# Results: +# Returns x and y co-ordinates, in Earth radii. +# Scale is true along the 40 deg 44 min parallels + +proc ::mapproj::toMollweide {lambda_0 lambda phi} { + variable degree + variable pi + variable halfpi + variable sqrt2 + variable sqrt8 + set lambda [expr {$lambda - $lambda_0 + 180.}] + if {$lambda < 0.0 || $lambda > 360.0} { + set lambda [expr {$lambda - 360. * floor($lambda / 360.)}] + } + set lambda [expr {($lambda - 180.) * $degree}] + set phi [expr {$phi * $degree}] + set theta [expr {2.0 * asin(2.0 * $phi / $pi)}] + set diff 1.0 + set pisinphi [expr {$pi * sin($phi)}] + while {abs($diff) >= 1.0e-4} { + set diff [expr {($theta + sin($theta) - $pisinphi) + / (1.0 + cos($theta))}] + set theta [expr {$theta - $diff}] + } + set theta [expr {0.5 * $theta}] + set x [expr {$sqrt8 * $lambda * cos($theta) / $pi}] + set y [expr {$sqrt2 * sin($theta)}] + return [list $x $y] +} + +# ::mapproj::fromMollweide -- +# +# Converts Mollweide projected map co-ordinates to latitude +# and longitude. +# +# Parameters: +# lambda_0 -- Longitude of the center of projection +# x,y -- normalized x and y co-ordinates of a point on the map +# +# Results: +# Returns a list consisting of the longitude and latitude in degrees. + +proc ::mapproj::fromMollweide {lambda_0 x y} { + variable pi + variable radian + variable degree + variable halfpi + variable sqrt2 + variable sqrt8 + set theta [expr {asin($y / $sqrt2)}] + set lambda [expr {$lambda_0 + $radian * $pi * $x / + ($sqrt8 * cos($theta)) + 180.}] + set phi [expr {asin((2.0 * $theta + sin(2.0 * $theta)) / $pi)}] + if {$lambda < 0.0 || $lambda > 360.0} { + set lambda [expr {$lambda - 360. * floor($lambda / 360.)}] + } + set lambda [expr {$lambda - 180.}] + return [list $lambda [expr {$phi * $radian}]] +} + +# ::mapproj::toEckertIV -- +# +# Project a latitude and longitude into the Eckert IV projection +# co-ordinates. +# +# Parameters: +# lambda_0 -- Longitude of the center of projection in degrees +# lambda -- Longitude of the point to be projected in degrees +# phi -- Latitude of the point to be projected in degrees +# +# Results: +# Returns x and y co-ordinates. Scale is true along the 40 deg 30 min +# parallels. + +proc ::mapproj::toEckertIV {lambda_0 lambda phi} { + variable degree + variable pi + variable halfpi + variable EckertIVK1 + variable EckertIVK2 + set lambda [expr {$lambda - $lambda_0 + 180.}] + if {$lambda < 0.0 || $lambda > 360.0} { + set lambda [expr {$lambda - 360. * floor($lambda / 360.)}] + } + set lambda [expr {($lambda - 180.) * $degree}] + set phi [expr {$phi * $degree}] + set theta [expr {$phi / 2}] + set diff 1.0 + set A [expr {(2.0 + $halfpi) * sin($phi)}] + while {abs($diff) >= 1.0e-4} { + set costheta [expr {cos($theta)}] + set sintheta [expr {sin($theta)}] + set diff \ + [expr {($theta + $sintheta * $costheta + 2.0 * sin($theta) - $A) + / (2.0 * $costheta * (1.0 + $costheta))}] + set theta [expr {$theta - $diff}] + } + set x [expr {$EckertIVK1 * $lambda * (1.0 + cos($theta))}] + set y [expr {$EckertIVK2 * sin($theta)}] + return [list $x $y] +} + +# ::mapproj::fromEckertIV -- +# +# Converts Eckert IV projected map co-ordinates to latitude +# and longitude. +# +# Parameters: +# lambda_0 -- Longitude of the center of projection +# x,y -- normalized x and y co-ordinates of a point on the map +# +# Results: +# Returns a list consisting of the longitude and latitude in degrees. + +proc ::mapproj::fromEckertIV {lambda_0 x y} { + variable pi + variable radian + variable degree + variable halfpi + variable sqrt2 + variable EckertIVK1 + variable EckertIVK2 + set sintheta [expr {$y / $EckertIVK2}] + set costheta [expr {sqrt(1.0 - $sintheta * $sintheta)}] + set theta [expr {atan2($sintheta, $costheta)}] + set phi [expr {asin(($theta + $sintheta*$costheta + 2.*$sintheta) + / (2. + $halfpi))}] + set lambda [expr {180.0 + $lambda_0 + + $radian / $EckertIVK1 * $x / (1.0 + $costheta)}] + if {$lambda < 0.0 || $lambda > 360.0} { + set lambda [expr {$lambda - 360. * floor($lambda / 360.)}] + } + set lambda [expr {$lambda - 180.}] + return [list $lambda [expr {$phi * $radian}]] +} + +# ::mapproj::toEckertVI -- +# +# Project a latitude and longitude into the Eckert IV projection +# co-ordinates. +# +# Parameters: +# lambda_0 -- Longitude of the center of projection in degrees +# lambda -- Longitude of the point to be projected in degrees +# phi -- Latitude of the point to be projected in degrees +# +# Results: +# Returns x and y co-ordinates. Scale is true along the 40 deg 30 min +# parallels. + +proc ::mapproj::toEckertVI {lambda_0 lambda phi} { + variable degree + variable pi + variable halfpi + variable EckertVIK1 + set lambda [expr {$lambda - $lambda_0 + 180.}] + if {$lambda < 0.0 || $lambda > 360.0} { + set lambda [expr {$lambda - 360. * floor($lambda / 360.)}] + } + set lambda [expr {($lambda - 180.) * $degree}] + set phi [expr {$phi * $degree}] + set theta [expr {$phi / 2}] + set diff 1.0 + set A [expr {(1.0 + $halfpi) * sin($phi)}] + while {abs($diff) >= 1.0e-4} { + set costheta [expr {cos($theta)}] + set sintheta [expr {sin($theta)}] + set diff \ + [expr {($theta + $sintheta - $A) + / (1.0 + $costheta)}] + set theta [expr {$theta - $diff}] + } + set x [expr {$lambda * (1.0 + cos($theta)) / $EckertVIK1}] + set y [expr {2.0 * $theta / $EckertVIK1}] + return [list $x $y] +} + +# ::mapproj::fromEckertVI -- +# +# Converts Eckert IV projected map co-ordinates to latitude +# and longitude. +# +# Parameters: +# lambda_0 -- Longitude of the center of projection +# x,y -- normalized x and y co-ordinates of a point on the map +# +# Results: +# Returns a list consisting of the longitude and latitude in degrees. + +proc ::mapproj::fromEckertVI {lambda_0 x y} { + variable pi + variable radian + variable degree + variable halfpi + variable sqrt2 + variable EckertVIK1 + puts [info level 0] + set theta [expr {0.5 * $EckertVIK1 * $y}] + puts [list theta = $theta] + set phi [expr {asin(($theta + sin($theta)) / (1.0 + $halfpi))}] + puts [list phi = $phi] + set lambda [expr {180.0 + $lambda_0 + $radian * $EckertVIK1 * $x + / (1 + cos($theta))}] + if {$lambda < 0.0 || $lambda > 360.0} { + set lambda [expr {$lambda - 360. * floor($lambda / 360.)}] + } + set lambda [expr {$lambda - 180.}] + return [list $lambda [expr {$phi * $radian}]] +} + +# ::mapproj::toRobinson -- +# +# Project a latitude and longitude into the Robinson projection. +# +# Parameters: +# lambda_0 -- Longitude of the center of projection in degrees +# lambda -- Longitude of the point to be projected in degrees +# phi -- Latitude of the point to be projected in degrees +# +# Results: +# Returns x and y co-ordinates, in Earth radii. +# Scale is true along the Equator. + +proc ::mapproj::toRobinson {lambda_0 lambda phi} { + variable RobinsonLatitude + variable RobinsonSplinePLEN + variable RobinsonSplinePDFE + variable RobinsonM + variable pi + variable degree + set lambda [expr {$lambda - $lambda_0 + 180.}] + if {$lambda < 0.0 || $lambda > 360.0} { + set lambda [expr {$lambda - 360. * floor($lambda / 360.)}] + } + set lambda [expr {$lambda - 180.}] + set y [math::interpolate::interp-cubic-splines $RobinsonSplinePDFE $phi] + set y [expr {$RobinsonM * $y}] + set s [math::interpolate::interp-cubic-splines $RobinsonSplinePLEN $phi] + set x [expr {$degree * $s * $lambda}] + return [list $x $y] +} + +# ::mapproj::fromRobinson -- +# +# Solve the Robinson projection for the +# latitude and longitude represented by a point on the map. +# +# Parameters: +# lambda_0 -- Longitude of the center of projection +# x,y -- normalized x and y co-ordinates of a point on the map +# +# Results: +# Returns a list consisting of the longitude and latitude in degrees. + +proc ::mapproj::fromRobinson {lambda_0 x y} { + variable RobinsonLatitude + variable RobinsonPDFE + variable RobinsonSplinePLEN + variable RobinsonSplinePDFE + variable RobinsonM + variable radian + + # We know that Robinson latitudes are equally spaced from [-90..90] + # at 5-degree intervals. Find the values for RobinsonPDFE that + # bracket the y co-ordinate. + + set y [expr {$y / $RobinsonM}] + set l 0 + set u [expr {[llength $RobinsonPDFE] - 1}] + while {$l < $u} { + set m [expr {($l + $u + 1) / 2}] + if {$y >= [lindex $RobinsonPDFE $m]} { + set l $m + } else { + set u [expr {$m - 1}] + } + } + set u [lindex $RobinsonLatitude [expr {$l+1}]] + set l [lindex $RobinsonLatitude $l] + for {set i 0} {$i < 12} {incr i} { + set m [expr {0.5 * ($u + $l)}] + set ystar [math::interpolate::interp-cubic-splines \ + $RobinsonSplinePDFE $m] + if {$ystar < $y} { + set l $m + } else { + set u $m + } + } + puts "latitude $m" + set s [math::interpolate::interp-cubic-splines $RobinsonSplinePLEN $m] + puts "parallel length $s" + set lambda [expr {180.0 + $lambda_0 + $radian * $x / $s}] + if {$lambda < 0.0 || $lambda > 360.0} { + set lambda [expr {$lambda - 360. * floor($lambda / 360.)}] + } + set lambda [expr {$lambda - 180.}] + return [list $lambda $m] +} + +# ::mapproj::toCassini -- +# +# Project a latitude and longitude into the Cassini projection. +# +# Parameters: +# lambda_0 -- Longitude of the center of projection in degrees +# phi_0 -- Latitude of the center of the sheet +# lambda -- Longitude of the point to be projected in degrees +# phi -- Latitude of the point to be projected in degrees +# +# Results: +# Returns x and y co-ordinates, in Earth radii. +# Scale is true along the central meridian. + +proc ::mapproj::toCassini {lambda_0 phi_0 lambda phi} { + variable degree + variable pi + variable twopi + variable quarterpi + set lambda [expr {$lambda - $lambda_0 + 180.}] + if {$lambda < 0.0 || $lambda > 360.0} { + set lambda [expr {$lambda - 360. * floor($lambda / 360.)}] + } + set lambda [expr {($lambda - 180.) * $degree}] + set phi [expr {$phi * $degree}] + set x [expr {asin(cos($phi) * sin($lambda))}] + set y [expr {atan2(tan($phi), cos($lambda)) - $degree * $phi_0}] + if {$y < -$pi} { + set y [expr {$y + $twopi}] + } elseif {$y > $pi} { + set y [expr {$y - $twopi}] + } + return [list $x $y] +} + +# ::mapproj::fromCassini -- +# +# Converts Cassini map co-ordinates to latitude and longitude. +# +# Parameters: +# lambda_0 -- Longitude of the center of projection +# phi_0 -- Latitude of the center of the sheet +# x,y -- normalized x and y co-ordinates of a point on the map +# +# Results: +# Returns a list consisting of the longitude and latitude in degrees. + +proc ::mapproj::fromCassini {lambda_0 phi_0 x y} { + variable degree + variable radian + set y [expr {$y + $degree * $phi_0}] + set phi [expr {$radian * asin(cos($x) * sin($y))}] + set lambda [expr {180. + $lambda_0 + $radian * atan2(tan($x), cos($y))}] + if {$lambda < 0.0 || $lambda > 360.0} { + set lambda [expr {$lambda - 360. * floor($lambda / 360.)}] + } + set lambda [expr {$lambda - 180.}] + return [list $lambda $phi] +} + +# ::mapproj::toPeirceQuincuncial +# +# Converts geodetic co-ordinates to the Peirce Quincuncial projection. +# +# Parameters: +# lambda_0 - Longitude of the central meridian. (Conventionally, 20.0). +# lambda - Longitude of the point to be projected in degrees +# phi - Latitude of the point to be projected in degrees. +# +# Results: +# Returns a list of the x and y co-ordinates. + +proc ::mapproj::toPeirceQuincuncial {lambda_0 lambda phi} { + variable degree + variable halfSqrt2 + variable pi + variable quarterpi + variable mquarterpi + variable threequarterpi + variable mthreequarterpi + variable PeirceQuincuncialScale + + # Convert latitude and longitude to radians relative to the + # central meridian + + set lambda [expr {$lambda - $lambda_0 + 180.}] + if {$lambda < 0.0 || $lambda > 360.0} { + set lambda [expr {$lambda - 360. * floor($lambda / 360.)}] + } + set lambda [expr {($lambda - 180.) * $degree}] + set phi [expr {$phi * $degree}] + + # Compute the auxiliary quantities 'm' and 'n'. Set 'm' to match + # the sign of 'lambda' and 'n' to be positive if |lambda| > pi/2 + + set cos_phiosqrt2 [expr {$halfSqrt2 * cos($phi)}] + set cos_lambda [expr {cos($lambda)}] + set sin_lambda [expr {sin($lambda)}] + set cos_a [expr {$cos_phiosqrt2 * ($sin_lambda + $cos_lambda)}] + set cos_b [expr {$cos_phiosqrt2 * ($sin_lambda - $cos_lambda)}] + set sin_a [expr {sqrt(1.0 - $cos_a * $cos_a)}] + set sin_b [expr {sqrt(1.0 - $cos_b * $cos_b)}] + set cos_a_cos_b [expr {$cos_a * $cos_b}] + set sin_a_sin_b [expr {$sin_a * $sin_b}] + set sin2_m [expr {1.0 + $cos_a_cos_b - $sin_a_sin_b}] + set sin2_n [expr {1.0 - $cos_a_cos_b - $sin_a_sin_b}] + if {$sin2_m < 0.0} {set sin2_m 0.0} + set sin_m [expr {sqrt($sin2_m)}] + if {$sin2_m > 1.0} { set sin2_m 1.0 } + set cos_m [expr {sqrt(1.0 - $sin2_m)}] + if {$sin_lambda < 0.0} { + set sin_m [expr {-$sin_m}] + } + if {$sin2_n < 0.0} { set sin2_n 0.0 } + set sin_n [expr {sqrt($sin2_n)}] + if {$sin2_n > 1.0} { set sin2_n 1.0 } + set cos_n [expr {sqrt(1.0 - $sin2_n)}] + if {$cos_lambda > 0.0} { + set sin_n [expr {-$sin_n}] + } + + # Compute elliptic integrals to map the disc to the square + + set x [ellFaux $cos_m $sin_m $halfSqrt2] + set y [ellFaux $cos_n $sin_n $halfSqrt2] + + # Reflect the Southern Hemisphere outward + + if {$phi < 0} { + if {$lambda < $mthreequarterpi} { + set y [expr {$PeirceQuincuncialScale - $y}] + } elseif {$lambda < $mquarterpi} { + set x [expr {-$PeirceQuincuncialScale - $x}] + } elseif {$lambda < $quarterpi} { + set y [expr {-$PeirceQuincuncialScale - $y}] + } elseif {$lambda < $threequarterpi} { + set x [expr {$PeirceQuincuncialScale - $x}] + } else { + set y [expr {$PeirceQuincuncialScale - $y}] + } + } + + # Rotate the square by 45 degrees to fit the screen better + + set X [expr {($x - $y) * $halfSqrt2}] + set Y [expr {($x + $y) * $halfSqrt2}] + + return [list $X $Y] +} + +# ::mapproj::fromPeirceQuincuncial -- +# +# Converts Peirce Quincuncial map co-ordinates to latitude and longitude. +# +# Parameters: +# lambda_0 -- Longitude of the center of projection +# x,y -- normalized x and y co-ordinates of a point on the map +# +# Results: +# Returns a list consisting of the longitude and latitude in degrees. + +proc ::mapproj::fromPeirceQuincuncial {lambda_0 x y} { + variable halfSqrt2 + variable radian + variable pi + variable halfpi + variable quarterpi + variable PeirceQuincuncialScale + variable PeirceQuincuncialLimit + + # Rotate x and y 45 degrees + + set X [expr {($x + $y) * $halfSqrt2}] + set Y [expr {($y - $x) * $halfSqrt2}] + + # Reflect Southern Hemisphere into the Northern + + set southern 0 + if {$X < -$PeirceQuincuncialLimit} { + set X [expr {-$PeirceQuincuncialScale - $X}] + set southern 1 + } elseif {$X > $PeirceQuincuncialLimit} { + set X [expr {$PeirceQuincuncialScale - $X}] + set southern 1 + } elseif {$Y < -$PeirceQuincuncialLimit} { + set Y [expr {-$PeirceQuincuncialScale - $Y}] + set southern 1 + } elseif {$Y > $PeirceQuincuncialLimit} { + set Y [expr {$PeirceQuincuncialScale - $Y}] + set southern 1 + } + + # Now we know that latitude will be positive. If X is negative, then + # longitude will be negative; reflect the Western Hemisphere into the + # Eastern. + + set western 0 + if {$X < 0.0} { + set western 1 + set X [expr {-$X}] + } + + # If Y is positive, the point is in the back hemisphere. Reflect + # it to the front. + + set back 0 + if {$Y > 0.0} { + set back 1 + set Y [expr {-$Y}] + } + + # Finally, constrain longitude to be less than pi/4, by reflecting across + # the 45 degree meridian. + + set complement 0 + if {$X > -$Y} { + set complement 1 + set t [expr {-$X}] + set X [expr {-$Y}] + set Y $t + } + + # Compute the elliptic functions to map the plane onto the sphere + + set cnx [cn $X $halfSqrt2] + set cny [cn $Y $halfSqrt2] + + # Undo the mapping to latitude and longitude + + set a1 [expr {acos(-$cnx * $cnx)}] + set a2 [expr {acos($cny * $cny)}] + set b [expr {0.5 * ($a1 + $a2)}] + set a [expr {0.5 * ($a1 - $a2)}] + set cos_a [expr {cos($a)}] + set cos_b [expr {-cos($b)}] + set lambda [expr {$quarterpi - atan2($cos_b, $cos_a)}] + set phi [expr {acos(hypot($cos_b, $cos_a))}] + + # Undo the reflections that were done above, to get correct latitude + # and longitude + + if {$complement} { + set lambda [expr {$halfpi - $lambda}] + } + if {$back} { + set lambda [expr {$pi - $lambda}] + } + if {$western} { + set lambda [expr {-$lambda}] + } + if {$southern} { + set phi [expr {-$phi}] + } + + # Convert latitude and longitude to degrees + + set lambda [expr {$lambda * $radian + 180. + $lambda_0}] + if {$lambda < 0.0 || $lambda > 360.0} { + set lambda [expr {$lambda - 360. * floor($lambda / 360.)}] + } + set lambda [expr {$lambda - 180.}] + set phi [expr {$phi * $radian}] + + return [list $lambda $phi] +} + +# ::mapproj::RotateCartesianY +# +# Rotates Cartesian co-ordinates about the y axis +# +# Parameters: +# phi - Angle (in degrees) about which to rotate. +# x,y,z - Cartesian co-ordinates +# +# Results: +# Returns a three-element list giving the rotated co-ordinates. + +proc ::mapproj::RotateCartesianY {phi x y z} { + variable degree + set phi [expr {$degree * $phi}] + set cos_phi [expr {cos($phi)}] + set sin_phi [expr {sin($phi)}] + return [list [expr {$x * $cos_phi - $z * $sin_phi}] \ + $y \ + [expr {$z * $cos_phi + $x * $sin_phi}]] +} + +# ::mapproj::ToCartesian -- +# +# Converts geodetic co-ordinates to Cartesian +# +# Parameters: +# lambda - Longitude of the point to be projected, in degrees +# phi - Latitude of the point to be projected, in degrees +# +# Results: +# Returns a three-element list, x, y, z where x is the component +# in the direction of longitude 0, latitude 0, y is the component +# in the direction of longitude 90 East, latitude 0, and +# z is the component in the direction of the North Pole +# +# Auxiliary procedure used in several projections to convert +# geodetic coordinates to Cartesian range and bearing. + +proc ::mapproj::ToCartesian {lambda phi} { + variable degree + set lambda [expr {$degree * $lambda}] + set phi [expr {$degree * $phi}] + set cos_phi [expr cos($phi)] + return [list [expr {$cos_phi * cos($lambda)}] \ + [expr {$cos_phi * sin($lambda)}] \ + [expr {sin($phi)}]] +} + +# ::mapproj::CartesianToRangeAndBearing +# +# Transforms view-relative Cartesian co-ordinates to range and +# bearing. +# +# Parameters: +# x,y,z - Cartesian co-ordinates relative to center of Earth; +# +x points to the viewer and +z to the "view-up" direction. +# +# Results: +# Returns a three-element list containing, in order, +# the cosine (easting) of the bearing, the sine (northing) of the +# bearing, and the range. + +proc ::mapproj::CartesianToRangeAndBearing {x y z} { + set c [expr {hypot($z, $y)}] + if {$c == 0} { + set cos_b 1.0 + set sin_b 0.0 + } else { + set cos_b [expr {$y / $c}] + set sin_b [expr {$z / $c}] + } + set range [expr {atan2($c, $x)}] + return [list $cos_b $sin_b $range] +} + +# ::mapproj::RangeAndBearingToCartesian -- +# +# Converts range and bearing to Cartesian co-ordinates. +# +# Parameters: +# cos_b, sin_b -- Cosine (easting) and sine (northing) of the bearing +# range - Range, in Earth radii +# +# Results: +# Returns Cartesian co-ordinates relative to center of Earth. +# x is toward the station, and z is "view up" + +proc ::mapproj::RangeAndBearingToCartesian {cos_b sin_b range} { + set c [expr {sin($range)}] + set x [expr {cos($range)}] + set y [expr {$cos_b * $c}] + set z [expr {$sin_b * $c}] + return [list $x $y $z] +} + + +# ::mapproj::CartesianToSpherical -- +# +# Transforms Cartesian x, y, z to spherical co-ordinates +# +# Parameters: +# x, y, z -- Coordinates of a point on the surface of the Earth, +# in Earth radii +# +# Results: +# Returns a two-element list comprising longitude and latitude +# in radians + +proc ::mapproj::CartesianToSpherical {x y z} { + return [list [expr {atan2($y, $x)}] [expr {atan2($z, hypot($y, $x))}]] +} + +# ::mapproj::toOrthographic -- +# +# Transforms latitude and longitude to x and y co-ordinates +# on an orthographic projection. Scale is true only at the +# point of projection. +# +# Parameters: +# lambda_0, phi_0 -- Longitude and latitude of the center of projection +# in degrees +# lambda, phi -- Longitude and latitude of the point to be projected +# in degrees +# +# Results: +# Returns map x and y co-ordinates, in Earth radii. + +proc ::mapproj::toOrthographic {lambda_0 phi_0 lambda phi} { + foreach {x y z} [ToCartesian [expr {$lambda-$lambda_0}] $phi] \ + break + foreach {x y z} [RotateCartesianY [expr {-$phi_0}] $x $y $z] \ + break + if {$x < 0} { + return {} + } else { + return [list $y $z] + } +} + +# ::mapproj::fromOrthographic -- +# +# Transforms x and y on an orthographic projection to latitude +# and longitude. +# +# Parameters: +# lambda_0, phi_0 -- Longitude and latitude of the center of projection +# in degrees +# x, y -- Co-ordinates of the projected point, in Earth radii +# +# Results: +# Returns a two element list containing longitude and latitude +# in degrees. + +proc ::mapproj::fromOrthographic {lambda_0 phi_0 x y} { + variable radian + set r [expr {hypot($x, $y)}] + set alpha [expr {asin($r)}] + set z [expr {sqrt(1.0 - $r*$r)}] + foreach {x y z} [RotateCartesianY $phi_0 $z $x $y] break + foreach {lambda phi} [CartesianToSpherical $x $y $z] break + + # Convert latitude and longitude to degrees + + set lambda [expr {$lambda * $radian + 180. + $lambda_0}] + if {$lambda < 0.0 || $lambda > 360.0} { + set lambda [expr {$lambda - 360. * floor($lambda / 360.)}] + } + set lambda [expr {$lambda - 180.}] + set phi [expr {$phi * $radian}] + + return [list $lambda $phi] +} + +# ::mapproj::toStereographic -- +# +# Transforms latitude and longitude to x and y co-ordinates +# on an orthographic projection. Scale is true only at the +# point of projection. +# +# Parameters: +# lambda_0, phi_0 -- Longitude and latitude of the center of projection +# in degrees +# lambda, phi -- Longitude and latitude of the point to be projected +# in degrees +# +# Results: +# Returns map x and y co-ordinates, in Earth radii. + +proc ::mapproj::toStereographic {lambda_0 phi_0 lambda phi} { + foreach {x y z} [ToCartesian [expr {$lambda-$lambda_0}] $phi] \ + break + foreach {x y z} [RotateCartesianY [expr {-$phi_0}] $x $y $z] \ + break + if {$x < -0.5} { + return {} + } else { + set y [expr {2. * $y / (1. + $x)}] + set z [expr {2. * $z / (1. + $x)}] + return [list $y $z] + } +} + +# ::mapproj::fromStereographic -- +# +# Transforms x and y on an orthographic projection to latitude +# and longitude. +# +# Parameters: +# lambda_0, phi_0 -- Longitude and latitude of the center of projection +# in degrees +# x, y -- Co-ordinates of the projected point, in Earth radii +# +# Results: +# Returns a two element list containing longitude and latitude +# in degrees. + +proc ::mapproj::fromStereographic {lambda_0 phi_0 x y} { + variable radian + variable halfpi + set denom [expr {4.0 + $x*$x + $y*$y}] + foreach {x y z} [list \ + [expr {(4.0 - $x*$x - $y*$y) / $denom}] \ + [expr {4. * $x / $denom}] \ + [expr {4. * $y / $denom}]] break + + foreach {x y z} [RotateCartesianY $phi_0 $x $y $z] break + foreach {lambda phi} [CartesianToSpherical $x $y $z] break + + # Convert latitude and longitude to degrees + + set lambda [expr {$lambda * $radian + 180. + $lambda_0}] + if {$lambda < 0.0 || $lambda > 360.0} { + set lambda [expr {$lambda - 360. * floor($lambda / 360.)}] + } + set lambda [expr {$lambda - 180.}] + set phi [expr {$phi * $radian}] + + return [list $lambda $phi] +} + +# ::mapproj::toGnomonic -- +# +# Transforms latitude and longitude to x and y co-ordinates +# on an orthographic projection. Scale is true only at the +# point of projection. +# +# Parameters: +# lambda_0, phi_0 -- Longitude and latitude of the center of projection +# in degrees +# lambda, phi -- Longitude and latitude of the point to be projected +# in degrees +# +# Results: +# Returns map x and y co-ordinates, in Earth radii. + +proc ::mapproj::toGnomonic {lambda_0 phi_0 lambda phi} { + foreach {x y z} [ToCartesian [expr {$lambda-$lambda_0}] $phi] \ + break + foreach {x y z} [RotateCartesianY [expr {-$phi_0}] $x $y $z] \ + break + if {$x < 0.01} { + return {} + } else { + set y [expr {$y / $x}] + set z [expr {$z / $x}] + return [list $y $z] + } +} + +# ::mapproj::fromGnomonic -- +# +# Transforms x and y on an orthographic projection to latitude +# and longitude. +# +# Parameters: +# lambda_0, phi_0 -- Longitude and latitude of the center of projection +# in degrees +# x, y -- Co-ordinates of the projected point, in Earth radii +# +# Results: +# Returns a two element list containing longitude and latitude +# in degrees. + +proc ::mapproj::fromGnomonic {lambda_0 phi_0 x y} { + variable radian + variable halfpi + set denom [expr {hypot(1.0, hypot($x, $y))}] + foreach {x y z} [list \ + [expr {1.0 / $denom}] \ + [expr {$x / $denom}] \ + [expr {$y / $denom}]] break + + foreach {x y z} [RotateCartesianY $phi_0 $x $y $z] break + foreach {lambda phi} [CartesianToSpherical $x $y $z] break + + # Convert latitude and longitude to degrees + + set lambda [expr {$lambda * $radian + 180. + $lambda_0}] + if {$lambda < 0.0 || $lambda > 360.0} { + set lambda [expr {$lambda - 360. * floor($lambda / 360.)}] + } + set lambda [expr {$lambda - 180.}] + set phi [expr {$phi * $radian}] + + return [list $lambda $phi] +} + +# ::mapproj::toAzimuthalEquidistant -- +# +# Transforms latitude and longitude to x and y co-ordinates +# on an orthographic projection. Scale is true only at the +# point of projection. +# +# Parameters: +# lambda_0, phi_0 -- Longitude and latitude of the center of projection +# in degrees +# lambda, phi -- Longitude and latitude of the point to be projected +# in degrees +# +# Results: +# Returns map x and y co-ordinates, in Earth radii. + +proc ::mapproj::toAzimuthalEquidistant {lambda_0 phi_0 lambda phi} { + foreach {x y z} [ToCartesian [expr {$lambda-$lambda_0}] $phi] \ + break + foreach {x y z} [RotateCartesianY [expr {-$phi_0}] $x $y $z] \ + break + foreach {cs sn range} [CartesianToRangeAndBearing $x $y $z] break + return [list [expr {$cs * $range}] [expr {$sn * $range}]] +} + +# ::mapproj::fromAzimuthalEquidistant -- +# +# Transforms x and y on an orthographic projection to latitude +# and longitude. +# +# Parameters: +# lambda_0, phi_0 -- Longitude and latitude of the center of projection +# in degrees +# x, y -- Co-ordinates of the projected point, in Earth radii +# +# Results: +# Returns a two element list containing longitude and latitude +# in degrees. + +proc ::mapproj::fromAzimuthalEquidistant {lambda_0 phi_0 x y} { + variable radian + variable halfpi + + set range [expr {hypot($y, $x)}] + set cos_b [expr {$x / $range}] + set sin_b [expr {$y / $range}] + foreach {x y z} [RangeAndBearingToCartesian $cos_b $sin_b $range] break + foreach {x y z} [RotateCartesianY $phi_0 $x $y $z] break + foreach {lambda phi} [CartesianToSpherical $x $y $z] break + + # Convert latitude and longitude to degrees + + set lambda [expr {$lambda * $radian + 180. + $lambda_0}] + if {$lambda < 0.0 || $lambda > 360.0} { + set lambda [expr {$lambda - 360. * floor($lambda / 360.)}] + } + set lambda [expr {$lambda - 180.}] + set phi [expr {$phi * $radian}] + + return [list $lambda $phi] +} + +# ::mapproj::toLambertAzimuthalEqualArea -- +# +# Transforms latitude and longitude to x and y co-ordinates +# on an orthographic projection. Scale is true only at the +# point of projection. +# +# Parameters: +# lambda_0, phi_0 -- Longitude and latitude of the center of projection +# in degrees +# lambda, phi -- Longitude and latitude of the point to be projected +# in degrees +# +# Results: +# Returns map x and y co-ordinates, in Earth radii. + +proc ::mapproj::toLambertAzimuthalEqualArea {lambda_0 phi_0 lambda phi} { + foreach {x y z} [ToCartesian [expr {$lambda-$lambda_0}] $phi] \ + break + foreach {x y z} [RotateCartesianY [expr {-$phi_0}] $x $y $z] \ + break + foreach {cs sn range} [CartesianToRangeAndBearing $x $y $z] break + set range [expr {2.0 * sin(0.5 * $range)}] + return [list [expr {$cs * $range}] [expr {$sn * $range}]] +} + +# ::mapproj::fromLambertAzimuthalEqualArea -- +# +# Transforms x and y on an orthographic projection to latitude +# and longitude. +# +# Parameters: +# lambda_0, phi_0 -- Longitude and latitude of the center of projection +# in degrees +# x, y -- Co-ordinates of the projected point, in Earth radii +# +# Results: +# Returns a two element list containing longitude and latitude +# in degrees. + +proc ::mapproj::fromLambertAzimuthalEqualArea {lambda_0 phi_0 x y} { + variable radian + variable halfpi + + set range [expr {hypot($y, $x)}] + set cos_b [expr {$x / $range}] + set sin_b [expr {$y / $range}] + set range [expr {2.0 * asin(0.5 * $range)}] + foreach {x y z} [RangeAndBearingToCartesian $cos_b $sin_b $range] break + foreach {x y z} [RotateCartesianY $phi_0 $x $y $z] break + foreach {lambda phi} [CartesianToSpherical $x $y $z] break + + # Convert latitude and longitude to degrees + + set lambda [expr {$lambda * $radian + 180. + $lambda_0}] + if {$lambda < 0.0 || $lambda > 360.0} { + set lambda [expr {$lambda - 360. * floor($lambda / 360.)}] + } + set lambda [expr {$lambda - 180.}] + set phi [expr {$phi * $radian}] + + return [list $lambda $phi] +} + +# ::mapproj::toHammer -- +# +# Transforms latitude and longitude to x and y co-ordinates +# on an orthographic projection. Scale is true only at the +# point of projection. +# +# Parameters: +# lambda_0-- Longitude of the center of projection +# in degrees +# lambda, phi -- Longitude and latitude of the point to be projected +# in degrees +# +# Results: +# Returns map x and y co-ordinates, in Earth radii. + +proc ::mapproj::toHammer {lambda_0 lambda phi} { + set lambda [expr {$lambda - $lambda_0 + 180.}] + if {$lambda < 0.0 || $lambda > 360.0} { + set lambda [expr {$lambda - 360. * floor($lambda / 360.)}] + } + set lambda [expr {$lambda - 180.0}] + foreach {x y z} [ToCartesian [expr {$lambda/2.}] $phi] \ + break + foreach {cs sn range} [CartesianToRangeAndBearing $x $y $z] break + set range [expr {2.0 * sin(0.5 * $range)}] + return [list [expr {2.0 * $cs * $range}] [expr {$sn * $range}]] +} + +# ::mapproj::fromHammer -- +# +# Transforms x and y on an orthographic projection to latitude +# and longitude. +# +# Parameters: +# lambda_0, phi_0 -- Longitude and latitude of the center of projection +# in degrees +# x, y -- Co-ordinates of the projected point, in Earth radii +# +# Results: +# Returns a two element list containing longitude and latitude +# in degrees. + +proc ::mapproj::fromHammer {lambda_0 x y} { + variable radian + variable halfpi + + set x [expr {0.5 * $x}] + set range [expr {hypot($y, $x)}] + set cos_b [expr {$x / $range}] + set sin_b [expr {$y / $range}] + set range [expr {2.0 * asin(0.5 * $range)}] + foreach {x y z} [RangeAndBearingToCartesian $cos_b $sin_b $range] break + foreach {lambda phi} [CartesianToSpherical $x $y $z] break + + # Convert latitude and longitude to degrees + + set lambda [expr {2.0 * $lambda * $radian + 180. + $lambda_0}] + if {$lambda < 0.0 || $lambda > 360.0} { + set lambda [expr {$lambda - 360. * floor($lambda / 360.)}] + } + set lambda [expr {$lambda - 180.}] + set phi [expr {$phi * $radian}] + + return [list $lambda $phi] +} + +# ::mapproj::toConicEquidistant +# +# Converts latitude and longitude to map co-ordinates on a +# conic equidistant projection. +# +# Parameters: +# lambda_0, phi_0 -- Longitude and latitude of the center of the sheet. +# phi_1, phi_2 -- Latitudes of the two standard parallels at which scale +# is true +# lambda, phi -- Longitude and latitude of the point to be projected +# +# Results: +# Returns a list of map x and y measured in Earth radii. + +proc ::mapproj::toConicEquidistant {lambda_0 phi_0 phi_1 phi_2 lambda phi} { + variable degree + set phi_0 [expr {$phi_0 * $degree}] + set phi_1 [expr {$phi_1 * $degree}] + set phi_2 [expr {$phi_2 * $degree}] + set lambda [expr {$lambda - $lambda_0 + 180.}] + if {$lambda < 0.0 || $lambda > 360.0} { + set lambda [expr {$lambda - 360. * floor($lambda / 360.)}] + } + set lambda [expr {($lambda - 180.0) * $degree}] + set phi [expr {$phi * $degree}] + set cos_phi_1 [expr {cos($phi_1)}] + set n [expr {($cos_phi_1 - cos($phi_2)) / ($phi_2 - $phi_1)}] + set G [expr {$cos_phi_1 / $n + $phi_1}] + set rho_0 [expr {$G - $phi_0}] + set theta [expr {$n * $lambda}] + set rho [expr {$G - $phi}] + set x [expr {$rho * sin($theta)}] + set y [expr {$rho_0 - $rho * cos($theta)}] + return [list $x $y] +} + +# ::mapproj::fromConicEquidistant -- +# +# Unprojects map x and y in a conic equidistant projection to +# latitude and longitude +# +# Parameters: +# lambda_0, phi_0 -- Longitude and latitude of the center of the sheet. +# phi_1, phi_2 -- Latitudes of the two standard parallels at which scale +# is true +# x, y -- Map co-ordinates in Earth radii +# +# Results: +# Returns a list of longitude and latitude in degrees. + +proc ::mapproj::fromConicEquidistant {lambda_0 phi_0 phi_1 phi_2 x y} { + variable degree + variable radian + set phi_0 [expr {$phi_0 * $degree}] + set phi_1 [expr {$phi_1 * $degree}] + set phi_2 [expr {$phi_2 * $degree}] + set cos_phi_1 [expr {cos($phi_1)}] + set n [expr {($cos_phi_1 - cos($phi_2)) / ($phi_2 - $phi_1)}] + set G [expr {$cos_phi_1 / $n + $phi_1}] + set rho_0 [expr {$G - $phi_0}] + set rho_0my [expr {$rho_0 - $y}] + set theta [expr {atan2($x, $rho_0my)}] + set rho [expr {sqrt($x*$x + $rho_0my * $rho_0my)}] + if {$n < 0.0} {set rho [expr {-$rho}]} + set phi [expr {($G - $rho) * $radian}] + set lambda [expr {($theta / $n * $radian) + $lambda_0 + 180.}] + if {$lambda < 0.0 || $lambda > 360.0} { + set lambda [expr {$lambda - 360. * floor($lambda / 360.)}] + } + set lambda [expr {$lambda - 180.}] + return [list $lambda $phi] +} + +# ::mapproj::toAlbersEqualAreaConic +# +# Converts latitude and longitude to map co-ordinates on a +# conic equal-area projection. +# +# Parameters: +# lambda_0, phi_0 -- Longitude and latitude of the center of the sheet. +# phi_1, phi_2 -- Latitudes of the two standard parallels at which scale +# is true +# lambda, phi -- Longitude and latitude of the point to be projected +# +# Results: +# Returns a list of map x and y measured in Earth radii. + +proc ::mapproj::toAlbersEqualAreaConic {lambda_0 phi_0 phi_1 phi_2 lambda phi} { + variable degree + set phi_0 [expr {$phi_0 * $degree}] + set phi_1 [expr {$phi_1 * $degree}] + set phi_2 [expr {$phi_2 * $degree}] + set lambda [expr {$lambda - $lambda_0 + 180.}] + if {$lambda < 0.0 || $lambda > 360.0} { + set lambda [expr {$lambda - 360. * floor($lambda / 360.)}] + } + set lambda [expr {($lambda - 180.0) * $degree}] + set phi [expr {$phi * $degree}] + set cos_phi_1 [expr {cos($phi_1)}] + set sin_phi_1 [expr {sin($phi_1)}] + set n [expr {0.5 * ($sin_phi_1 + sin($phi_2))}] + set theta [expr {$n * $lambda}] + set C [expr {$cos_phi_1 * $cos_phi_1 + 2.0 * $n * $sin_phi_1}] + set rho [expr {sqrt($C - 2.0 * $n * sin($phi)) / $n}] + set rho_0 [expr {sqrt($C - 2.0 * $n * sin($phi_0)) / $n}] + set x [expr {$rho * sin($theta)}] + set y [expr {$rho_0 - $rho * cos($theta)}] + return [list $x $y] +} + +# ::mapproj::fromAlbersEqualAreaConic -- +# +# Unprojects map x and y in a conic equal-area projection to +# latitude and longitude +# +# Parameters: +# lambda_0, phi_0 -- Longitude and latitude of the center of the sheet. +# phi_1, phi_2 -- Latitudes of the two standard parallels at which scale +# is true +# x, y -- Map co-ordinates in Earth radii +# +# Results: +# Returns a list of longitude and latitude in degrees. + +proc ::mapproj::fromAlbersEqualAreaConic {lambda_0 phi_0 phi_1 phi_2 x y} { + variable degree + variable radian + variable twopi + set phi_0 [expr {$phi_0 * $degree}] + set phi_1 [expr {$phi_1 * $degree}] + set phi_2 [expr {$phi_2 * $degree}] + set cos_phi_1 [expr {cos($phi_1)}] + set sin_phi_1 [expr {sin($phi_1)}] + set n [expr {0.5 * ($sin_phi_1 + sin($phi_2))}] + set C [expr {$cos_phi_1 * $cos_phi_1 + 2.0 * $n * $sin_phi_1}] + set rho_0 [expr {sqrt($C - 2.0 * $n * sin($phi_0)) / $n}] + set theta [expr {atan2($x, $rho_0 - $y)}] + set rho [expr {hypot($x, $rho_0 - $y)}] + set phi [expr {$radian * asin(($C - $rho*$rho*$n*$n) / (2.0 * $n))}] + set lambda [expr {($theta / $n * $radian) + $lambda_0 + 180.}] + if {$lambda < 0.0 || $lambda > 360.0} { + set lambda [expr {$lambda - 360. * floor($lambda / 360.)}] + } + set lambda [expr {$lambda - 180.}] + return [list $lambda $phi] +} + +# ::mapproj::toLambertConformalConic +# +# Converts latitude and longitude to map co-ordinates on a +# conformal conic projection. +# +# Parameters: +# lambda_0, phi_0 -- Longitude and latitude of the center of the sheet. +# phi_1, phi_2 -- Latitudes of the two standard parallels at which scale +# is true +# lambda, phi -- Longitude and latitude of the point to be projected +# +# Results: +# Returns a list of map x and y measured in Earth radii. + +proc ::mapproj::toLambertConformalConic {lambda_0 phi_0 phi_1 phi_2 lambda phi} { + variable degree + variable quarterpi + set phi_0 [expr {$phi_0 * $degree}] + set phi_1 [expr {$phi_1 * $degree}] + set phi_2 [expr {$phi_2 * $degree}] + set lambda [expr {$lambda - $lambda_0 + 180.}] + if {$lambda < 0.0 || $lambda > 360.0} { + set lambda [expr {$lambda - 360. * floor($lambda / 360.)}] + } + set lambda [expr {($lambda - 180.0) * $degree}] + set phi [expr {$phi * $degree}] + set cos_phi_1 [expr {cos($phi_1)}] + set sin_phi_1 [expr {sin($phi_1)}] + set tan1 [expr {tan($quarterpi + 0.5 * $phi_1)}] + set n [expr {log($cos_phi_1 / cos($phi_2)) + / log(tan($quarterpi + 0.5 * $phi_2) / $tan1)}] + set F [expr {$cos_phi_1 * pow($tan1, $n) / $n}] + set rho [expr {$F * pow(tan($quarterpi + 0.5 * $phi), -$n)}] + set rho_0 [expr {$F * pow(tan($quarterpi + 0.5 * $phi_0), -$n)}] + set x [expr {$rho * sin($n * $lambda)}] + set y [expr {$rho_0 - $rho * cos($n * $lambda)}] + return [list $x $y] +} + +# ::mapproj::fromLambertConformalConic -- +# +# Unprojects map x and y in a conformal conic projection to +# latitude and longitude +# +# Parameters: +# lambda_0, phi_0 -- Longitude and latitude of the center of the sheet. +# phi_1, phi_2 -- Latitudes of the two standard parallels at which scale +# is true +# x, y -- Map co-ordinates in Earth radii +# +# Results: +# Returns a list of longitude and latitude in degrees. + +proc ::mapproj::fromLambertConformalConic {lambda_0 phi_0 phi_1 phi_2 x y} { + variable degree + variable radian + variable quarterpi + set phi_0 [expr {$phi_0 * $degree}] + set phi_1 [expr {$phi_1 * $degree}] + set phi_2 [expr {$phi_2 * $degree}] + set cos_phi_1 [expr {cos($phi_1)}] + set sin_phi_1 [expr {sin($phi_1)}] + set tan1 [expr {tan($quarterpi + 0.5 * $phi_1)}] + set n [expr {log($cos_phi_1 / cos($phi_2)) + / log(tan($quarterpi + 0.5 * $phi_2) / $tan1)}] + set F [expr {$cos_phi_1 * pow($tan1, $n) / $n}] + set rho_0 [expr {$F * pow(tan($quarterpi + 0.5 * $phi_0), -$n)}] + set y [expr {$rho_0 - $y}] + set rho [expr {sqrt($x*$x + $y*$y)}] + if {$n < 0} { set rho [expr {-$rho}] } + set theta [expr {atan2($x, $y)}] + set phi [expr {$radian * 2 * atan(pow($F / $rho, 1.0 / $n)) - 90.}] + set lambda [expr {($theta / $n * $radian) + $lambda_0 + 180.}] + if {$lambda < 0.0 || $lambda > 360.0} { + set lambda [expr {$lambda - 360. * floor($lambda / 360.)}] + } + set lambda [expr {$lambda - 180.}] + return [list $lambda $phi] +} + +# Define commonly used cylindrical equal-area projections + +proc ::mapproj::toLambertCylindricalEqualArea {lambda_0 phi_0 lambda phi} { + toCylindricalEqualArea 0.0 $lambda_0 $phi_0 $lambda $phi +} +proc ::mapproj::fromLambertCylindricalEqualArea {lambda_0 phi_0 x y} { + fromCylindricalEqualArea 0.0 $lambda_0 $phi_0 $x $y +} +proc ::mapproj::toBehrmann {lambda_0 phi_0 lambda phi} { + toCylindricalEqualArea 30.0 $lambda_0 $phi_0 $lambda $phi +} +proc ::mapproj::fromBehrmann {lambda_0 phi_0 x y} { + fromCylindricalEqualArea 30.0 $lambda_0 $phi_0 $x $y +} +proc ::mapproj::toTrystanEdwards {lambda_0 phi_0 lambda phi} { + toCylindricalEqualArea 37.4 $lambda_0 $phi_0 $lambda $phi +} +proc ::mapproj::fromTrystanEdwards {lambda_0 phi_0 x y} { + fromCylindricalEqualArea 37.4 $lambda_0 $phi_0 $x $y +} +proc ::mapproj::toHoboDyer {lambda_0 phi_0 lambda phi} { + toCylindricalEqualArea 37.5 $lambda_0 $phi_0 $lambda $phi +} +proc ::mapproj::fromHoboDyer {lambda_0 phi_0 x y} { + fromCylindricalEqualArea 37.5 $lambda_0 $phi_0 $x $y +} +proc ::mapproj::toGallPeters {lambda_0 phi_0 lambda phi} { + toCylindricalEqualArea 45.0 $lambda_0 $phi_0 $lambda $phi +} +proc ::mapproj::fromGallPeters {lambda_0 phi_0 x y} { + fromCylindricalEqualArea 45.0 $lambda_0 $phi_0 $x $y +} +proc ::mapproj::toBalthasart {lambda_0 phi_0 lambda phi} { + toCylindricalEqualArea 50.0 $lambda_0 $phi_0 $lambda $phi +} +proc ::mapproj::fromBalthasart {lambda_0 phi_0 x y} { + fromCylindricalEqualArea 50.0 $lambda_0 $phi_0 $x $y +} diff --git a/tcllib/modules/mapproj/pkgIndex.tcl b/tcllib/modules/mapproj/pkgIndex.tcl new file mode 100755 index 0000000..c488e8f --- /dev/null +++ b/tcllib/modules/mapproj/pkgIndex.tcl @@ -0,0 +1,2 @@ +if {![package vsatisfies [package provide Tcl] 8.4]} {return} +package ifneeded mapproj 1.0 [list source [file join $dir mapproj.tcl]] |