Below a snippet of coordinate conversion between RD (EPGS:28992) and WGS84 (EPSG:4326) and vice versa.
Using osgeo.osr
osgeo can be used after installing GDAL (https://pypi.python.org/pypi/GDAL/).
from osgeo.osr import SpatialReference, CoordinateTransformation # Define the Rijksdriehoek projection system (EPSG 28992) epsg28992 = SpatialReference() epsg28992.ImportFromEPSG(28992) # correct the towgs84 epsg28992.SetTOWGS84(565.237,50.0087,465.658,-0.406857,0.350733,-1.87035,4.0812) # Define the wgs84 system (EPSG 4326) epsg4326 = SpatialReference() epsg4326.ImportFromEPSG(4326) rd2latlon = CoordinateTransformation(epsg28992, epsg4326) latlon2rd = CoordinateTransformation(epsg4326, epsg28992) # Check the transformation for a point close to the centre of the projected grid lonlatz = rd2latlon.TransformPoint(155000.0, 446000.0) print(lonlatz) # (5.387203018813555, 52.002375635973344, 43.614926571026444)
The towgs84 parameter is not always mentioned in a specific Coordinate System definition. If the towgs84 clause is contained in its definition, then this should be its "preferred" transformation. The clause towgs84 may not be included in the situations described here.
What is towgs84, actually? As written here, datum shifts can be approximated by 3 parameter spatial translations (in geocentric space), or 7 parameter shifts (translation + rotation + scaling). Those datum shifts can be described using the towgs84 parameter.
Without +towgs84 correction, the tranformation uses the official epsg:28992 definition. That gives a slightly different, erroneous result, which is in the order of ~1 cm for the point selected.
# Define the Rijksdriehoek projection system (EPSG 28992) epsg28992 = SpatialReference() epsg28992.ImportFromEPSG(28992) rd2latlon = CoordinateTransformation(epsg28992, epsg4326) latlon2rd = CoordinateTransformation(epsg4326, epsg28992) # Check the transformation for the same point, without using towgs84 correction lonlatz = rd2latlon.TransformPoint(155000.0, 446000.0) print(lonlatz) # (5.387202946158022, 52.00237563479786, 43.6057764403522)
Using pyproj
As GDAL/OGR uses a linkage to PROJ.4 for transforming between coordinate systems, you could also make use of pyproj, which is the python interface (wrapper) for proj4.
Alternatively to EPSG:28992, we could use another Amersfoort / RD New Coordinate System definition, which contains the proper towgs84 parameter, the SR-ORG:6781. Its Proj4 format is described here. As EPSG:28992, this is an oblique stereographic projection (see definition here).
from pyproj import Proj, transform
# copy the SR-ORG:6781 definition in Proj4 format from http://spatialreference.org/ref/sr-org/6781/
p1 = Proj("+proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +towgs84=565.237,50.0087,465.658,-0.406857,0.350733,-1.87035,4.0812 +units=m +no_defs")
p2 = Proj(proj='latlong',datum='WGS84')
# Transform point (155000.0, 446000.0) with SR-ORG:6781
lon, lat, z = transform(p1, p2, 155000.0, 446000.0, 0.0)
print(lon, lat, z) # 5.38720301881 52.002375636 43.614926571
See the difference with EPSG:28992 (without towgs84):
# define the projections p1 = Proj(init='epsg:28992') p2 = Proj(proj='latlong',datum='WGS84') # Transform point (155000.0, 446000.0) with EPSG:28992 lon, lat, z = transform(p1, p2, 155000.0, 446000.0, 0.0) print(lon, lat, z) # 5.38720294616 52.0023756348 43.6057764404
