# MGI / Balkans coordinate systems

T

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 topo-maps) 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.

## Setting the CRS parameters (R/ILWIS)

The "MGI / Balkans" coordinate system is the Gauss-Krueger-based coordinate systems, which has been in use since 1920's (Lapaine, 2001). The whole area of the ex-Yugoslavia was divided into three zones: Zone 5, Zone 6, Zone 7 (see figure), based on the central meridian.

 CRS_balkans.R : Some examples. gk_5.csy : MGI / Balkans zone 5 in ILWIS format. gk_6.csy : MGI / Balkans zone 6 in ILWIS format. gk_7.csy : MGI / Balkans zone 7 in ILWIS format. croatia.csy : HTRS96/TM coordinate system in ILWIS format.

 Coordinate systems used in the former Yugoslavia: division in zones (based on the central Meridian).

Table: Parameters of the MGI / Balkans coordinate systems
Coordinate system - Gauss Krueger Slovenia, Western Croatia (Z5) Eastern Croatia / BiH (Z6) Serbia, Montenegro, Macedonia (Z7)
Ellipsoid name Transverse Mercator
Ellipsoid name Bessel 1841
Semi-major 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 to the current referent system (WGS84) is non-systematic. The national land surveying agencies would typically set-up 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 country (see map of the errors for the HRG2000 in figure below). 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. For example, in Croatia, the global transformation parameters (HRG2000) have been estimated using 1780 referent points (the average 3D error is 0.7 meters) 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).

 Differences between GGM01C and HRG2000 geoid solutions (m). After Basic & Hecimovic (2005).

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:

```> T249 <- data.frame(Lat=as(char2dms("45d33'46.3998\"N"), "numeric"),
+    Lon=as(char2dms("18d41'47.80491\"E"), "numeric"))
> coordinates(T249) <- ~Lon+Lat
> proj4string(T249) <- CRS("+proj=longlat +datum=WGS84")
```

To derive the local coordinates we use the spTransform method available in the rgdal package:

```> T249.xy <- spTransform(T249, CRS("+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"))
> T249.xy
SpatialPoints:
Lon       Lat
[1,] 6554781.6 5046739.2
# should be: X=6554781.07,  Y=5046738.03
```

The difference between the true and transformed local coordinates is:

```> spDistsN1(matrix(c(6554781.07, 5046738.03), ncol=2), coordinates(T249.xy))
[1] 1.2709
```

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 re-estimate them locally using local referent points. This can be done for example by using 3-4 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:

```> system(paste('"c:/Program Files/firefox/firefox.exe"', ' -url http://maps.google.com/maps?&ll=',
+    coordinates(T249)[2],',',coordinates(T249)[1],'&t=h&z=19&q=',
+    coordinates(T249)[2],',',coordinates(T249)[1], sep=""), wait=F)
```

Or use the RgoogleMaps library to plot the maps in R:

```> library(RgoogleMaps)
> MyMap <- GetMap.bbox(center=rev(coordinates(T249)), zoom=18, destfile="T249.png", maptype ="satellite")
> PlotOnStaticMap(MyMap, lat=coordinates(T249)[,2], lon=coordinates(T249)[,1])
```

which will produce this plot:

 Referent point T249 as seen on the satellite image.

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

```> spTransform(T249, CRS("+proj=tmerc +lat_0=0 +lon_0=18 +k=0.9999 +x_0=6500000
+y_0=0 +ellps=bessel +towgs84=682,-199,480 +units=m"))
SpatialPoints:
Lon       Lat
[1,] 6554780.2 5046738.3```

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 set-up 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;
• 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;
• 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

A coordinate system used in Croatia to display the whole country is known as the HTRS96/TM, centered at 16°30 (Basic, 2001). Its parameters are:

```+proj=tmerc +lat_0=0 +lon_0=16.5 +k=0.9997 +x_0=2500000 +y_0=0 +ellps=bessel
+towgs84=550.499,164.116,475.142,5.80967,2.07902,-11.62386,0.99999445824
+units=m
```

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.

## 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., 11-22, Zagreb.
• Basic, T., Hecimovic, 2005. Latest geoid determinations for the Republic of Croatia. Newton's bulletin no3, pp. 82-91.
• Lapaine, M. 2001. Topografska kartografija u Hrvatskoj, Geografski horizont, 2, 1-53.