From 43533160c8cfad0f4d6cd7084f32fe950418e244 Mon Sep 17 00:00:00 2001 From: Guilhem Moulin Date: Sat, 17 Aug 2013 20:08:56 +0200 Subject: Abitility to use a different mapset. --- mkindex.pl | 155 ++++++++++++++++++++++++++++++++++++++++++------------------- 1 file 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 [B<-q>] [B<-d>] [B<--rename> I] [B<-s> I] [B<-r> I] [B<-b> I] [B<-f> I] [B<--s_srs> I] -[B<--t_srs> I] I [... I] I +[B<--t_srs> I] [-o I] [B<--crop>] I [... I] I B B<--man> @@ -95,6 +95,17 @@ If one of I or I 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, B<--overview=>I + +Use I 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 The source spatial reference set. If the input Is 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 and C files as overview maps for +the mapset C. 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] =~ /^(?\d*)(?:.*)?$/) { $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; -- cgit v1.2.3