#!/usr/local/bin/perl

#### geoconvert.pl -->  
####    general coordinate conversion utility for converting
####    between alvin xy, utm and lat/lon. Requires installation
####    of "proj".  Designed with UW alvin xy grid in mind.
####                  
####    In order to use in other locations, will need to reconfigure ORIGIN, 
####    UTMZONE, and HEMISPHERE environment variables. 
####                    
####    We assume that alvin xy grid is always located at 5000, 5000
####    Script also assumes northern hemisphere ... need to find out importance
####    of this for proj utility.
####
####    Created 08/06/00, FlowMow 2000 by Christian P. Sarason.Copyright is in 
####    spirit of GPL; free to use and modify; please e-mail improvements to:
####    cpsarason@oceaninquiry.net
####
####    Usage: 
####
####    geoconvert.pl [-xl<min> -lx<min> -lu<min> -ul<min> -xu -ux] <filename> 
####                      
####    arguments tell geoconvert how to convert.
####
####    Appeding "min" on argument says that preferred input/output is in 
####    deg/min/sec. Script expects stdin, other wise will slurp the filename 
####    provided.
####
####    Default input/output will be in degrees and dec. minutes for lat/lon 
####    To Do:
####    error checking --> sanity checks to make sure that input is the same
####        as expected for each kind of conversion (eg input is XY for -xl)
####    
####    interactive mode??
####    CGI mode??
####    Speed up???  Currently forks proj for every line.  Is there a faster 
####    way??


## $ENV{'ORIGIN'} = "47.08939 -129.1645167";
## $ENV{'ORIGIN_LAT'} = "47.08939";
## $ENV{'ORIGIN_LON'} = "-129.1645167";

#### UTM origin is in middle of "box": nominal coords for S&M structure.
#### UTM coords obtained by using xy2ll.kf. See output below:
####
####     [cpsarason@palau ctd]$ echo 4930 6000 | xy2ll.kf               
####     Lat/Lon = 	47.9478625	-129.0985884
####     [cpsarason@palau ctd]$ echo 47.9478625 -129.0985884 | proj +proj=utm +zone=9 +ellps=WGS84 -r
####     492638.34	5310510.09

$ENV{'UTMORIGIN_X'} = "492638.34";
$ENV{'UTMORIGIN_Y'} = "5310510.09";

$ENV{'UTMZONE'} = "9";
$ENV{'ELLPS'} = "WGS84";
$ENV{'HEMISPHERE'} = "north";

### Get flag to decide conversion scheme

$flag = shift(@ARGV);

### code below for XY to lat lon (dec. degrees) ###

if ($flag eq "-xl") {
    while (<>) {
       chop;
       ($alvinx, $alviny) = split(/\s+/);
       $utmx = $alvinx + $ENV{'UTMORIGIN_X'} - 4930;
       $utmy = $alviny + $ENV{'UTMORIGIN_Y'} - 6000;
       $command="echo $utmx $utmy | proj +proj=utm +zone=$ENV{'UTMZONE'} +ellps=$ENV{'ELLPS'}  -I -s -f \%.8f";
       $projout = `$command`;

       print "$projout"
  }
}


elsif ($flag eq "-xlmin") {
    while (<>) {
       chop;
       ($alvinx, $alviny) = split(/\s+/);
       $utmx = $alvinx + $ENV{'UTMORIGIN_X'} - 4930;
       $utmy = $alviny + $ENV{'UTMORIGIN_Y'} - 6000;
       $command="echo $utmx $utmy | proj +proj=utm +zone=$ENV{'UTMZONE'} +ellps=$ENV{'ELLPS'}  -I -s -f \%.8f";
       $projout = `$command`;
       ($dec_lat, $dec_lon) = split(/\s+/, $projout);

       print dec_to_degmin($dec_lat, $dec_lon);
  }
}

### code for converting LL to XY ---> 
###
###     This part of the code works by converting LL values to UTM and 
###     then subtracting the "origin" {which is at S&M, (4930, 6000)}
###     and adding up the appropriate coords to get to the Robigou XY grid.


elsif ($flag eq "-lx") { 
     while (<>) {
	chop;
 	$line = $_;
	$command = "echo $line | proj +proj=utm +zone=$ENV{'UTMZONE'} +ellps=$ENV{'ELLPS'}  -r -f \%.2f";
	$projout = `$command`; 
	($utmx, $utmy) = split(/\s+/, $projout);

	$alvinx = $utmx - $ENV{'UTMORIGIN_X'} + 4930;
	$alviny = $utmy - $ENV{'UTMORIGIN_Y'} + 6000;

	printf "%.2f %.2f\n", $alvinx, $alviny;
  }
}

elsif ($flag eq "-lxmin") { 
     while (<>) {
	chop;
	($deg_lat, $min_lat, $deg_lon, $min_lon) = split(/\s+/);
        $line = degmin_to_dec($deg_lat, $min_lat, $deg_lon, $min_lon);
	$command = "echo $line | proj +proj=utm +zone=$ENV{'UTMZONE'} +ellps=$ENV{'ELLPS'}  -r -f \%.2f";
	$projout = `$command`; 
	($utmx, $utmy) = split(/\s+/, $projout);
	$alvinx = $utmx - $ENV{'UTMORIGIN_X'} + 4930;
	$alviny = $utmy - $ENV{'UTMORIGIN_Y'} + 6000;

	printf "%.2f %.2f\n", $alvinx, $alviny;
  }
}

### code for converting between lat/lon and UTM below ###

elsif ($flag eq "-lu") {
 
 while (<>) {
	chop;
	$line = $_;
	$command = "echo $line | proj +proj=utm +zone=$ENV{'UTMZONE'} +ellps=$ENV{'ELLPS'}  -r -f \%.2f";
	$projout = `$command`; 

	print "$projout"
  }
}

elsif ($flag eq "-lumin") {
 
 while (<>) {

        chop;
	($deg_lat, $min_lat, $deg_lon, $min_lon) = split(/\s+/);
        $line = degmin_to_dec($deg_lat, $min_lat, $deg_lon, $min_lon);
	$command = "echo $line | proj +proj=utm +zone=$ENV{'UTMZONE'} +ellps=$ENV{'ELLPS'}  -r -f \%.2f";
	$projout = `$command`; 

	print "$projout"
  }
}

### code below for converting UTM to LL ###

elsif ($flag eq "-ul") {

 while (<>) {
	chop;
	$line = $_;
	$command = "echo $line | proj +proj=utm +zone=$ENV{'UTMZONE'} +ellps=$ENV{'ELLPS'} -I -s -f \%.5f";
	$projout = `$command`; 
	($dec_lat, $dec_lon) = split(/\s+/, $projout);

	print "$dec_lat $dec_lon\n";
  }
}
### code below for converting UTM to LL (deg/min) ###

elsif ($flag eq "-ulmin") {

 while (<>) {
	chop;
	$line = $_;
	$command = "echo $line | proj +proj=utm +zone=$ENV{'UTMZONE'} +ellps=$ENV{'ELLPS'} -I -s -f \%.5f";
	$projout = `$command`; 
	($dec_lat, $dec_lon) = split(/\s+/, $projout);

	print dec_to_degmin($dec_lat, $dec_lon);
  }
}


### code for XY to UTM ... calls xy2ll.kf, pipes to proj, outputs

elsif ($flag eq "-xu") { 

    while (<>) {
	chop;
	($alvinx, $alviny) = split(/\s+/);
        $utmx = $alvinx + $ENV{'UTMORIGIN_X'} - 4930;
        $utmy = $alviny + $ENV{'UTMORIGIN_Y'} - 6000;

	print "$utmx $utmy\n";
    }
}


### code for converting UTM to XY 

elsif ($flag eq "-ux") { 
    while (<>) {
	chop;
	($utmx, $utmy) = split(/\s+/);
	$alvinx = $utmx - $ENV{'UTMORIGIN_X'} + 4930;
	$alviny = $utmy - $ENV{'UTMORIGIN_Y'} + 6000;
	
	printf "%.2f %.2f\n", $alvinx, $alviny;
    }


}    

else { 

    print "Usage: geoconvert.pl [-xl<min> -lx<min> -lu<min> -ul<min> -xu -ux] <filename>\n";
    print "\n       The argument indicates direction of conversion, and\n";
    print "       min is appended to recieve or output in degrees/minutes.\n";
    print "\n       examples: echo 47d 8.067m -129d 5.912m | geoconvert.pl -lxmin\n";
    print "                 echo 5000 5000 | geoconvert.pl -xu\n\n";

}


sub dec_to_degmin {
    $dec_lat = $_[0];
    $dec_lon = $_[1];

    	($deg_lat, $fraction_lat) = split(/\./, $dec_lat);
	($deg_lon, $fraction_lon) = split(/\./, $dec_lon);

	$fraction_lat = "0.$fraction_lat";
	$fraction_lon = "0.$fraction_lon";

	$min_lat = $fraction_lat*60;
	$min_lon = $fraction_lon*60;
	
	$coord_string = sprintf "%2.0fd %2.3fm %2.0fd %2.3fm\n", $deg_lat,$min_lat, $deg_lon, $min_lon;
    return $coord_string;
}

sub degmin_to_dec {

    chop($_[0]); $deg_lat = $_[0];
    chop($_[1]); $min_lat = $_[1];
    chop($_[2]); $deg_lon = $_[2];
    chop($_[3]); $min_lon = $_[3];
 
    $frac_lat = $min_lat/60;
    $frac_lon = $min_lon/60;

    if ($deg_lat >=0 ) {
        $dec_lat = $deg_lat + $frac_lat;
    }
    else {
	$dec_lat = $dec_lat - $frac_lat;
    }

    if ($deg_lon >=0 ) {
	$dec_lon = $deg_lon + $frac_lat;
    }
    else { 
	$dec_lon = $deg_lon - $frac_lon;
    }

    return sprintf "%.5f %.5f", $dec_lat, $dec_lon;
}
    

  


