Skip to content

Instantly share code, notes, and snippets.

@BishopGIS
Created May 5, 2020 18:37
Show Gist options
  • Save BishopGIS/67d365e4d560f5197c20c1320a7f1d60 to your computer and use it in GitHub Desktop.
Save BishopGIS/67d365e4d560f5197c20c1320a7f1d60 to your computer and use it in GitHub Desktop.
Test coordinate transformation accuracy
from osgeo import gdal, osr
y = 7509137.5811
x = 4187591.89173
gsk2011 = '+proj=longlat +a=6378136.5 +rf=298.2564151 +towgs84=0.013,-0.092,-0.03,0.001738,0.003559,-0.004263,0.0074 +no_defs'
src = osr.SpatialReference()
dst = osr.SpatialReference()
src.ImportFromEPSG(3857)
dst.ImportFromProj4(gsk2011)
ct1 = osr.CoordinateTransformation(src, dst)
ct2 = osr.CoordinateTransformation(dst, src)
result = [x, y]
for i in xrange(0, 10000):
result = ct1.TransformPoint(result[0], result[1], 0.0)
result = ct2.TransformPoint(result[0], result[1], 0.0)
print('Orgignal: [{}:{}], result [{}:{}], delta [{}mm:{}mm]'.format(x, y, result[0], result[1], (x - result[0])*1000.0, (y - result[1])*1000.0))
@BishopGIS
Copy link
Author

Orgignal: [4187591.89173:7509137.5811], result [4187591.89173:7509137.5811], delta [9.31322574615e-07mm:1.86264514923e-06mm]

@BishopGIS
Copy link
Author

BishopGIS commented May 24, 2023

from osgeo import gdal, osr

y = 7509137.5811
x = 4187591.89173

src = osr.SpatialReference()
dst = osr.SpatialReference()
src.ImportFromEPSG(3857)
dst.ImportFromEPSG(7683)

ct1 = osr.CoordinateTransformation(src, dst)
ct2 = osr.CoordinateTransformation(dst, src)

result = [x, y]
for i in xrange(0, 1000000):
result = ct1.TransformPoint(result[0], result[1], 0.0)
result = ct2.TransformPoint(result[0], result[1], 0.0)

print('Original: [{}:{}], result [{}:{}], delta [{}mm:{}mm]'.format(x, y, result[0], result[1], (x - result[0])*1000.0, (y - result[1])*1000.0))

@BishopGIS
Copy link
Author

Original: [4187591.89173:7509137.5811], result [4187591.89173:7509137.5811], delta [9.31322574615e-07mm:1.86264514923e-06mm]

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment