[gdal-dev] How to use TransformPoints in python?

[Ari Jolma]

Frank Warmerdam kirjoitti:
> Hogendoorn, Rene wrote:
>>
>>     print proj.TransformPoints( len(polygon), x, y, z)
>>
>
> Rene,
>
> To the best of my knowledge TransformPoints() is not currently working 
> smoothly
> with Python.

But it works in Perl:

use Geo::GDAL;

$to = Geo::OSR::SpatialReference->create( EPSG=>28992 );
$to->SetTOWGS84(565.237, 50.0087, 465.658, -0.406857, 0.350733, 
-1.87035, 4.0812);

$t = Geo::OSR::CoordinateTransformation->new( $to,
    Geo::OSR::SpatialReference->create( EPSG=>4326 ) );

$polygon = [[132817.006604435708141, 550302.852720651309937, 0.],
        [131182.28895997320069, 558340.214472591876984, 0.],
        [132578.61028128489852, 558748.893883707583882, 0.],
        [136631.347774848196423, 553436.061539204441942, 0.],
        [136631.347774848196423, 553436.061539204441942, 0.],
        [132817.006604435708141, 550302.852720651309937, 0.]];

$t->TransformPoints($polygon);

for (@$polygon) { print "@$_\n"; }

this returns:

5.0572412099288 52.9393242597464 41.9509746224735
5.03233800529761 53.011479011305 41.8140474820234
5.05311376900775 53.0152112890499 41.8069483591898
5.11380123061856 52.967623034953 41.8972087137249
5.11380123061856 52.967623034953 41.8972087137249
5.0572412099288 52.9393242597464 41.9509746224735

It looks like you'd like to convert from 28992 to 4326 and not vice 
versa. Does the result look ok?

Cheers ;)

Ari

-- 
Prof. Ari Jolma
Environmental Management Information Technology
Teknillinen korkeakoulu / Helsinki University of Technology
tel: +358 9 451 3812 address: POBox 5300, 02015 TKK, Finland
Email: ari.jolma at tkk.fi 
http://geoinformatics.tkk.fi/bin/view/Main/AriJolmaHomePage