diff options
| -rwxr-xr-x | mkindex.pl | 269 |
1 files changed, 81 insertions, 188 deletions
@@ -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; } |