#!/usr/local/bin/perl
#
# geostate.pl Gerry Daumiller 2/8/99
# CGI program to convert from Latitude/Longitude to Montana State
# Plane Coordinates. This program assumes that the coordinate is
# in the Western Hemisphere.
# 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 ;
}
$Lat = $input{Lat} ;
$Lon = $input{Lon} ;
$dLat = $input{Dlat} ;
$mLat = $input{Mlat} ;
$sLat = $input{Slat} ;
$dLon = $input{Dlon} ;
$mLon = $input{Mlon} ;
$sLon = $input{Slon} ;
print "Content-type: text/html\n\n" ;
print "
" ;
# If decimal latitude was entered, convert it to degrees, minutes, and
# seconds. If latitude was entered in degrees, minutes, and seconds,
# convert to decimal degrees.
if ($Lat ne '')
{ $dLat = int($Lat) ;
$mLat = 60 * ($Lat - $dLat) ;
$sLat = 60 * ($mLat - int($mLat)) ;
$mLat = int($mLat) ;
}
elsif ($dLat ne '')
{ if ($mLat eq '') {$mLat = 0;}
if ($sLat eq '') {$sLat = 0;}
$Lat = $dLat + ($mLat / 60) + ($sLat / 3600) ;
}
else
{ print 'No latitude entered' ;
print "\n" ;
print "No value was entered for Latitude
\n" ;
print "" ;
exit ;
}
# If decimal longitude was entered, convert it to degrees, minutes, and
# seconds. If longitude was entered in degrees, minutes, and seconds,
# convert to decimal degrees.
if ($Lon ne '')
{ if ($Lon > 0) {$Lon = 0 - $Lon;}
$dLon = int($Lon) ;
$mLon = 60 * ($dLon - $Lon) ;
$sLon = 60 * ($mLon - int($mLon)) ;
$mLon = int($mLon) ;
}
elsif ($dLon ne '')
{ if ($dLon > 0) {$dLon = 0 - $dLon;}
if ($mLon eq '') {$mLon = 0;}
if ($sLon eq '') {$sLon = 0;}
$Lon = $dLon - ($mLon / 60) - ($sLon / 3600) ;
}
else
{ print 'No longitude entered' ;
print "\n" ;
print "No value was entered for Longitude
\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 latitude/longitude to a coordinate.
$Lat = $Lat * $AngRad ;
$Lon = $Lon * $AngRad ;
$t = sin($Pi4 - ($Lat / 2)) / cos($Pi4 - ($Lat / 2)) ;
$t = $t / (((1 - ($E * (sin($Lat)))) / (1 + ($E * (sin($Lat)))))**($E/2)) ;
$rho = $A * $F * ($t**$n) ;
$theta = $n * ($Lon - $M0) ;
$x = ($rho * sin($theta)) + $X0 ;
$y = $rho0 - ($rho * cos($theta)) ;
$Lat = $Lat / $AngRad ;
$Lon = $Lon / $AngRad ;
# Format the values for output.
$Lat = sprintf ("%.5f",$Lat) ;
$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 'Latitude/Longitude Converted to State Plane' ;
print "\n" ;
print "Latitude/Longitude Converted to State Plane
\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 " | \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 "