summaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
authorGuilhem Moulin <guilhem@fripost.org>2013-08-11 00:27:19 +0200
committerGuilhem Moulin <guilhem@fripost.org>2013-08-11 00:28:54 +0200
commit019fe72b12ca4bcaf074a9d79a3aa0a1d0916b7a (patch)
tree8d0ed81379fe41083ac4723fe1b5610c6abcf019
parentc50f7662900f583968ccb4ecd8d46150d47f0382 (diff)
Compute the zoom factor from the expected output size.
-rwxr-xr-xmkindex.pl148
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;
}