Python proj4 test Austrian datum shift to WGS84 from local coordinate system

test_datum_shift_at2wgs84

from pyproj import Proj, transform

## my test coordinates, location is on the end of a peer on lake Wörthersee
## if the tranformation is wrong our point will be in the water 

#   these coorinates below I believe I can trust
#   and are all TRUE  for reference ( source KAGIS Atlas 17.10.2014)

#       epsg        proj      coord
#       epsg:31252  GK M31	70534.8	165013.4
#  	epsg:31252  GK M34	-159211.2	166701.3
#  	epsg:31258  BMN M31	520534.8	165013.4
#  	epsg:4326   WGS84	14.253560	46.620721
#  	WGS84(GM)	14°15,21'	46°37,24'
#  	WGS84(GMS)	14°15'12,8''	46°37'14,6''
#  	epsg:3767  UTM 33N	442851.2	5163288.2
#   Geländehöhe(ALS): +439,97 müA	Objekthöhe(ALS): +0,02 m
#   EPSG 31253, 31252, 31251 correct +towgs84=577.326,90.129,463.919,5.137,1.474,5.297,2.4232

UTM_x = 442851.2; UTM_y = 5163288.2;
BMN_M31_x = 520534.8; BMN_M31_y = 165013.4;
GK_M31_x = 70534.8; GK_M31_y = 166701.3;
# GB_x = 2423346.99; GB_y = 5055619.87;
WGS84_lat = 46.6207211 # Decimal Degrees
WGS84_lon = 14.253560 # Decimal Degrees
# UTM_z = WGS84_z = 52.8 # Ellipsoidical height in meters

wgs84 = Proj(proj='latlong',datum='WGS84')
print 'proj4 library version = ',wgs84.proj_version
utm33 = Proj(proj='utm',zone='33')
bmn_m31 = Proj(init='epsg:31258', towgs84="577.326,90.129,463.919,5.137,1.474,5.297,2.4232")
# +towgs84=682,-203,480,0,0,0,0
# gaussb = Proj(init='epsg:26592',towgs84="-122.74,-34.27,-22.83,-1.884,-3.400,-3.030,-15.62")

# qgis epsg:31285
# +proj=tmerc +lat_0=0 +lon_0=13.33333333333333 +k=1 +x_0=450000 +y_0=0 +ellps=bessel +towgs84=577.326,90.129,463.919,5.137,1.474,5.297,2.4232 +units=m +no_defs

# qgis epsg: 31258
# +proj=tmerc +lat_0=0 +lon_0=13.33333333333333 +k=1 +x_0=450000 +y_0=-5000000 +ellps=bessel +towgs84=577.326,90.129,463.919,5.137,1.474,5.297,2.4232 +units=m +no_defs

print '\nBMN M31-->WGS84'
print 'Klagenfurt, Woerthersee BMN (converted):', xgb,ygb,zgb
back_lon, back_lat, back_z = transform(bmn_m31, wgs84, xgb, ygb, zgb)
print 'Klagenfurt, Woerthersee WGS84 (converted back):', back_lon,\
  back_lat, back_z
print 'Difference (seconds):', (back_lon-WGS84_lon)*3600,\
  (back_lat-WGS84_lat)*3600, back_z-WGS84_z, '(m)'

print '\nUTM-->WGS84'
print 'Klagenfurt, Woerthersee UTM33 (converted):', xutm33, yutm33
back_lon, back_lat, back_z = transform(utm33, wgs84, xutm33, yutm33)
print 'Klagenfurt, Woerthersee WGS84 (converted back):', back_lon,\
  back_lat, back_z
print 'Difference (seconds):', (back_lon-WGS84_lon)*3600,\
  (back_lat-WGS84_lat)*3600
  
print '\nWGS84-->BMN M31'
print 'Klagenfurt, Woerthersee BMN M31 (from KAGIS):', BMN_M31_x, BMN_M31_y
z = 0
xgb, ygb, zgb = transform(wgs84, bmn_m31,WGS84_lon, WGS84_lat, z)
print 'Klagenfurt, Woerthersee BMN M31  (converted):', xgb,ygb
print 'Difference (meters):', xgb-BMN_M31_x, ygb-BMN_M31_y