path: root/mkindex.pl

#!/usr/bin/perl

# Copyright © 2013 Guilhem Moulin <guilhem@fripost.org>
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# 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, see <http://www.gnu.org/licenses/>.

use 5.010_000;
use warnings;
use strict;

our $VERSION = "0.1, August 9, 2013";

use Getopt::Long qw/:config posix_default no_ignore_case gnu_compat
                            bundling auto_version auto_help/;
use Pod::Usage;
use File::Temp;

use Geo::GDAL;
use Geo::GDAL::Const; #qw/GRA_Bilinear/;
use Geo::OSR;
# http://geoinformatics.aalto.fi/doc/Geo-GDAL-1.9/html/

use Inline 'C';

my %config = ( progress   => Geo::GDAL->can('TermProgress_nocb') ?
                                sub { Geo::GDAL::TermProgress_nocb(@_[0,1]) } :
                                *Geo::GDAL::TermProgress
             , debug      => 0
             , resampling => &resamplings('near')
             , rename     => 's/\.[^.]+$//'
);


GetOptions( "q|quiet"  => sub { delete $config{progress} }
          , "d|debug"  => sub { $config{debug} = 1 }
          , "rename=s" => sub { $config{rename} = $_[1] }
          , "r|resampling=s"=> sub { $config{resampling} = &resamplings($_[1]) }
          , "s_srs=s" => sub { $config{s_srs} = Geo::OSR::SpatialReference::->new();
                               $config{s_srs}->SetFromUserInput( $_[1] ) }
          , "t_srs=s" => sub { $config{t_srs} = Geo::OSR::SpatialReference::->new();
                               $config{t_srs}->SetFromUserInput( $_[1] ) }
          , "man"     => sub { pod2usage(-exitstatus => 0, -verbose => 2) }
          )
    or pod2usage(2);
pod2usage(2) unless $#ARGV > 0;

my $index = pop;
my @mapset = @ARGV;

&debug( "Using mapset:", map {"  - ".$_} @mapset );
&debug( "Building index: $index" );


# Find the source and target SRS
unless ($config{s_srs}) {
    foreach my $map (@mapset) {
        &debug ("Opening '$map'");
        my $ds = Geo::GDAL::Open($map, 'ReadOnly') or exit 1;

        my $wkt = $ds->GetProjection();
        die "Error: Missing projection for '$map'.\n" unless $wkt;

        my $srs = Geo::OSR::SpatialReference::->new( $wkt );
        if ($config{s_srs}) {
            die "Error: Projection for '$map' differ from the previous input tiles'.\n"
                unless $config{s_srs}->IsSame( $srs );
        }
        else {
            $config{s_srs} = $srs;
        }
    }
}

$config{t_srs} //= $config{s_srs};
&debug ("Using source SRS: $config{s_srs}");
&debug ("Using target SRS: $config{t_srs}");


# Working directory
my $tmpdir = File::Temp::tempdir( CLEANUP => 1 );
&debug ("Using working directory: '$tmpdir'.");

# Ensure the existence of the internal driver.
my $driver = 'GTiff';
my $driver_opts = { INTERLEAVE => 'BAND'
                  , TILED      => 'YES'
                  , COMPRESS   => 'NONE'
                  , BIGTIFF    => 'IF_NEEDED'
                  , PROFILE    => 'BASELINE'
                  , TFW        => 'YES'
                  };
my $driver_type = 'Byte';
die "Error: Unknown driver: '$driver'.\n"
    unless Geo::GDAL::GetDriverByName($driver);

foreach my $map (@mapset) {
    &debug ("Opening '$map'");
    my $ds = Geo::GDAL::Open($map, 'ReadOnly') or exit 1;
    &expand_rgba( \$ds );
    &reproject(   \$ds );
}


#######################################################################

# Print debug messages.
sub debug {
    print STDERR (join '', map {"DEBUG: " .$_. "\n"} @_) if $config{debug};
}

# Create a temporary file. The file is created to defeat race
# conditions, but the file handle is closed right away.
sub tempfile {
    my ($fh, $name) = File::Temp::tempfile(DIR => $tmpdir, UNLINK => 1, @_);
    close $fh;
    &debug( "Creating file: '$name'" );
    return $name;
}

# The available resampling methods. The strings can be found
# with keys %Geo::GDAL::RESAMPLING_STRING2INT;
sub resamplings {
    my $r = shift;
    my %h = ( near => 'NearestNeighbour'
            , bilinear => 'Bilinear'
            , cubic => 'Cubic'
            , cubicspline => 'CubicSpline'
    );
    $h{$r};
}

# Expand a dataset to RGBA. This code is mostly a perl translation of pct2rgb(1).
sub expand_rgba {
    my $ds = shift;
    my $src_ds = $$ds;

    return unless $src_ds->{RasterCount} == 1;

    my $src_band = $src_ds->GetRasterBand(1);
    my $ct = $src_band->GetRasterColorTable() or return;

    &debug("Expansion to RGBA required.");

    my $dst_nb = 4;     # RGBA
    my $src_nc = 256;   # That's the best we can do, as we're using bytes

    my $lookup = [ map '', (0 .. $src_nc-1) ];
    $src_nc = $ct->GetCount() if $ct->GetCount() < $src_nc;
    for (my $i = 0; $i < $src_nc; $i++) {
        my @e = $ct->GetColorEntryAsRGB( $i );
        $lookup->[$_] .= chr $e[$_] for (0 .. $dst_nb-1);
    }

    my $dst_filename = tempfile(SUFFIX => '.tif');
    my ($xSize, $ySize) = $src_ds->Size();
    my $dst_ds = Geo::GDAL::GetDriverByName($driver)->Create( $dst_filename
                                                            , $xSize
                                                            , $ySize
                                                            , $dst_nb
                                                            , $driver_type
                                                            , $driver_opts
                                                            );

    $dst_ds->SetProjection( $src_ds->GetProjection() );
    $dst_ds->SetGeoTransform([ $src_ds->GetGeoTransform() ]);
    $dst_ds->SetGCPs( $src_ds->GetGCPs(), $src_ds->GetGCPProjection() )
        if $src_ds->GetGCPCount();

    &{$config{progress}}(0.0) if $config{progress};
    for (my $y = 0; $y < $ySize; $y++) {
        my $src_data = $src_band->ReadRaster(0, $y, $xSize, 1);

        for (my $b = 0; $b < $dst_nb; $b++) {
            local $_ = $src_data;
            # We use an inlined C function to speed up things a bit.
            # Essentially what we're doing is a splice through the
            # lookup table.
            # Using strings is much more efficients that huge perl
            # array of scalars. A pure perl version, not so much slower
            # than our C code, can be obtained using transliteration:
            #   my $src_str = join '', map chr, (0 .. $src_nc-1);
            #   eval "tr/\Q$src_str\E/\Q$lookup->[$b]\E/";
            &subst ( $_, $lookup->[$b], $xSize );
            $dst_ds->Band(1+$b)->WriteRaster(0, $y, $xSize, 1, $_);
        }

        &{$config{progress}}( (1+$y) / $ySize )
            if $config{progress};
    }

    # Erase (hence flush and close) the previous dataset
    $$ds = $dst_ds;
    return $dst_filename;
}


# Reproject the map. This code is inspired from
# http://jgomezdans.github.io/gdal_notes/reprojection.html
sub reproject {
    my $ds = shift;
    my $src_ds = $$ds;

    my ($src_xSize, $src_ySize) = $src_ds->Size();

    my $tx = Geo::OSR::CoordinateTransformation::->new( $config{s_srs}, $config{t_srs} );
    my @src_geo = $src_ds->GetGeoTransform();
    my ($ulx, $uly, $ulz) = $tx->TransformPoint( $src_geo[0], $src_geo[3] );
    my ($lrx, $lry, $lrz) = $tx->TransformPoint( $src_geo[0] + $src_geo[1]*$src_xSize,
                                               , $src_geo[3] + $src_geo[5]*$src_ySize );
    die "Error: Could not transform point.\n"
        if grep { $_ eq 'inf' } ($ulx, $uly, $ulz, $lrx, $lry, $lrz);

    my $ratio = 64; # TODO: should be computed
    my ($dst_xSize, $dst_ySize) = ( int( ($lrx-$ulx) / ($src_geo[1]*$ratio) )
                                  , int( ($lry-$uly) / ($src_geo[5]*$ratio) ) );
    &debug("Reducing map from ${src_xSize}x${src_ySize} to ${dst_xSize}x${dst_ySize}.");

    my $dst_filename = tempfile(SUFFIX => '.tif');
    my $dst_ds = Geo::GDAL::GetDriverByName($driver)->Create( $dst_filename
                                                            , $dst_xSize
                                                            , $dst_ySize
                                                            , $src_ds->{RasterCount}
                                                            , $driver_type
                                                            , $driver_opts
                                                            );

    my $t_srs_wkt = $config{t_srs}->ExportToWkt();
    $dst_ds->SetProjection( $t_srs_wkt );
    $dst_ds->SetGeoTransform([ $ulx, $src_geo[1]*$ratio, $src_geo[2]
                             , $uly, $src_geo[4],        $src_geo[5]*$ratio ]);

    my $res = Geo::GDAL::ReprojectImage( $src_ds, $dst_ds,
                                       , $config{s_srs}->ExportToWkt(), $t_srs_wkt
                                       , $config{resampling}, undef, undef
                                       , $config{progress} );
    die "Error: Failed to reproject map " . join (', ', map {"'".$_."'"}
                                                        $src_ds->GetFileList() )
                                          . ".\n"
        if $res;
    return $dst_filename;
}


__END__

__C__

void subst (unsigned char *data, unsigned char *lookup, int length) {
    int i;
    for (i = 0; i < length; i++)
        data[i] = lookup[ data[i] ];
}