#!/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 " \n" ; print " \n" ; print " \n" ; print " \n" ; print " \n" ; print " \n" ; print " \n" ; print " \n" ; print " \n" ; print " \n" ; print " \n" ; print " \n" ; print " \n" ; print " \n" ; print " \n" ; print " \n" ; print " \n" ; print " \n" ; print " \n" ; print " \n" ; print " \n" ; print " \n" ; print " \n" ; print " \n" ; print " \n" ; print " \n" ; print " \n" ; print "
State Plane Coordinate
NorthingEasting
Meters$y$x
Latitude:Longitude:
Decimal Degrees:$Lat$Lon
Degrees, Minutes, Seconds:$dLat  \n" ; print " $mLat  \n" ; print " $sLat  $dLon  \n" ; print " $mLon  \n" ; print " $sLon  
\n" ; print "\n" ; print "
\n" ; print "" ;