diff options
| -rwxr-xr-x | mkindex.pl | 73 |
1 files changed, 51 insertions, 22 deletions
@@ -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 {"'".$_."'"} |