Versions Compared

Old Version 4

changes.mady.by.user Giorgio Santinelli

Saved on

compared with

New Version 5

changes.mady.by.user Kees den Heijer

Saved on

Key

  • This line was added.
  • This line was removed.
  • Formatting was changed.

...

Code Block
languagepy
titleConvert 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
latlonzlonlatz = rd2latlon.TransformPoint(155000.0, 446000.0)
print latlonzlonlatz # (5.387203018813555, 52.002375635973344, 43.614926571026444)

...

Code Block
languagepy
titleConvert 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
latlonzlonlatz = rd2latlon.TransformPoint(155000.0, 446000.0)
print latlonzlonlatz # (5.387202946158022, 52.00237563479786, 43.6057764403522)

...

Code Block
languagepy
titleConvert 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
latlon, lonlat, z = transform(p1, p2, 155000.0, 446000.0, 0.0)
print latlon, lonlat, z # 5.38720301881 52.002375636 43.614926571

...

Code Block
languagepy
titleConvert 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
latlon, lonlat, z = transform(p1, p2, 155000.0, 446000.0, 0.0)
print latlon, lonlat, z # 5.38720294616 52.0023756348 43.6057764404