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

[Steve Gaffigan]

Rene,

You could also use the pyproj module (http://code.google.com/p/pyproj/).

> import os, osgeo, numpy
import pyproj
> from osgeo import ogr, osr
> 
> 
> def getProjection( EPSGcode ):
>     """ set a projection from ... to lat-long """
> 
>     srs = osgeo.osr.SpatialReference()
>     srs.ImportFromEPSG( EPSGcode )
>     # RD-NL lacks TOWGS84
>     if EPSGcode == 28992:
>         srs.SetTOWGS84( 565.237, 50.0087, 465.658, -0.406857, 0.350733, 
> -1.87035, 4.0812 )
> 
       return pyproj.Proj(srs.ExportToProj4())
> 
> 
> if __name__ == '__main__':
> 
>     proj = getProjection( 28992 )
> 
>     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.]]
> 
>     x = numpy.array( [p[0] for p in polygon] )
>     y = numpy.array( [p[1] for p in polygon] )
>     z = numpy.array( [p[2] for p in polygon] )
> 
       print proj(x, y, inverse=True)
> 
> returns:
> 
(array([ 5.05762415,  5.03271685,  5.05349572,  5.11419251,  5.11419251,
         5.05762415]), array([ 52.94040396,  53.01256714,  53.01629989, 
  52.96870616,
         52.96870616,  52.94040396]))

Steve