#!/usr/bin/perl -CSAL

use strict;
use warnings;


#################################################
# Command line parsing
#

if( (@ARGV != 1 && @ARGV != 2) || $ARGV[0] !~ /^(?:multipoint|polyline|polygon)$/i ) {
    print "Unknown shape type $ARGV[0].\n"
        if @ARGV && $ARGV[0] !~ /^(?:multipoint|polyline|polygon)$/i;
    print <<EOF;
usage: txt2shp.pl <type> [ <input file> ]
Generates a shapefile from a text file containing space-separated
longitude/latitude coordinates, one pair per line.  <type> may be multipoint,
polyline or polygon.  Empty lines in the input file separate the points of
different objects/records.  Lines starting with # are ignored.
EOF
    exit;
}

my $type= lc($ARGV[0]);

my $outfile= $ARGV[1] // "txt2shp";
$outfile =~ s/(?:\.[^.]*)?$/.shp/;

my $inhandle;

if( defined $ARGV[1] ) {
    open $inhandle, "<", $ARGV[1] or die "$!";
}
else {
    $inhandle= *STDIN{IO};
}


#################################################
# Input data parsing
#

my %overall;

sub addcoords
{
    my ($o, $x, $y)= @_;

    if( ! %$o ) {
        $$o{points}= [ [ $x, $y ] ];
        $$o{xmin}= $$o{xmax}= $x;
        $$o{ymin}= $$o{ymax}= $y;
    }
    else {
        push @{$$o{points}}, [ $x, $y ];
        $$o{xmin}= $x if $x < $$o{xmin};
        $$o{xmax}= $x if $x > $$o{xmax};
        $$o{ymin}= $y if $y < $$o{ymin};
        $$o{ymax}= $y if $y > $$o{ymax};
    }
    if( ! %overall ) {
        $overall{xmin}= $overall{xmax}= $x;
        $overall{ymin}= $overall{ymax}= $y;
    }
    else {
        $overall{xmin}= $x if $x < $overall{xmin};
        $overall{xmax}= $x if $x > $overall{xmax};
        $overall{ymin}= $y if $y < $overall{ymin};
        $overall{ymax}= $y if $y > $overall{ymax};
    }
    ++$overall{pointcount};
}

my $numre= qr/-?(?:\d+(?:\.\d*)?|\.\d+)(?:e[-+]?\d+)?/i ;

my @objects= ( {} );

while( defined($_ = <$inhandle>) ) {
    s/^\s+//;
    next if /^#/;
    chomp;
    if( /^$/ ) {
        push @objects, {} if %{$objects[-1]};
        next;
    }
    my @coords= split /\s+/;
    @coords >= 2 && $coords[0] =~ /$numre/ && $coords[1] =~ /$numre/
        or die "Non-numerical data on input line $.!";
    addcoords($objects[-1], @coords);
}

pop @objects unless %{$objects[-1]};

close $inhandle if defined $ARGV[1];


#################################################
# Constants and computations for shape file headers
#

# Formats for pack() and constants for header fields:
# Warning: pack supports floating point data in native endianness only!
# Doubles should be little endian.
my $shpheadfmt= "N7V2d8";
my $shpmagic= 9994;
my $shpversion= 1000;
my %shapetype= ( nullshape => 0, null => 0,
    point => 1, polyline => 3, polygon => 5, multipoint => 8,
    pointz => 11, polylinez => 13, polygonz => 15, multipointz => 18,
    pointm => 21, polylinem => 23, polygonm => 25, multipointm => 28,
    multipatch => 31 );
my $recheadfmt= "N2";
my $multiptheadfmt= "Vd4V";
my $multiptheadsz= 40;
my $polyheadfmt= "Vd4VV";
my $polyheadsz= 44;

# Lengths of parts of the shape file in 16-bit words:
# Length of shape header, i.e. record content without arrays
my $shapeheadlen= int(($type eq "multipoint" ? $multiptheadsz : $polyheadsz) / 2);
# Total length of part and point arrays of all records
my $totaldatasize= 8 * $overall{pointcount};
$totaldatasize += 2 * @objects if $type =~ /^poly/;
# Total shape file size
my $totalsize= 50 + (4 + $shapeheadlen) * @objects + $totaldatasize;


#################################################
# Write the main shape file:
#

open OUT, "> :raw :bytes", $outfile or die "Could not create $outfile: $!";

print OUT pack($shpheadfmt, $shpmagic, 0,0,0,0,0, $totalsize, $shpversion,
               $shapetype{$type}, $overall{xmin}, $overall{ymin},
               $overall{xmax}, $overall{ymax}, 0,0,0,0);

# Consecutive record number for record headers:
my $recnr= 1;

for (@objects) {
    $$_{recsize}= 4 + $shapeheadlen + 8 * @{$$_{points}};
    $$_{recsize} += 2 if $type =~ /^poly/;
    print OUT pack($recheadfmt, $recnr, $$_{recsize} - 4);
    if( $type eq "multipoint" ) {
        print OUT pack($multiptheadfmt, $shapetype{$type}, $$_{xmin}, $$_{ymin}, $$_{xmax}, $$_{ymax}, 0+@{$$_{points}});
    }
    else {
        print OUT pack($polyheadfmt . "V", $shapetype{$type}, $$_{xmin}, $$_{ymin}, $$_{xmax}, $$_{ymax}, 1, ($type eq "polyline"? 0:1)+@{$$_{points}}, 0);
    }
    for my $pt (@{$$_{points}}) {
        print OUT pack("d2", $$pt[0], $$pt[1]);
    }
    if( $type eq "polygon" ) {
        print OUT pack("d2", $$_{points}[0][0], $$_{points}[0][1]);
    }
    ++$recnr;
}

close OUT;


#################################################
# Write the index file:
#

$outfile =~ s/\.shp$/.shx/i;

open OUT, "> :raw :bytes", $outfile or die "Could not create $outfile: $!";

print OUT pack($shpheadfmt, $shpmagic, 0,0,0,0,0, 50+4*@objects, $shpversion,
               $shapetype{$type}, $overall{xmin}, $overall{ymin},
               $overall{xmax}, $overall{ymax}, 0,0,0,0);
my $offset= 50;
for (@objects) {
    print OUT pack("NN", $offset, $$_{recsize} - 4);
    $offset += $$_{recsize};
}

close OUT;


#################################################
# Write the dbase file:
#

$outfile =~ s/\.shx$/.dbf/i;

open OUT, "> :raw :bytes", $outfile or die "Could not create $outfile: $!";

my @time= localtime();
my $idwidth= 11;  # log($nrecords) / log(10) + 2;
# Only 11 seems to work for unknown reasons
print OUT pack("C4VvvV5", 3, $time[5], $time[4] + 1, $time[3],
                          0+@objects, 65, 12, 0, 0,0,0,0);
print OUT pack("a11AVCCvVVVC", "ID", "N", 0, $idwidth, 0, 0, 0,0,0, 13);
for (1..0+@objects) {
    printf OUT " \%*u", $idwidth, $_;
    # has to be one space longer than nominal length
}

close OUT;

