summaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
-rwxr-xr-xmkindex.pl269
1 files changed, 81 insertions, 188 deletions
diff --git a/mkindex.pl b/mkindex.pl
index 2766b9f..6502019 100755
--- a/mkindex.pl
+++ b/mkindex.pl
@@ -15,18 +15,18 @@
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
+our $VERSION = "0.2, August 12, 2013";
+
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 File::Basename;
-use File::Spec;
use Inline 'C';
# http://geoinformatics.aalto.fi/doc/Geo-GDAL-1.9/html/
@@ -41,8 +41,7 @@ my %config = ( progress => Geo::GDAL->can('TermProgress_nocb') ?
*Geo::GDAL::TermProgress
, debug => 0
, rename => 's/\.[^.]+$//'
- , zoom => '1024x1024'
- , resampling => &resamplings('near')
+ , outsize => [ 1024, 0 ]
, border => 10
, bordercolor=> 'red'
, fontsize => 25
@@ -55,17 +54,10 @@ GetOptions( "q|quiet" => sub { delete $config{progress} }
, "rename=s" => sub { pod2usage(2) unless $_[1] =~ /^s(.).+\1[ixgcadlu]*$/;
$config{rename} = $_[1]
}
- , "s|size=s" => sub { if ($_[1] =~ /^(\d*\.?\d+)(%?)$/) {
- $config{zoom} = $2 ? 100. / $1 : $1;
- }
- elsif ($_[1] =~ /^\d+x\d+$/) {
- $config{zoom} = $_[1];
- }
- else {
- pod2usage(2);
- }
+ , "s|size=s" => sub { pod2usage(2) unless $_[1] =~ /^(\d+)x(\d+)$/;
+ $config{outsize} = [ $1, $2 ];
}
- , "r|resampling=s"=> sub { $config{resampling} = &resamplings($_[1]) }
+ , "r|resampling=s"=> \$config{resampling}
, "b|border=s"=> sub { if ($_[1] =~ /^(?<b>\d*)(?<bc>:.*)?$/) {
$config{border} = $+{b} if $+{b};
$config{bordercolor} = $1 if $+{bc} =~ /^:(.+)/;
@@ -108,27 +100,27 @@ foreach my $map (@mapset) {
&debug ("Opening '$name'");
my $ds = Geo::GDAL::Open($name, 'ReadOnly') or exit 1;
- $map->{wkt} = $ds->GetProjection() || $config{s_srs}->ExportToWkt();
+ $map->{wkt} = $ds->GetProjection();
+ my $srs;
unless ($config{s_srs}) {
- die "Error: Missing projection for '$name'.\n" unless $map->{wkt};
+ die "Error: Missing projection for '$name'.\n"
+ unless $map->{wkt} or $config{s_srs};
+ $map->{wkt} ||= $config{s_srs}->ExportToWkt();
- $map->{srs} = Geo::OSR::SpatialReference::->new( $map->{wkt} );
+ $srs = Geo::OSR::SpatialReference::->new( $map->{wkt} );
if ($s_srs) {
- unless ($s_srs->IsSame($map->{srs})) {
- my $msg = "Projection for '$name' differs from the previous input tiles'.";
- die "Error: " .$msg. "\n" unless $config{t_srs};
- warn "WARN: " .$msg. "\n";
- }
+ die "Error: Projection for '$name' differs from the previous input tiles'\n"
+ unless $config{t_srs} or $s_srs->IsSame($srs);
}
else {
- $s_srs = $map->{srs};
+ $s_srs = $srs;
}
}
- &debug( "Using source SRS: ". ($config{s_srs} // $map->{srs}) );
+ &debug( "Using source SRS: ". ($config{s_srs} // $srs)->ExportToProj4() );
my $tx = Geo::OSR::CoordinateTransformation::->new(
- $config{s_srs} // $map->{srs},
+ $config{s_srs} // $srs,
$config{t_srs} // $config{s_srs} // $s_srs
);
die "Error: Could compute transformation between source and target SRS.\n"
@@ -151,53 +143,7 @@ foreach my $map (@mapset) {
}
$config{t_srs} //= $config{s_srs};
-$config{s_wkt} = $config{s_srs}->ExportToWkt() if $config{s_srs};
-$config{t_wkt} = $config{t_srs}->ExportToWkt();
-&debug ("Using target SRS: $config{t_srs}");
-
-
-# Compute the optimal zoom factor if the user is expecting a size.
-if ($config{zoom} =~ /^(\d+)x(\d+)$/) {
- my ($xSize, $ySize) = ($1,$2);
-
- # Get the maximum resolution, and corners of the final tile;
- my ($xRes, $yRes) = (0,0);
- my ($ulx, $uly, $lrx, $lry);
- foreach (@mapset) {
- # TODO: What to do with multiple tiles?
- &debug ("Opening '$_->{filename}'");
- my $src_ds = Geo::GDAL::Open($_->{filename}, 'ReadOnly') or exit 1;
- my $dst_ds = Geo::GDAL::AutoCreateWarpedVRT( $src_ds
- , $config{s_wkt} // $_->{wkt}
- , $config{t_wkt} );
- $_->{t_transform} = [ $dst_ds->GetGeoTransform() ];
- $config{ratio} = (0. + $dst_ds->{RasterXSize}) / $dst_ds->{RasterYSize};
- &debug( "Using ratio: $config{ratio}" );
- undef $src_ds;
- undef $dst_ds;
-
- $xRes = $_->{s_transform}->[1] if abs $xRes < abs $_->{s_transform}->[1];
- $yRes = $_->{s_transform}->[5] if abs $yRes < abs $_->{s_transform}->[5];
-
- $ulx = $_->{t_corners}->[0]->[0]
- unless defined $ulx and $ulx < $_->{t_corners}->[0]->[0];
- $uly = $_->{t_corners}->[0]->[1]
- unless defined $uly and $uly > $_->{t_corners}->[0]->[1];
- $lrx = $_->{t_corners}->[3]->[0]
- unless defined $lrx and $lrx > $_->{t_corners}->[3]->[0];
- $lry = $_->{t_corners}->[3]->[1]
- unless defined $lry and $lry < $_->{t_corners}->[3]->[1];
- }
-
- my $xRatio = (0.+$lrx-$ulx) / ($xRes*$xSize);
- my $yRatio = (0.+$lry-$uly) / ($yRes*$ySize);
-
- $config{zoom} = $xRatio < $yRatio ? $yRatio : $xRatio;
- map { $_->{t_transform}->[1] = $_->{s_transform}->[1] * $config{zoom}
- ; $_->{t_transform}->[5] = $_->{s_transform}->[5] * $config{zoom}
- } @mapset;
-}
-&debug( "Using zoom: $config{zoom}" );
+&debug ("Using target SRS: ". $config{t_srs}->ExportToProj4()) if $config{t_srs};
# Working directory
@@ -217,13 +163,8 @@ my $driver_type = 'Byte';
die "Error: Unknown driver: '$driver'.\n"
unless Geo::GDAL::GetDriverByName($driver);
-foreach my $map (@mapset) {
- &debug ("Opening '$map->{filename}'");
- my $ds = Geo::GDAL::Open($map->{filename}, 'ReadOnly') or exit 1;
- &expand_rgba( \$ds, $map );
- my $f = &reproject( \$ds, $map );
- &annotate( $f, $map );
-}
+map &expand_rgba($_), @mapset;
+&annotate( &merge(), $index );
#######################################################################
@@ -242,23 +183,12 @@ sub tempfile {
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 $map = shift;
- my $src_ds = $$ds;
+
+ &debug ("Opening '$map->{filename}'");
+ my $src_ds = Geo::GDAL::Open($map->{filename}, 'ReadOnly') or exit 1;
return unless $src_ds->{RasterCount} == 1;
@@ -277,9 +207,9 @@ sub expand_rgba {
$lookup->[$_] .= chr $e[$_] for (0 .. $dst_nb-1);
}
- my $dst_filename = &tempfile(SUFFIX => '.tif');
+ $map->{filename} = &tempfile(SUFFIX => '.tif');
my ($xSize, $ySize) = $src_ds->Size();
- my $dst_ds = Geo::GDAL::GetDriverByName($driver)->Create( $dst_filename
+ my $dst_ds = Geo::GDAL::GetDriverByName($driver)->Create( $map->{filename}
, $xSize
, $ySize
, $dst_nb
@@ -316,129 +246,92 @@ sub expand_rgba {
# Erase (hence flush and close) the previous dataset
undef $src_ds;
- $$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 $map = shift;
- my $src_ds = $$ds;
-
- my ($ul,$lr) = @{$map->{t_corners}}[0,3];
-
- # TODO: it might be xSize that's determined from ySize
- my $dst_xSize = int( ($lr->[0]-$ul->[0]) / $map->{t_transform}->[1] );
- my $dst_ySize = int( $dst_xSize / $config{ratio} );
-
- my ($src_xSize, $src_ySize) = $src_ds->Size();
- &debug("Reducing map from ${src_xSize}x${src_ySize} to ${dst_xSize}x${dst_ySize}.");
-
+# Build the index. We use gdalwarp(1) because it does some clever
+# computation to estimate the size of the reprojected tiles, and the
+# GDALSuggestedWarpOutput() is not available in the Perl (or Python)
+# API.
+sub merge {
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
- );
- $dst_ds->SetProjection( $config{t_wkt} );
- $dst_ds->SetGeoTransform( $map->{t_transform} );
-
- my $res = Geo::GDAL::ReprojectImage( $src_ds, $dst_ds,
- , $map->{wkt}, $config{t_wkt}
- , $config{resampling}, undef, undef
- , $config{progress} );
- die "Error: Failed to reproject map " . join (', ', map {"'".$_."'"}
- $src_ds->GetFileList() )
- . ".\n"
- if $res;
-
- # Erase (hence flush and close) the previous dataset
- undef $src_ds;
+ &debug("Wrapping into file '$dst_filename'...");
+
+ my @opts;
+ push @opts, '-s_srs', $config{s_srs}->ExportToProj4() if $config{s_srs};
+ push @opts, '-t_srs', $config{t_srs}->ExportToProj4() if $config{t_srs};
+ push @opts, '-r', $config{resampling} if $config{resampling};
+ push @opts, '-ts', $config{outsize}->[0], $config{outsize}->[1];
+ push @opts, '-q' unless $config{progress};
+ push @opts, '-of', $driver, map { ('-co', $_.'='.$driver_opts->{$_}) }
+ (keys %$driver_opts);
+ my @src_filenames = map {$_->{filename}} @mapset;
+
+ system 'gdalwarp', @opts, @src_filenames, $dst_filename;
+ $? == 0 or die "gdalwarp failed.\n";
return $dst_filename;
}
-# Annotate the map (add border & caption).
+# Annotate the maps (add border & caption).
sub annotate {
- my $src_filename = shift;
- my $map = shift;
-
- my $img = Image::Magick::->new();
+ my $filename = shift;
- $img->Read($src_filename);
- my ($xSize, $ySize) = $img->Get(qw/width height/);
+ my $ds = Geo::GDAL::Open($filename, 'ReadOnly') or exit 1;
+ my ($ok,@invt) = Geo::GDAL::InvGeoTransform([ $ds->GetGeoTransform() ]);
+ undef $ds;
+ die "Error: Couldn't inverse affine transformation\n" unless $ok;
- &debug( "Border: $config{border}:$config{bordercolor}" );
- &debug( "Font: $config{fontsize}:". ($config{fontcolor} // $config{bordercolor}).
- ":$config{font}" );
- &debug( "Text: $map->{text}");
- # To get the list of fonts:
- # &debug( "Supported fonts:", map {"\t".$_} $img->QueryFont() );
+ my $img = Image::Magick::->new();
+ $img->Read( $filename );
+ my ($width, $height) = $img->Get(qw/width height/);
my @draw = ( fill => 'none'
, stroke => $config{bordercolor}
, strokewidth => $config{border}
+ , primitive => 'polygon'
);
my @ann = ( fill => $config{fontcolor} // $config{bordercolor}
, strokewidth => 0
, font => $config{font}
, pointsize => $config{fontsize}
- , gravity => 'Center'
- , text => $map->{text}
);
- my @points;
+ foreach my $map (@mapset) {
+ # In case the projection have changed, the map may have been rotated,
+ # so we draw a polygon instead of a rectangle. Also, we're using
+ # a mask to avoid the text to be out of bounds.
- if ($config{t_srs}->IsSame($config{s_srs} // $map->{srs})) {
- push @draw, primitive => 'rectangle';
- push @points, [ 0, 0 ], [ $xSize-1, $ySize-1 ];
- }
- else {
- # The projection changed. Since the map might have been rotated,
- # we draw a polygon instead of a rectangle. Also, instead of
- # adding the border *then* reproject we use a mask in order to
- # have nice anti-aliasing and better looking fonts.
+ my @points = map {[ Geo::GDAL::ApplyGeoTransform(\@invt, @$_[0,1] ) ]}
+ @{$map->{t_corners}}[0,1,3,2];
+ my $points = join (' ', map {join ',', @$_} @points);
+ &debug( "Polygon: " . $points );
+ $img->Draw( @draw, points => $points );
- my ($success,@invgeo) = Geo::GDAL::InvGeoTransform($map->{t_transform});
- die "Error: Couldn't inverse affine transformation\n" unless $success;
+ # Find the middle of the polygon; We add a Y-offset to be
+ # vertically centered.
+ my ($x,$y) = (0,0);
+ map { $x+= $_->[0] / (1.+$#points); $y+= $_->[1] / (1.+$#points); }
+ @points;
+ push @ann, x => $x, y => $y + $config{fontsize}/2., align=> 'Center';
- push @draw, primitive => 'polygon';
- push @points, map {[ Geo::GDAL::ApplyGeoTransform(\@invgeo, @$_[0,1] ) ]}
- @{$map->{t_corners}}[0,1,3,2];
- my $mask = Image::Magick::->new( size => "${xSize}x${ySize}" );
+ my $text = Image::Magick::->new( size => "${width}x${height}" );
+ $text->ReadImage('xc:none');
+
+ my $mask = Image::Magick::->new( size => "${width}x${height}" );
$mask->ReadImage('xc:black');
$mask->Draw( primitive => 'polygon'
- , fill => 'white'
- , points => join (' ', map {join ',', @$_} @points) );
- $img->Mask( mask => $mask );
+ , strokewidth => 0
+ , fill => 'white'
+ , points => $points );
+ $text->Mask( mask => $mask );
undef $mask;
- }
- &debug( ucfirst {@draw}->{primitive} .": ".
- join (' ', map {join ',', @$_} @points) );
- $img->Draw( @draw, points => join (' ', map {join ',', @$_} @points) );
- $img->Annotate( @ann );
- my $dst_filename = &tempfile(SUFFIX => '.tif');
- $img->Write($dst_filename);
-
- # The georeference didn't change when annotating. Hence we simply
- # link the former world file.
- my ($src_tfw, $dst_tfw) = map { my ($name,$path) = fileparse($_, qr/\.tiff?$/);
- File::Spec->catfile($path, $name.'.tfw');
- }
- ($src_filename, $dst_filename);
-
- if (-f $src_tfw and not -f $dst_tfw) {
- link $src_tfw, $dst_tfw;
- &debug( "Linking worldfiles: $src_tfw -> $dst_tfw" );
+ $text->Annotate( @ann, text => $map->{text} );
+ $img->Composite( image => $text, compose => 'Over' );
+ undef $text;
}
-
+ $img->Write($index);
undef $img;
- $map->{thumbnail} = $dst_filename;
}