Test coordinate transformation accuracy

test_osr.py

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

The same but with fiona that is easy to be installed in Python venv:

from fiona.transform import transform

x, y = 4187591.89173, 7509137.5811

crs_merc = "epsg:3857"
crs_gsk2011 = dict(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=True)

x_ = [x]
y_ = [y]
for _ in range(10000):
    x_, y_ = transform(crs_merc, crs_gsk2011, x_, y_)
    x_, y_ = transform(crs_gsk2011, crs_merc, x_, y_)

print(
    "Orgignal: [{}:{}], result [{}:{}], delta [{}mm:{}mm]".format(
        x, y, x_[0], y_[0], (x - x_[0]) * 1000.0, (y - y_[0]) * 1000.0
    )
)

Orgignal: [4187591.89173:7509137.5811], result [4187591.89161:7509137.58102], delta [0.117402058095mm:0.0778120011091mm]

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, 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))

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

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

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