#!/usr/bin/perl

use Geo::Coordinates::UTM;
use GD::Polyline;
use LWP::UserAgent;

if (scalar(@ARGV) == 0) {
 print "You must provide the OSM relation code for the polygon you want to calculate:\n";
 print " http://wiki.openstreetmap.org/wiki/WikiProject_Germany/Grenzen\n";
 exit 1;
}

my $polyline = new GD::Polyline;

print "Fetching data from OSM ...\n";

my $ua = LWP::UserAgent->new;
my $res = $ua->get("http://api.openstreetmap.org/api/0.6/relation/".$ARGV[0]."/full");
my $xml = $res->content;

my $zone = "";
my $i = 0;

print "Parsing XML ...\n";

while ($xml =~ m/<node id="\d*" lat="(\d\d\.\d*)" lon="(\d\.\d*)"/g) {
 my $lat = $1;
 my $lon = $2;

 my ($tzone, $east, $north) = latlon_to_utm(23, $lat, $lon);
 $polyline->addPt($east, $north);
 $i++;

 if ($zone eq "") {
  $zone = $tzone;
 }

 if ($zone ne $tzone) {
  print "ERROR: Polygon crosses zone border!\n";
  exit 1;
 }
}

print "Polygon contains ".$i." points.\n";
print "Calculating centroid ...\n";

my ($cx, $cy) = $polyline->centroid();
my ($clat, $clon) = utm_to_latlon(23, $zone, $cx, $cy);

$clat_deg = int($clat);
$clat_min = ($clat - int($clat)) * 60;

$clon_deg = int($clon);
$clon_min = ($clon - int($clon)) * 60;

print $zone." ".int($cx + 0.5)." ".int($cy + 0.5)."\n";
print "N ".$clat_deg." ".sprintf("%.3f", $clat_min)." E ".$clon_deg." ".sprintf("%.3f", $clon_min)."\n";

exit 0

