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
_______________________________________________
gdal-dev mailing list
[email protected]
http://lists.osgeo.org/mailman/listinfo/gdal-dev
