summaryrefslogtreecommitdiffstats
path: root/mkindex.pl
diff options
context:
space:
mode:
Diffstat (limited to 'mkindex.pl')
-rwxr-xr-xmkindex.pl263
1 files changed, 263 insertions, 0 deletions
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 <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] ];
+}
generated by cgit v1.2.3 (git 2.25.1) at 2024-04-20 11:05:53 +0000