summaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
-rwxr-xr-xmkindex.pl73
1 files changed, 51 insertions, 22 deletions
diff --git a/mkindex.pl b/mkindex.pl
index dcaf35a..3ed4293 100755
--- a/mkindex.pl
+++ b/mkindex.pl
@@ -82,7 +82,7 @@ foreach my $map (@mapset) {
&debug ("Opening '$name'");
my $ds = Geo::GDAL::Open($name, 'ReadOnly') or exit 1;
- $map->{wkt} = $ds->GetProjection();
+ $map->{wkt} = $ds->GetProjection() || $config{s_srs}->ExportToWkt();
unless ($config{s_srs}) {
die "Error: Missing projection for '$name'.\n" unless $map->{wkt};
@@ -110,18 +110,31 @@ foreach my $map (@mapset) {
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} = [ &pix2proj(\@geo,
+ ( [ 0, 0]
+ , [ 0, $ds->{RasterYSize}-1 ]
+ , [ $ds->{RasterXSize}-1, 0 ]
+ , [ $ds->{RasterXSize}-1, $ds->{RasterYSize}-1 ]
+ ))
];
- # $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}});
+
+
+ { local $_ = basename $name; eval $config{rename}; $map->{text} = $_; }
+}
+
+sub pix2proj {
+ my $geo = shift;
+ map [ $geo->[0] + $_->[0]*$geo->[1] + $_->[1]*$geo->[2],
+ $geo->[3] + $_->[0]*$geo->[4] + $_->[1]*$geo->[5] ],
+ @_
}
$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}");
@@ -133,6 +146,18 @@ if ($config{zoom} =~ /^(\d+)x(\d+)$/) {
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} );
+ $_->{transform2} = [ $dst_ds->GetGeoTransform() ];
+ $config{ratio} = (0. + $dst_ds->{RasterXSize}) / $dst_ds->{RasterYSize};
+ &debug( "Using ratio: $config{ratio}" );
+ undef $src_ds;
+ undef $dst_ds;
+
$xRes = $_->{transform}->[1] if abs $xRes < abs $_->{transform}->[1];
$yRes = $_->{transform}->[5] if abs $yRes < abs $_->{transform}->[5];
@@ -140,17 +165,21 @@ if ($config{zoom} =~ /^(\d+)x(\d+)$/) {
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];
+ $lrx = $_->{corners}->[3]->[0]
+ unless defined $lrx and $lrx > $_->{corners}->[3]->[0];
+ $lry = $_->{corners}->[3]->[1]
+ unless defined $lry and $lry < $_->{corners}->[3]->[1];
}
my $xRatio = (0.+$lrx-$ulx) / ($xRes*$xSize);
my $yRatio = (0.+$lry-$uly) / ($yRes*$ySize);
+
$config{zoom} = $xRatio < $yRatio ? $yRatio : $xRatio;
+ map { $_->{xRes} = $_->{transform}->[1] * $config{zoom}
+ ; $_->{yRes} = $_->{transform}->[5] * $config{zoom}
+ } @mapset;
}
-&debug( "Using zoom: $config{zoom}." );
+&debug( "Using zoom: $config{zoom}" );
# Working directory
@@ -229,7 +258,7 @@ sub expand_rgba {
$lookup->[$_] .= chr $e[$_] for (0 .. $dst_nb-1);
}
- my $dst_filename = tempfile(SUFFIX => '.tif');
+ my $dst_filename = &tempfile(SUFFIX => '.tif');
my ($xSize, $ySize) = $src_ds->Size();
my $dst_ds = Geo::GDAL::GetDriverByName($driver)->Create( $dst_filename
, $xSize
@@ -281,15 +310,16 @@ sub reproject {
my $src_ds = $$ds;
my @geo = @{$map->{transform}};
- my ($ul,$lr) = @{$map->{corners}};
+ my ($ul,$lr) = @{$map->{corners}}[0,3];
- my $dst_xSize = int( ($lr->[0]-$ul->[0]) / ($geo[1]*$config{zoom}) );
- my $dst_ySize = int( ($lr->[1]-$ul->[1]) / ($geo[5]*$config{zoom}) );
+ # TODO: it might be xSize that's determined from ySize
+ my $dst_xSize = int( ($lr->[0]-$ul->[0]) / $map->{xRes} );
+ 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}.");
- my $dst_filename = tempfile(SUFFIX => '.tif');
+ my $dst_filename = $index; #&tempfile(SUFFIX => '.tif');
my $dst_ds = Geo::GDAL::GetDriverByName($driver)->Create( $dst_filename
, $dst_xSize
, $dst_ySize
@@ -297,14 +327,13 @@ sub reproject {
, $driver_type
, $driver_opts
);
-
- my $t_srs_wkt = $config{t_srs}->ExportToWkt();
- $dst_ds->SetProjection( $t_srs_wkt );
- $dst_ds->SetGeoTransform([ $ul->[0], $geo[1]*$config{zoom}, $geo[2]
- , $ul->[1], $geo[4], $geo[5]*$config{zoom} ]);
+ $dst_ds->SetProjection( $config{t_wkt} );
+ $map->{transform2}->[1] = $map->{xRes};
+ $map->{transform2}->[5] = $map->{yRes};
+ $dst_ds->SetGeoTransform( $map->{transform2} );
my $res = Geo::GDAL::ReprojectImage( $src_ds, $dst_ds,
- , $map->{wkt}, $t_srs_wkt
+ , $map->{wkt}, $config{t_wkt}
, $config{resampling}, undef, undef
, $config{progress} );
die "Error: Failed to reproject map " . join (', ', map {"'".$_."'"}