MGI / Balkans coordinate systems
his article describes the coordinate systems used in the former federal republic of Yugoslavia, which are still valid in most of the countries in the region, and explains the ways to set correctly the coordinate system parameters. The focus is put on defining the correct proj4 string, so that one can convert local coordinates (e.g. from the topomaps) to other coordinate systems, and especially to geographic WGS84 coordinates. Here you can also download the ILWIS coordinate system files that carry the correct definition of the coordinate system and the geodetic datum.
Contents 
Setting the CRS parameters (R/ILWIS)
The Militärgeographisches Institut "MGI / Balkans" coordinate system is the GaussKruegerbased coordinate systems, which has been in use since 1920's (Lapaine, 2001). The whole area of the exYugoslavia was divided into three zones: Zone 5, Zone 6, Zone 7 (see figure), based on the central meridian.

Coordinate system  Gauss Krueger  Slovenia, Western Croatia (Z5)  Eastern Croatia / BiH (Z6)  Serbia, Montenegro, Macedonia (Z7) 

Projection type  Transverse Mercator  
Ellipsoid name  Bessel 1841  
Semimajor axis  6377483.865 m  
Flattening  299.257223563  
dX  550.499 m  550.499 m  574.027 m 
dY  164.116 m  164.116 m  170.175 m 
dZ  475.142 m  475.142 m  401.545 m 
Rx  5.80967  5.80967  4.88786 
Ry  2.07902  2.07902  0.66524 
Rz  11.62386  11.62386  13.24673 
M_BF  5.541764 ppm  5.541764 ppm  6.88933 ppm 
Central meridian  15° E  18° E  21° E 
Central parallel  0° N  
Scale at central meridian  0.9999  
False X  5500000  6500000  7500000 
False Y  0 
The difference from the original Bessel's ellipsoid (socalled Hermannskogel datum) to the current referent system (WGS84) is nonsystematic. The national land surveying agencies would typically setup a network of referent points and then calculate the mean transformation parameters (local geoid) that are in average accurate enough for the most of the area of interest. There are in fact several sources for the global transformation parameters that can be used to link the local datum with the WGS84 system. The most accurate global transformation parameters are the one provided by the national land survey agency (see for example the transformation parameters for each county in Croatia). There are also global transformation parameters that can be used for a wider area, but then the mean positional accuracy can exceed 1 m. For example, In Croatia, transformation parameters for HRG2000 have been estimated using 1780 referent points (the average 3D error is 0.7 meters; see map of the errors for the HRG2000 in figure below) by the State Geodetic Authority. In Serbia, the global transformation parameters were determined based on 1217 referent points (the average 3D error is 0.4 meters).
The MGI / Balkans Zone 6 in the PROJ.4 should be set as:
+proj=tmerc +lat_0=0 +lon_0=18 +k=0.9999 +x_0=6500000 +y_0=0 +ellps=bessel +towgs84=550.499,164.116,475.142,5.80967,2.07902,11.62386,0.99999445824 +units=m
where tmerc is the Transverse Mercator projection system, ellps is the Bessel 1841 ellipsoid, and +towgs84 carries the 7 datum transformation parameters (translation + rotation + scaling) to WGS84; also known as the Bursa Wolf method.
This coordinate system corresponds to the EPSG 31276:
+init=epsg:31276
Or in ILWIS format (note that there are small differences in the way the datum parameters are expressed):
[Ilwis] Description=Coordinate System Projection "MGI / Balkans zone 6" Class=Coordinate System Projection Type=CoordSystem [Domain] Type=DomainCoord [CoordSystem] CoordBounds=1e+308 1e+308 1e+308 1e+308 Width=28 Decimals=2 UnitSize=1.000000 Type=Projection Projection=transverse mercator Datum=User Defined Ellipsoid=Bessel 1841 [Projection] False Easting=6500000.000000 False Northing=0.000000 Central Meridian=18.000000 Central Parallel=0.000000 Scale Factor=0.9999000000 [Datum] dx=550.499 dy=164.116 dz=475.142 Type=User Defined BursaWolf rotX=0.000028166075 rotY=0.000010079373 rotZ=0.000056354064 dS=0.00000554176
Validation of CRS parameters
To validate these parameters, we can use a referent point (T249) for which we should have both WGS84 geographic coordinates and the local coordinates:
UNIQ60c374027135c4d9geshi00000003QINU
To derive the local coordinates we use the spTransform method available in the rgdal package:
UNIQ60c374027135c4d9geshi00000004QINU
The difference between the true and transformed local coordinates is:
UNIQ60c374027135c4d9geshi00000005QINU
Hence the difference is in fact significant (but can be ignored for more general scales e.g. <1:25000). The only way to increase the precision of the datum shift parameters would be to reestimate them locally using local referent points. This can be done for example by using 34 local reference points and the find datum transformation parameters function in ILWIS.
Next, we can check if the point falls where it should in Google maps: UNIQ60c374027135c4d9geshi00000006QINU
Which shows the correct location in the Google Maps:
<widget type="googlemap" width="300">
<marker lat="45.56288883" lon="18.696612475"> T249 </marker>
</widget>
Or use the RgoogleMaps library to plot the maps in R:
UNIQ60c374027135c4d9geshi00000007QINU
which will produce this plot:
We can also try reprojecting/transforming the coordinates in ILWIS:
?TRANSFORM(COORD(6554781.07,5046738.03,gk_6), LatLonWGS84) = 45°33'46.4''N, 18°41'47.8''E
which again matches the known coordinates.
A simplified (but also less precise) solution is to use the 3 parameter datum shift (so called Molodensky method):
UNIQ60c374027135c4d9geshi00000009QINU
Note that the definition of the MGI / Balkans coordinate systems (EPSG:31275, EPSG:31276, EPSG:31277) available via http://spatialreference.org is incomplete because this database does not contain any definition of the geodetic datum. If the geodetic datum is missing, the difference from the true location will be even few hundreds of meters. In fact, if only one datum shift parameter is missing or is in error, transformation of coordinates will lead to frustrations. On the other hand, even if you correctly setup the (global) coordinate system parameters, you will not be able to achieve precision better than 0.5 meters (see above).
Datum shifts by various sources
Datum shifts Bessel ellipsoid to WGS84 as specified by various sources:
 HRG2000 geoid: dX, dY, dZ, Rx, Ry, Rz, M_BF = 550.499, 164.116, 475.142, 5.80967, 2.07902, 11.62386, 0.99999445824;
 HTRS96 geoid: dX, dY, dZ, Rx, Ry, Rz, M_BF = 550.5670, 164.6118, 474.1386, 5.976766, 2.099773, 11.495481, 0.999994552075;
 estimated by the Serbian Geodetic Authority: dX, dY, dZ, Rx, Ry, Rz, M_BF = 574.027, 170.175, 401.545, 4.88786, 0.66524, 13.24673, 0.99999311067;
 estimated by the Austrian Militärgeographisches Institut (MGI) and valid for the whole Western Balkans region: dX, dY, dZ, Rx, Ry, Rz, M_BF = 678.059, 179.019, 585.545, 14.43, 0.42, 18.02, 0.999997490;
 calculated by J. Hendrikse (from three points in Slavonia) and used in ILWIS 3.0: dX, dY, dZ = 682, 199, 480;
 according to GARMIN in GPS receivers: dX, dY, dZ = 675, 205, 475;
 according to NIME: dX, dY, dZ = 682, 203, 480;
Other coordinate systems of interest
The most recent coordinate system used in Croatia to display the whole country is known as the Croatian Terrestrial Reference System (HTRS96/TM) centered at 16°30. Its parameters are (Lapaine, 2007):
+proj=tmerc +lat_0=0 +lon_0=16.5 +k=0.9999 +x_0=500000 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
which is aligned on the European Terrestrial Reference System (ETRS89). The most recent estimate of the global transformation parameters from the old geodetic datum (from year 1901) to HRTS96/TM for Croatian territory is:
+towgs84=550.5670,164.6118,474.1386,5.976766,2.099773,11.495481,0.999994552075
Another coordinate system used by Military and similar organizations is the Universal Transverse Mercator system (Balkans falls in two zones: 33N and 34N). The parameters for Croatia are:
+proj=utm +zone=33 +ellps=WGS84 +datum=WGS84 +units=m +no_defs
The bounding box coordinates for the UTM33N coordinate system for Croatia are: Xmin=360100, Ymin=4670200, Xmax=849400, Ymax=5156700.
To generate the grid cell id in the European Reference Grid System (EPSG:3035), we can run:
UNIQ60c374027135c4d9geshi0000000DQINU
This will create a shape file with grid nodes placed exactly on the reference grid. Each grid cell has an unique name, which is abbreviation of the XY coordinates of the origin of the grid. E.g. 4572_2625 corresponds to the grid node X=4572500, Y=2625500.
Literature:
 Basic, T. 2001. Detailed Geoid Model of the Republic of Croatia (in Croatian), Reports of the State Geodetic Administration on Scientific and Professional Projects in the Year 2000, Ed. Landek I., 1122, Zagreb.
 Basic, T., Hecimovic, 2005. Latest geoid determinations for the Republic of Croatia. Newton's bulletin no3, pp. 8291.
 Annoni, A. (Edt) 2003. European Reference Grids  EUR Report 21494 EN. Proceedings of the "European Reference Grids" workshop, Ispra, 2729 October 2003.
 Lapaine, M. 2001. Topografska kartografija u Hrvatskoj, Geografski horizont, 2, 153.
 Lapaine, M. 2007. New Official Map Projection of Croatia  HTRS96/TM. Kartografija i geoinformacije, Vol. 6 No. izv. / spec.
<Rating> Rate this article: 1 Poor and in error 2 Far from ready 3 Needs improvements 4 Minor corrections 5 Very useful!! </Rating>