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