Guest User

Untitled

a guest
Jun 16th, 2012
223
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. #!/usr/bin/perl
  2.  
  3. use warnings;
  4. use strict;
  5.  
  6. package LEST97;
  7.  
  8. # A package to convert L-EST 97 coordinates to WSG84 lat / lon
  9. # and vice versa. See [1] for algorithm details and [2] for
  10. # specifics on L-EST 97.
  11. #
  12. # [1] http://www.linz.govt.nz/geodetic/conversion-coordinates/
  13. #     projection-conversions/lambert-conformal-conic/index.aspx
  14. # [2] http://www.geo.ut.ee/~raivo/ESTCOORD.HTML
  15. #
  16. # Copyright (c) 2012 University of Tartu
  17. #
  18. # Permission is hereby granted, free of charge, to any person obtaining a copy of
  19. # this software and associated documentation files (the "Software"), to deal in
  20. # the Software without restriction, including without limitation the rights to
  21. # use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies
  22. # of the Software, and to permit persons to whom the Software is furnished to do
  23. # so, subject to the following conditions:
  24. #
  25. # The above copyright notice and this permission notice shall be included in all
  26. # copies or substantial portions of the Software.
  27. #
  28. # THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
  29. # IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  30. # FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
  31. # AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  32. # LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
  33. # OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
  34. # SOFTWARE.
  35.  
  36. use Exporter 'import';
  37. use Math::Trig;
  38.  
  39. our @EXPORT = qw(
  40.     lest97_to_wsg84
  41.     wsg84_to_lest97
  42. );
  43.  
  44. # Projection constants
  45.  
  46. use constant {
  47.    
  48.     # GRS-80 ellipsoid parameters
  49.    
  50.     R       => 6378137,
  51.     ECC     => 0.0033528106811823188,
  52.    
  53.     # L-EST97 projection parameters
  54.    
  55.     PHI1    => 1.01229096615671, # 1st std parallel, 58d00' N
  56.     PHI2    => 1.03556202279179, # 2nd std parallel, 59d20' N
  57.     PHI0    => 1.00387069499363, # Origin parallel, 57d31'03.19415" N
  58.     LAMBDA0 => 0.41887902047864, # Origin latitude, 24d00' E
  59.     N0      => 6375000,          # False Northing
  60.     E0      => 500000,           # False Easting
  61. };
  62.  
  63. # Derived projection constants
  64.  
  65. my $e    = sqrt(2*ECC - ECC**2);
  66. my $m1   = _const_m(PHI1);
  67. my $m2   = _const_m(PHI2);
  68. my $t0   = _const_t(PHI0);
  69. my $t1   = _const_t(PHI1);
  70. my $t2   = _const_t(PHI2);
  71. my $n    = (log($m1)-log($m2)) / (log($t1)-log($t2));
  72. my $F    = $m1 / ($n*($t1**$n));
  73. my $rho0 = _const_rho($t0);
  74.  
  75. sub lest97_to_wsg84 {
  76.    
  77.     # Convert L-EST 97 coordinates to
  78.     # WSG84 latitude and longitude (degrees).
  79.    
  80.     my ($x, $y) = @_;
  81.  
  82.     my $lat;
  83.     my $lon;
  84.     my $Nprim = $y - N0;
  85.     my $Eprim = $x - E0;
  86.     my $rhoprim = sqrt($Eprim**2 + ($rho0-$Nprim)**2); 
  87.     my $tprim = ($rhoprim / (R * $F))**(1/$n);
  88.     my $gammaprim = atan($Eprim / ($rho0 - $Nprim));
  89.    
  90.     $lat = pi/2 - 2*atan($tprim);
  91.     $lon = ($gammaprim/$n) + LAMBDA0;
  92.    
  93.     my $old_lat = 0; while (abs($old_lat - $lat) > 1e-6) {
  94.         my $ax;    
  95.         $old_lat = $lat;
  96.         $ax = $tprim * ((1-$e*sin($lat)) / (1+$e*sin($lat)))**($e/2);
  97.         $lat = pi/2 - 2*atan($ax);
  98.     }
  99.    
  100.     (rad2deg($lat), rad2deg($lon));
  101. }
  102.  
  103. sub wsg84_to_lest97 {
  104.    
  105.     # Convert WSG84 latitude and longitude (degrees)
  106.     # to L-EST 97 coordinates.
  107.    
  108.     my ($lat, $lon) = @_;
  109.    
  110.     $lat = deg2rad($lat);
  111.     $lon = deg2rad($lon);
  112.    
  113.     my $t = _const_t($lat);
  114.     my $rho = _const_rho($t);
  115.     my $gamma = $n*($lon - LAMBDA0);
  116.     my $y = int(N0 + $rho0 - $rho*cos($gamma));
  117.     my $x = int(E0 + $rho*sin($gamma));
  118.    
  119.     ($x, $y);
  120. }
  121.  
  122. sub _const_m {
  123.     my ($phi) = @_;
  124.     cos($phi) / sqrt(1 - $e**2 * sin($phi)**2);
  125. }
  126.  
  127. sub _const_t {
  128.     my ($phi) = @_;
  129.     my $ax = tan(pi/4 - $phi/2);
  130.     my $bx = (
  131.         (1 - $e*sin($phi)) /
  132.         (1 + $e*sin($phi))
  133.     ) ** ($e / 2);
  134.    
  135.     $ax / $bx;
  136. }
  137.  
  138. sub _const_rho {
  139.     my ($t) = @_;  
  140.     R * $F*($t**$n);   
  141. }
  142.  
  143. 1;
RAW Paste Data