##### Child pages
• Python convert coordinates
Go to start of banner

# Python convert coordinates

Go to start of metadata

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/).

Convert coordinates using towgs84 parameter
```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.

Convert coordinates without using towgs84 parameter
```# 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).

Convert coordinates using SR-ORG:6781 (with towgs84 parameter)
```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):

Convert coordinates using EPSG:28992 (without towgs84 parameter)
```# 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
```
• No labels