summaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
authorGuilhem Moulin <guilhem@fripost.org>2013-08-11 20:00:58 +0200
committerGuilhem Moulin <guilhem@fripost.org>2013-08-11 20:00:58 +0200
commit71ff18a52e5a63ceef5ff83aa39b1d7e23e93ee6 (patch)
tree7f745011a391ee4bc3a5fbdb8e8c81c3c70e60d5
parentc0ba2b1fb8a0128fdba50dcae7d873a45f45854e (diff)
Use the builtin Geo::GDAL::ApplyGeoTransform.
-rwxr-xr-xmkindex.pl96
1 files changed, 42 insertions, 54 deletions
diff --git a/mkindex.pl b/mkindex.pl
index 661655a..2766b9f 100755
--- a/mkindex.pl
+++ b/mkindex.pl
@@ -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 );