Georeferencing of Multi-Sheet Maps Based on Least Squares with Constraints—First Military Mapping Survey Maps in the Area of Czechia

Abstract

The article deals with the possibility of georeferencing old multi-sheet map works. Various approaches to problem solving and a workable method for using the least squares method with the conditions of the adjacency of map sheets are discussed. To increase reliability, the IRLS robust statistical method is used, which uses iterative weighting of individual measurements based on Huber’s M-estimate. The method is applied to the First Military Mapping of the Habsburg monarchy as a typical representative of old topographic maps, which are not easy to georeference due to unknown parameters of the used cartographic projection. A georeferenced layer of the above mentioned mapping is available on the Mapire.eu portal as well. A basic analysis of the comparison of georeferencing results using our method and the mentioned portal is performed.

1. Introduction

If we want to learn what our landscape looked like 200 or more years ago, we have the opportunity to work with the first topographic maps, which record individual topographic elements in sufficient detail. These topographic maps serve as a memory of the landscape and show the state of our environment in the period before the Industrial Revolution. Cartographers can process and make these maps available in the current reference coordinate systems to a number of experts whose work is based on working with spatial data. These can be historians, geographers, environmentalists, demographers, ecologists, and others. Unfortunately, the topographic maps of that time were created with no geodetic basis, thus, it is not possible to use the standard georeferencing procedures that we know from use in newer mappings.

In general, it is possible to consider the use of the map sheet index and image transformation using a projective transformation and map sheet corner coordinates. This method can be applied if we know the parameters of the used cartographic projection, the reference projection surface and the planar coordinate system [1]. This method cannot be applied to the First Military Survey (1763–1787) of the Habsburg Empire. The Austrian literature agrees that the survey has no geodetic and, therefore, projection background [2,3]. The second option is to use the coordinates of ground control points (GCPs) with different types of transformations. Since the maps transformed in this way do not fit to each other after georeferencing due to inaccuracies of GCPs, it is possible to apply a dual approach. It is either possible to create a continuous mosaic of map sheets before georeferencing and then transform the whole image at once (used, e.g., in [1,4]), or use a type of transformation that preserves the continuity of map sheets even after georeferencing (e.g., in [5,6]). Our considerations on the issue follow the latter option, because we believe that, in this way, it is possible to eliminate some problematic aspects of the transformation of the entire mosaic at once. The aim of this article is to present a method for georeferencing multi-sheet map works, including its evaluation and analysis. The area of Bohemia was chosen as the test area on the maps of the First Military Mapping.

In the second half of the 18th century (1763–1787), the so-called First Military Mapping was carried out on the territory of the former Habsburg monarchy. It arose on the basis of the needs of the Austrian army after the Seven Years’ War (1756–1763) [2]. Unfortunately, these maps were created without geodetic networks, using the à la vue method, and their accuracy is not sufficient from today’s point of view [3]. Individual map sections, brouillons, have rectangular shapes and dimensions of 61.8 cm by 40.8 cm. From this, we can deduce the use of a cylindrical type of projection. In some provinces (Bohemia, Moravia, Silesia, Hungary), it was possible to use existing Müller maps [7] as a basis, probably also made in a cylindrical projection [8]. Molnár and Timár use Cassini’s projection as a probable projection of the map [9]. As each map sheet has been processed separately, the map drawing does not fit exactly on the adjacent edges. Despite this, for a number of historical or environmental studies, there is need to work with mosaics of larger territories, preferably directly-mapped provinces.

The ideal situation for the georeferencing of old maps is the knowledge of the used cartographic projection, including all parameters. If this information is not known, it is possible to try to estimate the projection and its parameters. This issue is thoroughly discussed in [10,11], nevertheless, it is successfully applied only to small-scale maps. The procedure is actually based on the display of meridians and parallels on individual maps. For large and medium scale maps, to which topographic maps belong, it is practically impossible to reveal the used cartographic projection. The differences between projections are, in fact, at the graphic level of erroneous GCPs which are utilized for the calculation. Therefore, if we do not know the projection and cannot estimate it, it is possible to work with the most probable variant. Thus, for example, Cassini’s projection was used to georeference a map of France [12]. If this most probable variant still does not correspond exactly to the processed maps, it is possible to add a second step by some of the interpolation methods, e.g., grid shift binary (GSB) as in [9].

Georeferencing can also be used without the knowledge of the projection used. Then, it is necessary to work with GCPs coordinates and use one of the transformation methods. A large group of global methods (similarity, affine, polynomial transformations) and/or local transformations using interpolation mechanisms (spline, IDW) is offered. The division of transformation methods is described in [13].

Among the global methods, an affine transformation is very often used. However, if the cartographic projection of the old map and the reference coordinate system are significantly different, it is necessary to use more complex polynomial transformations. They also successfully reduce errors at GCPs after the transformation. The testing of the degrees of polynomial transformations for georeferencing is performed in [14]. It is important to note that the least squares method (LSM) is usually used to calculate global transformation methods [15]. To calculate the adjustment correctly, there is a need to know how individual GCPs and their surroundings behave after the transformation. It is especially important to evenly and continuously cover the transformed area with GCPs [16]. When balancing using the least squares method, it is also possible to consider assigning weights to GCPs based on their characteristics. This is addressed in [17]. A standard adjustment using LSM has the disadvantage of not being able to reveal erroneous or poor-quality measurements. The GCP which is subject to an error and is not detected affects the adjustment results. Therefore, it is good to consider statistical methods that are able to filter out these measurements or suppress their influence. There are many statistical methods that can test residuals at GCPs. We can mention McKay–Nair’s method of testing of maximal residual based on known standard deviation. Another possibility can be Pearson–Chandrasekhar’s test or testing using mean residual. The use of statistical tests for GCPs or the use of a robust adjustment based on the maximum likelihood method, the so-called M-estimator [18], can be recommended.

For local transformation methods, transformations based on a radial function are used. Spline transformations are the most used, namely the thin plate spline (TPS) method [19]. However, it is possible to use other interpolations. Their common feature is non-residuality on GCPs and thus image distortion. Therefore, it is important to have a relatively large number of GCPs in the map drawing so that excessive deformations do not occur in areas not covered by GCPs. An alternative to using interpolation methods is to transform the image in parts (triangles, quadrilaterals) [20]. The division of the old map image into triangles and their separate transformations is called rubber-sheeting [21]. However, the disadvantage of these procedures may be visible edges between individual parts. For example, line elements are bent on them. Choosing a particular transformation method is a very difficult task. The comparison of the standard polynomial transformation and the method of rubber-sheeting can be found in [22]. A comparative study of transformations when applied to historical map of France is presented in [23]. The distribution and the number of GCPs and the use of kernel-based methods are discussed. The reliability of individual GCPs and thus the assignment of different weights is also used.

Our approach applies global transformation methods, affine or polynomial transformations, with the conditions of the adjacency of neighboring map sheets after georeferencing. The result is a continuous mosaic without gaps or overlaps. The iterative reweighted least squares (IRLS) method [24] is used to calculate LSM, using the reassignment of weights to individual measurements within iterations based on Huber’s M-estimate [25]. In the article, we will test the quality of the results of the mosaic of the First Military Mapping in the area of Bohemia and compare the results with the mosaic created on the Mapire.eu portal [26].

2. Data and Methods

The First Military Mapping survey of the Habsburg monarchy was carried out in Bohemia in the years 1764–1767. This incredibly short time was made possible by the fact that the mapping took place from the eye, only by officers on horseback. Enlargements of the Müller map of Bohemia may have been used, which had been made in 1720 at a scale of approximately 1:132,000. The measuring table was probably used in part, the lengths were estimated or stepped. In particular, militarily important objects were mapped: roads, bridges, brick buildings, meadows, pastures, forests, or watercourses. The elevations shown on these maps are made by hachures. In winter, clean drawings of map sheets measuring 61.8 cm by 40.8 cm at a scale of 1:28,800 were created from the original maps created in the field. In addition to the map sheets, textual comments and military-geographical descriptions of the area were supplemented. In addition to 273 map sheets in Bohemia, 19 volumes of topographic descriptions of the landscape were made. Due to the imperfections of the maps, which were manifested in the wars with Prussia in 1778 and 1779, the rectification of the most important (northern) map sheets was ordered. Thus, 141 map sheets were newly mapped and two map sheets were corrected. The differences between the older and rectified maps are not very obvious, for example, the drawing of forest areas differs slightly. Of course, all the resulting maps were confidential. The original map sheets are stored in the Österreichisches Staatsarchiv—Kriegsarchiv. The J. E. Purkyně University in Ústí nad Labem was involved in making these maps accessible [27]. We have a total of 273 raster images of maps of the first military mapping in a resolution of 400 dpi and a color depth of 24 bits (TIFF files without compression). The data volume of each sheet is more than 300 MB. The quality of this digital data is sufficient for most of the intended analyzes. A resolution of 600 dpi would be more appropriate for automatic image pattern searching, automatic text recognition, and other graphical operations.

A good identification of the corners of map sheets and GCPs is crucial for georeferencing maps and creating continuous mosaics. However, identifying corners is very problematic. The map sheets have a marked map frame on the north and east sides only (see Figure 1). The northeast corner is, thus, clearly identifiable. For others, finding the exact point is problematic, most, of course, in the southwest corner, where it is only possible to infer its position according to the map drawing. Identifying GCPs is easier. In particular, permanently identifiable objects (churches, chapels, farmsteads etc.) were selected for georeferencing; if such objects were missing, then road intersections, sacral objects, or other map elements were selected. The current location of these objects was determined according to the state map work of the Czech Republic, or according to the already georeferenced Second Military Mapping. Out of 273 map sheets, only 250 were selected to create the mosaic. The remaining 23 map sheets are marginal areas where it is not possible to read the corners of the sheets well and the numbers of GCPs are insufficient on them. Therefore, they were excluded from these calculations. Their addition to the resulting mosaic is possible after individual georeferencing using a polynomial transformation. A total of 6849 GCPs were collected on 250 map sheets, averaging more than 27 points on the map sheet. In addition, a regular distribution within the maps was important (see Figure 2).

For georeferencing, we chose a global method using the least squares method. In order for individual sheets to fit together, it was necessary to include conditions in the adjustment. The calculation of intermediate measurements with conditions is used for the calculation. Standard correction equations include a matrix of plan A, a vector of unknowns x, and a measurement vector of l:

$$\mathit{v}=\mathit{A}\mathit{x}-\mathit{l}$$

(1)

As all unknowns in the equations are in a linear form, there is no need to linearize the equations. The derivatives of transformation equations according to individual variables appear in the matrix A. For the affine transformation of one map sheet we have a total of six unknowns, for the polynomial transformation of the second degree we have 12 unknowns.

The equations of conditions must be added to these equations. In our case, the conditions define the identity of the corresponding edges. In the case of an affine transformation, the identity of the end points of the edges after the adjustment is ensured, in the case of a polynomial transformation, the identity of the points in the middle of the edges after the adjustment must also be used:

$$\mathit{B}\mathit{x}=\mathbf{0}.$$

(2)

The solution then has the form:

$$\left(\begin{array}{cc}{\mathit{A}}^{\mathit{T}}\mathit{A}& \mathit{B}\\ {\mathit{B}}^{\mathit{T}}& \mathbf{0}\end{array}\right)\left(\begin{array}{c}\mathit{x}\\ \mathit{k}\end{array}\right)=-\left(\begin{array}{c}{\mathit{A}}^{\mathit{T}}\mathit{l}\\ \mathit{b}\end{array}\right).$$

(3)

For the total calculation of the joint adjustment it is necessary to prepare the matrices A, B and the vector l. In the joint adjustment of 250 map sheets by means of an affine transformation, there will be a total of 1500 unknowns. The matrix A will have 1500 columns. The number of rows in the matrix will correspond to the number of measurements, i.e., the number of GCPs multiplied by two. In our case, it will be 13,698 lines. The number of conditions, and thus the number of rows in the matrix B, is given by the layout of the map sheets in the mosaic. It is important that the conditions are independent of each other. The calculation would collapse if one of the conditions were a linear combination of the others. In our case, the layout of the map sheets is obvious from Figure 2. For a polynomial transformation of the second degree we will have a total of 3000 unknowns, the number of measurements remains the same, the number of conditions increases due to the identity of the center points of the edges.

The result of the adjustment are transformation parameters for individual map sheets. After applying them to the image data, the adjacent edges fit exactly together. The drawing behind the map frames remains a problem. This must be masked within the mosaic using defined masks. After the affine transformation, the masks are parallelograms, after the polynomial transformation of the second degree, they are formations composed of four parabolic arcs.

In order to increase the reliability of our results, it is also possible to use the robust IRLS method. It is based on the repeated calculation of the adjustment, with the weights being reassigned to individual measurements in the calculation on the basis of an estimate of their remoteness. To adjust the weights, it is crucial to choose the correct score function. The parabolic function is used for LSM. To obtain a robust estimate, we use Huber’s solution, where the edge parts of the parabola are replaced by straight lines. Then, it is possible to derive an influence function and a weight function that assigns specific weights. As a parameter c for Huber’s estimate, we chose c = 1.5, i.e., the case where the expected contamination of measurements is 4%. More about this statistical approach can be found in [25]. After each iteration, it is possible to calculate the ratio of the resulting mean standard errors (RMSE). The end of the calculation is set with a difference of less than 0.1%.

The whole calculation is quite demanding and, therefore, we decided to create the MultiGeoref software that solves this problem [28]. The first prototype of the software was written in Matlab, but it turned out that it would be better to use the possibilities of C++ for more demanding calculations. The software uses the Qt graphics library and works as a Windows application. It was also used for georeferencing Müller maps [29].

Another possibility to refine the results seems to be the introduction of weights into the calculation of adjustment. The weights should reflect different expected accuracy of individual types of GCPs (e.g., churches, yards, crossings, etc.) and their implementation requires additional knowledge of conditions in the area and especially the correct determination of weights based on experimentally determined a posteriori mean errors of individual adjusted coordinates. The scales then enter the adjustment as a diagonal matrix P, which contains inverse values of the mean square errors of the individual previously adjusted measurements.

This can be done individually for all coordinates, but deriving weights only from a posteriori mean errors can be impacted with errors similarly to the measurements themselves. However, it is possible to use the categorization of GCPs into similar types and assess these individual categories in aggregate, as a similar degree of accuracy may be assumed to determine the coordinates of objects of similar categories while maintaining the same form of mapping throughout the province, which was generally met within the First Military Survey.

In the experimentally tested area in northern Bohemia, mean errors achieved on individual types of GCPs were evaluated. It was categorized as follows: churches, farmsteads/courtyards/castles, mills, road crossings, other types of objects (bridging of watercourses, crossing of land cultures, etc.)

Several foreign projects deal with the georeferencing of old topographic maps. The Mapire.eu website is a map portal showing available maps of selected European countries in a seamless view. In particular, detailed topographic maps made during the 18th and 19th centuries are included, while seamless georeferencing makes it possible to display an extensive part of Europe as a large summary map, which has no analogues between contemporary map projects.

One of the layers available on the portal is the First Military Survey of the Habsburg monarchy. Georeferenced layers connected into a seamless mosaic are stored in the form of TLS map services and available within the Mapire.eu portal as separate map layers according to individual periods.

Already on the basis of a visual comparison of the map layer published on Mapire.eu with the original sheets, it is evident that the territory was probably transformed in parts (triangular networks defined by GCPs or other per partes method), which is clearly visible in breakages of line objects, which appear abundantly within the maps (see Figure 3). Together with the knowledge of not as many GCPs used for the whole adjustment, this raises the question of whether georeferencing conceived in this way is sufficiently accurate and whether higher geometric accuracy could not be achieved by using global transformation with edge connection conditions and applying robust statistical methods to eliminate gross deviations, as outlined above.

Thus, it is hypothesized that the adjustment of the IRLS method with the conditions of continuity of the edges of the map sheets leads to a lower total load of the resulting mosaic by geometric errors, which may stem mainly from a significantly higher number of identical points used. To verify this hypothesis, another group of control points (different from the original GCPs used for the initial adjustment) was collected, on which deviations of these points from the actual positions within both adjusted results were tested, i.e., both Mapire.eu maps and IRLS-adjusted mosaics.

3. Results

The developed MultiGeoref software [28] was used for georeferencing our area of interest, and 250 map sheets were adjusted in 4 variants. The second-order polynomial and affine transformations were tested, both in the classical version with an LSM adjustment and in the IRLS version with Huber’s estimation. The results of the resulting mean errors and the number of iterations are shown in

Table 1. Results of georeferencing.

Transformation	Affine LSM	Affine IRLS	Polynomial 1 LSM	Polynomial 1 IRLS
RMSE	1335 m	1043 m	390 m	280 m
No. of iterations	1	13	1	10

¹ Second order polynomial.

.

It is obvious that the use of IRLS leads to the reduction of the resulting RMSE, in the case of the polynomial transformation even quite significantly. The resulting error of 280 m is at a very good level according to the parameters of the map creation. The results differ from those published in [28], mainly due to different contamination settings of the input measurements. Furthermore, a different number of map sheets were used (two were discarded) and with a different set of GCPs (some were also discarded or modified after the inspection). The distributions of the resulting deviations (corrections) on GCPs are shown in Figure 4. The maximum frequencies reach deviations in the range from 100 m to 300 m.

The IRLS method used, together with the robust Huber’s estimate, provides an estimate of the adjusted characteristic, which is, compared to the classical LSM, not unbiased but is less affected by completely incorrectly determined GCPs. The disadvantage is also the higher computational complexity of the method. The comparison of the computational time of georeferencing using classical LSM and IRLS is described in [28].

The resulting image data in the mosaic show proper characteristics and usually copy the rectangular shape of map sheets. Some edge sheets are significantly deformed because they contain only a small area covered by GCPs. The rest of the sheet is, thus, significantly affected, especially with the polynomial transformation used, by the need to extrapolate outside the area covered by the data (see Figure 5). Problems can also be caused by a sparse population and generally not a very rich drawing in some peripheral parts of Bohemia (especially in northern, northeastern and southwestern mountainous areas), which makes it difficult to find GCPs; their lack has consequently a secondary effect on the result of georeferencing.

As part of testing the possibility of introducing weights by analyzing the behavior of mean errors of individual adjusted coordinates in test areas in northern Bohemia and their classification into individual defined categories, the existence of a correlation between this categorization and the achieved mean errors was not reliably proven. A higher level of accuracy of the representation of churches and other fixed, important objects (farmsteads, châteaux, yards, mills) has been demonstrated, but assessing the accuracy of road crossings and other “lower” types of GCPs is difficult due to their great diversity and, most likely, a lower priority of the geometric accuracy of their rendition (they could often have been a schematic “background”, mere iconographic symbols, complementing the substantial content of maps). Examples of collected GCPs with their position on the OpenStreetMap can be seen on Figure 6.

The accuracy of the quantities entering within individual categories of GCPs can be roughly estimated according to the following

Table 2. Weight estimates of particular GCP categories.

GCP Category	Weight
church	1
château, yard, farmstead, monastery	1.4
mill	2
path crossing	more than 3

(the weight of churches represents the value of 1 and the other weights are related to it).

To give an idea of the degree of accuracy of the mosaic alignment (especially in relation to the data layer on the Mapire.eu website), a new control layer of more than 50 GCPs was created, which were approximately evenly distributed over the area of Bohemia. When designing these points, those that had already served as GCPs of the main IRLS adjustment method or points in their vicinity were excluded in order to avoid influencing the results of the comparison. The GCPs used for Mapire.eu are not known, however, given the small number of them, it can be assumed that the probability of using one of them in the control comparison is very low.

All types of usable objects were chosen as points in equal numbers (with the exception of churches, which practically had been completely used in the main adjustment before) and, in addition to the current maps, the maps of the Second Military Survey were also used to identify the objects. The coordinate differences between the ideal position and the position within the adjusted map layer were tested.

From the results it is not possible to unambiguously decide which dataset is better adjusted. Fifty GCPs are not able to convincingly cover the territory of Bohemia and it is especially the Mapire.eu data—which seem to show relatively large local deformations—that would need more GCPs for quality assessment. The median of the position deviation is around 250 m in both cases. The mean value is slightly weaker to the detriment of Mapire.eu, but the difference is negligible. Additionally, while Mapire.eu shows a larger amount of minor deviations, the mosaic using IRLS does not show errors greater than 1300 m, which value is exceeded by Mapire.eu in four cases, and the deviations end beyond the value of 2 km here.

We visualized the resulting mosaic in the ArcGIS Pro environment as a mosaic dataset. This data type allows storing raster data along with their footprints, which define the mask of the image in the mosaic. The footprints are also the output of our software, which automatically generates shapefiles from transformed map sheet corner coordinates, creating a new map sheet definition. These generated shapefiles are then paired with the footprints and used to crop the out-of-frame image of individual map sheets.

When adjusted map frames are available for individual map sheets, it is possible to create a seamless map layer using a mosaic dataset (we used the mosaic dataset of the ArcGIS software, ESRI platform) similar to that available on the mapire.eu website. The Mosaic Dataset allows a number of settings for the resulting visual display, rendering order, color corrections, rendering speed, etc. Using ESRI tools, it is possible to publish a raster mosaic as a web map service, both OGC-typed and proprietary.

The Image service, which is available in the ArcGIS Server extension named Image Server, was used to publish the First Military Survey mosaic. This type of service provides fast tiled data, and, at the same time, this map layer can be published as OGC standard WMTS.

4. Discussion

In the article, we introduced the method of joint georeferencing of a multi-sheet map work with the conditions of continuity of the adjacent edges of map sheets. It can be used where the parameters of the used cartographic projection are not known. The methods of georeferencing old maps using simple geometric transformations and the adjustment using LSM are generally known and used. Alternatively, interpolation algorithms for rubber-sheeting methods are applied. Examples of georeferencing single-sheet early maps can be found in [30,31]. In our case, we work with multi-sheet map works, which is significantly more complicated. Mathematically, it is possible to use the conditions of a sheet continuity during the adjustment [5,6], georeferencing the entire mosaic merged in the original coordinate system [1,4], or use a local transformation method based on, e.g., grid shift binary (GSB) as a second step [9]. For the georeferencing of the First Military Survey, the Arcanum project used the method of georeferencing the whole mosaic of each province using 30–50 GCPs and a polynomial transformation in as early as 2011 [32]. The accuracy of such georeferencing was not great. This method was later supplemented by a local type of transformation. The transformation of the whole mosaic using the local GSB method after adjusting the parameters of the estimated Cassini’s projection was also used for georeferencing Moldavia and Walachia [33]. The accuracy of local methods reaches average deviations on GCPs of around 200 m. Roughly corresponding accuracy is achieved by our method, which, however, does not use Cassini’s projection and GSB, but only a polynomial transformation with the conditions of the continuity of map sheets. Therefore, the comparison of the accuracy of both approaches on a set of checkpoints will be described below.

A group of 50 control GCPs that did not match the points originally used for the global adjustment was used to assess accuracy. As presented above, the results do not clearly show that one of the methods would be significantly geometrically more accurate, or rather the set of GCPs used is not able to prove this. The data show the predominance of smaller deviations (sometimes in the order of meters or tens of meters, i.e., well below the total mean square error value) in the case of Mapire.eu, but they are balanced by relatively large positional differences, in some cases exceeding 1 km (which practically represents a borderline value for IRLS) and extending up to 2 km. This coincides with the visual observation where some parts of Bohemia show significant positional differences and are strongly deformed. When looking at adjusted map sheets using IRLS, it is clear that substantial deformations of the sheets, except for several peripheral positions, practically do not occur.

Certainly, within the research a large amount of GCPs was collected and the used mathematical method is more complex than the GSB method available in common GIS software. Anyway, the results of our adjustment are more resistant to gross errors and the deformation of the map sheets after the transformation is not as noticeable. Additionally, the transformation outside the convex hull of GCPs is much better using polynomial (or affine) transformation than using interpolation based on GSB. That is why we considered our approach to be beneficial.

The medians of positional deviations are practically the same in both cases, which may indicate a similar explanatory value of both approaches and different errors at individual points due to random fluctuations of both methods rather than due to a systemically defective approach of some of them. Multiple GCPs would be needed to make a good assessment of the method used on Mapire.eu, which will be done in further research. The difficulty of this approach is a considerable exhaustion of some sites in terms of identifiable points where their compliance can be guaranteed. Unfortunately, due to the Industrial Revolution, industrialization and associated transformations of the landscape, the course of many roads changed, settlements expanded and the rural landscape structure transformed. Fortunately, the time-proximity of the Second Military Survey maps, which are geometrically very accurate (at least in the order of magnitude), makes it possible to better capture the image rendered on the sheets of the First Military Survey and compare it with reality with a higher dose of certainty.

In general, it can be stated that since the beginning of the 19th century, new cartographic works began to emerge, the creation of which is based on geodetic networks. Their accuracy is then significantly higher, which can be seen in [34,35]. The approaches to their georeferencing are then much simpler, especially with the use of the knowledge of the parameters of cartographic projections and the layouts of map sheets. The results of our article are, therefore, usable especially for maps from the second half of the 18th century, when a great expansion of detailed topographic maps began. Their use seems problematic mainly due to the difficulty of georeferencing. If it is possible to publish the mosaics of entire provinces of the First Military Survey in sufficient quality using standard web map services, the usability of these maps will be significantly higher. Whether we are researching industrial areas that have changed significantly since the Industrial Revolution in the 19th century [36], or we are interested in researching, for example, historical hydrological data [37], historical topographic maps are an excellent source of information.

5. Conclusions

The study brings new possibilities for the georeferencing, publication, viewing, and analysis of an old cartographic work with a great informative value, which the maps of the First Military Survey undoubtedly represent. The aim was to present the possibilities of the connection and global coordinate location of a large cartographic work which can be modernized by digital knowledge and using the appropriate mathematical and statistical apparatus, although it has no geometric basis and was never intended as a global seamless map. Moreover, these maps are spatially placed so that they can be viewed by means of geographic information systems or, completely intuitively, used as an underlying map layer.

With regard to the state of development of cartographic methods of their time, the First Military Survey maps represent a nicely rendered, but geometrically harder to capture depiction of Central Europe of the second half of the 18th century, which is valuable especially for the unity of the design, detail and its capability to visualize the landscape before the great transformations caused by the Industrial Revolution and associated societal and landscape changes. In this respect, the respective maps and similar map works published at the appropriate time in other European countries represent a unique and valuable probe into the late-Baroque to early-modern landscape and, in conjunction with other, later mappings of the same area, can serve as a suitable basis for studying the landscape development and comparison.

The new georeferencing method will be applied to all map works that have more than one map sheet and, at the same time, the parameters of the used cartographic projection are not known with sufficient accuracy. The results of our work can be compared with the results published within the Mapire.eu. Here, it was not possible to convincingly show whether the new method is better or worse in terms of mean errors, but at least it is comparable to the map layer published on the Mapire.eu portal. The fact remains that, using our transformation method, the map sheets are not deformed inside and they merely follow global transformation equations. Using the GSB methods may lead to distorted map sheets with, e.g., geographic names that cannot be read. Our method can also better control the quality of GCPs when outlying measurements are suppressed by the assignment of weights during the IRLS process.