diff options
| -rwxr-xr-x | mkindex.pl | 148 |
1 files changed, 107 insertions, 41 deletions
@@ -26,10 +26,9 @@ use Getopt::Long qw/:config posix_default no_ignore_case gnu_compat use Pod::Usage; use File::Temp; +# http://geoinformatics.aalto.fi/doc/Geo-GDAL-1.9/html/ 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'; @@ -37,14 +36,26 @@ 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/\.[^.]+$//' + , zoom => '1024x1024' + , resampling => &resamplings('near') ); GetOptions( "q|quiet" => sub { delete $config{progress} } - , "d|debug" => sub { $config{debug} = 1 } - , "rename=s" => sub { $config{rename} = $_[1] } + , "d|debug" => \$config{debug} + , "rename=s" => sub { $config{rename} = $_[1] # TODO: sanitize + } + , "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); + } + } , "r|resampling=s"=> sub { $config{resampling} = &resamplings($_[1]) } , "s_srs=s" => sub { $config{s_srs} = Geo::OSR::SpatialReference::->new(); $config{s_srs}->SetFromUserInput( $_[1] ) } @@ -56,37 +67,92 @@ GetOptions( "q|quiet" => sub { delete $config{progress} } pod2usage(2) unless $#ARGV > 0; my $index = pop; -my @mapset = @ARGV; +my @mapset = map {{ filename => $_ }} @ARGV; -&debug( "Using mapset:", map {" - ".$_} @mapset ); +&debug( "Using mapset:", map {" - ".$_->{filename}} @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 ); +# We also store the corders points and the affine transformation, to +# avoid to compute it again later on. +my $s_srs; +foreach my $map (@mapset) { + my $name = $map->{filename}; + + &debug ("Opening '$name'"); + my $ds = Geo::GDAL::Open($name, 'ReadOnly') or exit 1; + $map->{wkt} = $ds->GetProjection(); + + unless ($config{s_srs}) { + die "Error: Missing projection for '$name'.\n" unless $map->{wkt}; + + $map->{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"; + } } else { - $config{s_srs} = $srs; + $s_srs = $map->{srs}; } } + + &debug( "Using source SRS: ". ($config{s_srs} // $map->{srs}) ); + my $tx = Geo::OSR::CoordinateTransformation::->new( + $config{s_srs} // $map->{srs}, + $config{t_srs} // $config{s_srs} // $s_srs + ); + die "Error: Could compute transformation between source and target SRS.\n" + unless $tx; + + my @geo = $ds->GetGeoTransform(); + $map->{transform} = \@geo; + $map->{corners} = [ [ $geo[0], $geo[3] ] + , [ $geo[0] + $geo[1]*$ds->{RasterXSize} + , $geo[3] + $geo[5]*$ds->{RasterYSize} ] + ]; + # $map->{corners}->[0]: [ $ulx, $uly, $ulz ] + # $map->{corners}->[1]: [ $lrx, $lry, $lrz ] + $tx->TransformPoints( $map->{corners} ); + die "Error: Could not transform point.\n" + if grep !/^[-+]?\d*\.?\d+$/, (map @$_, @{$map->{corners}}); } $config{t_srs} //= $config{s_srs}; -&debug ("Using source SRS: $config{s_srs}"); &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) { + $xRes = $_->{transform}->[1] if abs $xRes < abs $_->{transform}->[1]; + $yRes = $_->{transform}->[5] if abs $yRes < abs $_->{transform}->[5]; + + $ulx = $_->{corners}->[0]->[0] + unless defined $ulx and $ulx < $_->{corners}->[0]->[0]; + $uly = $_->{corners}->[0]->[1] + unless defined $uly and $uly > $_->{corners}->[0]->[1]; + $lrx = $_->{corners}->[1]->[0] + unless defined $lrx and $lrx > $_->{corners}->[1]->[0]; + $lry = $_->{corners}->[1]->[1] + unless defined $lry and $lry < $_->{corners}->[1]->[1]; + } + + my $xRatio = (0.+$lrx-$ulx) / ($xRes*$xSize); + my $yRatio = (0.+$lry-$uly) / ($yRes*$ySize); + $config{zoom} = $xRatio < $yRatio ? $yRatio : $xRatio; +} +&debug( "Using zoom: $config{zoom}." ); + + # Working directory my $tmpdir = File::Temp::tempdir( CLEANUP => 1 ); &debug ("Using working directory: '$tmpdir'."); @@ -105,10 +171,10 @@ 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 ); + &debug ("Opening '$map->{filename}'"); + my $ds = Geo::GDAL::Open($map->{filename}, 'ReadOnly') or exit 1; + &expand_rgba( \$ds, $map ); + &reproject( \$ds, $map ); } @@ -143,6 +209,7 @@ sub resamplings { # 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; return unless $src_ds->{RasterCount} == 1; @@ -172,8 +239,8 @@ sub expand_rgba { , $driver_opts ); - $dst_ds->SetProjection( $src_ds->GetProjection() ); - $dst_ds->SetGeoTransform([ $src_ds->GetGeoTransform() ]); + $dst_ds->SetProjection( $map->{wkt} ); + $dst_ds->SetGeoTransform( $map->{transform} ); $dst_ds->SetGCPs( $src_ds->GetGCPs(), $src_ds->GetGCPProjection() ) if $src_ds->GetGCPCount(); @@ -200,6 +267,7 @@ sub expand_rgba { } # Erase (hence flush and close) the previous dataset + undef $src_ds; $$ds = $dst_ds; return $dst_filename; } @@ -209,21 +277,16 @@ sub expand_rgba { # http://jgomezdans.github.io/gdal_notes/reprojection.html sub reproject { my $ds = shift; + my $map = shift; my $src_ds = $$ds; - my ($src_xSize, $src_ySize) = $src_ds->Size(); + my @geo = @{$map->{transform}}; + my ($ul,$lr) = @{$map->{corners}}; - 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 $dst_xSize = int( ($lr->[0]-$ul->[0]) / ($geo[1]*$config{zoom}) ); + my $dst_ySize = int( ($lr->[1]-$ul->[1]) / ($geo[5]*$config{zoom}) ); - 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) ) ); + my ($src_xSize, $src_ySize) = $src_ds->Size(); &debug("Reducing map from ${src_xSize}x${src_ySize} to ${dst_xSize}x${dst_ySize}."); my $dst_filename = tempfile(SUFFIX => '.tif'); @@ -237,17 +300,20 @@ sub reproject { 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 ]); + $dst_ds->SetGeoTransform([ $ul->[0], $geo[1]*$config{zoom}, $geo[2] + , $ul->[1], $geo[4], $geo[5]*$config{zoom} ]); my $res = Geo::GDAL::ReprojectImage( $src_ds, $dst_ds, - , $config{s_srs}->ExportToWkt(), $t_srs_wkt + , $map->{wkt}, $t_srs_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; return $dst_filename; } |