Implementation of Local Reference Systems without Length Distortion in GIS Software

. In geoinformatics, there is a considerable demand for dealing with different layers given in different geodetic datums and map projections within the same project. In the majority of the software, the computations are completed on-the-fly using the open-source program PROJ. In recent years, the definition of national reference systems for a wide range of countries has already been developed and published; therefore, users have the opportunity to set them in GIS applications. Another significant challenge is to set local reference systems in individual mapping products, developed and maintained for a limited area, typically for industrial zones. These local reference systems may have special requirements as regards avoiding length distortions caused by map projections. This paper presents a solution for defining such kinds of local coordinate systems and implementing it in any mapping system using the PROJ library. The applied method can also be imported into different GNSS units at the same time.


Introduction
For simplicity, the Earth's shape is approximated by an ellipsoid of revolution while maps are developed traditionally on a plane. Thus, a geometric or mathematical connection is needed to flatten the globe's surface into a plane. The projection can either be performed to the plane directly or to a geometric body (cylinder or cone), which can be developed into the plane (Hofmann-Wellenhof et al., 2003). An ellipsoid of revolution cannot be mapped without distortion of distances onto a plane (Maling, 1992). In geodesy, primarily conformal projections with no distortions in angles are preferred while the distortion of lengths, the scale factor, is optimized. In national map projections, the scale factor often exceeds a few hundred parts per million (ppm), which can be problematic for various engineering projects.
Efforts have been invested in developing low distortion map projections, especially for engineering surveying tasks (Dennis, 2016). In this particular case, the scale factor might remain below a negligible level (e.g. 1 ppm) (Pickford and Gibbings, 2009). Since engineering surveying projects are typically carried out within a few kilometres, it is acceptable to approximate the figure of the Earth by plane in practice. Consequently, traditional engineering surveying tasks are carried out using a two-dimensional rectangular coordinate system, which is usually defined at the average reduced level of the area to avoid applying height correction to distance measurements. A horizontal shift of the origin is often applied to introduce false easting and northing to separate x and y coordinates by their magnitude. The system can be oriented either to the North or to a specific direction, e.g. the principal axis of an engineering object. The local reference system is connected to the national one, using common points which have coordinates in both the local and the national systems. The connection is typically a horizontal similarity transformation; its accuracy depends mainly on the true position errors of the reference points.
Recently, a demand has arisen to also use the geographic coordinate system in engineering surveying. There are at least two main reasons. Firstly, in many cases, the measurement technique based on the Global Positioning System (GPS) is an efficient alternative to traditional measurement techniques. Secondly, geoinformatics systems are widely used in parallel with different layers in various coordinate systems. Therefore, an adequate map projection needs to be defined. Granted that this is used only in a limited area, the map projection can be defined so that its scale factor remains less than an acceptable limit. Moreover, the developed map projection is supposed to fit a network of existing reference points whose coordinates were determined by classical surveying techniques, from various historical periods.
In Geographic Informatic Systems (GIS), it is quite common to use different layers given in different geodetic datums and map projections, within the same project. In many cases, transformation computations are completed using the open-source library, known as PROJ (PROJ, 2020). PROJ is a generic coordinate transformation software that transforms geospatial coordinates from one coordinate reference system to another. This includes a wide range of cartographic projections as well as geodetic transformations. In addition to most open-source software (QGIS, GRASS, MapServer, PostGIS for instance), some commercial applications also use PROJ in the background. Developers can apply PROJ relatively flexibly due to its command-line applications and application programming interface (API). Transformation on a set of input points from one coordinate reference system to a second can be performed using its cs2cs commandline application. The parameters of coordinate reference systems, including map projection parameters, can be given using either control parameters to the command line applications or referring to predefined settings stored in databases.

Oblique Mercator map projection
The Oblique Mercator for the sphere or the ellipsoid of revolution is equivalent to a regular Mercator projection that has been altered by wrapping a cylinder around the sphere so that it touches the surface along the great circle path chosen for the central line, instead of along the Earth's Equator (Snyder, 1987). The tilt angle (azimuth) of the central line can be given in two different ways. In the first case, the azimuth is given directly. In the second case, the azimuth is given indirectly by specifying two central line points (PROJ, 2020). Equations for both variants to compute easting and northing from latitude and longitude (forward direction) as well as latitude and longitude from easting and northing (the reverse direction) were developed by Martin Hotine (Hotine, 1947) using hyperbolic functions, which were reformed into a closed form by Snyder (Snyder, 1987).
The Oblique Mercator projection is useful for mapping areas having their greatest extent along a direction that is neither north-south nor east-west (PROJ, 2020). Several map projections can be derived from the original Oblique Mercator, with a wide range of countries using one of them in their national mapping systems. One of the most classical applications is map projection for the Alaskan panhandle. Switzerland uses a different implementation of Oblique Mercator by Rosenmund, while Madagascar uses the Laborde version (Kennedy and Kopp, 2001). Oblique Mercator projection was proposed for Moldova to decrease length distortion in large scale mapping (Chiriac and Vlascenco, 2015). In addition to the map projections of national mapping systems, a further variation of Oblique Mercator has been applied in remote sensing since the 1970s (Colvocoreses, 1974). Taking advantage of this kind of map projection has little scale distortion within the sensing range of an orbiting mapping satellite (Engels and Grafarend, 1995).
Normally, for Mercator projection, the scale factor depends only on latitude : The scale factor reaches 1 ppm at 0.081°, which corresponds to about 9 km on Earth's surface. Consequently, the Mercator projection can be applied with negligible scale factor within a distance of 9 km to the line where the sphere (or the cylinder) of the projection touches the planet regardless of its orientation (normal, transversal or oblique). For comparison, the function of distance and scale factor on stereographic projection is also plotted in Fig. 1.   Fig. 1. The scale factor of Mercator and stereographic projections for short distances The scale factor of Oblique Mercator projection can be computed using the formula (Snyder, 1987) or (Hooijberg, 2008): where , are ellipsoidal latitude and longitude, 0 is the longitude of the centre of the projection, = 6384076.910 , = 1.000752020 are calculated constants of the Oblique Mercator projection, , are the ellipsoid parameters, and is the coordinate along the centre line of the projection, which can be calculated either from the geographic or the rectified ( , ) coordinates. When investigating the scale factor, the inverse method is preferred since the conversion from geographic to rectified is available using specific PROJ commands (see listings later in this paper): where 0 , 0 are false easting and northing at the natural origin, while is the rectified bearing of the central line.

Study area
Our method to implement an Oblique Mercator map projection with almost negligible scale factor on a limited area has been tested in the area of a large industrial zone, located in Hungary. Its approximate size is about 4 km 2 . During its construction, even during the operation phase, there has been a wide range of engineering surveying tasks with enhanced accuracy requirements. Consequently, networks of both horizontal and vertical reference points were developed back in the 1980s. The majority of the horizontal reference points were established by reinforced concrete pillars (Fig. 2). While the vertical network has been measured regularly, every two years recently, the horizontal network's latest measurement was performed more than 20 years ago. The horizontal coordinates were determined with mm accuracy by adjusting precise angle and distance measurements taken by a Kern DKM2 theodolite and a Mekometer electronic distance measuring instrument. The local coordinate system was defined at the average reduced level of the area (96.5 m), with the system oriented to the North. The system has a false origin to separate easting and northing by their magnitude and have only positive coordinates.

Fig. 2. Horizontal reference points in the network measured by GPS
In 2019, a GPS based measurement campaign was launched to determine any possible crustal movements in the region. Within the horizontal network, 14 points were chosen and measured by dual-frequency receivers in static mode. The 6-hour sessions and the observations were post-processed by the Bernese scientific post-processing software (Dach et al. 2015). The standard deviation of the horizontal coordinates was typically ± 1 mm or even less, while the standard deviation of the reduced levels was ± 1 or 2 mm. One of the campaign's main results were the geographic coordinates of the 14 points in the ETRS89 system. Therefore, these points can be considered common points to define the necessary map projection and coordinate transformation between ETRS89 and the local coordinate system.

Results
At first, an Oblique Mercator map projection of the WGS84 ellipsoid was defined. The centre point of the map projection is located in the geometrical centre of the area of interest. In the first step, neither an alignment nor a scale was applied (Listing 1).

Listing 1. The PROJ command to project points given with geographic coordinates to an Oblique
Mercator map projection defined only in its centre point. No transformation is performed in this initial step.
cs2cs +proj=latlong +datum=WGS84 +to +proj=omerc +lat_0=46.57668421 +lonc=18.85441407 +x_0=0.0 +y_0=0.0 +k_0=1.0 +alpha=0 +gamma=0 +ellp=WGS84 +units=m The projected coordinates were then transformed into the local system of reference points using a two-dimensional similarity transformation. Parameters were estimated using an open-source project, GeoEasy (Siki, 2018) applying a traditional least-squares solution with ± 7 mm root mean square error (Table 1). The shift values correspond to the x and y coordinates of the centre point, while the scale factor comes from height correction. It is also worth mentioning that the elevation above the ellipsoid (150 m on average) has to be taken into account for the height correction of distances. A small rotation angle is also detected since the local network's orientation is slightly tilted to the North in the ETRS89 system. The maximum residual of the reference points is 18 mm, while the majority are only 2 mm (Fig 3.). At first, the most probable explanation of the residuals is the significant time difference between the terrestrial and satellite-based determination of the points; hence their slight movement cannot be excluded. On the other hand, a few millimetres residual might occur due to the accuracy of both determinations. Nevertheless, the accuracy of the transformation is acceptable from a practical point of view.  . 3. Residuals of transformation in the common points of the transformation. The arrows' length is proportional to its absolute value, which is also given in millimetres as labels. Arrows are rotated to present the direction of residuals. (Base map and data from OpenStreetMap and OpenStreetMap Foundation) With the pipeline option of PROJ, it is possible to perform a set of operations chained together on the same input data (PROJ, 2020); hence, the geographic coordinates can be converted into the local system in a single PROJ command (Listing 2). On the other hand, the pipeline option is available only in version 5 or above, and the majority of GIS applications currently still use version 4.
Listing 2. The PROJ pipeline command. In addition to the projection, a two-dimensional Helmert transformation is also performed. The transformation is given by control parameters.
cs2cs +proj=latlong +datum=WGS84 +to +proj=pipeline +step +proj= omerc +lat_0=46.57668421 +lon_0=18.85441407 +ellps=WGS84 +step +proj=helmert +convention=coordinate_frame +x=1395.698 +y=11825.906 +s=1.000023273 +theta=521 Involving that transformation into the parameters of the map projection can fix this issue. PROJ offers the false origin and scale factor for the majority of projections. Moreover, rotation is also involved in Oblique Mercator projection by defining the azimuth and the rectified bearing of the centre line (Listing 3).
Listing 3. The PROJ command, which is equivalent to the pipeline version in Listing 2 but transformation is included in the map projection parameters cs2cs +proj=latlong +datum=WGS84 +to +proj=omerc +lat_0=46.57668421 +lonc=18.85441407 +x_0=1395.698 +y_0=11825.906 +k_0=1.000023273 +alpha=-0.14472222 +gamma=0 +ellp=WGS84 +units=m As described previously, the Oblique Mercator projection can be used with a maximum of 1 ppm scale factor within approximately a 9 km distance from its centre line. The statement has been confirmed by computing the scale factor applying the precise formula (2) in a grid in the area of interest. It is proved that the scale factor reaches 1 ppm approximately 9 km from the centre line (Fig 4.). The 9 km figure meets all the requirements since specialist engineering jobs are typically carried out within a few kilometres. As the Oblique Mercator projection is an option in a wide variety of GNSS receivers when setting the coordinate system, the developed map projection can also be adopted in most GPS units (Fig 5.). The settings were checked by measuring some of the reference points with a Leica GPS receiver operating in RTK mode. The coordinates determined by the GPS receiver and projected using the Oblique Mercator were directly compared to the reference points' known values. Residuals remained below the expected accuracy.

Conclusion
In this paper, seamless integration of engineering surveying and global reference systems in QGIS (QGIS, 2020.) and other PROJ based GIS software, as well as geodetic GNSS receivers, have been presented. The system can be considered scale-factor-less within an approximately 9 km range of the central line, which exceeds the specific area of interest in engineering surveying. The method is based upon an Oblique Mercator projection of WGS84 ellipsoid and a horizontal similarity transformation. However, the transformation can be skipped by setting up adequate parameters of the map projection. While the majority of map projections support involving shift and scale of transformation, there is no option for the rotation. The main advantage of applying Oblique Mercator projection is that its centre line's azimuth is equivalent to a rotation of the coordinate system. The proposed solution offers an excellent opportunity to maintain engineering surveying accuracy while simultaneously using external (national or global) map sources.