diff options
| -rwxr-xr-x | mkindex.pl | 96 |
1 files changed, 42 insertions, 54 deletions
@@ -134,30 +134,22 @@ foreach my $map (@mapset) { die "Error: Could compute transformation between source and target SRS.\n" unless $tx; - my @geo = $ds->GetGeoTransform(); - $map->{transform} = \@geo; - $map->{corners} = [ &pix2proj(\@geo, - ( [ 0, 0] - , [ 0, $ds->{RasterYSize}-1 ] - , [ $ds->{RasterXSize}-1, 0 ] - , [ $ds->{RasterXSize}-1, $ds->{RasterYSize}-1 ] - )) - ]; - $tx->TransformPoints( $map->{corners} ); + $map->{s_transform} = $ds->GetGeoTransform(); + $map->{t_corners} = [ map {[ Geo::GDAL::ApplyGeoTransform($map->{s_transform}, @$_) ]} + ( [ 0, 0] + , [ 0, $ds->{RasterYSize}-1 ] + , [ $ds->{RasterXSize}-1, 0 ] + , [ $ds->{RasterXSize}-1, $ds->{RasterYSize}-1 ] + ) + ]; + $tx->TransformPoints( $map->{t_corners} ); die "Error: Could not transform point.\n" - if grep !/^[-+]?\d*\.?\d+$/, (map @$_, @{$map->{corners}}); + if grep !/^[-+]?\d*\.?\d+$/, (map @$_, @{$map->{t_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(); @@ -178,31 +170,31 @@ if ($config{zoom} =~ /^(\d+)x(\d+)$/) { my $dst_ds = Geo::GDAL::AutoCreateWarpedVRT( $src_ds , $config{s_wkt} // $_->{wkt} , $config{t_wkt} ); - $_->{transform2} = [ $dst_ds->GetGeoTransform() ]; + $_->{t_transform} = [ $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]; - - $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}->[3]->[0] - unless defined $lrx and $lrx > $_->{corners}->[3]->[0]; - $lry = $_->{corners}->[3]->[1] - unless defined $lry and $lry < $_->{corners}->[3]->[1]; + $xRes = $_->{s_transform}->[1] if abs $xRes < abs $_->{s_transform}->[1]; + $yRes = $_->{s_transform}->[5] if abs $yRes < abs $_->{s_transform}->[5]; + + $ulx = $_->{t_corners}->[0]->[0] + unless defined $ulx and $ulx < $_->{t_corners}->[0]->[0]; + $uly = $_->{t_corners}->[0]->[1] + unless defined $uly and $uly > $_->{t_corners}->[0]->[1]; + $lrx = $_->{t_corners}->[3]->[0] + unless defined $lrx and $lrx > $_->{t_corners}->[3]->[0]; + $lry = $_->{t_corners}->[3]->[1] + unless defined $lry and $lry < $_->{t_corners}->[3]->[1]; } my $xRatio = (0.+$lrx-$ulx) / ($xRes*$xSize); my $yRatio = (0.+$lry-$uly) / ($yRes*$ySize); $config{zoom} = $xRatio < $yRatio ? $yRatio : $xRatio; - map { $_->{transform2}->[1] = $_->{transform}->[1] * $config{zoom} - ; $_->{transform2}->[5] = $_->{transform}->[5] * $config{zoom} + map { $_->{t_transform}->[1] = $_->{s_transform}->[1] * $config{zoom} + ; $_->{t_transform}->[5] = $_->{s_transform}->[5] * $config{zoom} } @mapset; } &debug( "Using zoom: $config{zoom}" ); @@ -296,7 +288,7 @@ sub expand_rgba { ); $dst_ds->SetProjection( $map->{wkt} ); - $dst_ds->SetGeoTransform( $map->{transform} ); + $dst_ds->SetGeoTransform( $map->{s_transform} ); $dst_ds->SetGCPs( $src_ds->GetGCPs(), $src_ds->GetGCPProjection() ) if $src_ds->GetGCPCount(); @@ -336,11 +328,10 @@ sub reproject { my $map = shift; my $src_ds = $$ds; - my @geo = @{$map->{transform}}; - my ($ul,$lr) = @{$map->{corners}}[0,3]; + my ($ul,$lr) = @{$map->{t_corners}}[0,3]; # TODO: it might be xSize that's determined from ySize - my $dst_xSize = int( ($lr->[0]-$ul->[0]) / $map->{transform2}->[1] ); + my $dst_xSize = int( ($lr->[0]-$ul->[0]) / $map->{t_transform}->[1] ); my $dst_ySize = int( $dst_xSize / $config{ratio} ); my ($src_xSize, $src_ySize) = $src_ds->Size(); @@ -355,7 +346,7 @@ sub reproject { , $driver_opts ); $dst_ds->SetProjection( $config{t_wkt} ); - $dst_ds->SetGeoTransform( $map->{transform2} ); + $dst_ds->SetGeoTransform( $map->{t_transform} ); my $res = Geo::GDAL::ReprojectImage( $src_ds, $dst_ds, , $map->{wkt}, $config{t_wkt} @@ -412,25 +403,22 @@ sub annotate { # adding the border *then* reproject we use a mask in order to # have nice anti-aliasing and better looking fonts. - my ($ul1,$ll1,$ur1,$lr1) = @{$map->{corners}}; - my ($rx,$ry) = @{$map->{transform2}}[1,5]; - my ($ul2,$ll2,$ur2,$lr2) = &pix2proj( $map->{transform2}, - ( [ 0, 0] - , [ 0, $ySize-1 ] - , [ $xSize-1, 0 ] - , [ $xSize-1, $ySize-1 ] - )); + my ($success,@invgeo) = Geo::GDAL::InvGeoTransform($map->{t_transform}); + die "Error: Couldn't inverse affine transformation\n" unless $success; + push @draw, primitive => 'polygon'; - push @points, [ int(($ul1->[0]-$ul2->[0]) / $rx) - , int(($ul1->[1]-$ul2->[1]) / $ry) ] - , [ int(($ll1->[0]-$ll2->[0]) / $rx) - , int(($ll1->[1]-$ll2->[1]) / $ry) + $ySize ] - , [ int(($lr1->[0]-$lr2->[0]) / $rx) + $xSize - , int(($lr1->[1]-$lr2->[1]) / $ry) + $ySize ] - , [ int(($ur1->[0]-$ur2->[0]) / $rx) + $xSize - , int(($ur1->[1]-$ur2->[1]) / $ry) ] - ; + push @points, map {[ Geo::GDAL::ApplyGeoTransform(\@invgeo, @$_[0,1] ) ]} + @{$map->{t_corners}}[0,1,3,2]; + my $mask = Image::Magick::->new( size => "${xSize}x${ySize}" ); + $mask->ReadImage('xc:black'); + $mask->Draw( primitive => 'polygon' + , fill => 'white' + , points => join (' ', map {join ',', @$_} @points) ); + $img->Mask( mask => $mask ); + undef $mask; } + &debug( ucfirst {@draw}->{primitive} .": ". + join (' ', map {join ',', @$_} @points) ); $img->Draw( @draw, points => join (' ', map {join ',', @$_} @points) ); $img->Annotate( @ann ); |