#!/usr/local/bin/perl
#
# stategeo.pl Gerry Daumiller 2/8/99
# CGI program to convert from Montana State Plane Coordinates
# to Latitude/Longitude
# The formulae this program is based on are from "Map Projections,
# A Working Manual" by John P. Snyder, U.S. Geological Survey
# Professional Paper 1395, 1987, pages 295-298
# Resolve and unencode name/value pairs into %input
foreach (split('&', $ENV{'QUERY_STRING'})) {
s/\+/ /g ;
($name, $value)= split('=', $_, 2) ;
$name=~ s/%(..)/chr(hex($1))/ge ;
$value=~ s/%(..)/chr(hex($1))/ge ;
$input{$name}.= $value ;
}
$y = $input{y} ;
$x = $input{x} ;
print "Content-type: text/html\n\n" ;
print "
" ;
if (($y eq '') or ($x eq ''))
{ print 'No value entered' ;
print "\n" ;
print "No value was entered for Northing or Easting
\n" ;
print "" ;
exit ;
}
# Set up the coordinate system parameters.
# $A = 6378206.4 # major radius of ellipsoid, map units (NAD27)
# $E = 0.08227185422 # eccentricity of ellipsoid (NAD27)
$A = 6378137 ; # major radius of ellipsoid, map units (NAD83)
$E = 0.08181922146 ; # eccentricity of ellipsoid (NAD83)
$AngRad = 0.01745329252 ; # number of radians in a degree
$Pi4 = 3.141592653582 / 4 ; # Pi / 4
$P1 = 45 * $AngRad ; # latitude of first standard parallel
$P2 = 49 * $AngRad ; # latitude of second standard parallel
$P0 = 44.25 * $AngRad ; # latitude of origin
$M0 = -109.5 * $AngRad ; # central meridian
$X0 = 600000 ; # False easting of central meridian, map units
# Calculate the coordinate system constants.
$m1 = cos($P1) / sqrt(1 - (($E**2) * ((sin($P1))**2))) ;
$m2 = cos($P2) / sqrt(1 - (($E**2) * ((sin($P2))**2))) ;
$t1 = sin($Pi4 - ($P1 / 2)) / cos($Pi4 - ($P1 / 2)) ;
$t1 = $t1 / (((1 - ($E * (sin($P1)))) / (1 + ($E * (sin($P1)))))**($E/2)) ;
$t2 = sin($Pi4 - ($P2 / 2)) / cos($Pi4 - ($P2 / 2)) ;
$t2 = $t2 / (((1 - ($E * (sin($P2)))) / (1 + ($E * (sin($P2)))))**($E/2)) ;
$t0 = sin($Pi4 - ($P0 / 2)) / cos($Pi4 - ($P0 / 2)) ;
$t0 = $t0 / (((1 - ($E * (sin($P0)))) / (1 + ($E * (sin($P0)))))**($E/2)) ;
$n = log($m1 / $m2) / log($t1 / $t2) ;
$F = $m1 / ($n * ($t1**$n)) ;
$rho0 = $A * $F * ($t0**$n) ;
# Convert the coordinate to Latitude/Longitude.
# Calculate the Longitude.
$x = $x - $X0 ;
$Pi2 = $Pi4 * 2 ;
$rho = sqrt(($x**2) + (($rho0 - $y)**2)) ;
$theta = atan2($x,($rho0 - $y)) ;
$t = ($rho / ($A * $F))**(1 / $n) ;
$Lon = $theta / $n + $M0 ;
$x = $x + $X0 ;
# Estimate the Latitude
$Lat0 = $Pi2 - (2 * atan2($t,1)) ;
# Substitute the estimate into the iterative calculation that
# converges on the correct Latitude value.
$part1 = (1 - ($E * sin($Lat0))) / (1 + ($E * sin($Lat0))) ;
$Lat1 = $Pi2 - (2 * atan2(($t * ($part1**($E/2))),1)) ;
do {
$Lat0 = $Lat1 ;
$part1 = (1 - ($E * sin($Lat0))) / (1 + ($E * sin($Lat0))) ;
$Lat1 = $Pi2 - (2 * atan2(($t * ($part1**($E/2))),1)) ;
} until abs($Lat1 - $Lat0) < 0.000000002 ;
# Convert from radians to degrees.
$Lat1 = $Lat1 / $AngRad ;
$Lon = $Lon / $AngRad ;
# Calcuate degrees, minutes, and seconds.
$dLat = int($Lat1) ;
$mLat = 60 * ($Lat1 - $dLat) ;
$sLat = 60 * ($mLat - int($mLat)) ;
$mLat = int($mLat) ;
$dLon = int($Lon) ;
$mLon = 60 * ($dLon - $Lon) ;
$sLon = 60 * ($mLon - int($mLon)) ;
$mLon = int($mLon) ;
# Format the values for output.
$Lat = sprintf ("%.5f",$Lat1) ;
$Lon = sprintf ("%.5f",$Lon) ;
$sLat = sprintf ("%.2f",$sLat) ;
$sLon = sprintf ("%.2f",$sLon) ;
$x = sprintf ("%.0f",$x) ;
$y = sprintf ("%.0f",$y) ;
# Adjust minutes and seconds, in case they have rounded to 60.
if ($sLat == 60.00)
{ $sLat = 0 ;
$mLat = $mLat + 1;
}
if ($mLat ge 60)
{ $mLat = $mLat - 60 ;
$dLat = $dLat + 1 ;
}
if ($sLon == 60.00)
{ $sLon = 0 ;
$mLon = $mLon + 1;
}
if ($mLon ge 60)
{ $mLon = $mLon - 60 ;
$dLon = $dLon - 1 ;
}
#Send the HTML with the answers.
print 'State Plane Converted to Latitude/Longitude' ;
print "\n" ;
print "State Plane Converted to Latitude/Longitude
\n" ;
print "\n" ;
print "\n" ;
print " \n" ;
print " | \n" ;
print " State Plane Coordinate |
\n" ;
print " \n" ;
print " | \n" ;
print " Northing | \n" ;
print " Easting |
\n" ;
print " \n" ;
print " Meters | \n" ;
print " $y | \n" ;
print " $x |
\n" ;
print " |
\n" ;
print " |
\n" ;
print " \n" ;
print " | \n" ;
print " Latitude: | \n" ;
print " Longitude: |
\n" ;
print " \n" ;
print " Decimal Degrees: | \n" ;
print " $Lat | \n" ;
print " $Lon |
\n" ;
print " \n" ;
print " | \n" ;
print " | \n" ;
print " |
\n" ;
print " \n" ;
print " Degrees, Minutes, Seconds: | \n" ;
print " $dLat \n" ;
print " $mLat \n" ;
print " $sLat | \n" ;
print " $dLon \n" ;
print " $mLon \n" ;
print " $sLon |
\n" ;
print "
\n" ;
print "
\n" ;
print "
\n" ;
print "" ;