Spatial reference systems and transformations in Austria

Transforming data from Austrian local coordinate system to WGS84 should be a no brainer but this is GIS stuff we are talking about and therefor it isn't :)

First I am posting some great references for my own sanity that I find it again and that other GIS, Geoinformatik students find some help.

Link to epsg codes for Austria ESRI Synergis: http://www.esri-austria.at/downloads/coords_at.html

• Für Datumsübergänge von MGI nach WGS 1984 sollte, falls keine genauen Werte für die 7 Parameter Transformation vorliegen, jene der ÖK 50 vom BEV verwendet werden - dies entspricht MGI_To_WGS_1984_3 (EPSG-code 1618).

Link 1: We in Austria DO have a Datum shift in order to move from a Projected Coordinate System to Geographic Coordinate system and the two magic code numbers are in this file below (1618(the rest of us in Austria and 1194(for Steirmark))

http://svn.osgeo.org/metacrs/geotiff/trunk/libgeotiff/csv/datum_shift.csv

BIG NOTE on PROJ4 taken from the homepage:

OTHER PROGRAMS
       The  proj  program  is  limited  to   converting
       between  geographic  and  projection coordinates
       within one datum.

       The cs2cs program operates similarly, but allows
       translation  between any pair of definable coor‐
       dinate  systems,  including  support  for  datum
       translation.

How to call this from Python using the python implementatino of Proj4 library (code taken from the source code of pyProj4 library) and simplified for your reading pleasure.

	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

view raw test_datum_shift_at2wgs84 hosted with ❤ by GitHub

This github gist is here https://gist.github.com/da5a72f894aa72a28032.git

Have fun playing with projections!



http://epsg.io/
http://www.gomogi.com