From 549390f42b7cce09ad59c6882258b73f593cbbd6 Mon Sep 17 00:00:00 2001 From: Guilhem Moulin Date: Sat, 10 Aug 2013 18:05:20 +0200 Subject: Initial commit --- mkindex.pl | 263 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 263 insertions(+) create mode 100755 mkindex.pl (limited to 'mkindex.pl') diff --git a/mkindex.pl b/mkindex.pl new file mode 100755 index 0000000..cfc6ac3 --- /dev/null +++ b/mkindex.pl @@ -0,0 +1,263 @@ +#!/usr/bin/perl + +# Copyright © 2013 Guilhem Moulin +# +# 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 . + +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] ]; +} -- cgit v1.2.3