#!/usr/bin/perl -l -w
#------------------------------------------------------------------------
#  $Id: creategmt.pl 1273 2009-09-04 14:12:18Z bzfraack $
#------------------------------------------------------------------------

use Encode;
use Getopt::Std;
use POSIX;

# expected input format:
# Nodefile :
#  nodename;logitude;latitude;size[;color]
#     color denotes the entry number in the nodecolor table
#
#  The node radius is caculated: rad := size/maxsize*maxrad
#
# Linkfile :
#  linkname;source;target;flow;capacity
#
#  Source and target are given as nodenames.
#  Capacity is used to calculate the width of a link: width := cap/maxcap*maxwidth;
#  Flow and capacity are used to calculate the color of a link
#  from the 2 defined colorvectors Cfull and Cempty:
#  usage := flow/cap ;
#  if usage < 0.5 then color := Cempty + 2 * usage * Cfull;
#  else (usage >= 0.5) color := Cfull + 2 * (1-usage) * Cempty ;

# TODO:
#  - features:
#     --> draw parallel links as arcs
#     --> plot more than one layer in parallel in a perspective 3d view

sub usage(){
   print 'usage: creategmt.pl [options] <nodefile> <linkfile> ';
   print '  options: ';
   print '     -c <country> (us|germany|austria|poland|eu)';
   print '                                         default = germany ';
   print '                                                           ';
   print '     -r <maximal node radius>            default = 50      ';
   print '     -R 0|1  use node radius             default = 1       ';
   print '     -s <std. node radius>               default = 4       ';
   print '     -n <nodecolor as r/g/b>             default = 0/0/100 ';
   print '                                                           ';
   print '     -N 0|1  write node names            default = 1       ';
   print '     -E 0|1  write edge names            default = 0       ';
   print '     -t textcolor as r/g/b               default = 0/0/0   ';
   print '     -g textbgcolor as r/g/b             default = 255/255/255';
   print '     -z textsize                         default = 7       ';
   print '                                                           ';
   print '     -w <maximal edge width>             default = 8       ';
   print '     -W 0|1  use edge width              default = 1       ';
   print '     -b <std. link width>                default = 1       ';
   print '     -e <color of empty link as r/g/b >  default = 0/240/0 ';
   print '     -f <color of full link as r/g/b>    default = 255/0/0 ';
   print '     -C 0|1  use link color              default = 1       ';
   print '     -l <std. link color as r/g/b>       default = 0/0/100 ';
   print '     -U show unused links                default = 1       ';
   print '                                                           ';
   print '     -v <elevation 0-90>                 default = 90      ';
   print '     -a <azimuth 0-360>                  default = 180     ';
   print '                                                           ';
   print '     -o <output file>                    default = out.eps ';
   print '     -q <title>                          default = \$country';
   print '                                                           ';
   print '     -h show this usage information and exit               ';

   exit(1);
}

$TMP_SH_SCRIPT = POSIX::tmpnam();
$outfile       = "out.eps";
$title         = "My title";

$z         = "0";
$elevation = 90;
$azimuth   = 180;

# maximal radius for demand node
$maxrad = 50;
$userad = 1; # node radius depends on emanating demand?
$stdrad = 4; # radius for all nodes if not demanddep.
$minrad = 3; # minimal radius for nodes


# maximal line width
$maxwidth    = 8;
$usewidth    = 1;   # linewidth depends on capacity
$stdwidth    = 1;   # linewidth to use for all nodes if not capdep.
$showUnused  = 1;   # show links with cap = 0 as thin $unusedcolor lines?
$minwidth    = 1;   # minimal width for lines
$unusedwidth = 1;   # minimal width for unused links
#$unusedwidth = 0.2; # minimal width for unused links

# set used colors
$usecolor    = 1;          # use capacity usage for color of links
$stdcolor    = "0/0/100";  # color for all links if not usage depending
#$unusedcolor = "10";      # color of unused links
$unusedcolor = "50";      # color of unused links
$nodecolor   = "0 0 100";
$nodeborderc = "0/0/100";
$colorempty = "0/240/0";
$colorfull = "255/0/0";

#define nodecolor table:
# entry 0 = default nodecolor
@nodecarr = ( $nodecolor, "200 0 0", "0 200 0");

# set text properties
$textcolor = "0/0/0";
$textbgc = "255/255/255";
$textsize = 9;
$nodetextprop = "0 5 ML";   # angle fontno justify( (B)ottom (M)iddle (T)op / (L)eft (C)enter (R)ight )
$edgetextprop = "0 5 MC";
$textBoxSize = "0.05/0.05"; # size of textbox for node and edge names

# write node and edge names
$writeNodeNames = 1;
$writeEdgeNames = 0;

# set country;
#$country = "us";
$country = "germany";

my $options = ();
# set default values
$options{'r'} = $maxrad;
$options{'R'} = $userad;
$options{'s'} = $stdrad;
$options{'w'} = $maxwidth;
$options{'W'} = $usewidth;
$options{'b'} = $stdwidth;
$options{'l'} = $stdcolor;
$options{'c'} = $country;
$options{'e'} = $colorempty;
$options{'f'} = $colorfull;
$options{'N'} = $writeNodeNames;
$options{'E'} = $writeEdgeNames;
$options{'n'} = $nodecolor;
$options{'C'} = $usecolor;
$options{'U'} = $showUnused;
$options{'g'} = $textbgc;
$options{'t'} = $textcolor;
$options{'z'} = $textsize;
$options{'o'} = $outfile;
$options{'v'} = $elevation;
$options{'a'} = $azimuth;
$options{'h'} = 0;
$options{'q'} = $country;

getopts( "r:R:w:W:c:e:f:N:E:n:C:o:hU:l:b:s:t:g:z:v:a:q:", \%options );

#foreach $opt (keys %options){
#   print "$opt : $options{$opt}";
#q}

$maxrad          =  $options{'r'};
$userad          =  $options{'R'};
$stdrad          =  $options{'s'};
$maxwidth        =  $options{'w'};
$usewidth        =  $options{'W'};
$stdwidth        =  $options{'b'};
$stdcolor        =  $options{'l'};
$country         =  $options{'c'};
$colorempty      =  $options{'e'};
$colorfull       =  $options{'f'};
$writeNodeNames  =  $options{'N'};
$writeEdgeNames  =  $options{'E'};
$nodecolor       =  $options{'n'};
my $tmpncol      =  $nodecolor;
$tmpncol         =~ s/\ /\//g;
$nodeborderc     =  $tmpncol; #use same color for node and nodeborder;
$usecolor        =  $options{'C'};
$showUnused      =  $options{'U'};
$textcolor       =  $options{'t'};
$textbgc         =  $options{'g'};
$textsize        =  $options{'z'};
$outfile         =  $options{'o'};
$elevation       =  $options{'v'};
$azimuth         =  $options{'a'};
$title           =  $options{'q'};

$nodecolor =~ s/\//\ /g;
$nodecarr[0] = $nodecolor;

if( @ARGV < 2 || $options{'h'} ){
   usage();
}

$nodefile = $ARGV[0];
$linkfile = $ARGV[1];

$tmplinks   = POSIX::tmpnam();
$tmpnodes   = POSIX::tmpnam();
$tmptext    = POSIX::tmpnam();
$tmpcolor   = POSIX::tmpnam();
$tmpoutfile = POSIX::tmpnam();

sub min($$)
{
   if( $_[0] < $_[1] ){
      return $_[0];
   } else {
      return $_[1];
   }
}

sub substitute ($){
   my $text = $_[0];
   $t = encode( "utf8", $text );
   $t =~ s/Ä/\\276/g;
   $t =~ s/Ö/\\331/g;
   $t =~ s/Ü/\\335/g;
   $t =~ s/ä/\\342/g;
   $t =~ s/ö/\\363/g;
   $t =~ s/ü/\\370/g;
   $t =~ s/ß/\\373/g;
   return $t;
}

sub getColor($$$) {

   my ( $cempt, $cfull, $load ) = ( $_[0], $_[1], $_[2] );
   @cvempty  = split( /\//, $cempt );
   @cvfull   = split( /\//, $cfull );

   if ( $load < 0.5 ) {
      $red   = $cvempty[0] + 2 * $load * $cvfull[0];
      $green = $cvempty[1] + 2 * $load * $cvfull[1];
      $blue  = $cvempty[2] + 2 * $load * $cvfull[2];
   }
   else {
      $red   = $cvfull[0] + 2 * (1-$load) * $cvempty[0];
      $green = $cvfull[1] + 2 * (1-$load) * $cvempty[1];
      $blue  = $cvfull[2] + 2 * (1-$load) * $cvempty[2];
   }

   $red   = sprintf( "%.0f", min( $red,   255 ) );
   $green = sprintf( "%.0f", min( $green, 255 ) );
   $blue  = sprintf( "%.0f", min( $blue,  255 ) );
   return "$red/$green/$blue";
}

if( $country eq 'nocountry' ) {
   $AREA="-R5.6/15.4/47.2/55.12";
   $PROJ="-Jm10.5/51/3.0";
   $LEGEND="";
   $WATER=""; #-W2/70/70/70" # -S200/200/255" #" # -S0/0/255" #
   $BOUNDARY=""; #"-Bafg"; #geswn" #"-Ba5f1geswn"
#   COLR="relief.cpt"
   $PORTRAITMODE="-P";
   $NATBOUND="-N1/0pt/255/255/255";
#   $RES="-Di";
   $RES="-Dl";
   # textlabel offsets
   $nodeXoffset = 0; # offset for nodename x position
   $nodeYoffset = .1 # offset for nodename y position
}

if( $country eq 'germany' ) {
   $AREA="-R5.6/15.4/47.2/55.12";
   $PROJ="-Jm10.5/51/3.0";
   $LEGEND="-L6.6/54.8/51/100";
   $WATER="-S100/100/155"; #-W2/70/70/70" # -S200/200/255" #" # -S0/0/255" #
   $BOUNDARY=""; #"-Bafg"; #geswn" #"-Ba5f1geswn"
#   COLR="relief.cpt"
   $PORTRAITMODE="-P";
   $NATBOUND="-N1/2/70/70/70";
#   $RES="-Di";
   $RES="-Dl";
   # textlabel offsets
   $nodeXoffset = 0; # offset for nodename x position
   $nodeYoffset = .1 # offset for nodename y position
}

if( $country eq "austria" ) {
   $AREA="-R9.4/17.2/46.3/49.1";
   $PROJ="-Jm13.3/47.7/5.0";
   $LEGEND="-L10.4/49/47.7/100";
   $WATER="";  #"-S200/200/255";
   $BOUNDARY=""; #"-Ba3f2gESwn";
#   COLR="relief.cpt"  #$GMTHOME/share/GMT_sealand.cpt"
   $NATBOUND="-N1/2/70/70/70";
   $RES="-Di";
   $PORTRAITMODE="";
   # textlabel offsets
   $nodeXoffset = 0; # offset for nodename x position
   $nodeYoffset = .1 # offset for nodename y position
}

if( $country eq "poland" ) {
   $AREA="-R14.0/26.0/49.0/54.9";
   $PROJ="-Jm19.0/51.5/2.8";
   $LEGEND="-L15/49.3/52/100";
   $WATER="-S100/100/155";
   $BOUNDARY=""; #"-Ba3f2gESwn";
#   COLR="relief.cpt"  #$GMTHOME/share/GMT_sealand.cpt"
   $NATBOUND="-N1/2/70/70/70";
   $RES="-Di";
   $PORTRAITMODE="-P";
   # textlabel offsets
   $nodeXoffset = 0; # offset for nodename x position
   $nodeYoffset = .1 # offset for nodename y position
}

if( $country eq "eu" ) {
   $AREA="-R-12.0/29.0/36.0/62.0";
   $PROJ="-Jm17.0/47.0/0.7";
   $LEGEND="-L-4.0/38/45.0/500";
   $WATER="-S100/100/155";
   $BOUNDARY=""; #"-Ba3f2gESwn";
#   COLR="relief.cpt"  #$GMTHOME/share/GMT_sealand.cpt"
   $NATBOUND="-N1/2/70/70/70";
   $RES="-Dl";
   $PORTRAITMODE="-P";
   # textlabel offsets
   $nodeXoffset = 0; # offset for nodename x position
   $nodeYoffset = .1 # offset for nodename y position
}

if( $country eq 'us' ){
   # west/east/south/north boundary
   $AREA="-R-125.0/-68.0/25.0/48.0";
   # mercator projection, long/lat/scale
   $PROJ="-Jm-103.0/37.5/0.42";
   # print map scale at position long/lat/scale-lat/length
   $LEGEND="-L-120/27/37.5/500";
   # paint wet areas (r/g/b)
   $WATER="-S100/100/155";
   # boundary tickmarks: see psbasemap
   $BOUNDARY=""; #"-Ba3f2gESwn"
   # COLR="relief.cpt"
   # portrait mode (default: landscape)?
   $PORTRAITMODE="-P"; # -P
   # national boundaries: US (1) and state (2) boundaries, width/r/g/b
   $NATBOUND="-N1/2/50/50/50 -N2/1/200/200/200";
   # resolution: intermediate
   $RES="-Dl";
   # textlabel offsets
   $nodeXoffset = 0; # offset for nodename x position
   $nodeYoffset = .3 # offset for nodename y position
}

if( $country eq 'us-ca' ){
   $AREA="-R-126.0/-68.0/25.0/52.0";
   # mercator projection, long/lat/scale
   $PROJ="-Jm-103.0/40.5/0.42";
   # print map scale at position long/lat/scale-lat/length
   $LEGEND="-L-120/27/37.5/500";
   # paint wet areas (r/g/b)
   $WATER="-S100/100/155";
   # boundary tickmarks: see psbasemap
   $BOUNDARY=""; #"-Ba3f2gESwn"
   # COLR="relief.cpt"
   # portrait mode (default: landscape)?
   $PORTRAITMODE="-P"; # -P
   # national boundaries: US (1) and state (2) boundaries, width/r/g/b
   $NATBOUND="-N1/2/50/50/50 -N2/1/200/200/200";
   # resolution: intermediate
   $RES="-Dl";
   # textlabel offsets
   $nodeXoffset = 0; # offset for nodename x position
   $nodeYoffset = .3 # offset for nodename y position
}

my %nodedata = ();
my @linkdata;
my $maxdem = 0;
my $maxcap = 0;


# $nodeid = 0;

# read nodefile
open inNodes, "< $nodefile" or die "can't open node file: $!\n";
while ( my $line = <inNodes> )
{
   next if ( $line =~ /^\s*$/ || $line =~ /^\s*#/ );  # ignore empty and comment lines
   chomp $line;
   my @linelist = split( ';', $line ); # name;y;x;emdemand
   $node = $linelist[0];
   $nodedata{$node} = \@linelist;
   if ( $linelist[3] > $maxdem ) {
      $maxdem = $linelist[3]; # remember maximal emanating demand for radius computation
   }
}
close inNodes;

# read linkfile
open inLinks, "< $linkfile" or die "can't open link file: $!\n";
while ( my $line = <inLinks> )
{
   next if ( $line =~ /^\s*$/ || $line =~ /^\s*#/ );  # ignore empty and comment lines
   chomp $line;
   my @linelist = split( ';', $line ); # name;source;target;flow;cap
   push( @linkdata, \@linelist );
   if( $linelist[4] > $maxcap ){
      $maxcap = $linelist[4]; # remember maximal capacity for linewidth computation
   }
}
close inLinks;

# Check if all capacities are 0.
# If so probably no solution was given and all links should be drawn,
# and get a color other than $unusedcolor.
my $nosolution = 0;
if( $maxcap == 0 ){
   # set maxcap to 1 to avoid computational problems from division by 0
   $maxcap     = 1;
   $nosolution = 1;
}

# draw map bg with pascoast
open TMPSH, "> $TMP_SH_SCRIPT" or die "Can't open $TMP_SH_SCRIPT for writing: $!\n";
print TMPSH "gmtset BASEMAP_TYPE plain ANOT_FONT_SIZE 8 WANT_EURO_FONT TRUE DOTS_PR_INCH 600";
print TMPSH "pscoast -E$azimuth/$elevation $AREA $PROJ $LEGEND $WATER $BOUNDARY $PORTRAITMODE $NATBOUND $RES -A20 -K >> $tmpoutfile";

my @allLinkRecs = ("");
my $linkcount = 0;
my $linkRecordCount = 0;

# draw all links with cap = 0 first
if ( $showUnused || $nosolution )
{
   $color = ( $nosolution ? $stdcolor : $unusedcolor );
   $width = ( $nosolution ? $stdwidth : $unusedwidth );

   foreach my $linkref ( @linkdata )
   {
      my @linkinfo = @$linkref;
      my $cap      = $linkinfo[4];
      unless( $cap > 0 )
      {
         @source = @{$nodedata{$linkinfo[1]}};
         @target = @{$nodedata{$linkinfo[2]}};
         $linkrecord = "> -W$width"."p/$color \n$source[1] $source[2] $z \n$target[1] $target[2] $z \n";
         $allLinkRecs[$linkRecordCount] .= $linkrecord;

         if ( ++$linkcount >= 1000 )
         {
            push( @allLinkRecs, "");
            $linkcount = 0;
            ++$linkRecordCount;
         }
      }
   }
}

foreach my $linkref ( @linkdata )
{
   my @linkinfo = @$linkref;
   my $cap   = $linkinfo[4];
   my $width = $stdwidth;
   if ( $usewidth )
   {
      $width = $maxwidth/$maxcap * $cap;
      if( $width < $minwidth ) {
         $width = $minwidth;
      }
   }

   ( $red, $green, $blue ) = ( 0, 0, 0 );

   # set load representing color if usecolor is set
   if ( $cap > 0 )
   {
      $color = $stdcolor;
      if ( $usecolor )
      {
         $flow  = $linkinfo[3];
         $load  = 1.0 * $flow/$cap;
         $color = getColor( $colorempty, $colorfull, $load );
      }

      @source = @{$nodedata{$linkinfo[1]}};
      @target = @{$nodedata{$linkinfo[2]}};
      $linkrecord = "> -W$width"."p/$color \n$source[1] $source[2] $z \n$target[1] $target[2] $z \n";
      $allLinkRecs[$linkRecordCount] .= $linkrecord;
      if( ++$linkcount >= 1000 )
      {
         push( @allLinkRecs, "");
         $linkcount = 0;
         $linkRecordCount += 1;
      }
   }
}

# write linkdata into tmpfile and call psxyz
open TMPFILELINKS, "> $tmplinks" or die "can't open tmp link file $tmplinks: $!\n";
foreach $rec ( @allLinkRecs ) {
   $rec =~ s/\n$//;
   print TMPFILELINKS "$rec";
}
print TMPSH "psxyz $tmplinks -E$azimuth/$elevation $AREA $PROJ -O -M -K >> $tmpoutfile";
#$linkpsxy = "psxy $AREA $PROJ -W0.4p/0/0/100  -O -M -K -W1p/0/0/100 >> $tmpoutfile";
#$linkpsxy = "psxy $AREA $PROJ -W0.4p/0/0/100  -O -M -K >> $tmpoutfile";

close TMPFILELINKS;

my $allTextRecs = "";

# write edge names for each edge if writeEdgeNames is set
if ( $writeEdgeNames )
{
   foreach my $linkref ( @linkdata )
   {
      my @linkinfo = @$linkref;
      my $cap      = $linkinfo[4];
      my $cost     = $linkinfo[5];

      if ( $showUnused || $cap > 0 ) { # if unused are shown or link is used

         @source = @{$nodedata{$linkinfo[1]}};
         @target = @{$nodedata{$linkinfo[2]}};
         # create pstext parameter string
         # calc middle x/yy position between source and target node
         $midx = ( $source[2] + $target[2] ) / 2.0;
         $midy = ( $source[1] + $target[1] ) / 2.0;
         $text = substitute( $linkinfo[5] );
         # echo textdata and pipe to pstext
         $allTextRecs .= "$midy $midx $textsize $edgetextprop $text\n";
      }
   }
}

my @allNodeRecs = ("");
my $nodecount = 0;
my $nodeRecordCount = 0;

# draw nodes with psxy
while ( my ( $node, $nodeinforef ) = each %nodedata )
{
   my @nodeinfo = @$nodeinforef;
   my $rad   = $stdrad;

   if ( $userad )
   {
      $rad = $maxrad * sqrt( 1.0 * $nodeinfo[3] / $maxdem );
      # if radius < minimal radius: use minimal rad
      if ( $rad < $minrad ) {
         $rad = $minrad ;
      }
   }

   my $ncolor = ( scalar(@nodeinfo) > 4 ? $nodeinfo[4] : 0 );

   # my $noderecord = "> -G$nodecolor \n$nodeinfo[1] $nodeinfo[2] $z $rad\n";
   my $noderecord = "$nodeinfo[1] $nodeinfo[2] $z $ncolor $rad\n";
   $allNodeRecs[$nodeRecordCount] .= $noderecord;
   if ( ++$nodecount >= 1000 )
   {
      push( @allNodeRecs, "" );
      $nodecount = 0;
      ++$nodeRecordCount;
   }
}

# write a tmp colorfile
open TMPFILECOLOR, "> $tmpcolor" or die "can't open tmp color file $tmpcolor: $!\n";
my $i = 0;
foreach $cf ( @nodecarr ) {
   print TMPFILECOLOR "$i\t$cf\t", 1+$i, "\t$cf";
   ++$i;
}
close TMPFILECOLOR;

# write nodedata into tmpfile and call psxyz
open TMPFILENODES, "> $tmpnodes" or die "can't open tmp node file $tmpnodes: $!\n";
foreach $rec ( @allNodeRecs ) {
   $rec =~ s/\n$//;
   print TMPFILENODES $rec;
}
close TMPFILENODES;
# print TMPSH "psxyz $tmpnodes -E$azimuth/$elevation $AREA $PROJ -Scp -O -K -M -C$tmpcolor -W0p/$nodeborderc  >> $tmpoutfile"
   print TMPSH "psxyz $tmpnodes -E$azimuth/$elevation $AREA $PROJ -Shp -O -K -M -C$tmpcolor  >> $tmpoutfile";

# write node names for each node if writeNodeNames is set
if ( $writeNodeNames )
{
   while ( my ( $node, $nodeinforef ) = each %nodedata )
   {
      my @nodeinfo  = @$nodeinforef;
      my $nodex     = $nodeinfo[2] + $nodeXoffset ;
      my $nodey     = $nodeinfo[1] + $nodeYoffset ;
      my $text      = substitute( $nodeinfo[0] );
      $allTextRecs .= "$nodey $nodex $textsize $nodetextprop $text\n";
   }
}

# remove final newline
$allTextRecs =~ s/\n$//;

# write textdata into tmpfile and call pstext
open  TMPFILETEXT, "> $tmptext" or die "can't open tmp text file $tmptext: $!\n";
print TMPFILETEXT $allTextRecs;
close TMPFILETEXT;
print TMPSH "pstext $tmptext -Z$z -E$azimuth/$elevation $AREA $PROJ -C$textBoxSize -G$textcolor -W${textbgc}o -O -K >> $tmpoutfile";

# postprocess ps file to set title
$tmpoutfile2 = POSIX::tmpnam();
print TMPSH "awk '/^\%\%Title/ { print \$1, \"$title\"; next; } { print; }' $tmpoutfile > $tmpoutfile2";

# convert ps to eps
print TMPSH "ps2eps -q -f < $tmpoutfile2 > $outfile";
print TMPSH "rm -f $tmpoutfile $tmpoutfile2";

# call the shell with all the commands
close TMPSH;
system( "sh $TMP_SH_SCRIPT" );

# clean up
unlink $TMP_SH_SCRIPT;
unlink $tmpcolor;
unlink $tmpnodes;
unlink $tmplinks;
unlink $tmptext;

# End.