Pareto-based evolutionary algorithms for the calculation of transformation parameters and accuracy assessment of historical maps

Computers & Geosciences 57 (2013) 124–132

Contents lists available at SciVerse ScienceDirect

Computers & Geosciences journal homepage: www.elsevier.com/locate/cageo

Pareto-based evolutionary algorithms for the calculation of transformation parameters and accuracy assessment of historical maps F. Manzano-Agugliaro a,n, C. San-Antonio-Gómez b, S. López c, F.G. Montoya a, C. Gil c a

Department of Engineering, University of Almería, 04120 Almería, Spain Department of Cartographic Engineering, Geodesy and Photogrammetry, Polytechnic University of Madrid, 28040 Madrid, Spain c Department of Computer Science, University of Almería, 04120 Almería, Spain b

art ic l e i nf o

a b s t r a c t

Article history: Received 10 October 2012 Received in revised form 18 February 2013 Accepted 9 April 2013 Available online 29 April 2013

When historical map data are compared with modern cartography, the old map coordinates must be transformed to the current system. However, historical data often exhibit heterogeneous quality. In calculating the transformation parameters between the historical and modern maps, it is often necessary to discard highly uncertain data. An optimal balance between the objectives of minimising the transformation error and eliminating as few points as possible can be achieved by generating a Pareto front of solutions using evolutionary genetic algorithms. The aim of this paper is to assess the performance of evolutionary algorithms in determining the accuracy of historical maps in regard to modern cartography. When applied to the 1787 Tomas Lopez map, the use of evolutionary algorithms reduces the linear error by 40% while eliminating only 2% of the data points. The main conclusion of this paper is that evolutionary algorithms provide a promising alternative for the transformation of historical map coordinates and determining the accuracy of historical maps in regard to modern cartography, particularly when the positional quality of the data points used cannot be assured. & 2013 Elsevier Ltd. All rights reserved.

Keywords: Pareto front Evolutionary algorithm Historical cartography Accuracy assessment Transformation parameters Coordinate transformation

1. Introduction Historical maps are an important part of our cultural heritage (Jenny and Hurni, 2011). These maps not only represent valuable physical artefacts but also provide an important information source for historians and geographers, who frequently incorporate historical data into Geographical Information System (GIS) (Weir, 1997; Audisio et al., 2009). The scales, coordinate systems, projections, and surveying and mapping techniques used in historical maps vary widely (Podobnikar, 2009). The different reference systems employed in different historical maps often require coordinate transformations between them (Tierra et al., 2008). As historical maps typically exhibit higher degrees of inaccuracy and uncertainty compared to contemporary cartographic databases, it is not surprising that these two issues are of particular concern in historical cartography studies and historical GIS applications (Plewe, 2003). Accuracy analysis of early maps is therefore an important topic in historical cartography (Harley, 1968). The positional accuracy of a point on a map is defined as the difference between its recorded location on the map and its actual location on the ground or location on a source of known higher accuracy (Tucci and Giordano, 2011).

n

Corresponding author. Tel.:+34 950015396; fax: +34 950015491. E-mail address: [email protected] (F. Manzano-Agugliaro).

0098-3004/$ - see front matter & 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.cageo.2013.04.010

The coordinate method calculates the correlation between two sets of map coordinates of points identified by modern latitudes and longitudes (Tobler, 1966). The deviation between each computer-generated point (digitised from an early map) and the corresponding point on the modern map can be displayed as a vector indicating the direction and magnitude of the error (Ravenhill and Gilg, 1974). A related problem is the transformation of coordinates between two different geodetic reference systems. A set of points with known coordinates in both reference systems is used to obtain the transformation parameters. The Helmert transformation is one example of a commonly used method for geodesic transformations between two reference frames (Vanícek and Krakiwsky, 1986; Vanícek and Steeves, 1996). The quality of a transformation between two sets of coordinates depends on the positional quality of the points used to calculate the transformation parameters. However, a problem arises when there is substantial spatial uncertainty in the historical data points. In this case, certain points with gross errors will not be used in the calculation of the accuracy of the entire map. When the inaccuracy of a measurement is not objectively known, as is frequently the case for historical maps, the measured feature is defined as uncertain (Hunter and Goodchild, 1993). Although uncertainty and inaccuracy are ontologically distinct concepts, it is often difficult to measure the two separately in practice, particularly in the context of historical maps (Hu, 2010).

F. Manzano-Agugliaro et al. / Computers & Geosciences 57 (2013) 124–132

In previous papers (Manzano-Agugliaro et al., 2012), the coordinates obtained from the historical map were displaced by the average latitude and longitude error to correct the absolute displacement error in the historical map, i.e., to correct the georeferencing error. Any points whose displacements from the corresponding points in the modern map exceeded a specified distance were then eliminated from the analysis. These points were considered to contain gross errors, either because of incorrect identification with the current points or due to a gross error in the historical map. The final estimate of the map accuracy is obtained based on the root mean square (RMS) displacement after the removal of the points with gross errors. In this historical cartography, we cannot be certain that any given point is more accurate than another. To displace the coordinates on a historical map to positions such that the final (RMS) errors are minimal, many combinations must be performed, each time discarding the points that exceed a fixed maximum RMS for a gross or unacceptable error and then calculating the new errors. A historical map may exhibit rotation (due to projection effects) as well as horizontal and vertical scale errors in addition to latitudinal and longitudinal displacement, making the number of possible combinations that must be considered even higher. If we also aim to discard as few points as possible from the historical map, then the problem becomes multi-objective. The multi-objective nature of these mapping problems makes the decision-making process complex. Fortunately, the increase in computational resources in recent years has allowed researchers to develop efficient computational algorithms for handling complex optimisation problems. In particular, multi-objective evolutionary algorithms (MOEAs) are known for their ability to optimise several objective functions simultaneously to provide a representative Pareto front, which is a set of problem solutions representing a trade-off between the objectives (Márquez et al., 2011). The aim of this paper is to assess the performance of evolutionary algorithms in determining the accuracy of historical maps in regard to modern cartography.

2. Concepts in multi-objective optimisation Many of the problems faced in engineering and other disciplines are optimisation problems. Many optimisation problems are difficult to solve because of features such as non-linear formulations, constraints, and NP-hard complexity. The techniques for solving optimisation problems generally fall into two categories. Exact techniques provide the optimal solution to a given problem but are impractical for handling NP-hard problems because of prohibitive computation time and/or memory requirements. Non-exact techniques, such as metaheuristic methods (Glover and Kochenberger, 2003), provide satisfactory (though not necessarily optimal) solutions to complex problems in a reasonable amount of time. Most computational optimisation research has been focused on solving single-objective problems, including constraints in some cases. Nevertheless, many real-world problems require the simultaneous optimisation of several competing objectives. Several authors have proposed multi-objective algorithms based on Pareto optimisation (Baños et al., 2009, 2011) to solve these multiobjective optimisation problems (MOPs). In contrast to single-objective optimisation problems, the solution to a MOP consists of a set of non-dominated solutions known as the Pareto optimal set rather than a single solution. A solution that belongs to this set is said to be a Pareto optimum and when the solutions of this set are plotted in objective space, they are collectively known as the Pareto front. Obtaining the Pareto front is the main goal in multi-objective optimisation.

125

Evolutionary algorithms are particularly desirable in the solution of multi-objective optimisation problems because they simultaneously handle a set of possible solutions, yielding an entire set of Pareto optimal solutions in a single run of the algorithm (Coello, 1999). While a single-objective optimisation problem may have only one optimal solution, a multi-objective optimisation problem may have an uncountable set of solutions. When evaluated, these solutions produce vectors whose components represent a tradeoff between the various objectives of the problem. At this point, an expert in the problem, referred to as the “decision maker” (DM), must implicitly choose one (or several) solutions by selecting one or more of these vectors. For the sake of generality, we can assume that although in the following definitions (Talbi, 2009) we use the term minimisation for all of the objectives, there are problems in which several (or all) of the objectives must be maximised to generate the optimum solutions. A multi-objective optimisation problem can be informally defined (Osyczka, 1985) as the problem of finding: “a vector of decision variables that satisfies the constraints and optimises a vector function whose elements represent the objective functions. These functions form a mathematical description of the performance criteria, which are usually in conflict with each other. Hence, the term ‘optimise’ means finding a solution that would give the values of all of the objective functions that are acceptable to the decision maker.” Definition 1. Multi-objective optimisation problem. A MOP is defined as MOP≜min FðxÞ ¼ ðf 1 ðxÞ; f 2 ðxÞ; …; f M ðxÞÞ; x∈S where M≥2, is the number of objectives, x ¼ ðx1 ; …; xk Þ is the vector representing the decision variables and S represents the set of feasible solutions associated with equality and inequality constraints and explicit bounds. FðxÞ ¼ ðf 1 ðxÞ; f 2 ðxÞ; …; f M ðxÞÞ is the vector of objectives to be optimised. Since in real-world MOPs the criteria are usually in conflict, there is a need to establish other concepts to consider optimality. In that sense, a partial order relation, known as Pareto-dominance relation can be defined. Definition 2. Pareto dominance. An objective function vector is defined as a dominating vector (denoted by FðxÞ≺F′ðxÞ) if and only if no component of it is smaller than the corresponding component of and at least one component of it is strictly smaller, that is, ′

′

∀i∈f1; …; M g : f i ðxÞ ≤f i ðxÞ∧∃i∈f1; …; M g : f i ðxÞ of i ðxÞ Definition 3. Pareto Optimality. A solution xn ∈S is Pareto optimal if for every, does not dominate, that is FðxÞ≺Fðxn Þ. The concept of Pareto optimality is directly related to the dominance concept and was initially proposed by Edgeworth (1881) and extended by Pareto (1896). Definition 4. Pareto optimal set. For a given MOP, the Pareto optimal set is defined as ℘n ¼ xn . Definition 5. Pareto front. For a given MOP and its Pareto optimal   set ℘n , the Pareto front is defined as ℘ℑn ¼ ℑðxn Þ ℘n ¼ xn . Thus the Pareto front is the image of the Pareto optimal set in the objective space, and as such it constitutes the main goal of multi-objective optimisation. As a consequence, all the solutions belonging to the Pareto front are indifferent because any change in any objective xi that would improve it degrades any amount of other objectives. The concept of indifference is illustrated in Fig. 1.

126

F. Manzano-Agugliaro et al. / Computers & Geosciences 57 (2013) 124–132

Fig. 1. Pareto-dominance relations in a two-objective problem.

3. Data: historical cartography analysed The map used in this paper is taken from the first edition of the Atlas Geográfico de España (AGE) by Tomás López (1730–1802), published in 1804. The historical context of the cartographic work of Tomás López has its origins in the arrival of Philip V to the Spanish throne. His prime minister, the marquis of La Ensenada (1702–1781) devised a set of reforms to improve the administration of country's mainland territories and called for a national cartographic base on which to surface and quantify taxes (SanAntonio-Gómez et al., 2011). Tomas Lopez was one of two men chosen for this important task. Between 1752 and 1760, Lopez (a pensioner of the king) served as one of the 18th century's most prestigious cartographers. He lived in Paris and followed the teachings of Jean Baptiste Bourgignon d’Anville (1697–1782). The cartographic method used by Lopez, known as desk cartography, was based on a questionnaire sent to each village priest. In addition, each priest was required to provide a sketch of his village and 2 leagues surrounding it. Lopez composed a new sketch of the region based on the sketches provided by the village priests, other existing maps from the 16th, 17th and 18th centuries, and vast quantities of other documents (answers to questionnaires, local histories, geographical descriptions, and cartographic sketches), which have been preserved in the National Library of Spain (San-Antonio-Gómez et al., 2011). Although he never performed any fieldwork, Lopez was aware of the latest advances in scientific cartography and was acknowledged for his cartographic skill. The AGE was the first comprehensive and detailed map of Spain. The accuracy of the maps by Lopez is generally considered to be low (Chias and Abad, 2012), although few studies have quantitatively addressed the accuracy of his maps compared to current maps. The AGE sheet used to for this study is number 64, corresponding to the “Reyno de Jaen” [Kingdom of Jaen], Fig. 2. This map was provided in digital format by the maps library at the Spanish Geographical Information Centre (CENIG); the page was scanned from the original map in the Atlas. The map is an engraving whose original dimensions were 39  42 cm; its frame is marked by integer coordinates centred on 381 latitude and 131 longitude, and the minutes are graduated in steps of five. The referenced longitude source is the Teide peak, as indicated in the top and central parts of the map. Tomas Lopez employed a spherical rather than an ellipsoidal Earth model, with an assumed Earth radius of 6340 km (Manzano-Agugliaro et al., 2012). He was aware that the Earth more closely followed an ellipsoidal model but argued that no major discrepancies were introduced by the use of the spherical model (Lopez, 1795). The sheet georeferencing was performed using the ArcGis v9.3 software. The 4 central points of the outer frame of the map were

used as input because their coordinates are labelled. The RMS obtained in this operation was 0.003561, corresponding to 0.396 km. Following the georeferencing, all of the settlements on the map were digitised (see the example in Fig. 3). A modern reference ellipsoid (European Datum 1950) was used to obtain a list of settlements with their original names and latitude and longitude coordinates on the original map. In total, 233 settlements of varying hierarchical size were digitised, including 163 ciudades (cities and large towns), villas (towns) and lugares (locations, places or sites). These settlements were compared with current cities; 85 of the settlements were identified by name, while no correspondence could be found for the other 75 settlements. Only 53% of the ciudades indicated by Tomás López coincide with current towns and cities. Bear in mind that comparing maps made 210 years apart is no easy task; however, we find that approximately half of the historical settlements have either disappeared or had their names changed.

4. Methodology 4.1. Procedure The procedure can be divided into two stages (see Fig. 4): Stage 1: Data filtering. The locations on which the algorithm is to be run are filtered as follows. First, the coordinates of the present-day locations and Tomas Lopez's towns are read in, and the data are stored in memory. We then calculate the distances between the locations as given by the two data sources. Once these distances have been obtained, we can calculate the mean error of the set of historical towns by taking the arithmetic mean. This error is used in the final step to select the locations for which the calculations will be performed; eliminating any points such that the distance between the current and historical (as provided by Lopez) coordinates exceeds 20 km. Through this selection, we eliminate those input values that are too far removed from the true location and classify them as data collection errors. Stage 2: Optimisation using evolutionary algorithms. The evolutionary algorithm is run to obtain the data to be analysed. The following steps are performed. 2-1: Generation of the initial population. The first generation of individuals to be evolved in each iteration (generation) are specified. The initial population is normally generated by assigning random genes (values) to the individuals. 2-2: Evaluation of the initial population. The evaluation consists of calculating the distances between the coordinates

F. Manzano-Agugliaro et al. / Computers & Geosciences 57 (2013) 124–132

127

Fig. 2. Kingdom of Jaén (Tomás López 1787).

of the present-day settlements and those obtained from the maps of Tomas Lopez, recalculating the mean error in each solution and subsequently arranging them in ascending order. 2-3: Selection. Natural selection is simulated in this step; the individuals who adapt best to the situation are chosen. The most well-adapted individuals comprise the solutions that exhibit the smallest mean errors. 2-4: Reproduction. This phase involves generating new individuals (solutions) from the individuals selected in the previous stage. The individuals generated inherit characteristics from their parents and potentially yield offspring that can improve the results of their forebears by combining the best genes (values) from each parent. 2-5: Evaluation of the new individuals. Once the offspring have been obtained, we have a new population composed of the previous generation plus the offspring generated in the previous step. The individuals are reassessed by following the same procedure as applied previously (in step 2-2). 2-6: Elitism. This is the final operation in the generational cycle. The evaluation of the previous step is considered, and the

solutions (individuals) with the least favourable results are eliminated. A more natural expression of this process would be “the survival of the fittest”, with each successive generation attaining a higher likelihood of achieving fit offspring by maintaining only the most well-adapted individuals. 2-7: Finally, a check is performed to assess whether the required optimisation criteria are satisfied. The satisfaction of the criteria is usually limited by the number of generations used. If the objectives have not been met, then we return to the selection step. 2-8: Once the objectives have been met, the individuals who have survived are returned. These individuals form our new set of solutions.

4.2. Problem formulation The programme solves a multi-objective optimisation problem using parallel evolutionary algorithms. The evolutionary algorithm used in this case is msPESA (Gil et al., 2007), with 50,000 iterations

128

F. Manzano-Agugliaro et al. / Computers & Geosciences 57 (2013) 124–132

Fig. 3. Example of digitisation of AGE, page 64.

(hence 50K). This algorithm has been successfully employed in several engineering applications (Alcayde et al., 2011; Baños et al., 2009; Montoya et al., 2010). It is a hybrid algorithm that combines a variant of the archive update strategy of PESA (Corne et al., 2000) and adopts several concepts from the selection mechanism in NSGA-II (Deb et al., 2002). The msPESA algorithm stores the nondominated solutions of each iteration in the archive of candidate solutions, without removing the solutions that are dominated by new solutions, to improve the genetic variability. Once the archive is full, a squeeze factor (Gil et al., 2007) is used to remove excess solutions from the archive. There are two objectives to be optimised. The mean error is minimised when calculating the distances between the presentday urban areas and the towns of Tomas Lopez, and the number of locations included in the calculations is maximised. Software objectives:

 Minimise the average distance between actual coordinates and Tomás López's coordinates. The distance between two points is calculated using the haversine formula, which gives great-circle distances between two points on a sphere from their longitudes and latitudes. It assumes a spherical Earth and ignores ellipsoidal effects but remains particularly well-conditioned for numerical computation even at small distances (Montavont and Noel, 2006):   n n D ¼ ∑ haversin dri ¼ ∑ haversinðϕi2 −ϕi1 Þ i¼1

i¼1

þcosðϕi1 Þ cos ðϕi2 Þhaversinðψ i2 −ψ i1 Þ

where, D is the average distance (the mean error in km between the coordinates of Tomás López's towns and the equivalent for present-day locations), haversin is the haversine function:   θ 1−cos θ haversinðθÞ ¼ sin2 ¼ 2 2 di is the distance between actual coordinates and Thomas López's coordinates, r is the radius of Earth (6371 km), ϕi2 ; ϕi1 : Latitude of point i1, Latitude of point i2, and ψ i2 ; ψ i1 : Longitude of point i1, Longitude of point i2. Let h denote the haversin (d/r). One can then solve for d either by simply applying the inverse haversine or by using the arcsin (inverse of sine) function: d ¼ r haversin−1 ðhÞ ¼ 2r arcsinðhÞ

 Minimise the number of excluded locations ðli Þ. n

N ¼ ∑ li i¼1

where N is the number of locations that have not been excluded from analysis. Calculation parameters (genotype) The parameters change in the evolutionary algorithm until the solution is arrived at (what is known as the genotype). These are the parameters that are used to perform the three transformations

F. Manzano-Agugliaro et al. / Computers & Geosciences 57 (2013) 124–132

on the map coordinates, see Fig. 5, and obtain the final coordinates (Lat′, Lon′) – Rotation origin (Lat0, Lon0): the centre coordinates where the points are rotated, calculated when beginning the solution. This is the centre of mass of the towns: Lon0 ¼

1 n ∑ Loni ; ni¼1

Lat 0 ¼

1 n ∑ Lat i ni¼1

Using the centre of mass, we obtain coordinates that will be approximately at the centre of the nucleus of towns, so the variation is small. If we use as the rotation origin the

129

coordinates of 01 latitude and 01 longitude, we would then have to increase the latitudinal and longitudinal displacements. – Rotation: the rotation applied, in degrees (α). In tests of programming, no rotation had greater than 0.21 for the obtained solutions. Since the evolutionary algorithm does not discriminate any solution a priori, a limitation to 0.51 was made in the software, this limitation narrows the search space solution, and the 50,000 iterations are used within a space solution with best possible solutions. Lat′ ¼ Lat  cosðαÞ−Lon  sinðαÞ Lon′ ¼ Lon  sinðαÞ þ Lon  cosðαÞ

Fig. 4. Methodology flowchart.

130

F. Manzano-Agugliaro et al. / Computers & Geosciences 57 (2013) 124–132

– Latitudinal displacement in sexagesimal degrees (ΔLat). Sum of the original value plus that calculated by the algorithm. Lat″ ¼ Lat′ þ ΔLat – Longitudinal displacement in sexagesimal degrees (ΔLon). Lon″ ¼ Lon′ þ ΔLon – Horizontal scale factor (SLon): the scale based on the coordinates, on the longitudinal axis. Limited to 20% in order for the calculations not to be shifted too much (between 0.8 and 1.2). Lon′′′ ¼ Lon″ þ SLon – Vertical scale factor (SLat): the scale based on the coordinates on the latitudinal axis. Lat′′′ ¼ Lat″ þ SLat

In the evaluation process we must calculate both objectives, by applying the changes described in the genotype for each individual to each of the coordinates on the Tomás López map (TL) for each location in the individual (solution) for the town. This gives us transformed coordinates, for each of which we calculate the distance to the equivalent town in the present-day nucleus. Calculate the mean error with the arithmetical mean: Error ¼

1 n ∑ D ni¼1 i

Di ¼distance between locations (TL map and current). Maximise the number of locations. This is done by minimising the locations that are eliminated, i.e.: Number of locations eliminated ¼ total number of locations −number of locations used

When we run the programme based on the msPesa 50K algorithm on the 83 towns in the Kingdom of Jaen, we obtain the solutions that are partially shown in Table 1. Those solutions that eliminate the smallest number of towns are included. The Pareto front obtained consists of the best results from among 50,000 combinations of possible solutions, each based on 80 crosses (crossovers) and 10 mutations (mutation). The algorithm performs evaluations each time it runs, and the best solution is maintained for the Pareto front. The Pareto front is therefore constructed from the best solution for 1 excluded town, the best solution for 2 excluded towns, and so forth. A total of 50,000 possible solutions have been tried when the procedure terminates. In Table 1, the algorithm has also been compared with a classical solution in which possible systematic error is eliminated from the historical map by applying a latitudinal and longitudinal displacement equal to the mean error between the two maps. This classical solution is superior only to the evolutionary solution in which one point is removed. All of the solutions obtained by the evolutionary algorithm with removal of more than one point provide superior results compared to the classical algorithm. Fig. 6 shows the Pareto front for the solutions obtained. Observe that by eliminating only two locations from the original map, the error is reduced to 5.4 km. The error decreases very slowly as 4, 5 and 6 settlements are eliminated and falls somewhat more rapidly (to 4.8 km) when 8 locations are eliminated. It is then necessary to eliminate many more locations (19) to achieve a slight further reduction of the error to 4.4 km. In general, we find that the solutions that use the algorithm require a negative latitudinal displacement and positive longitudinal displacement, as for the classical solution, although the required displacements are smaller when rotation is allowed. The longitudinal displacement is greater than the latitudinal displacement, confirming the past lack of resources to measure the geodesic longitude between two points. The latitude was easier to

SLat

Lat

Δ Lon

(Lon0, Lat0)

Δ Lat

Average Error (Km)

Solution map

SLon

5. Results

α SLon

Original map SLat

8.5 8 7.5 7 6.5 6 5.5 5 4.5 4

Pareto front

0

Lon

2

4

6

8

10

12

Classical solution

14

16

18

20

Excluded Towns

Fig. 5. Transformations applied to the original map.

Fig. 6. Pareto front with evolutionary solutions and classical solution.

Table 1 Adjusted results obtained with evolutionary algorithm msPesa 50K. Classical Solution

Error (km) N (excluded) ΔLat (deg) ΔLon (deg) Lat0 Lon0 α (deg) SLon SLat

7.5427 0 −0.0229 0.1062 0 0 0 0 0

Solution with evolutionary algorithm 1

2

3

4

5

6

8.0131 1 −0.0199 0.0574 37.98101 −3.63599 −0.08998 0.8 1.0733

5.4328 2 −0.0140 0.0869 37.98101 −3.63599 −0.09230 1.0304 1.1234

5.3742 4 −0.0143 0.0863 37.98101 −3.63599 −0.09290 1.0714 1.1131

5.3230 5 −0.0118 0.0887 37.98101 −3.63599 −0.10052 0.9847 1.1126

5.1039 6 −0.0118 0.0887 37.98101 −3.63599 −0.06780 1.0541 1.1391

4.7941 8 −0.0158 0.0812 37.98101 −3.63599 −0.05900 1.0458 1.1420

F. Manzano-Agugliaro et al. / Computers & Geosciences 57 (2013) 124–132

measure in the past. The solutions exhibiting the best behaviour are those with negative rotation, although in all cases this rotation is very small. We see that that the longitudinal scale factor (X axis) can be either positive or negative, whereas the Y scale factor (latitude) always remains positive. This may suggest that the mapmaker reduced this magnitude too much when transferring it to the map. Of the solutions obtained with the evolutionary msPesa 50K algorithm, the first three eliminate approximately the same number of towns. The advantage of having a solutions front is that it enables

131

us to choose the most appropriate solution based on the error tolerance. For the historical map under study, we have succeeded in reducing the error by 40% by eliminating only 2% of the points. We find that in all of the solutions obtained from the evolutionary algorithm, the same towns have been eliminated, although some runs eliminated more towns than others to reduce the error. This indicates that the eliminated towns exhibit the lowest cartographic precisions and highest uncertainties on the historical map from among all of the combinations of parameters over the 50,000 runs.

Fig. 7. Classical solution.

Fig. 8. Evolutionary solution 2 (two towns excluded).

132

F. Manzano-Agugliaro et al. / Computers & Geosciences 57 (2013) 124–132

Fig. 7 shows the cartographic precision of the TL map, considering only the classical solution, for the mean error displacement. The errors with this solution have been compared with the present-day map. We see that the southern part of the map is the area with the greatest error of up to 17.49 km. Fig. 8 shows solution 2 obtained with the evolutionary msPesa 50K algorithm, which eliminated two points (towns) to calculate the coordinate transformation parameters (displacements, rotation and scale factors). The towns that do not belong to solution 2 – Toya in the east and Santa Ana in the south – are marked with black triangles. The multi-objective optimisation shown in Fig. 8 enables superior solutions compared the classical solution, with greater homogeneity in the results. 6. Conclusions The multi-objective nature of the historical cartography optimisation problem allows the problem to be tackled using Paretobased optimisation techniques, such as the msPesa 50K algorithm. The advantage of using Pareto fronts is that a solutions front, rather than a single solution, is obtained. The most suitable solution can therefore be adopted according to the specific requirements of the historical map in question. Multi-objective optimisation allows any point or town with a high cartographic uncertainty to be discounted in the calculation of the transformation parameters. Using genetic algorithms, the linear error has been reduced by 40% for the map under study by eliminating only 2% of the data points. Evolutionary algorithms provide a promising alternative for the transformation of historical map coordinates and determining the accuracy of historical maps in regard to modern cartography, particularly when the quality of the position information cannot be assured for the set of points used. Acknowledgements This study was funded by project HAR2009-12937 (GIS Systematic Analysis for the Planimetric Accuracy of the Geographic Atlas of Spain of Tomas Lopez, 1804), a project of the Spanish Ministry of Science and Innovation. References Alcayde, A., Baños, R., Gil, C., Montoya, F.G., Moreno-Garcia, J., Gómez, J., 2011. Annealing-tabu PAES: a multi-objective hybrid meta-heuristic. Optimization 60 (12), 1473–1491. Audisio, C., Nigrelli, G., Lollino, G., 2009. A GIS tool for historical instability processes data entry: an approach to hazard management in two Italian Alpine river basins. Computers and Geosciences 35 (8), 1735–1747. Baños, R., Manzano-Agugliaro, F., Montoya, F.G., Gil, C., Alcayde, A., Gómez, J., 2011. Optimization methods applied to renewable and sustainable energy: a review. Renewable and Sustainable Energy Reviews 15 (4), 1753–1766. Baños, R., Gil, C., Reca, J., Martinez, J., 2009. Implementation of scatter search for multi-objective optimization: a comparative study. Computational Optimization and Applications 42, 421–441.

Coello, C., 1999. A comprehensive survey of evolutionary based multi-objective optimization techniques. Knowledge and Information Systems 1 (3), 269–308. Corne, D., Knowles, J.D., Oates, M.J., 2000. The Pareto envelope-based selection algorithm for multiobjective optimization. Parallel problem solving from nature PPSN VI. Lecture Notes in Computer Science 1917, 839–848. Chias, P., Abad, T., 2012. The art of describing the territory: historic maps and plans of the bridge of Alcántara (Cáceres, España). Informes de la Construción 64 (extra), 121–134. Deb, K., Pratab, A., Agrawal, S., Meyarivan, T., 2002. A fast elitist non-dominated sorting genetic algorithm for multi-objective optimization: NSGA-II. IEEE Transactions on Evolutionary Computation 6, 182–197. Edgeworth, F., 1881. Mathematical psychics: an essay on the application of mathematics to the moral sciences. C. Kegan Paul and Co., London. Gil, C., Márquez, A., Baños, R., Montoya, M.G., Gómez, J., 2007. A hybrid method for solving multi-objective global optimization problems. Journal of Global Optimization 38 (2), 265–281. Glover, F.W., Kochenberger, G.A., 2003. Handbook of Metaheuristics. Kluwer, Dordrecht. Harley, J.B., 1968. The evaluation of early maps: towards a methodology. Imago Mundi 22, 62–74. Hu, B., 2010. Application of geographic information systems (GIS) in the history of cartography. Proceedings of World Academy of Science. Engineering and Technology 66, 1562–1565. Hunter, G., Goodchild, M.F., 1993. Managing uncertainty in spatial databases: putting theory into practice. Journal of Urban and Regional Information Systems Association 5, 55–62. Jenny, B., Hurni, L., 2011. Studying cartographic heritage: analysis and visualization of geometric distortions. Computers and Graphics 35 (2), 402–411. Lopez, T. 1795. Principios geográficos aplicados al uso de los mapas. Ed. Don Benito Cano. Manzano-Agugliaro, F., Martinez-García, J., San-Antonio-Gómez, C., 2012. GIS analysis of the accuracy of Tomas Lopez´s historical cartography in the Canary Islands (1742–1746). Scientific Research and Essays 7 (2), 199–210. Márquez, A.L., Baños, R., Gil, C., Montoya, M.G., Manzano-Agugliaro, F., Montoya, F.G., 2011. Multi-objective crop planning using Pareto-based evolutionary algorithms. Agricultural Economics 42 (6), 649–656. Montoya, F.G., Baños, R., Gil, C., Espín, A., Alcayde, A., Gómez, J., 2010. Minimization of voltage deviation and power losses in power networks using Pareto optimization methods. Engineering Applications of Artificial Intelligence 23 (5), 695–703. Montavont, J., Noel, T., 2006. IEEE 802.11 handovers assisted by GPS information. IEEE International Conference on Wireless and Mobile Computing, Networking and Communications 2006, WiMob 2006, Art. no. 1696358, pp. 166–172. Osyczka, A., 1985. Multicriteria optimization for engineering design. In: Gero, John S. (Ed.), Design Optimization. Academic Press, pp. 193–227. Pareto, V., 1896. Cours d'economie politique. Rouge Lausanne, Switzerland1896. Plewe, B.S., 2003. Representing datum-level uncertainty in historical GIS. Cartography and Geographic Information Science 30, 319–334. Podobnikar, T., 2009. Georeferencing and quality assessment of Josephine survey maps for the mountainous region in the Triglav National Park. Acta Geodaetica et Geophysica Hungarica 44 (1), 49–66. Ravenhill, W., Gilg, A., 1974. The accuracy of early maps? Towards a computer-aided method. Cartographic Journal 11 (1), 48–52. San-Antonio-Gómez, C., Velilla, C., Manzano-Agugliaro, F., 2011. Tomas Lopez's geographic atlas of Spain in the Peninsular war: a methodology for determining errors. Survey Review 43 (319), 30–44. Talbi, E., 2009. Metaheuristics: from design to implementation. John Wiley & Sons, Inc., New York. Tierra, A., Dalazoana, R., De Freitas, S., 2008. Using an artificial neural network to improve the transformation of coordinates between classical geodetic reference frames. Computers and Geosciences 34 (3), 181–189. Tobler, W.R., 1966. Medieval distortions: the projections of ancient maps. Annals of the Association of American Geographers 56 (2), 60–351. Tucci, M., Giordano, A., 2011. Positional accuracy, positional uncertainty, and feature change detection in historical maps: results of an experiment. Computers, Environment and Urban Systems 35 (6), 452–463. Vanícek, P., Krakiwsky, E., 1986. Geodesy: The Concepts, second ed. Elsevier Science B.V., Amsterdam697. Vanícek, P., Steeves, R., 1996. Transformation of coordinates between two horizontal geodetic datums. Journal of Geodesy 70 (11), 740–745. Weir, M.J.C., 1997. A century of forest management mapping. Cartographic Journal 34 (1), 5–12.