#!/usr/bin/perl

# Copyright (C) 2010 Thorsten Kukuk
#
# This program is free software; you can redistribute it and/or
# modify it under the terms of the GNU General Public License
# in Version 2 as published by the Free Software Foundation.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
# MA  02110-1301, USA.

=head1 NAME

geo-centroid - Calculate area and centroid of a polygon

=head1 SYNOPSIS

geo-centroid [options] <input.txt>

=head1 DESCRIPTION

geo-centroid calculates the area and centroid of a polygon.
There are two ways to calculate the centroid: one based on
the area and one as average of all corners.

=head1 OPTIONS

  --degdec            Print only decimal format
  --degdec1           Print only decimal format with prefix
  --mindec            Print in DDD° MM.MMM format
  --mindec1           Shorter variant
  --mindec2           Print in DDD.MM.MMM format
  --utm               Print in UTM format
  --version           Print version number and exit
  --man               Print manual page
  --usage             Display usage
  -h|-?|--help        Help

=head1 INPUT

geo-centroid expects as input a list of latitude/longitude pairs,
line by line for every corner of the polygon.

=head1 COPYRIGHT

Copyright (c) 2010 by Thorsten Kukuk.  All rights reserved.

This package is free software; you can redistribute it and/or modify
it under the GPL version <http://gnu.org/licenses/gpl2.html>.
There is NO WARRANTY, to the extent permitted by law.

=cut


use strict;
use warnings;
use Pod::Usage;
use Getopt::Long;
use GEO::Coords;
use Geo::Coordinates::UTM;

my $Version = '(geo-tools) 1.23';

my $help = 0;
my $man = 0;
my $usage = 0;
my $version = 0;
my $print_degdec = 0;
my $print_degdec1 = 0;
my $print_mindec = 0;
my $print_mindec1 = 0;
my $print_mindec2 = 0;
my $print_dms = 0;
my $print_utm = 0;
my $cnt = 0;
my @wpts;

GetOptions('degdec' => \$print_degdec,
           'degdec1' => \$print_degdec1,
           'mindec' => \$print_mindec,
           'mindec1' => \$print_mindec1,
           'mindec2' => \$print_mindec2,
           'dms' => \$print_dms,
           'utm' => \$print_utm,
	   'version' => \$version,
           'man' => \$man,
	   'usage|u' => \$usage,
           'help|h|?' => \$help) or pod2usage(2);
pod2usage(0) if $help;
pod2usage(-exitstatus => 0, -verbose => 2) if $man;
pod2usage(-exitstatus => 0, -verbose => 0) if $usage;

if ($version) {
  print "geo-centroid $Version\n";
  exit;
}

$print_mindec = !($print_degdec || $print_degdec1 || $print_mindec ||
		  $print_mindec1 ||$print_mindec2 || $print_dms ||
		  $print_utm);

while (<>) {
  my $line = $_;
  chomp ($line);

  ($wpts[$cnt][0], $wpts[$cnt][1]) = GEO::Coords::parse2($line);
  $cnt++;
}

my $zone;
my $area = 0;
for (my $i = 0; $i < $cnt; $i++) {
  my $j = ($i + 1) % $cnt;
  my ($xi, $xj, $yi, $yj);

  if (!$zone) {
    ($zone,$xi,$yi) = latlon_to_utm('wgs84',$wpts[$i][0],$wpts[$i][1]);
  } else {
    ($zone,$xi,$yi) = latlon_to_utm_force_zone ('wgs84', $zone, $wpts[$i][0],$wpts[$i][1]);
  }
  ($zone,$xj,$yj) = latlon_to_utm_force_zone('wgs84', $zone, $wpts[$j][0],$wpts[$j][1]);

  $area += ($xi * $yj - $xj * $yi);
}
$area = abs ($area/2);
print "Area: $area m^2\n";

my $cx = 0;
my $cy = 0;
my $lats = 0;
my $lons = 0;
my $centroid = 0;
for (my $i = 0; $i < $cnt; $i++) {
  my $j = ($i + 1) % $cnt;
  my ($zone1,$xi,$yi) = latlon_to_utm('wgs84',$wpts[$i][0],$wpts[$i][1]);
  my ($zone2,$xj,$yj) = latlon_to_utm('wgs84',$wpts[$j][0],$wpts[$j][1]);

  $cx += ($xi+$xj)*($xi*$yj - $xj*$yi);
  $cy += ($yi+$yj)*($xi*$yj - $xj*$yi);

  $lats += $wpts[$i][0];
  $lons += $wpts[$i][1];
}

$cx /= (6*$area);
$cy /= (6*$area);

my ($x, $y) = utm_to_latlon ('wgs84', $zone, $cx, $cy);

print "Centroid (via area):\n";
GEO::Coords::print_decdeg($x, $y) if ($print_degdec);
GEO::Coords::print_decdeg($x, $y, 1) if ($print_degdec1);
GEO::Coords::print_mindec($x, $y, 0) if ($print_mindec);
GEO::Coords::print_mindec($x, $y, 1) if ($print_mindec1);
GEO::Coords::print_mindec($x, $y, 2) if ($print_mindec2);
GEO::Coords::print_dms($x, $y, 0) if ($print_dms);
GEO::Coords::print_utm($x, $y) if ($print_utm);

$lats /= $cnt;
$lons /= $cnt;
printf "Centroid (via sum):\n";
GEO::Coords::print_decdeg($lats, $lons) if ($print_degdec);
GEO::Coords::print_decdeg($lats, $lons, 1) if ($print_degdec1);
GEO::Coords::print_mindec($lats, $lons, 0) if ($print_mindec);
GEO::Coords::print_mindec($lats, $lons, 1) if ($print_mindec1);
GEO::Coords::print_mindec($lats, $lons, 2) if ($print_mindec2);
GEO::Coords::print_dms($lats, $lons, 0) if ($print_dms);
GEO::Coords::print_utm($lats, $lons) if ($print_utm);

1;
