Compute the zoom factor from the expected output size.

1 files changed, 107 insertions, 41 deletions

diff --git a/mkindex.pl b/mkindex.pl
index cfc6ac3..dcaf35a 100755
--- a/mkindex.pl
+++ b/mkindex.pl
@@ -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;
 }