Received 5 January 2016; accepted 14 May 2016; published 17 May 2016
Interchange and mutual comparison between geo-referenced data sets require the precise knowledge of both the cartographic and the geodetic systems on which they are based. The lack of this type of information is a major obstacle in the integration of information in geographic information system (GIS), precluding its analysis and hence the benefit of mutual comparisons and relationship. The broad (and fully justified) use of Google Earth imagery in surveying requires that external data are brought in the geodetic standard of Google Earth (namely: either geographic or UTM coordinates, with WGS84 as datum).
During the last fifty years, the east-central African republic of Burundi has benefited from systematic and detailed surveys in numerous fields as geophysics (ground or airborne), geology, agronomy/soil science, land-use, biology, socio-economics, politics, resulting in maps (local, general), or database records. These data have been collected, and the maps produced using different standards as the latter have been changing with time and the habits of the authors.
The topographical map of Burundi at scale 1/50k  , made available as from 1981 has been the base for a number of these surveys. Moreover, older surveys, based on maps of a lesser quality or simply older, have often been transferred onto the 1981 map (as for instance the geological maps  ). The 1981 map is thus a milestone for a wide community of users.
The topographical map of Burundi at scale 1/50,000 in 38 sheets published in 1981 by the IGEBU was established using the Transverse Mercator projection (Central Meridian 30˚E/0˚N, Xo = 500.000 m, Yo = 10.000.000 m, f = not given, but probably 0.9999, as is the usual practice) and the Clarke 1880 modified ellipsoid  , altitudes being given relatively to sea level (MSL).
The position of the reference geodetic points used for this map have not been published, nor the geodetic datum, i.e. the position and orientation of the ellipsoid used relatively to the mass center of the Earth. This latter information became indeed pertinent only after the emergence of satellite-based geodesy and the subsequent use of the mass center of the Earth as geometric center (as does datum WGS84). As long as this information is not available, the 1981 topographical map and other maps derived from it cannot be used in a GIS for comparison with data or maps based on another system.
Burundi is located in the domain of Arc 1950 geodetic network. The only indication published on relations between Arc 1950 and the WGS84 datum is a local datum (“Arc 1950-Burundi, ARF-H Code”) based on 3 directly measured geodetic points, with the values ∆x = −153 ± 20 m, .∆Y = −5 ± 20 m, .∆z = −292 ± 20 m (NGA 1991 in  . It is not known whether these translation parameters are specific to the 1981 map or if they are anterior.
The recent realization  of a digital air cover in natural colors with 0.5 m resolution, complemented by a DEM with 10 m resolution, based on the ellipsoid and WGS84 datum, but with a UTM projection (zone 35), ellipsoidal heights (h 1 in the table), tied to a novel geodetic network of 17 points, can provide an adequate estimate of the datum to attribute to the topographic map of 1981.
This information would allow the use the all-important 1981 map in a GIS and engage in comparisons with other geo-referenced data or maps using other datum.
The determination of the transformation parameters required for passing from a projection using ellipsoid A to another projection using ellipsoid B necessitates the knowledge of not only the shape (length of axes) of ellipsoids A and B, but also the position of their center and angles of rotation around their 3 axes.
By making the simplifying assumption that the difference between the two ellipsoids (Clarke 1880 mod. and WGS84 in this case) results only from a translation (i.e. without rotations), this difference can be estimated from the difference in cartesian (=geocentric) X, Y, Z coordinates of a point common to the two standards and calculated in each of them (see  or  ).
h: elevation above ellipsoid;
e: eccentricity =;
a, b: semi-axes of the ellipsoid;
N: vertical curvature =.
In order to calculate X, Y, Z, we need, besides the shape of the ellipsoid (a and b, as above), the geographic coordinates (longitude, latitude) and elevation of the point. Usually, one makes use of geodetic points of the national geodetic network, of known (i.e. assumed) coordinates in this (old) system, and re-determines the position of the same points in the (new) WGS84-based system. This latter information is readily obtained from GPS measurements. The transformation parameters are obtained from the difference between the old and the new system, as shown below.
Because the (old) geodetic points are no longer visible in the field (nor displayed on the 1981 map), we will use this map itself instead of the geodetic points on which the map is based and compare featured elements of the map with the same objects on the WGS84 orthorectified imagery.
In the case of the 2013 imagery, the coordinates are derived readily from the image and associated DEM.
The 1981 topographic map, duly geo-referenced in its own system (see above), provides the required geographic (latitude/longitude) coordinates as well as elevation via the presence of leveled points (NB MSL altitude assimilated to the geoid: H in Table 1). These leveled points, frequently placed at road intersections, are practically the only ones unequivocally recognizable on both documents.
The topographic map at 1/50 k has a resolution that can be estimated to 5 m, thus differing by an order of magnitude of that of the 2013 image; the map is further symbolized widely, especially as regards buildings (it is quite obsolete in this respect); these are source of inaccuracy. It is also established from multiple aerial photographs, which is a possible source of local and random distortions.
Being random, these errors can be compensated by averaging a number of points spread over a large area (38 in this case). Several points were considered for each site, placed in the vicinity (more or less 1 km) of the 2013 geodetic points, in the hope of minimizing possible distortions in the 2013 image. It would obviously have been more practical to work directly with the 17 geodetic points created by the producer of the 2013 image (readily visible) if these items had been recognizable on the 1981 map, which is not the case.
Regarding the vertical datum, that is to say the distance between the ellipsoid and the geoid used in the 1981 map, it can be assumed that the 1981 geodesist wished to “paste” the reference ellipsoid to the geoid (as this is the norm, resulting in a local datum), and that this distance is consequently zero. We therefore calculate the geocentric coordinates considering that the altitude refers to the ellipsoid.
The arithmetic difference between the geocentric coordinates (x, y and z) in both systems will thus constitute a measure of the value of the parameter for passing from the ellipsoid used for the 1981 map and to the one used for the 2013 image.
The geocentric coordinates can be calculated, with a geodetic calculator, for instance  (as we did) or from the formulas above in a spread-sheet.
3. Step-by-Step Procedure
1) The reference map is geo-coded; the latitude, longitude of any point can then be read with an appropriate GIS;
2) A point is selected on the map, preferably a leveled point, thus providing the elevation. The latitude, longitude and elevation are now available;
3) The latitude, longitude of the same point is observed on the WGS84-referenced image and the elevation read on the associated digital elevation model;
4a) The geocentric coordinates of 2 are calculated, NB: The map elevation is considered as being relative to the ellipsoid, not MSL;
4b) The geocentric coordinates of 3 are calculated;
5) The differences in X, Y and Z between the two sets of coordinates are computed. These differences represent the transformation (translation) parameters we are searching for;
6) The working GIS in completed with a new reference system/datum (=shape of ellipsoid and translation parameters) proper to the map.
Steps 2 - 5 are repeated for a number of points, if desired, and the average value is provided to step 6. Step 3 observations can be provided by a Google Earth image.
The procedure described above has been repeated for 38 points and the difference calculated. The average and standard deviation for the latter is (Table 1):
Table 1.Geographic and geocentric coordinates of points common to the 1981 map and the 2013 image,and difference between the geocentric coordinates.The latter(d),and their average(m),provide the requested transformation parameters.
The dispersion of the values is more important than what is obtained from direct measurements on geodetic points (see for instance  , where error margins are less by an order of magnitude). This is probably the result of combined inexactitudes in the process of map production (exaggerated wideness of roads for instance) and the resolution of the map itself (5 m at best for a 1/50 k map).
Regarding the absolute values, these results are more precise but indeed very close to those proposed by the NGA  . Either set of values can thus be used in the process of integrating 1981-based maps in a GIS.
Note that the average value of N (the difference between h and H in WGS84) is −9.97 ± 5.33 m in Table 1. When calculating this value from h and H provided by a geodetic computer of  for 17 geodetic reference points, −10.39 ± 0.39 m is obtained.
An alternative, but nevertheless very appropriate data set of WGS84 imagery is provided by Google Earth especially in areas where high resolution imagery is available notwithstanding the poor resolution in elevation.
An early version of this paper received constructive comments by Pierre Voet (NGI, Belgium).