summaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
Diffstat
-rwxr-xr-xmkindex.pl155
1 files changed, 107 insertions, 48 deletions
diff --git a/mkindex.pl b/mkindex.pl
index b787ead..467a296 100755
--- a/mkindex.pl
+++ b/mkindex.pl
@@ -29,7 +29,7 @@ mkindex.pl - Create an index file from a mapset
B<mkindex.pl> [B<-q>] [B<-d>] [B<--rename> I<perlexpr>] [B<-s> I<geometry>]
[B<-r> I<method>] [B<-b> I<border>] [B<-f> I<font>] [B<--s_srs> I<SRS>]
-[B<--t_srs> I<SRS>] I<map> [... I<map>] I<index>
+[B<--t_srs> I<SRS>] [-o I<map>] [B<--crop>] I<map> [... I<map>] I<index>
B<mkindex.pl> B<--man>
@@ -95,6 +95,17 @@ If one of I<width> or I<font> is left empty, the corresponding default
is used (but the other can be specified regardless, for instance in
C<::DejaVu-Sans-Book>).
+=item B<-o> I<map>, B<--overview=>I<map>
+
+Use I<map> as the overview map. This option can be used serveral times
+to use multiple (warped) overview maps.
+
+=item B<--crop>
+
+Alternative overview map may be much bigger than the actual mapset. Use
+B<--crop> to trim away the useless part. Note that the resulting index file
+may be smaller than wanted.
+
=item B<--s_srs=>I<SRS>
The source spatial reference set. If the input I<map>s include a SRS, it is
@@ -143,6 +154,12 @@ opaque blue.
Use DejaVu Sans with a pointsize of 75 for the captions. Also, display
only the non-redundant part of the filenames.
+=item mkindex.pl -o overview1.tif -o overview2.tif --crop maps*.tif index.tif
+
+Use the C<overview1.tif> and C<overview2.tif> files as overview maps for
+the mapset C<maps*.tif>. Also, crop the parts of the index that are not
+in the mapset.
+
=back
=head1 ENVIRONMENT
@@ -181,6 +198,8 @@ use Getopt::Long qw/:config posix_default no_ignore_case gnu_compat
use Pod::Usage;
use File::Temp;
use File::Basename;
+use File::Spec;
+use POSIX qw/floor ceil/;
use Inline 'C';
# http://geoinformatics.aalto.fi/doc/Geo-GDAL-1.9/html/
@@ -193,9 +212,9 @@ use Image::Magick;
my %config = ( progress => sub {
# Geo::GDAL::TermProgress_nocb writes to STDOUT
open my $x, '>&', \*STDOUT or die "Error: Can't dup STDOUT: $!";
- open STDOUT, '>&', STDERR or die "Error: Can't dup STDERR: $!";
+ open *STDOUT, '>&', \*STDERR or die "Error: Can't dup STDERR: $!";
Geo::GDAL::TermProgress_nocb(@_[0,1]);
- open STDOUT, '>&', $x or die "Error: Can't dup X: $!";
+ open *STDOUT, '>&', $x or die "Error: Can't dup X: $!";
close $x;
}
, debug => 0
@@ -221,6 +240,8 @@ GetOptions( "q|quiet" => sub { delete $config{progress} }
$config{outsize} = [ $1, $2 ];
}
, "r|resampling=s"=> \$config{resampling}
+ , "o|overview=s@" => \$config{overviews}
+ , "crop" => \$config{crop}
, "b|border=s"=> sub { if ($_[1] =~ /^(?<b>\d*)(?<bc>:.*)?$/) {
$config{border} = $+{b} if ($+{b} // '');
$config{bordercolor} = $1 if ($+{bc} // '') =~ /^:(.+)/;
@@ -246,7 +267,6 @@ GetOptions( "q|quiet" => sub { delete $config{progress} }
)
or pod2usage(2);
# TODO: -n|--numeric: put numbers instead of text
-# TODO: -o|--overview: use a different overview mapset
# TODO: when the input is RGB (not RGBA), we want add an alpha channel
# but keep the nodata value as matte within the map bounds.
# This seems to be possible with a clever use of
@@ -259,14 +279,20 @@ my @mapset = map {{ filename => $_ }} @ARGV;
$config{fontcolor} //= $config{bordercolor};
&debug( "Using mapset:", map {" - ".$_->{filename}} @mapset );
+if ($config{overviews}) {
+ $config{overviews} = [ map {{ filename => $_ }} @{$config{overviews}} ];
+ &debug( "Using overview:", map {" - ".$_->{filename}} @{$config{overviews}} );
+}
+else {
+ $config{overviews} = \@mapset;
+}
&debug( "Building index: $index" );
-# Find the source and target SRS
-# We also store the corders points and the affine transformation, to
-# avoid to compute it again later on.
+# Find the source and target SRS, and ensure the source files are in a
+# decent format.
my $s_srs;
-foreach my $map (@mapset) {
+foreach my $map (@{$config{overviews}}) {
my $name = $map->{filename};
&debug ("Opening '$name'");
@@ -308,50 +334,69 @@ foreach my $map (@mapset) {
). "\n";
}
- my $srs;
- $map->{wkt} = $ds->GetProjection();
- unless ($config{s_srs}) {
- die "Error: Missing projection for '$name'.\n"
- unless $map->{wkt} or $config{s_srs};
- $map->{wkt} ||= $config{s_srs}->ExportToWkt();
-
- $srs = Geo::OSR::SpatialReference::->new( $map->{wkt} );
+ if ($config{s_srs}) {
+ # The 's_srs' option takes precedence.
+ $map->{wkt} = $config{s_srs}->ExportToWkt();
+ }
+ else {
+ # Croak if the source files don't have a valid projection.
+ $map->{wkt} = $ds->GetProjection()
+ or die "Error: Missing projection for '$name'.\n";
+
+ # Croak unless either 't_srs' is set, or all the map are in the
+ # same projection.
+ my $srs = Geo::OSR::SpatialReference::->new( $map->{wkt} );
if ($s_srs) {
- die "Error: Projection for '$name' differs from the previous input tiles'\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 = $srs;
}
}
+}
+
+$config{t_srs} //= $config{s_srs} // $s_srs;
+&debug ("Using target SRS: ". $config{t_srs}->ExportToProj4());
+
- &debug( "Using source SRS: ". ($config{s_srs} // $srs)->ExportToProj4() );
- my $tx = Geo::OSR::CoordinateTransformation::->new(
- $config{s_srs} // $srs,
- $config{t_srs} // $config{s_srs} // $s_srs
- );
- die "Error: Could compute transformation between source and target SRS.\n"
- unless $tx;
+# Store the corners of each tile, an the text to display.
+foreach my $map (@mapset) {
+ my $name = $map->{filename};
+
+ &debug ("Opening '$name'");
+ my $ds = Geo::GDAL::Open($name, 'ReadOnly') or exit 1;
- $map->{s_transform} = $ds->GetGeoTransform();
- $map->{t_corners} = [ map {[ Geo::GDAL::ApplyGeoTransform($map->{s_transform}, @$_) ]}
+ my $geot = [ $ds->GetGeoTransform() ];
+ $map->{t_corners} = [ map [ Geo::GDAL::ApplyGeoTransform($geot, @$_) ],
( [ 0, 0]
, [ 0, $ds->{RasterYSize}-1 ]
, [ $ds->{RasterXSize}-1, 0 ]
, [ $ds->{RasterXSize}-1, $ds->{RasterYSize}-1 ]
)
];
- $tx->TransformPoints( $map->{t_corners} );
- die "Error: Could not transform point.\n"
- if grep !/^[-+]?\d*\.?\d+$/, (map @$_, @{$map->{t_corners}});
+ my $srs = $config{s_srs};
+ unless ($srs) {
+ my $wkt = $ds->GetProjection()
+ or die "Error: Missing projection for '$name'.\n";
+ $srs = Geo::OSR::SpatialReference::->new( $wkt );
+ }
+ &debug( "Using source SRS: ". $srs->ExportToProj4() );
+
+ unless ($config{t_srs}->IsSame($srs)) {
+ my $tx = Geo::OSR::CoordinateTransformation::->new( $srs, $config{t_srs} );
+ die "Error: Could compute transformation between source and target SRS.\n"
+ unless $tx;
+
+ $tx->TransformPoints( $map->{t_corners} );
+ die "Error: Could not transform point.\n"
+ if grep !/^[-+]?\d*\.?\d+$/, (map @$_, @{$map->{t_corners}});
+ }
{ local $_ = basename $name; eval $config{rename}; $map->{text} = $_; }
}
-$config{t_srs} //= $config{s_srs};
-&debug ("Using target SRS: ". $config{t_srs}->ExportToProj4()) if $config{t_srs};
-
# Working directory
my $tmpdir = File::Temp::tempdir( CLEANUP => 1 );
@@ -366,7 +411,7 @@ die "Error: Unknown driver: '$driver'.\n"
$driver_opts->{$_} //= $gdal_opts->{$_} for (keys %$gdal_opts);
my $driver_ext = Geo::GDAL::GetDriverByName($driver)->Extension();
-map &expand_rgba($_), @mapset;
+map &expand_rgba($_), @{$config{overviews}};
my $mosaic = &merge();
&annotate( $mosaic );
@@ -390,20 +435,22 @@ sub tempfile {
# Expand a dataset to RGBA. This code is nearly a perl clone of pct2rgb(1).
sub expand_rgba {
my $map = shift;
+ my $name = $map->{filename};
- &debug ("Opening '$map->{filename}'");
- my $src_ds = Geo::GDAL::Open($map->{filename}, 'ReadOnly') or exit 1;
+ &debug ("Opening '$name'");
+ my $src_ds = Geo::GDAL::Open($name, 'ReadOnly') or exit 1;
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.");
+ $map->{filename} = &tempfile(SUFFIX => '.'.$driver_ext);
if ($config{progress}) {
- print STDERR "Processing input file $map->{filename}.\n";
+ print STDERR "Processing input file $name.\n";
&{$config{progress}}(0.0);
}
- &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
@@ -415,7 +462,6 @@ sub expand_rgba {
$lookup->[$_] .= chr $e[$_] for (0 .. $dst_nb-1);
}
- $map->{filename} = &tempfile(SUFFIX => '.'.$driver_ext);
my ($xSize, $ySize) = $src_ds->Size();
my $dst_ds = Geo::GDAL::GetDriverByName($driver)->Create( $map->{filename}
, $xSize
@@ -424,9 +470,8 @@ sub expand_rgba {
, $driver_type
, $driver_opts
);
-
$dst_ds->SetProjection( $map->{wkt} );
- $dst_ds->SetGeoTransform( $map->{s_transform} );
+ $dst_ds->SetGeoTransform([ $src_ds->GetGeoTransform() ]);
$dst_ds->SetGCPs( $src_ds->GetGCPs(), $src_ds->GetGCPProjection() )
if $src_ds->GetGCPCount();
@@ -481,15 +526,14 @@ sub merge {
delete $driver_opts->{COMPRESS} if $driver_opts->{COMPRESS};
push @opts, '-of', $driver, map { ('-co', $_.'='.$driver_opts->{$_}) }
(keys %$driver_opts);
- my @src_filenames = map {$_->{filename}} @mapset;
-
+ my @src_filenames = map $_->{filename}, @{$config{overviews}};
{
# gdalwarp writes to STDOUT; we add >&2
open my $x, '>&', \*STDOUT or die "Error: Can't dup STDOUT: $!";
- open STDOUT, '>&', STDERR or die "Error: Can't dup STDERR: $!";
+ open *STDOUT, '>&', \*STDERR or die "Error: Can't dup STDERR: $!";
system 'gdalwarp', @opts, @src_filenames, $dst_filename;
$? == 0 or die "Error: gdalwarp failed.\n";
- open STDOUT, '>&', $x or die "Error: Can't dup X: $!";
+ open *STDOUT, '>&', $x or die "Error: Can't dup X: $!";
close $x;
}
return $dst_filename;
@@ -554,6 +598,7 @@ sub annotate {
$mask->Set( colorspace => 'RGB' );
}
+ my ($lx0,$uy0,$rx0,$ly0);
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.
@@ -571,10 +616,10 @@ sub annotate {
# Get the extremes of this polygon.
my ($lx,$uy,$rx,$ly);
foreach (@points) {
- $lx = int $_->[0] unless defined $lx and $lx < int $_->[0];
- $uy = int $_->[1] unless defined $uy and $uy < int $_->[1];
- $rx = int $_->[0] unless defined $rx and $rx > int $_->[0];
- $ly = int $_->[1] unless defined $ly and $ly > int $_->[1];
+ $lx = floor $_->[0] unless defined $lx and $lx < floor $_->[0];
+ $uy = floor $_->[1] unless defined $uy and $uy < floor $_->[1];
+ $rx = ceil $_->[0] unless defined $rx and $rx > ceil $_->[0];
+ $ly = ceil $_->[1] unless defined $ly and $ly > ceil $_->[1];
};
# The region where to confine the caption.
my $cell = ($rx-$lx).'x'.($ly-$uy).'+'.$lx.'+'.$uy;
@@ -586,6 +631,13 @@ sub annotate {
, @ann
, geometry => $cell
, text => $map->{text} );
+
+ if ($config{crop}) {
+ $lx0 = $lx unless defined $lx0 and $lx0 < $lx;
+ $uy0 = $uy unless defined $uy0 and $uy0 < $uy;
+ $rx0 = $rx unless defined $rx0 and $rx0 > $rx;
+ $ly0 = $ly unless defined $ly0 and $ly0 > $ly;
+ }
}
# Change the colorspace (now to non-linear gray), and turn off
@@ -601,6 +653,13 @@ sub annotate {
$raster->Composite( image => $rgb, compose => 'Over' );
undef $rgb;
+ # Throw away the outside.
+ if ($config{crop}) {
+ my $crop = ($rx0-$lx0).'x'.($ly0-$uy0).'+'.$lx0.'+'.$uy0;
+ &debug( "Crop: $crop" );
+ $raster->Crop( geometry => $crop );
+ }
+
binmode STDOUT if $index =~ /^(?:[[:alnum:]]*:)?-$/; # Useful for pipes
$raster->Write($index);
undef $raster;
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-21 12:21:12 +0000