#!/usr/bin/perl -w
# $Id: lonlat2glims_id,v 1.5 2002-07-24 15:21:31-06 braup Exp braup $
# Program to take a list of lon/lat values as input (filename on the
# command line or as standard input) and generate GLIMS glacier IDs,
# printed to standard output.

($progname = $0)        =~ s!^.*/!!;  # get basename of program

use Getopt::Std;
use POSIX;

$version = '$Revision: 1.5 $';
($version) = $version =~ /^\$Revision:\s*(\d+\.\d*)/;

$usage = "$progname version $version
Usage:
  $progname  -h            (prints this help message and exits)
OR
  $progname  [-v] [-l] [lon_lat_filename]
where
  -l    Output will be 3 columns:  lon, lat, and ID.  Otherwise, only
        the ID is output.

  -v    specifies verbose mode

  This program reads a list of lon/lat coordinates and creates a GLIMS
  glacier ID for each one, and prints the results to standard output.  If
  no files are specified on the command line, then standard input is read.

  If the longitude is not in the range 0-360, then a warning is printed to
  standard error, and it is fixed.  If a latitude is outside the range
  [-90,90], then processing is aborted.

  The input is expected to have two columns, in lon lat order.  Comments
  are allowed, and are preceded by a '#' character.  Blank lines are also
  allowed.
";

$opt_h = $opt_v = $opt_l = '';

if (! getopts('lhv') ) {
  die "$usage\n";
}

print STDERR "$progname version $version\n" if $opt_v;

die "$usage\n" if $opt_h;

while (<>) {
  next if /^\s*#/;      # skip comments
  next if /^\s*$/;      # skip blank lines
  s/\s*#.*$//;          # strip comments off of data lines
  ($lon,$lat) = split;
  $id = &make_glims_id( $lon, $lat );
  if ($opt_l) {
    write;
  } else {
    print "$id\n";
  }
}

sub make_glims_id {
  my ($lon, $lat) = @_;
  $lon = &fix_lon( $lon );
  &check_lat( $lat );
  my $SN = $lat < 0 ? "S" : "N";
  $lat = abs($lat);
  my $id = sprintf( "G%07.3fE%06.3f$SN", $lon, $lat );
  $id =~ s/\.//g;
  return $id;
}

sub fix_lon {
  my $lon = $_[0];
  my $absolute_lower_limit = -360;
  my $absolute_upper_limit = 360;

  if ($lon >= 0 && $lon < 360) {
    return $lon;
  } else {
    warn "Input longitude not in range [0,360).  Fixing.\n";
    # This next check isn't necessary for the fix, but it might be a good
    # idea for sanity.
    if ($lon < $absolute_lower_limit || $lon > $absolute_upper_limit) {
      die "Input longitude is outside the range [$absolute_lower_limit,$absolute_upper_limit].
      Maybe something is really wrong!  Aborting.\n";
    }
    $lon = fmod($lon, 360.0);
    $lon += 360.0 if $lon < 0;
  }
  return $lon;
}

sub check_lat {
  my $lat = $_[0];

  if ($lat < -90 || $lat > 90) {
    die "Input latitude is outside the range [-90,90].  Aborting.\n";
  }
}

format STDOUT =
@<<<<<<<<<<<<<<<<<<  @<<<<<<<<<<<<<<<<<<  @<<<<<<<<<<<<<<<<<<<<<<
$lon,$lat,$id
.