Cartometric analysis of selected works by Leonardo da Vinci

. This research presents the results of the cartometric analysis of two Leonardo da Vinci’s works, i.e. the Mapa mundi and A map of Imola . The analysis of the Mapa mundi was conducted by employing various distortion measures and creating maps to show their distribution by distortion isograms. A map of Imola was analysed utilising MapAnalyst, taking advantage of the software’s tools, such as distortion grids, isolines of scale and local map rotation angles.


Introduction
Leonardo da Vinci (1452-1519) lived and worked in the Renaissance, which was a period of great discoveries and progress in science, technology, and art.This period was characterised by a significant progress in cartography, to which Leonardo also made a contribution.The most well-known cartographic works by Leonardo da Vinci include: A map of Imola, the map of the Chiana Valley, map of western Tuscany, the map of the Pontine Marshes and the Mapa mundi.Apart from these works, he also created numerous other maps and topographic sketches, mainly for military, hydrological, and engineering purposes.The cartographic works of Leonardo da Vinci were described in several publications, including Isaacson (2017), Kemp (2007), Puceković (2013), Ristujczina (2020), Snyder (1993), andTyler (2017).
The paper presents the results of the cartometric analysis of the Leonardo da Vinci's Mapa mundi and A map of Imola.These works were selected due to the availability of modern analytical methodologies and software that support such analysis.Although Leonardo da Vinci did not provide mathematical functions describing the projection he used for the Mappa mundi, they were introduced later, e.g. in Bower (2012).The projection functions provide a basis for conducting a detailed analysis of a projection taking advantage of modern methods and computer software.On the other hand, A map of Imola has a geometric basis.It was created according to detailed field measurements, so its cartometric properties can be defined.The analysis of the map is limited only to the area inside the city walls.As the applied methods and software require the identification of the corresponding points on the historical map and on an equivalent contemporary map, it would be difficult to identify such points outside this area.Other maps created by Leonardo da Vinci cover large areas drawn from a bird's eye view, which were not based on such precise measurements as the map of Imola.
A map of Imola (Figure 1) was drawn in 1502.As opposed to most town maps from this era, which were made using oblique projection, A map of Imola employed the orthogonal projection.The author measured the length of streets using steps.The directions were measured The Mappa mundi (Figure 2), drawn in eight parts (octants), uses the Leonardo da Vinci's projection that has an original structure (Royal Collection Online, n.d.b).Its original surface is a sphere divided into octants, creating eight equilateral spherical triangles whose borders are the meridians and the equator (Figure 3).
In the Leonardo da Vinci's projection, the octant of the sphere is the Reuleaux triangle.It is a curved triangle whose sides are the arcs of circles with their centres and ends at the vertices of an equilateral triangle.A Reuleaux triangle is constructed by drawing an equilateral triangle and then drawing circles in the vertices of the triangle, of the radii equal to the sides of the triangle.The intersection of those circles forms the Reuleaux triangle (Figure 4).

Methodology for the analysis of the Leonardo da Vinci's projection
This section analyses the metric properties of the projection used by Leonardo da Vinci in his Mappa mundi.The shape of the graticule and the projection distortions were analysed.The metric properties of the projection used by Leonardo da Vinci were analysed with the use of projection functions according to Bower (2012): (1) where:  For the reverse projection, the formulas take the following form: The parametric scales in the Leonardo da Vinci's projection are as follows: The area distortion scale p is expressed by the same formula as μ λ .
The m and n extreme scales and angular distortions are calculated based on the parametric scales and the area scale with the use of formulas presented in relevant literature, e.g. in Pędzich (2014).The distortions of the whole area are determined by the so-called integral measures (Biernacki, 1949).They allow for the comparison of the distortions of the given area in several different cartographic projections.
For the purposes of this research, the mean square Airy distortion was applied.It was calculated in the nodes of the cartographic grid, with the following formula: The analysis of the metric properties of the projection was limited to a single octant , as all other octants are projected in a similar way.

The analysis of the graticule in the Leonardo da Vinci's projection
The images of parallels in this projection are The graticule (ϕ = 5°, ∆λ = 5°) has the form presented in Figure 5.
The extension of the projection to the whole globe results in a rather unusual representation (Figure 6).One may notice, that the projection is not regular for the whole globe.The parallel of the is projected in form of a point.Parts of the images of the western and eastern hemi-π 2 spheres overlap, so this projection may be used only in a limited area between the parallels of latitude ϕ = 90° (1 -√3) and ϕ = 90° and the meridians λ = -135° and λ = 135°.For example, the graticule for ¼ of the globe is shown in Figure 7.

Local distortions in the Leonrdo da Vinci's projection
Calculating the measures of local distortions within the borders of the projected triangle does not pose any major difficulties, with the exception of the geographic pole ϕ = π 2 , where the value of the scale of linear distortion towards the parallel and the scale of area distortions are obtained in the form of indeterminacy, i.e. μ λ = p = 0 0 .Thus, the limit of the function that describes the scale at ϕ → π 2 was calculated.The scales and the distortions take the extreme values in the points of the coordinates: The relevant formulas were applied to calculate the values of distortions in the nodes of the regular grid and maps were developed that present the distribution of local distortions in the projected triangle -i.e. the representation of ⅛ of the surface of the sphere (Figures 8-12).
As a result, the following distortion distributions were obtained.The isolines of the scales of distortion for areas and the scales of length distortions towards the parallels are the same in both cases and have the shape of arcs of circles that are parallel to the images of parallels, while the isolines of length distortion scales towards meridians and extreme n scales are similar.The isolines of extreme scales m and angular distortions also have similar shapes.These maps allow us to notice the places where the smallest and largest distortions occur and the rate at which the projection distortions change.

Modifications of the projection
In order to obtain slightly different properties, the projection functions may be multiplied by the constant coefficient m 0. Such transformation will not cause changes to the shapes of meridians and parallels, i.e. ⅛ of the sphere will still be projected as the Reuleaux triangle, but it will result in a different distribution of projection distortions.
For example, in the projection discussed here, the length of the sides of the Reuleaux triangle equals where R is the radius of the Earth.They are representations of ¼ of the equator and two halves of meridians of the lengths of π 2 R. In order to obtain a projection that would maintain the length of these three lines, the projection functions should be multiplied by the coefficient m 0 = 3 π R.This results from simple calculations.The length of ½ of the meridian, i.e. the side of the triangle being ⅛ of the sphere equals: and the length of the corresponding arc of the circle -the image of this meridian (the side of the Reuleaux triangle): Hence, the ratio of these lengths is: Thus, in order to increase the length of the image of the meridian to the same value as in the original, r should be multiplied by the reverse of this quotient.
The projection may also be modified so that the surface areas of ⅛ of the sphere and its image are equal (this does not mean an equal area projection), then the projection functions should be multiplied by the coefficient , which results from the following calculations: The area of the Reuleaux triangle equals: where d is equal to the length of the side of the triangle, which, for this projection before modification, is equal to the length of half meridian d = π 2 R. At the same time, the area of ⅛ of the sphere: Thus, the proportion of the areas of the triangles in the original and in the image plane equals: Thus, the coefficient that should be applied is equal to the square root from the reverse of this quotient.

Comparison of the Leonardo da Vinci's projection with other projections with similar shapes of graticules
For the purpose of our comparison, we use projections that have similar shapes of the graticules, i.e. the Bonne projection, with the parallel ϕ 0 = 45° that is projected without distortions, oblique azimuthal projections: equal area, conformal, and equidistant towards the meridians, with the main point of the coordinates ϕ 0 = 45°, λ 0 = 0°.Graticules for ⅛ of the sphere were created for the projections mentioned above.The graticules are presented in Figure 13.
These projections were also compared in terms of distortions, where mean length distor-tions were measured in a grid of regularly placed points in ⅛ of the sphere, according to Airy's measure, based on formula (4), where m and n are extreme scales of length distortions.The results are presented in Table 1.In spite of the fact that lower distortions were obtained in the azimuthal projections and the Bonne projection, visually the graticule in the Leonardo da Vinci's projection performs much better.

Methodology for the analysis of A map of Imola
The analysis presented in this study was conducted based on a scanned map published on the Royal Collection Trust website (Royal Collection Online, n.d. a).The map's physical dimensions are 44.0 cm × 60.2 cm, while the scan is 2000 by 1483 pixels.As mentioned above, the analysis was performed only for the area contained within the city walls.
The cartometric analysis of A map of Imola was performed using commonly applied tools for analysing historical maps and MapAnalyst software.
MapAnalyst was developed at the Institute of Cartography in Zurich.The software was written in Java and it may be downloaded free of charge from the website www.mapanalyst.org(Jenny, 2020).It enables effective identification and management of control points on historical maps and modern reference maps.Apart from that, it calculates and displays grids of distortion and vectors of error on historical maps, as well as isolines of scale and local rotation angles.It offers a wide range of parameters to adjust the created graphics.The software also calculates the scale of historical maps and their rotation angle, provides statistical indicators and offers an interactive tool to analyse the local changes in the shift, scale and rotation angle of the map (Jenny, 2006).
The analysis of the cartometric properties in MapAnalyst software begins with uploading the historical and modern maps.The user identifies pairs of corresponding points and places them on both maps in an interactive way.Then, the appropriate method of transformation is selected.The next step is the selection of one of the three visualisation methods.The displacement vector is the simplest method.Each line starts at the previously identified point on the historical map and ends in the location where this point would be situated if the historical map was as accurate as the reference map.This end point is the result of the application of the transformation method between the sets of corresponding points.Long vectors identify large errors on the map, while short vectorssmall errors.Another visualisation method is a distortion grid, which displays local distortions and rotation angles on the historical map.This algorithm is based on the multi-quadratic interpolation method.Finally, the isoline method The interpretation of the results should take into consideration the fact that the applied methods could reveal inaccuracies of historical maps, which are caused by two factors: − the cartographer might have made mistakes at various stages of creating the map (e.g. during measurements and data integration or while drawing and reproducing the map), − the material on which the map was drawn is subject to deformations with the passing of time and conditions in which it is stored (Jenny, 2006).
In this research, the central part of A map of Imola was analysed using the contemporary OSM map as a reference material.70 pairs of corresponding points were selected on both maps (Figure 14).The analysis was conducted utilising Helmert transformation and visualisation methods, including the distortion grid and isolines of scale and rotation angles.

Analysis of the cartometric properties of A map of Imola
According to the obtained calculation results, the map's scale is 1:4,640, and the rotation angle in comparison to the current map is 13 degrees.There are some local changes in the scale and rotation angles on the map, which are presented in Figures 15-17.The differences in distance between individual points were also calculated.The positional accuracy obtained using root mean square error was 12.800 m.
Visualisation of a distortion grid reveals that the largest distortions occur in the south-western and south-eastern parts of the map (see Figure 15, left).Figure 16 presents the isolines of the local rotation angles of the map in reference to the OSM map.The interval of 1 degree and a 300 m radius of the circle were adopted.The map shows that the largest rotation angles of the analysed map compared to the reference map are 18 and 17 degrees.Finally, Figure 17 represents the changes in scale of the analysed map.The scale varies from approx.1:4,400 in the bottom left corner of the map to approx.1:4,900 in the right bottom corner.This means that 1 cm on the map corresponds to approx.44 m and 49 m respectively in these two locations.
The differences between A map of Imola and the modern map may be caused by numerous factors, including: − changes to the physical dimension of the map over the last 500 years, − the content of the map was supplemented by using cartographic materials that had been developed earlier, by other authors, − the use of less accurate measurement methods.

Conclusion
Leonardo da Vinci was the author of numerous cartographic studies, including maps and cartographic projection.A map of Imola is the earliest surviving town map in the orthogonal projection entirely created using detailed field measurements and principles of geometry.As such, it has a cartometric basis, which has been analysed in this article.Another analysed work by Leonardo da Vinci is his Mappa mundi.It presents the image of the whole globe on a plane, divided into octants represented by the Releaux triangle.The map was created in a pseudo-conical projection, which is neither a conformal nor an equal area projection.The central meridian of each octant is projected as a straight line, with its length maintained.The other meridians are projected as symmetrical curves in reference to the central meridian, with meridians of the longitudes λ = ± π 4 projected as arcs of circles.The parallels are projected as arcs of concentric circles.The advantage of this type of projection is the relatively low distortion that results from the fact that the sphere is divided into several sections, each represented separately.When the parts are appropriately arranged, they create a map that gives an impression of continuity of the represented seas and land masses.
from the Palazzo Comunale to the centres of street junctions (Royal Collection Online, n.d.a).The map is currently part of the Royal Collection in Windsor and its dimensions are 44 × 60.2 cm.

Figure 1 .
Figure 1.A map of Imola (Royal Collection Trust / © His Majesty King Charles III 2023)

Figure 2 .Figure 3 .
Figure 2. Mappa mundi in the Leonardo da Vinci's projection (northern hemisphere) (Royal Collection Trust / © His Majesty King Charles III 2023) arcs of circles of the radii r = R� π 2 -ϕ� for the specific latitude ϕ = ϕ.The image of each parallel starts and ends in the points located on the images of meridians that are the borders of the Reuleaux triangle.The images of meridians in the da Vinci's projection have the form of certain curves, provided that two meridians of the longitudes λ = ± π 4 are represented as the arcs of circles of the radius r = R and the centres in the points of the coordinates x 1 = -√3 4 Rπ, y 1 = R π 4 and x 2 = -√3 4 Rπ, y 2 = -R π 4 .The central meridian λ = 0 is projected to a section of a straight line of the length R π 2 .The ends of the section have the coordinates x = 0, y = 0 and x = -R π 2 , y = 0.

Figure 5 .
Figure 5.The graticule in the projection used by Leonardo da Vinci

Figure 7 .
Figure 7.The graticule of ¼ of the globe

Figure 13 .
Figure 13.Presentation of graticules in the following projections (left to right): da Vinci, Bonne, oblique azimuthal conformal, equal area, and equidistant projections

Figure 14 .
Figure 14.Sets of corresponding points on Leonardo da Vinci's map (Royal Collection Trust / © His Majesty King Charles III 2023) and on the modern OSM map

Figure 17 .
Figure 17.Isolines of scale (interval 1:100, the radius of the search circle 300 m) on Leonardo da Vinci's map (Royal Collection Trust / © His Majesty King Charles III 2023) and on the modern OSM map