[Dutch] Waarom vanuit Python wel RD => WGS84 maar andersom niet?

[Stefan de Konink]

Om op sommige plekken niet afhankelijk te zijn van een PostGIS installatie
had ik een simpel Python scriptje gemaakt die van RD naar WGS84 kan. Nu
dacht ik slim te zijn en het scriptje uit te breiden naar de omgekeerde
variant maar nu komt er totale bagger uit.

"""
select X(loc), Y(loc) from (select transform(setsrid(makepoint(5.96, 52.2),
4326), 28992) AS loc) as x;
       x         |        y
------------------+------------------
194159.178375368 | 468142.275781605
(1 row)
"""

from rd import wgs84rd, rdwgs84
wgs84rd(5.96, 52.2)
(128410, 445807) # crap dus

rdwgs84(194159, 468142)
(5.959997, 52.199998) # komt aardig overeen

Wat ik gebruik; het enige andere bij de tweede functie is dus de omdraaing
van src/dst.

from osgeo import osr, ogr
src_string = '+proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889
+k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +units=m
+towgs84=565.2369,50.0087,465.658,-0.406857330322398,0.350732676542563,-1.8703473836068,4.0812
+no_defs no_defs <>'
dst_string = '+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs'

src = osr.SpatialReference()
src.ImportFromProj4(src_string)
dst = osr.SpatialReference()
dst.ImportFromProj4(dst_string)

def rdwgs84(x, y):
    ogr_geom = ogr.CreateGeometryFromWkt('POINT(%d %d)' % (x, y))
    ogr_geom.AssignSpatialReference(src)
    ogr_geom.TransformTo(dst)

    return (round(ogr_geom.GetX(), 6), round(ogr_geom.GetY(), 6))

def wgs84rd(lon, lat):
    ogr_geom = ogr.CreateGeometryFromWkt('POINT(%d %d)' % (lon, lat))
    ogr_geom.AssignSpatialReference(dst)
    ogr_geom.TransformTo(src)

    return (int(round(ogr_geom.GetX())), int(round(ogr_geom.GetY())))



Wat gaat er fout?

Stefan
------------- volgend deel ------------
Een HTML-bijlage is gescrubt...
URL: <http://lists.osgeo.org/pipermail/dutch/attachments/20120527/86076cb0/attachment.html>