# **Geographic Information Science (GIScience) and Geospatial Approaches for the Analysis of Historical Visual Sources and Cartographic Material**

Edited by Motti Zohar Printed Edition of the Special Issue Published in *ISPRS International Journal of Geo-Information*

www.mdpi.com/journal/ijgi

**Geographic Information Science (GIScience) and Geospatial Approaches for the Analysis of Historical Visual Sources and Cartographic Material**

## **Geographic Information Science (GIScience) and Geospatial Approaches for the Analysis of Historical Visual Sources and Cartographic Material**

Editor

**Motti Zohar**

MDPI ' Basel ' Beijing ' Wuhan ' Barcelona ' Belgrade ' Manchester ' Tokyo ' Cluj ' Tianjin

*Editor* Motti Zohar Geography and Environmental Studies University of Haifa Haifa Israel

*Editorial Office* MDPI St. Alban-Anlage 66 4052 Basel, Switzerland

This is a reprint of articles from the Special Issue published online in the open access journal *ISPRS International Journal of Geo-Information* (ISSN 2220-9964) (available at: www.mdpi.com/ journal/ijgi/special issues/historical visual sources).

For citation purposes, cite each article independently as indicated on the article page online and as indicated below:

LastName, A.A.; LastName, B.B.; LastName, C.C. Article Title. *Journal Name* **Year**, *Volume Number*, Page Range.

**ISBN 978-3-0365-4124-2 (Hbk) ISBN 978-3-0365-4123-5 (PDF)**

© 2022 by the authors. Articles in this book are Open Access and distributed under the Creative Commons Attribution (CC BY) license, which allows users to download, copy and build upon published articles, as long as the author and publisher are properly credited, which ensures maximum dissemination and a wider impact of our publications.

The book as a whole is distributed by MDPI under the terms and conditions of the Creative Commons license CC BY-NC-ND.

## **Contents**


## **About the Editor**

#### **Motti Zohar**

Motti Zohar is a geographer studying the natural and human environments using Geo-Information Science. In his studies, he integrates GIScience in conjunction with various scientific disciplines to resolve spatial and temporal scenarios.

## *Editorial* **GIScience and Historical Visual Sources: A Promising Look at Past Scenarios and Sceneries**

**Motti Zohar**

Department of Geography and Environmental Studies, University of Haifa, Haifa 3498838, Israel; motti.zohar@univ.haifa.ac.il; Tel.: +972-4-8240444

**Keywords:** historical geography; GIS; GIScience; HGIS; visual sources; spatial approaches; cartography

The discipline of historical geography evolved rapidly in the 20th century [1,2] and has been applied in many studies worldwide. The main difference between historical geography and historical studies is that the former investigates past scenarios from a geographical point of view [3]. For this task, it utilizes spatial approaches to interrogate textual and written sources as well as various cartographic material. The prominent themes of historical geography include landscape and cityscape reconstructions, rural and urban development, past regional demographic examinations, and historical trends associated with geographic locations. Over the years, the discipline has seen great popularity, but lately, it seems that its growth has somewhat slowed down in parallel to the developing discipline of environmental history [4,5]. Generally speaking, one can claim that historical geography is seeking new conceptualizations and theories in order to rise again [6,7].

The introduction of GIS (geographic information systems) and GIScience (geographical information science) have revolutionized the way scientists analyze spatial data. To date, GIS serves as a framework for a wide range of spatial applications [8]. Using GIS, one can perform cartographic analyses [9–11]; create deep maps [12]; leverage geographic information using geodesign approaches [13]; inspect narratives using map stories [14]; and conduct place-based studies using platial GIS [15]. GIS is also used to create new forms of virtual knowledge [16]; perform 2D or 3D landscapes reconstructions [17–22]; and resolve complex scenarios of past phenomena [23–25]. The GIScience sub-discipline defining historical analyses using GIS is referred to as HGIS (historical GIS) [26]. Similar to other GIScience sub-disciplines, HGIS is evolving and being integrated significantly with historical geography [27,28].

The study of historical scenarios benefits from written accounts but also from visual sources such as drawings, maps, early photographs, and air-photos. Until GIScience was introduced, much of the interrogation of these sources was qualitative in nature. Moreover, occasionally, it overlooked vast amounts of information and geographic features one might extract from these sources. GIScience can bridge this gap and offers a wide range of approaches to interpret and analyze historical visual sources that are of relevance to historical geography as well as other social sciences and humanities. In my opinion, the potential of using GIScience in historical geography studies is enormous and can be the "game changer" that historical geography is desperately looking for [29].

The aim of this Special Issue entitled "Geographic Information Science (GIScience) and Geospatial Approaches for the Analysis of Historical Visual Sources and Cartographic Material" is to address the exploitation of historical visual sources for resolving past scenarios and sceneries. With the theoretical and methodological contributions made in this Special Issue, it is possible to shed some light on the potential benefits of using GIScience in conjunction with such visual sources, leading to promising perspectives of the past.

The Special Issue includes eight published papers. The publications are equally divided into methodological and theoretical contributions. The first group of publications

**Citation:** Zohar, M. GIScience and Historical Visual Sources: A Promising Look at Past Scenarios and Sceneries. *ISPRS Int. J. Geo-Inf.* **2022**, *11*, 286. https://doi.org/10.3390/ ijgi11050286

Received: 23 April 2022 Accepted: 26 April 2022 Published: 28 April 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

contains tools, techniques, and approaches for a quantitative interpretation of various visual sources as well as deriving spatial insights of their creators. This corpus includes a visual graph-based framework for matching geographical areas through time [30]; assessing the impact of the cartographer's position and topographic accessibility on the accuracy of the extracted historical land use information [31]; measurement analysis of urban spatial layouts using the square grid method [32]; and coupling historical maps and LiDAR data to identify man-made landforms [33]. The second group includes analyses and examinations of past scenarios: GIScience and historical cartography for evaluating land use changes and resulting effects on carbon balance [34]; geospatial analysis of the non-surveyed (estimated) coastlines in Inoh's Map, 1821 [35]; agricultural land-use changes in the Judean region from the end of the Ottoman Empire to the end of the British Mandate [36]; and the historical Vltava River Valley–various historical sources within web mapping environment [37].

These studies emphasize the enormous amount of information one can extract from old visual sources using GIScience. Seemingly, such studies are being conducted for ages, but the road is still long, and there is much more to be discovered as the interrogation of the available visual sources is far from being complete. In a broader context, GIScience has a crucial part in fostering the historical geography discipline as well as developing other digital humanities domains by facilitating research, education, public outreach, and data dissemination. I hope that this SI and the included papers have contributed a little to achieve that goal.

I would like to acknowledge the authors of the SI (as well as those eventually not included in it) for their valuable efforts, the reviewers who contributed significant comments and remarks, the editorial office, and the *ISPRS International Journal of Geo-Information* for their great assistance and support.

**Funding:** This research received no external funding.

**Acknowledgments:** I wish to thank the Israel Science Foundation for their support and contribution in fostering the integration of GIScience with historical geography studies.

**Conflicts of Interest:** The author declares no conflict of interest.

#### **References**


## *Article* **Agricultural Land-Use Changes in the Judean Region from the End of the Ottoman Empire to the End of the British Mandate: A Spatial Analysis**

**Gad Schaffer**

Multidisciplinary Studies Department, Tel-Hai Academic College, Upper Galilee 1220800, Israel; schaffergad@telhai.ac.il; Tel.: +972-55-2236800

**Abstract:** Vines and olives are two important and widespread traditional agricultural crops that are also connected to the Judeo–Christian–Muslim tradition. The goal of the research was to demonstrate the importance of using cartographical sources to obtain a more accurate and complete view of the past. To this end, the aims were: (1) to reconstruct the former agricultural land-use in three periods, 1873–1874, 1917, and 1943–1945; (2) to analyze the different spatial physical factors that could explain the spatial distribution of traditional agricultural landscapes; (3) to identify the changes which took place between the three reconstructed timestamps. The research employed different cartographic sources and the implemented analyses were conducted using GIS tools and methods. The results show that, in the past, the distribution of vines and olive groves greatly depended on several physical geographic factors (climate, slopes, direction). Nonetheless, human factors such as political instability, cultural and religious beliefs contributed as well. Moreover, this research showed how GIS has advanced historical geography research. Lastly, the research demonstrated that obtaining the most complete view of the past can be achieved by a combination of sources together with the use of GIS tools and methods.

**Keywords:** land use/land cover (LULC); landscapes; historical maps; Geographic Information System (GIS); agriculture; vineyards; olive groves; Ein Karem; Bethlehem; Hebron

#### **1. Introduction**

The interdisciplinary field of Historical Geographic Information Systems (HGIS) was formed towards the end of the 1990s [1]. The novelty of HGIS lies in the use of Geographic Information System (GIS) software programs for the study of history, and later on for geography as well [2–4]. GIS software programs allow for the extraction, analysis, and synthesis of spatial components found in cartographic sources, a task that was previously almost impossible or very time-consuming [2]. The introduction of GIS to the field of History and Geography has increased the use of historical cartographical sources (maps, drawings, photographs, and aerial photos) for research [5,6].

In the last two decades, much research has been performed on the reconstruction of past landscapes and the examination of changes in land use/land cover (LULC) over time [7–10]. LULC research has increased awareness of LULC changes and their links to current environmental changes [11]. Moreover, landscapes comprise cultural meanings, values, beliefs, and ideologies, making each one of them unique [12,13]. Landscapes, or more precisely, cultural landscapes, were recognized by UNESCO in 1992 as important "combined works of nature and man" to be acknowledged and protected [14] (p. 18). In 2014, UNESCO declared the Battir cultural landscape as a World Heritage Site. Battir is a small Palestinian village in the West Bank, located south-west of Jerusalem, in the Judean region (Figure 1). The nomination of Battir as a World Heritage Site was due to its special farmed valleys composed of old stone terraces, which was the traditional way of irrigating

**Citation:** Schaffer, G. Agricultural Land-Use Changes in the Judean Region from the End of the Ottoman Empire to the End of the British Mandate: A Spatial Analysis. *ISPRS Int. J. Geo-Inf.* **2021**, *10*, 319. https:// doi.org/10.3390/ijgi10050319

Academic Editors: Wolfgang Kainz and Motti Zohar

Received: 23 February 2021 Accepted: 6 May 2021 Published: 8 May 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

agricultural land in this part of the world. The title that UNESCO awarded to the Battir Cultural Landscape was Land of Olives and Vines [15].

Following the agricultural revolution about 10,000 years ago, humans started to settle down where good physical conditions such as fertile soil, water, and moderate temperatures could increase the production of good and abundant crops, thus providing plenty of food [16]. Gradually, agricultural crops became important for the survival of the increased population [16]. Agriculture was also important since its products could be traded for other products or for money. The agricultural crops in each geographical area were different, especially because of the unique physical conditions of each region. According to various archaeological findings and written historical sources, the main crops in Palestine were wheat, barley, legumes, and grapes, as well as fruit trees, mainly olive, date and pomegranate [17]. Often, these crops were given additional cultural symbolism and religious importance. This was the case with the olives and vines [18].

In many parts of the Mediterranean region and in the eastern regions of the Middle East, vineyards and olive groves have been a part of the landscape for centuries [19–21]. Olive trees originated in the Mediterranean region and were domesticated around 6000–7000 years ago [22]. Vines (Vitis vinifera) are climbing plants that originated in different places around the world [17]. Viticulture as a source of wine has always been one of the most important and widespread agricultural land-uses due to its significant presence in the Judeo–Christian tradition [23,24]. Likewise, the olive tree and its products are deeply rooted in the Judeo– Christian–Muslim tradition [25]. In the Bible, vines are mentioned hundreds of times, reflecting their substantial role in agriculture. Furthermore, the names of many historical settlements in Palestine derived from the word vines [24]. For example, *Gat,* which in Hebrew means a winepress device, and *Ein Karem,* which in Hebrew means the spring of vineyards. While olives and vines were important crops in Palestine, the number of vines and olives fluctuated throughout history. Some causes of these fluctuations were climate-related, for example, years of drought destroyed vine plantations [26,27]. Another reason was that throughout history, various pests attacked and damaged crops, such as the grape phylloxera insect, which attacks vines [23]. Other reasons were anthropogenic, such as political instability generated by conflicts and wars in the area that caused the locals to migrate and abandon their crops [25]. Additionally, the fluctuations were connected to trends in chain supply and demand, which changed over the years. Changes in culture and religious beliefs also affected crops, such as the prohibition of wine production, to varying degrees, due to the application of Islamic laws in several periods [23,28]. For example, during the Muslim conquest (636–1099 AD), the Caliphate ordered the destruction of all wine-producing vineyards and only allowed the growth of a type of vineyard which produced table grapes [23]. Additionally, in the Ottoman period, the prohibition of growing vines for wine was mostly enforced until the middle of the nineteenth century [23].

At the beginning of the nineteenth century, the power of the Ottoman empire that ruled Palestine from 1516 began to diminish, and Palestine was left to fend for itself [29]. In that period, Palestine suffered from political instability, security problems, and a shrinking population [29,30]. With the continuing decline of the Ottoman empire towards the midnineteenth century, European interests in the region began to increase together with the beginning of Jewish immigration towards the end of the century. This was apparent from the capitulation treaties, which facilitated foreign investment, development, and exploration of Palestine [29,31]. The late nineteenth century is considered as the beginning of Palestine's modernization, which consisted of the construction of new buildings and infrastructures such as the Hejaz railroad, the start of foreign trade and a significant growth in population that enhanced further development [32]. During World War I, the Ottoman Empire collapsed, and Palestine came under the rule of the British Mandate for the next 31 years. Under the British Mandate, Palestine was classified as an agricultural land, as was the case for most British colonies [33]. During the last years of Ottoman rule, olive groves were present all across Palestine, but vineyards could only be found in two locations in the Galilee (around Kefr Birim, Sasa, and el Fish, and between Beit Jenn and Rameh), in a

*ISPRS Int. J. Geo-Inf.* **2021**, *10*, x FOR PEER REVIEW 3 of 16

single location in Samaria (around Mezrah esh Sherkiyeh), and in three locations in Judea (around Ein Karim, Bethlehem, and Hebron). locations in the Galilee (around Kefr Birim, Sasa, and el Fish, and between Beit Jenn and Rameh), in a single location in Samaria (around Mezrah esh Sherkiyeh), and in three locations in Judea (around Ein Karim, Bethlehem, and Hebron).

olive groves were present all across Palestine, but vineyards could only be found in two

**Figure 1.** The three research areas—Ein Karem, Bethlehem, and Hebron. Each research area is circled in red, which delineates a 5 km radius from the main settlement outwards. The watershed ridge line of the Judaean and Samarian mountains and the phytogeographical regions are also shown (the phytogeographical regions were made on the basis of the 1984 Floristic Regions of the Land of Israel map [34]). **Figure 1.** The three research areas—Ein Karem, Bethlehem, and Hebron. Each research area is circled in red, which delineates a 5 km radius from the main settlement outwards. The watershed ridge line of the Judaean and Samarian mountains and the phytogeographical regions are also shown (the phytogeographical regions were made on the basis of the 1984 Floristic Regions of the Land of Israel map [34]).

The goal of the research was to demonstrate the importance of using cartographical sources to obtain a more accurate and complete view of the past. The research spans 72 years that consist of over three timestamps, between 1873 and 1945. The reasons for choosing this period were twofold. First is the rapid change including substantial geo-political shifts. The second reason for ending the research period with 1945 was practically technical. Between 1948 and 1967, two of the research areas located in the West Bank were under Jordanian control. The research only found reprints of older maps with no new additions. Large-scale maps of other periods were not available, and small-scale maps did not represent the examined categories. Moreover, in the years since 1945, the categories found on the maps have changed, and in some maps the categories examined (i.e., vineyards and olive groves) became more general (i.e., agricultural area), which made it impossible to compare them with previous years. The goal of the research was to demonstrate the importance of using cartographical sources to obtain a more accurate and complete view of the past. The research spans 72 years that consist of over three timestamps, between 1873 and 1945. The reasons for choosing this period were twofold. First is the rapid change including substantial geopolitical shifts. The second reason for ending the research period with 1945 was practically technical. Between 1948 and 1967, two of the research areas located in the West Bank were under Jordanian control. The research only found reprints of older maps with no new additions. Large-scale maps of other periods were not available, and small-scale maps did not represent the examined categories. Moreover, in the years since 1945, the categories found on the maps have changed, and in some maps the categories examined (i.e., vineyards and olive groves) became more general (i.e., agricultural area), which made it impossible to compare them with previous years.

To achieve this goal, the research examined the agricultural landscapes found in the Judean region, focusing on vineyards and olive groves. More specifically, three objectives were set towards achieving the goal of this research: (1) to reconstruct the former agricultural land-use over three periods, 1873–1874, 1917, and 1943–1945; (2) to analyze the different spatial physical factors that could explain the spatial distribution of agricultural landscapes; and (3) to identify the changes which took place between the three reconstructed timestamps. To achieve this goal, the research examined the agricultural landscapes found in the Judean region, focusing on vineyards and olive groves. More specifically, three objectives were set towards achieving the goal of this research: (1) to reconstruct the former agricultural land-use over three periods, 1873–1874, 1917, and 1943–1945; (2) to analyze the different spatial physical factors that could explain the spatial distribution of agricultural landscapes; and (3) to identify the changes which took place between the three reconstructed timestamps.

#### **2. Materials and Methods**

#### *2.1. Research Area*

The geographic areas chosen for the research were three 5 km radius areas around three settlements in the Judean region—Ein Karem, Bethlehem, and Hebron (Figure 1). Each radius consisted of an area of 78.5 km<sup>2</sup> . In the preliminary stage of the research, the 1881 Palestine Exploration Fund (PEF) map, which was the first source to be used in this research, was quickly scanned in order to identify areas which had both vineyards and orchards in the same location. Due to the presence of more areas in the Judea region than in the Galilee and Samaria, it was decided to focus on the three areas there, as the available sample was broader and could thus paint a clearer picture. The choice of a 5 km radius around these settlements was made because the distance from the center of these settlements to the boundaries of the agricultural crops found on the PEF map was approximately 5 km.

#### *2.2. Research Materials*

To reconstruct the changes in land-use over time, three sets of British maps were used. The first was the PEF map at a scale of 1:63,360, published in 1881 [35] (Table 1). The PEF map is considered the most accurate and precise map of nineteenth-century Palestine [36–38]. The survey was carried out between 1871–1878 by two British army officers, Claude Conder and Horatio Kitchener [39]. Specifically, the areas studied in this research were surveyed between 1873–1874. The second set of maps were British military maps at a scale of 1:25,000, published in 1917 (Table 1); these maps were copies of the original 1881 PEF map, with updated information. The third set of maps used were British Mandatory maps at a scale of 1:20,000, published between 1943–1945 (Table 1). Another source used was a world digital elevation model (DEM) at a resolution of 1-ARC (Table 1). A physiographic map of Israel from 1995 and a 1984 Floristic Region map of Israel were also used (Table 1) [34,40]. A written source that was widely used in this research was the Survey of Western Palestine: Memoirs of the Topography, Orography, Hydrography, and Archaeology—Judea [41]. The memoirs were written along with the PEF survey and describe the surveyed areas in great detail.


**Table 1.** The cartographical sources used in the research. For each cartographical source, the following details were provided: the agency that created the map, the name of the source, the year the source was published, the scale of the source and, if the source was georeferenced in this research, and the total Root Mean Square Error (RMSE).

#### *2.3. Research Methods*

The georeferencing, digitization, analysis, and examination of the different datasets were all performed using ArcGIS (10.5.1) software (Figure 2). The first step of this research was to geo-reference the cartographical sources. The GIS software allows us to examine the accuracy of the maps by comparing them to more accurate present-day sources. During the georeferencing process, the GIS software also calculates the Root Mean Squared Error (RMSE) [52]. The higher the RMSE number, the less spatially accurate the cartographical source. The 1881 PEF map, as well as the 2014 World DEM image, were already georeferenced and available in GCS WGS 1984, a world coordinate system [36,53]. The 1917 maps and the 1943–1945 maps were georeferenced using the graticule intersections (the coordinate grid points on the maps). The 1917 and 1943–1945 maps had a low RMSE, meaning they were quite accurate and could be used (Table 1). Since the coordinates in these sources were in the old local geographic coordinate system, the Israeli Cassini Soldner, the sources were first georeferenced to Israel's new geographic coordinate system, the Israeli Transverse Mercator, and were then projected into the world coordinate system, GCS WGS 1984, using the Project tool in ArcGIS. The physiographic map of Israel and the Floristic Regions of the Land of Israel were georeferenced using ten recognizable shared control points (i.e., Haifa Bay, Rosh HaNikra, Cape Costigan in the Dead Sea, etc.) found between the maps to an already georeferenced National Geographic map.

The second step of the research was to determine the agricultural categorization of the land-use. In the PEF and 1917 maps, the only two categories of agricultural land-use were vineyards and orchards. In the 1943–1945 maps, agricultural land-use was divided into five crop categories—vineyard, olive grove, banana grove, citrus grove, and orchards. In this research, the four land-use categories chosen for digitization were vineyards, olive groves, orchards, and built-up areas. The research focused on traditional agricultural crops used in this region, particularly vineyards and olive groves. For a better understanding of the settlement's linkage with local agricultural land-use, the built-up areas were also digitized. In the third step of the research, the land-use categories were digitized in all three time periods—1873–1874, 1917, and 1943–1945. The digitization was conducted on a scale of 1:5000–1:15,000.

The fourth step of the research was to compare and analyze the newly created digitized layers with different spatial factors. To quantitatively identify the changes in land-use categories over time, the Summarize tool was used in ArcGIS. In addition, the Tabulate Area (Spatial Analyst) tool was used to examine whether the land-use area categories of 1873–1874 remained the same in 1943–1945.

The fifth and last step was to examine whether physical spatial factors could explain the distribution and changes in land-use categories, and accordingly, the digitized layer was compared with elevation, slope, and aspect to the sun, using the 2014 world DEM. To compare the land-use categories to the elevation, the 2014 DEM raster was transferred into a polygon using the Raster to Polygon (Conversion) tool. Following that, the Intersect tool was used to perform a calculation between the land-use digitized layer and the polygon DEM layer. Finally, the compared data were extracted with a focus on the average area using the Summarize tool. To examine whether a link exists between the land-use layers and slopes, the 2014 DEM was used to measure the gradient of each area using the Slope (Spatial Analyst) tool. The new created slope layer was compared to the land-use layers using the Zonal Statistics as Table tool, from which data were extracted. Lastly, to examine in which direction the land-use categories found were facing with respect to the sun, the 2014 DEM was used to create a layer of directions to the sun using the Aspect (Spatial Analyst) tool. The created polygon layer was composed of ten different directions. To make the results easier to comprehend, the Aspect layer was first transformed into a raster using the Polygon to Raster (Conversion) tool. Later, using the Reclassify (Spatial Analyst) tool, the classes of the original Aspect layer were reclassified into four major directions: north, from 0◦–45◦ and from 315◦–0◦ , east from 45◦–135◦ , south from 135◦–225◦ , and west from 225◦–315◦ . Once the reclassification was complete and a new raster was produced, it was

compared with the land-use layers using the Tabulate Area (Spatial Analyst) tool, from which the data were extracted to an Excel sheet. shared control points (i.e., Haifa Bay, Rosh HaNikra, Cape Costigan in the Dead Sea, etc.) found between the maps to an already georeferenced National Geographic map.

The georeferencing, digitization, analysis, and examination of the different datasets were all performed using ArcGIS (10.5.1) software (Figure 2). The first step of this research was to geo-reference the cartographical sources. The GIS software allows us to examine the accuracy of the maps by comparing them to more accurate present-day sources. During the georeferencing process, the GIS software also calculates the Root Mean Squared Error (RMSE) [52]. The higher the RMSE number, the less spatially accurate the cartographical source. The 1881 PEF map, as well as the 2014 World DEM image, were already georeferenced and available in GCS WGS 1984, a world coordinate system [36,53]. The 1917 maps and the 1943–1945 maps were georeferenced using the graticule intersections (the coordinate grid points on the maps). The 1917 and 1943–1945 maps had a low RMSE, meaning they were quite accurate and could be used (Table 1). Since the coordinates in these sources were in the old local geographic coordinate system, the Israeli Cassini Soldner, the sources were first georeferenced to Israel's new geographic coordinate system, the Israeli Transverse Mercator, and were then projected into the world coordinate system, GCS WGS 1984*,* using the Project tool in ArcGIS. The physiographic map of Israel and the Floristic Regions of the Land of Israel were georeferenced using ten recognizable

*ISPRS Int. J. Geo-Inf.* **2021**, *10*, x FOR PEER REVIEW 5 of 16

*2.3. Research Methods*

**Figure 2.** Scheme showing the steps conducted in this research.

#### **3. Results**

#### *3.1. Overall Changes in Land-Use Categories over Time*

The agricultural land-use areas (orchards, vineyards, and olive groves) amount to 12–23% of the total research areas in 1873–1874, 10–23.5% in 1917, and 27.5–41% in 1945 (Table 2).

As shown in Table 2, in 1873–1874, Bethlehem had the largest area of orchards (20.73% of the total research area), followed by Ein Karem and Hebron. In contrast, Hebron had the largest area of vineyards (15.5% of the total research area), with Ein Karem and Bethlehem far behind. In 1917, as well as between 1873–1874, a similar trend was found, with Bethlehem with the largest area of orchards and Hebron with the largest area of vineyards. Between 1873–1874 and 1917, the results show a decrease in all rations of agricultural land areas with the exception of an increased percentage of total area of vineyards in Hebron (+18%). In 1943–1945, unlike previous examined years, a new separate land-use category emerged-olive groves—which, in both the 1873–1874 and 1917 maps, was included in the orchard category. Interestingly, in 1943–1945, Hebron still had the largest percentage of total area consisting of vineyards (26.11%) with an addition of orchards (21.91%). Between 1917 and 1943–1945, there was a noticeable increase in all three research areas in the total

area of all agricultural lands. As depicted in Figure 3, the dominant land-use in Ein Karem and Bethlehem was orchards, while in Hebron it was vineyards. Moreover, in Hebron to a larger extent, and in Bethlehem to a smaller extent, the orchards in the 1873–1874 and 1917 maps do not expand much to the east. While in Hebron in the 1943–1945 maps, the agricultural land has expanded across the entire research area and to the east, in Bethlehem the expansion of agricultural land was mostly to the north, west and south, leaving the east almost empty of this type of land. research areas in the total area of all agricultural lands. As depicted in Figure 3, the dominant land-use in Ein Karem and Bethlehem was orchards, while in Hebron it was vineyards. Moreover, in Hebron to a larger extent, and in Bethlehem to a smaller extent, the orchards in the 1873–1874 and 1917 maps do not expand much to the east. While in Hebron in the 1943–1945 maps, the agricultural land has expanded across the entire research area and to the east, in Bethlehem the expansion of agricultural land was mostly to the north, west and south, leaving the east almost empty of this type of land.

largest percentage of total area consisting of vineyards (26.11%) with an addition of orchards (21.91%). Between 1917 and 1943–1945, there was a noticeable increase in all three

*ISPRS Int. J. Geo-Inf.* **2021**, *10*, x FOR PEER REVIEW 7 of 16

**Figure 3.** Changes in agricultural land-use in the research areas. The figure is divided into the three research areas, Ein Karem, Bethlehem and Hebron, and each one is divided into the three time periods, 1873–1874, 1917, and 1943–1945. Land uses are divided into four categories—builtup areas, olive groves, orchards, and vineyards. **Figure 3.** Changes in agricultural land-use in the research areas. The figure is divided into the three research areas, Ein Karem, Bethlehem and Hebron, and each one is divided into the three time periods, 1873–1874, 1917, and 1943–1945. Land uses are divided into four categories—built-up areas, olive groves, orchards, and vineyards.

**Table 2.** Agricultural land-use changes in the research areas (Ein Karem, Bethlehem, and Hebron) for each period examined, 1873–1874, 1917, and 1943–1945. The total area for each category is presented both in square kilometers and percentage of the total research area. Additionally, the dynamic percentage of land use changes in each category is presented for the periods between 1873–1874–1917 and 1917–1943–1945. **Table 2.** Agricultural land-use changes in the research areas (Ein Karem, Bethlehem, and Hebron) for each period examined, 1873–1874, 1917, and 1943–1945. The total area for each category is presented both in square kilometers and percentage of the total research area. Additionally, the dynamic percentage of land use changes in each category is presented for the periods between 1873–1874–1917 and 1917–1943–1945.


An analysis of the continuity of area types found that only a small part of agricultural land in 1873–1874 was still agricultural in 1943–1945 (Table 3). For example, in Hebron in 1943–1945, only 2.9 km<sup>2</sup> of vineyards were located in the same area that had vineyards in 1873–1874. In Ein Karem and Bethlehem, the numbers were minimal (Table 3). It is interesting to note that many olive-grove areas in 1943–1945 occupy the same areas as orchards in 1873–1874. For example, the total amount of olive groves in Bethlehem in 1943–1945 was 13.2 km<sup>2</sup> (Table 2) and more than half of this area, 7.7 km<sup>2</sup> , was found to be orchard areas in 1873–1874 (Table 3). This could mean that more than half of the previous orchards in Bethlehem were olive groves. The same might be true of Ein Karem and Hebron (Tables 2 and 3).

**Table 3.** The continuity of location conversion between the agricultural land-use categories in 1873–1874 with the land cover in 1943–1945. The table shows the amount of change in square kilometers.


#### *3.2. Physical Factors Explaining Spatial Distribution of Land-Use Categories*

#### 3.2.1. Elevation and Slope

Results of the correlation between land-cover categories, elevation, and slopes are presented in Table 4. With respect to elevation, in all the three research areas, vineyards were found on average in higher areas than orchards. For example, in Hebron, vineyards could be found at an elevation of 938–939 m while orchards were found at an elevation of 890–923 m. A similar trend can be seen in the other two research areas. With regard to olive groves, in Ein Karem, the average elevation of olive groves (624 m) in 1943–1945 resembled the average elevation of orchards (626 m), a difference of 2 m (Table 4).

The slope results show that the land-use categories could be found in areas between 26 and 35 degrees. In Hebron, the slope degree was higher than in Bethlehem and Ein Karem. Moreover, in Hebron, vineyards were found on higher slopes (31 degrees) than orchards (26–29 degrees) in 1873–1874 and 1917, unlike Ein Karem and Bethlehem. In 1943–1945, olive groves could be found across a large range of slopes (Table 4).

**Table 4.** The correlation between land-cover categories to elevation and slopes for each period examined.


#### 3.2.2. Aspect

The results of the *Aspect* examination are presented in Table 5. Most of the vineyards in Ein Karem and Bethlehem across all periods were facing eastwards, while in Hebron they mostly faced southwards. With regard to orchards, in Ein Karem and Hebron most of them faced southwards, while in Bethlehem they mostly faced eastwards. The olive-grove aspect strongly resembled that of orchards. Most olive groves in Ein Karem and Hebron were southwards-facing, like the orchards in those areas, and in Bethlehem most olive groves faced eastwards like the rest of the orchards in this area.

**Table 5.** The downhill direction in relation to the sun. Each research area, shown on the left side, is subdivided into one of the three periods (1873–1874, 1917, and 1943–1945) and then further subdivided into land-use categories (vineyards, orchards and olive groves). The upper columns preset the percentage out of the total area that shows the division into four downhill directions in relation to the sun—north, east, south, and west.


#### **4. Discussion**

*4.1. Vineyards—Land-Use Change*

According to the PEF map, in the period of the survey of Palestine between 1871–1878, there were only six large areas of vineyards. Two areas were found in the Galilee region, one in the Samaria region, and three in the Judean region. Indeed, in general terms these three regions have good conditions for growing vines [17,54–56]. These regions are located at a latitude of 31◦ , the average temperature is around 10–20 ◦C and usually does not drop below −4 ◦C nor exceed 30 ◦C. Additionally, the average annual rainfall is between 550 and 700 mm. Despite the good conditions in many places in Palestine, the Ottoman occupation and the Islamic ban on drinking alcohol led to a decline in wine production

and the vineyards that produced it [17,57]. Notwithstanding the ban on wine, it was possible to assume that vineyards would be found in cities where Jews and Christians Arabs lived, but that was not the case. In the Judea region at the end of the nineteenth century, Bethlehem had a Christian majority (the population was 5000 people of whom only 300 were Muslims) and Ein Karem was an all Christian village of 600 people, and yet the area size of the vineyards there was very small compared to the size of the vineyards in Hebron, which had a Muslim majority (the population of Hebron numbered 10,000 of whom 1000–12,000 were Jews) [41]. The reason vineyard areas were very extensive in Hebron, a predominantly Muslim settlement, was that these vineyards produced other nonalcoholic vine products [57]. Vines were grown for table grapes, raisins, and the production of a concentrated grape syrup [17]. With the slow decline of the Ottoman Empire in the late nineteenth century, there was a gradual increase in vine plantations for wine production. Initially, this was carried out mainly by Christian Arabs who, for religious reasons, were allowed to plant and produce wine for personal use and for ceremonies. Nonetheless, these plantations were small in size. The first extensive vineyards were planted by new Jewish settlers in Rishon LeZion in 1885, followed by other vineyards in Rosh Pina and Zikhron Ya'akov [23]. The new Jewish settlements brought vines for wine production from France. Nonetheless, this did not last long, since in 1890 vines were attacked by the grape phylloxera insect pest and many vines had to be uprooted [23,58]. During the First World War, other agricultural crops were more crucial to the war effort, and the vine industry declined further. It was only from 1936 on that vines started to be planted in increasing amounts [17]. The results of this research show that, during World War I (1917), there was a reduction of 9–25% in the total area of orchards and vineyards (Table 2).

This research illustrates that the largest area of vineyards in all three periods examined herein was in Hebron, as also described in the PEF memoirs: "the neighborhood of the Hebron hills is one of the principal vine-growing districts" [41] (p. 319). Hebron has better physical factors for growing vines than the other two examined areas. The annual rainfall average in Hebron is around 520–550 mm/year, its altitude (927 m) results in lower temperatures (yearly average temperature in winter is 7.5–10 ◦C and 22 ◦C in the summer). Interestingly, if we examine the soil type, it is the same in all three research areas—Terra Rossa and Rendzina [34]. Nonetheless, the PEF memoirs mention that Hebron's valleys "have good soil in them" [41] (p. 297). Ilan (1984) argues that, in the nineteenth century, most of the vineyards in the Hebron area were planted in good, deep, fertile Terra Rossa soil, which has good drainage [33]. Indeed, for vines, the most important aspect of soil is the draining of water [55]. Lastly, most vineyards in Hebron were found on south-facing slopes, a direction that is an advantage in the Northern Hemisphere, because it attracts the warmth of the afternoon sun, which improves ripeness and can protect against frost [55].

It is interesting to note that the PEF memoirs describe Hebron as "surrounded with fine vineyards" (p. 309) and maintain that the vineyards "extended over 6 miles", that is, 9.6 km. Although, as the 1873–1874 map shows, the vineyard areas were extensive, the study did not find vineyards which extended beyond 5.5 km from Hebron.

Another interesting point is that written sources tell us that, at the end of the nineteenth century, new vineyards were planted by Christian Arab and Jewish immigrants [23]. However, the second period examined in this research (1917) shows a strong decrease in vineyard areas in the research areas (except in Hebron). The reason for that might be that World War I had taken a toll on the vineyard industry.

If we take a look at the changes in vineyards throughout the years, the results show a marked increase in vineyards from 1917 to 1943–1945 in all three research areas, and especially in Hebron and Ein Karem. The first aid that the British Mandate gave to the agricultural sector was during its military regime (1917–1920), following the extensive damage caused to the agricultural sector in Palestine during World War I. A large amount of agricultural produce in Palestine was intended for export, and during the war there was a drastic decline in exports and many vineyards and orchard areas shrunk [17,57]. Under the direction of Colonel E.R. Sawer, the British decided on several actions, among them the renewal of traditional agricultural crops, such as vines and olives. They also examined the options for developing new agricultural crops, chiefly, but not only, citrus [33].

The increase in vineyards can especially be seen in the 1943–1945 period, even though the precipitation in these years was slightly below average, which could have potentially reduced the vineyard areas [26]. One explanation for the increase in vineyard areas was that, during the first years of the Mandate in Palestine, the British encouraged the promotion of native crops such as vines [33]. A second reason for the increase in vineyards was the introduction of wine varieties by Jewish settlers and the British [23]. Lastly, the British decided to promote the development of additional water sources for agriculture [33]. Under the civil regime (1920–1948), the British promoted training services to encourage the cultivation of olives and vines, and provided irrigation services [33]. This could explain the planting of vineyards in areas that, before 1943–1945, would have been less hospitable for these plantations. For example, in Hebron in 1943–1945, vineyards were found in the south-east side that were previously uncultivated (Figure 3). Another interesting point to note is that while in the 1873–1874 and 1917 periods vineyards were located in one or two distinct areas, in 1943–1945, vineyards are scattered all over, in large and small areas. This is especially visible in the Ein Karem area (Figure 3). One explanation might be that with the encouragement of agriculture, guidance, and help from the British, and the realization that a larger market for vine products existed, many small settlements decided to plant vineyards as one of their many crops.

While in Palestine until 1917, vineyards were mostly grown for products other than wine, in the Mediterranean region, wine was their main product. Until the middle of the nineteenth century, when French vineyards were attacked by the phylloxera insect, global wine production was dominated by France. From 1880 onwards, France started to import wine from Spain, Italy, and its North African colonies which all saw an increase in the plantation of vines [59]. Moreover, Greece, Spain and Turkey were the main producers of raisins, which were used mostly for wine production as well as in other industries [60]. In 1890, France began to recover and reproduce wine, dictate tariffs and regulations, which continued into the twentieth century, and as a result wine and raisin production in other countries dropped, which in turn decreased the cultivation of vines [59,61,62]. While these changes occurred in the Mediterranean and European regions, in Palestine vineyard growth was not heavily affected by these events and was mainly a result of an increase in domestic consumption. The export of wine from Palestine only began in 1899, and in that year only 366 kilos were exported [63]. In the beginning of the twentieth century, the numbers started to soar in 1936, during the British Mandate, when the export of wine really rose [64].

#### *4.2. Olive Groves—Land-Use Change*

Olive groves only appear as a separate crop on the 1943–1945 maps. In the 1873–1874 and 1917 maps, the only categories are vineyards and orchards. Two sources can help us fill the gap regarding olive groves in previous periods. The first source is the PEF map memoirs, which describe the surveyed land in great detail. The memoirs give us hints as to where olive trees could be found. For example, in Hebron it is stated that "the hill above Hebron is terraced with stone walls and olive plantations" [41] (p. 307) (Figure 4A). Moreover, "vineyards of Hebron extend over above 6 miles, and olives are also grown there" [41] (p. 297). The description of Bethlehem mentions that "the valleys on the north and south are deep, the sides carefully terraced, vines and olives, figs and either trees are grown along the slopes" [41] (p. 28) (Figure 4B). These descriptions of the PEF memoirs, indicate that other crops existed apart from olive groves, "pomegranates, figs, quinces, and apricots" could be found around Hebron [41] (p. 308). While these sources clearly show that orchards on the maps also included olive groves, they do not specify the amount or size of these areas. Another way to tackle this issue is by examining where olive groves appear on the 1943–1945 maps and whether these exact areas were once orchard areas. There is a strong probability that olive trees that were planted and found in an area in the past could

still be found in that same area after many years. Planting fruit trees is usually done with the thought of long-term profit and especially with regard to olive trees. The results of this research allowed us to estimate the percentages of olive groves (Tables 2 and 3). The assumption according to the results is that Bethlehem had the largest olive groves in the past and still has the largest area in 1943–1945. in the past could still be found in that same area after many years. Planting fruit trees is usually done with the thought of long-term profit and especially with regard to olive trees. The results of this research allowed us to estimate the percentages of olive groves (Tables 2 and 3). The assumption according to the results is that Bethlehem had the largest olive groves in the past and still has the largest area in 1943–1945.

indicate that other crops existed apart from olive groves, "pomegranates, figs, quinces, and apricots" could be found around Hebron [41] (p. 308). While these sources clearly show that orchards on the maps also included olive groves, they do not specify the amount or size of these areas. Another way to tackle this issue is by examining where olive groves appear on the 1943–1945 maps and whether these exact areas were once orchard areas. There is a strong probability that olive trees that were planted and found in an area

*ISPRS Int. J. Geo-Inf.* **2021**, *10*, x FOR PEER REVIEW 12 of 16

**Figure 4.** Olive groves surroundings the city of Hebron in a photo from 1890–1900 (**A**) [65]. Olive trees in the front of the photo showing the city of Bethlehem from 1890–1900 (**B**) [66]. **Figure 4.** Olive groves surroundings the city of Hebron in a photo from 1890–1900 (**A**) [65]. Olive trees in the front of the photo showing the city of Bethlehem from 1890–1900 (**B**) [66].

Among the three research areas, Bethlehem and Ein Karem had the best physical conditions for olive groves in respect to elevation, precipitation, and soil [67,68]. For example, both Bethlehem and Ein Karem are located in the optimal elevation for olive groves (400–800 m), unlike Hebron, which exceeds the optimal elevation (927 m) [68]. Indeed the PEF memoirs also argue that in Hebron "the olive does not flourish well in any part within this Sheet, but the villages in the low hills have a few" [41] (p. 319). While both Ein Karem and Bethlehem had optimal physical conditions for growing olive groves, Bethlehem also had other reasons to have more olive groves. In 1881, Bethlehem was a large settlement with a population of 5000, while Ein Karem was a very small village with a population of only 600 [41]. This could partly explain the smaller amount of olive groves in Ein Karem. However, the stronger explanation is that, from the seventeenth century onwards, Bethlehem's olive trees, unlike other parts in Palestine, were mainly used to produce handcrafted religious souvenirs which would be later sold in Jerusalem or shipped abroad [69]. The Bethlehem orphanage, which was established in 1864, trained the orphanage children to carve in olive wood [69]. This is another indication that a large number of the orchards appearing on the 1873–1974 and 1917 maps were olive groves. Moreover, during the British Mandate, Bethlehem olive-grove areas expanded as olive wood craft continued to flourish in the city [70]. This goes hand in hand with the research results that show that in Bethlehem the size of the olive groves was the largest of all three research areas. Among the three research areas, Bethlehem and Ein Karem had the best physical conditions for olive groves in respect to elevation, precipitation, and soil [67,68]. For example, both Bethlehem and Ein Karem are located in the optimal elevation for olive groves (400–800 m), unlike Hebron, which exceeds the optimal elevation (927 m) [68]. Indeed the PEF memoirs also argue that in Hebron "the olive does not flourish well in any part within this Sheet, but the villages in the low hills have a few" [41] (p. 319). While both Ein Karem and Bethlehem had optimal physical conditions for growing olive groves, Bethlehem also had other reasons to have more olive groves. In 1881, Bethlehem was a large settlement with a population of 5000, while Ein Karem was a very small village with a population of only 600 [41]. This could partly explain the smaller amount of olive groves in Ein Karem. However, the stronger explanation is that, from the seventeenth century onwards, Bethlehem's olive trees, unlike other parts in Palestine, were mainly used to produce handcrafted religious souvenirs which would be later sold in Jerusalem or shipped abroad [69]. The Bethlehem orphanage, which was established in 1864, trained the orphanage children to carve in olive wood [69]. This is another indication that a large number of the orchards appearing on the 1873–1974 and 1917 maps were olive groves. Moreover, during the British Mandate, Bethlehem olive-grove areas expanded as olive wood craft continued to flourish in the city [70]. This goes hand in hand with the research results that show that in Bethlehem the size of the olive groves was the largest of all three research areas.

There are two other notably interesting points. The first, has to do with the lack of olive groves east of Bethlehem found across all three research periods. One very logical reason is the fact that the Mediterranean region ends to the east of Bethlehem, and the Irano-Turanian region begins (Figure 1). The Irano-Turanian region does not possess the ideal physical characteristics that both olive trees and vineyards need. Moreover, the east side of the valleys deepen rapidly which make it hard to grow agricultural crops [41]. The second interesting point is that, in all three research areas across all the periods examined, the orchards and olive groves were closer to built-up areas, while the vineyards were located further away from them (Figure 3). The results of the research failed to find a direct There are two other notably interesting points. The first, has to do with the lack of olive groves east of Bethlehem found across all three research periods. One very logical reason is the fact that the Mediterranean region ends to the east of Bethlehem, and the Irano-Turanian region begins (Figure 1). The Irano-Turanian region does not possess the ideal physical characteristics that both olive trees and vineyards need. Moreover, the east side of the valleys deepen rapidly which make it hard to grow agricultural crops [41]. The second interesting point is that, in all three research areas across all the periods examined, the orchards and olive groves were closer to built-up areas, while the vineyards were located further away from them (Figure 3). The results of the research failed to find a direct link between the distance of vineyards and orchards from built-up areas to the different examined physical factors.

The increase in olive grove areas between 1873 and 1945 found in this research was not a unique phenomenon, but the reasons for it were different. From the beginning of the nineteenth century, all across the Mediterranean region, a movement of Europeans began to encourage the plantation of olive trees in the coastal areas of Northern Africa and the

Levant. This was boosted by the increase in demand for olive oil in France and Italy [71]. In Spain, for example, the increase in olive groves was mostly due to the relative low capital investment needed by the farmers to grow olive trees and the fact that the trees could produces a variety of staple products (olives as food, olive oil, wood for heating, etc.) [72]. In Palestine, olive trees were mainly used in the production of soap with only a small part used to produce edible olive oil [71,73]. In 1885, the export of olive oil reached 916 liters and had peaked in 1890, with 2639 liters of oil but then reduced to a minimum of 32 liters in 1904 [63]. From the end of the nineteenth and the beginning of the twentieth centuries, the increase in olive grove plantations in Palestine was due to the increase in demand for soap—an important source of employment, mostly among the Arab population [71]. In 1885, 0.43 tons of soap were exported, with this export peaking at 4.4 and 5.2 tons in 1895 and 1899, respectively [63]. Moreover, as was previously mentioned, in Bethlehem, olive trees were used to produce souvenirs. From the end of the nineteenth century, with the opening of Palestine to foreign European nations, the tourism industry also increased, as well as the demand for religious handicraft [69,70].

#### **5. Conclusions**

The goal of the research was to demonstrate the importance of using cartographical sources to obtain a more accurate and complete view of the past. The research examined agricultural land-use changes in three research areas in the Judean region, from the end of the Ottoman Empire to the end of the British Mandate. The research employed different cartographic sources using GIS tools and methods.

The results of this research showed that, in the past, the distribution of vines and olive groves in the Judean region were greatly dependent on physical geographic factors (climate, slopes, direction). The research has shown that Hebron had the optimal physical environmental conditions for vines, and Bethlehem and Ein Karem had the optimal physical conditions for olive trees. Nonetheless, the research also demonstrated that human-related factors also had a large influence on the area size of these land-use categories. For example, although the main olive tree product at that time was designated for the soap industry, in Bethlehem, the hand-crafted religious souvenirs that were made of the wood of olive trees necessitated the growing of olive trees in this area. Despite the ban on wine production during the Ottoman period, Hebron was one of the centers for vines that were grown for products other than wine. Additionally, at the start of the twentieth century, the cultivation of vineyards in Palestine was mainly a result of an increase in domestic consumption, not local-international trade. The research has also demonstrated that, with the passage of time, agricultural land-use area increased and was found in less optimal physical environmental areas, such as the agricultural areas east of Hebron from 1917 onwards. This was brought about by political changes in Palestine, which led to developments and investments in the agricultural sector, mostly during the British Mandate.

While not an aim in and of itself—since this research is founded on historical cartographical sources—it was also critical to verify the level of accuracy and completeness of these sources. One method of examining completeness is by comparing the cartographical source to other cartographical sources from the same place and period [38]. However, we rarely have two different historical sources from the same place and time. A second method is to use a combination of historical sources such as diaries, memoirs, route notes, and travel maps from that same period, or different cartographical sources from the same place but from other periods and examine whether the features appearing on them make sense. This research used a combination of sources to strengthen the level of completeness of the 1873–1874 PEF map. The PEF memoirs, which describe the different regions and provide details on the area's demography, geography, and landscape, were used to check the accuracy of the map. These descriptions helped to complete the overall picture illustrated by the map. Furthermore, two scenery photographs of two of the research areas from 1890–1900 were added to this research to add further reference, albeit general, to the data on the map. Lastly, this research also examined other maps from later periods, 1917

and 1943–1945, to further strengthen our understanding of the past (1873–1874) and of long-term land-use changes.

The introduction of GIS into historical-geographic research has shed light both on the importance of cartographical material and on the large amount of information found in these sources. Lastly, this research has proved that one good way of painting a clearer and broader picture of the past is through a combination of different sources.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Data available on request.

**Acknowledgments:** This research could not have been done without the kind help and the open and free accessibility to use cartographic materials from the following institutions: The Tel-Hai Collage Library Cartographic Collection and the Israel National Library—Eran Laor Cartographic Collection. Thank you.

**Conflicts of Interest:** The author declares no conflict of interest.

#### **References**


## *Article* **Geospatial Analysis of the Non-Surveyed (Estimated) Coastlines in Inoh's Map, 1821**

**Yuki Iwai 1,2,\* and Yuji Murayama <sup>3</sup>**


**Abstract:** The history of modern maps in Japan began with Inoh's map that was made by surveying the whole of Japan on foot 200 years ago. Inoh's team investigated coastlines, major roads, and geographical features such as rivers, lakes, temples, forts, village names, etc. The survey was successively conducted ten times from 1800 to 1816. Inoh's map is known as the first scientific map in Japan using a systematic method. However, the actual survey was conducted only for 75% of the coastlines in Japan and the remaining 25% was drawn by Inoh's estimation (observation). This study investigated how the non-surveyed (estimated) coastlines were distributed in the map and why the actual survey was not conducted in these non-surveyed coastlines. Using GIS, we overlaid the geometrically corrected Inoh's map (Digital Inoh's Map Professional Edition) with the current map published by the Geospatial Information Authority (GSI) of Japan for examining the spatial difference. We found that the non-surveyed coastlines were in places where the practice of actual surveying was topographically difficult because of the limited surveying technology of those days. The analytical result shows that 38.6% of the non-surveyed coastlines were cliffs, 25.7% were rocky beaches, and 6.2% were wetlands and tidal lands (including rice fields and tidal flats).

**Keywords:** Inoh's map; historical maps; coastlines; terrain; land use; historical GIS

#### **1. Introduction**

Inoh's map is known as Japan's first scientific map with high accuracy [1]. In the later period of the Edo era (1603–1868), mapping techniques had developed incredibly, and the Shogunate (government) had started editing the Japanese archipelago's precise map to protect national land from foreign invaders [2]. At first, the national map had been compiled by connecting local maps of different regions based on different surveying standards. Finally, in 1821, the national map of Japan, titled "Dai Nihon Enkai Yochi Zenzu" (commonly known as "Inoh's map"), was completed by surveyor Inoh Tadataka. He surveyed the coastlines in the whole of Japan with a uniform standard unlike previous national maps of Japan.

Inoh's map consists of several different scales to represent the whole of Japan (3 smallscale maps about 1:432,000, 8 medium-scale maps about 1:216,000, and 214 large-scale maps about 1:36,000). To create a precise national map, Inoh Tadataka carried out 10 survey trips between 1800 and 1816 (Figure 1). Hokkaido was surveyed during the first expedition launched in 1800. Surveying in western Japan began in 1805. From the fifth survey in the Kinki and Chugoku regions, the project had gained enough financial support and direct management by the Shogunate government [3]. After that, the size of Inoh's team had been upgraded drastically with rich equipment. The combination of the open traversing and forward intersection functioned efficiently for accurate surveying at that time [4]. Using the latitude astronomical observation, Inoh Tadataka attempted to correct the map errors caused by open traversing [5]. Furthermore, he added the measurement of coastlines

**Citation:** Iwai, Y.; Murayama, Y. Geospatial Analysis of the Non-Surveyed (Estimated) Coastlines in Inoh's Map, 1821. *ISPRS Int. J. Geo-Inf.* **2021**, *10*, 580. https:// doi.org/10.3390/ijgi10090580

Academic Editors: Motti Zohar and Wolfgang Kainz

Received: 30 June 2021 Accepted: 24 August 2021 Published: 27 August 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

from the sea by boat to overcome topographical obstacles (in this paper, this method is considered estimation). coastlines from the sea by boat to overcome topographical obstacles (in this paper, this method is considered estimation). Many studies have demonstrated that Inoh's map is very accurate [6–8]. However,

[4]. Using the latitude astronomical observation, Inoh Tadataka attempted to correct the map errors caused by open traversing [5]. Furthermore, he added the measurement of

*ISPRS Int. J. Geo-Inf.* **2021**, *10*, x FOR PEER REVIEW 2 of 13

Many studies have demonstrated that Inoh's map is very accurate [6–8]. However, the coastlines surveyed by Inoh accounted for only 75% of the total in the country (Figure 2). The remaining 25% were drawn by his estimation (including indirect surveys from the sea and visual observation) [9], which is defined as the non-surveyed coastlines in this study. For the scientific discussion on the accuracy of Inoh's map, it is important to specify the places that were not surveyed, identify their geographical characteristics, and find out why these coastlines were not directly surveyed. A survey at Inoh's time might be complex, especially for the rias (jagged coastlines) in northern Japan, where steep terrains were dominant in the coastal region. It is also essential to examine the accuracy of the observed non-surveyed coastlines to discuss the limitations of surveying in that period. Based on the above background, this study investigates the geographical characteristics of the non-surveyed coastlines by analyzing their distributions, physical conditions, and differences between the estimated and actual coastlines at that time. the coastlines surveyed by Inoh accounted for only 75% of the total in the country (Figure 2). The remaining 25% were drawn by his estimation (including indirect surveys from the sea and visual observation) [9], which is defined as the non-surveyed coastlines in this study. For the scientific discussion on the accuracy of Inoh's map, it is important to specify the places that were not surveyed, identify their geographical characteristics, and find out why these coastlines were not directly surveyed. A survey at Inoh's time might be complex, especially for the rias (jagged coastlines) in northern Japan, where steep terrains were dominant in the coastal region. It is also essential to examine the accuracy of the observed non-surveyed coastlines to discuss the limitations of surveying in that period. Based on the above background, this study investigates the geographical characteristics of the non-surveyed coastlines by analyzing their distributions, physical conditions, and differences between the estimated and actual coastlines at that time.

**Figure 1.** Roads and coastlines in Inoh's map. The date of the survey is based on [10]. Mamiya Rinzo learned surveying from Inoh Tadataka, and he was in charge of surveying Hokkaido (Ezo) [11]. Source: Digital Inoh's Map. **Figure 1.** Roads and coastlines in Inoh's map. The date of the survey is based on [10]. Mamiya Rinzo learned surveying from Inoh Tadataka, and he was in charge of surveying Hokkaido (Ezo) [11]. Source: Digital Inoh's Map.

GIS-based analysis on historical maps has been conducted for restoring landscapes and geographical features to their state at the map's creation. Previous studies have examined historical maps by employing residual components for georeferencing (e.g., X, Y, etc.) [12–17]. For instance, Hirai [18] introduced a sophisticated method to analyze the distortion based on residual information from georeferenced historical maps. For the

**Figure 2.** Coastlines in Inoh's large-scale maps (1821). Only major islands of the Japanese archipelago (Hokkaido, Honshu, Shikoku, and Kyushu) are displayed. The image is superimposed on the georeferenced Inoh's large-scale map. Source: Digital Inoh's Map. **Figure 2.** Coastlines in Inoh's large-scale maps (1821). Only major islands of the Japanese archipelago (Hokkaido, Honshu, Shikoku, and Kyushu) are displayed. The image is superimposed on the georeferenced Inoh's large-scale map. Source: Digital Inoh's Map.

location errors of geographical features in historical maps, many studies have demonstrated the usefulness of overlay with other accurate maps [19–21]. Iwai and Murayama [22] discussed the distortion factors in Inoh's Tokyo map regarding the projection aspects and the surrounding geographical environment. Further, there were extensive studies on the changes in the landscape between the past and present [23–30]. The spatial transformation of the landscape and regional structure over a long period was the main target of these studies. However, their use of GIS and spatial analysis was not much concerned with investigating landscape features concerning how the historical maps were made (or where the surveying was conducted). **2. Materials and Methods** Unfortunately, to the best of our knowledge, accurate historical maps showing nonsurveyed coasts do not exist today. So that, we cannot obtain the actual (correct) shape of GIS-based analysis on historical maps has been conducted for restoring landscapes and geographical features to their state at the map's creation. Previous studies have examined historical maps by employing residual components for georeferencing (e.g., X, Y, etc.) [12–17]. For instance, Hirai [18] introduced a sophisticated method to analyze the distortion based on residual information from georeferenced historical maps. For the location errors of geographical features in historical maps, many studies have demonstrated the usefulness of overlay with other accurate maps [19–21]. Iwai and Murayama [22] discussed the distortion factors in Inoh's Tokyo map regarding the projection aspects and the surrounding geographical environment. Further, there were extensive studies on the changes in the landscape between the past and present [23–30]. The spatial transformation of the landscape and regional structure over a long period was the main target of these studies. However, their use of GIS and spatial analysis was not much concerned with investigating landscape features concerning how the historical maps were made (or where the surveying was conducted).

#### coastlines have retained the same condition until now. Hence, we use the current map **2. Materials and Methods**

published by the Geospatial Information Authority (GSI) of Japan to elucidate the geographical features of those days (According to the GSI website, the GSI map defines the Unfortunately, to the best of our knowledge, accurate historical maps showing nonsurveyed coasts do not exist today. So that, we cannot obtain the actual (correct) shape of the

the coastlines in the early 19th century. However, it is reasonable to assume that most

analysis.

other factors are discussed later.

coastlines in the early 19th century. However, it is reasonable to assume that most coastlines have retained the same condition until now. Hence, we use the current map published by the Geospatial Information Authority (GSI) of Japan to elucidate the geographical features of those days (According to the GSI website, the GSI map defines the coastline as the boundary between the land and the sea when the sea level reaches its highest level. This coastline definition is the same as Inoh's map [3]). In addition, we exclude the coastlines changed by anthropogenic activities (terrain modification) in this analysis. Inoh's Map"), which is a georeferenced Inoh's large-scale map (scale 1: 36,000) [31]. In the Digital Inoh's Map, georeferencing was conducted to collate the curved sections and branching points of the existing roads and coasts on the GSI map. Control points were set up at shrines and temples, whose locations have been fixed until the present. At least 50 control points (at most 200 control points) were acquired for a single Inoh's large-scale map [31]. Figure 4 shows the non-surveyed coastlines in Inoh's map and the actual coastlines of those days (substituted by GSI map). The shape of each non-surveyed coastline is

By overlaying Inoh's map with the GSI map, we can compare the non-surveyed (estimated) coastlines with the actual coastlines. The Inoh's map used in this study is the Digital Inoh's Map, Professional Edition (from now on referred to simply as "Digital

coastline as the boundary between the land and the sea when the sea level reaches its highest level. This coastline definition is the same as Inoh's map [3]). In addition, we exclude the coastlines changed by anthropogenic activities (terrain modification) in this

It is noteworthy that the remarkable transformation in coastlines in Japan started after the Meiji period (1868–1912) due to rapid modernization, such as constructing ports, reclaiming wetlands, etc. In this study, such changed coastlines are designated as "unknown" because we cannot visualize their exact shape in the early 19th century. We focus on the non-surveyed (estimated) coastlines of "cliff," "rocky beaches," and "wetlands" (including rice fields and tidal flats) that were supposed to be unchanged during the last 200 years (Figure 3) (According to the GSI website, a cliff is defined as a steep slope at least 3 m high and 75 m long. The rocky beach is a beach made up of rocks that are scattered around. Wetlands are a minimum of 75 m by 75 m or 50 m by 125 m. A tidal flat is defined as an area of sand, mud, etc., at least 50 m by 50 m). The effects of erosion and

It is noteworthy that the remarkable transformation in coastlines in Japan started after the Meiji period (1868–1912) due to rapid modernization, such as constructing ports, reclaiming wetlands, etc. In this study, such changed coastlines are designated as "unknown" because we cannot visualize their exact shape in the early 19th century. We focus on the non-surveyed (estimated) coastlines of "cliff," "rocky beaches," and "wetlands" (including rice fields and tidal flats) that were supposed to be unchanged during the last 200 years (Figure 3) (According to the GSI website, a cliff is defined as a steep slope at least 3 m high and 75 m long. The rocky beach is a beach made up of rocks that are scattered around. Wetlands are a minimum of 75 m by 75 m or 50 m by 125 m. A tidal flat is defined as an area of sand, mud, etc., at least 50 m by 50 m). The effects of erosion and other factors are discussed later. different from the shape of the actual coastline. This is because the non-surveyed coastline was drawn from the tip of the surveyed coastline by estimation (observation). Therefore, it is considered that the difference is larger further away from the surveyed coastline. It is preferable to use the topographic map of the same era, but there are no national maps showing land use or terrain in the late Edo period. Therefore, the GSI map was used in this study to supplement national scale topographic and land use information. We constructed the non-surveyed sections using the coastline data of the GSI map and combined the topographic and land use information with Inoh's map. The criterion for determining the topography and land use of the non-surveyed coast was based on the nearest adjacent map symbol in the GSI map (Figure 5).

*ISPRS Int. J. Geo-Inf.* **2021**, *10*, x FOR PEER REVIEW 4 of 13

**Figure 3.** Topographical/land conditions: (**a**) cliffs, (**b**) rocky beaches, and (**c**) wetlands and mudflats. Source of photo: Google Earth Pro. **Figure 3.** Topographical/land conditions: (**a**) cliffs, (**b**) rocky beaches, and (**c**) wetlands and mudflats. Source of photo: Google Earth Pro.

By overlaying Inoh's map with the GSI map, we can compare the non-surveyed (estimated) coastlines with the actual coastlines. The Inoh's map used in this study is the Digital Inoh's Map, Professional Edition (from now on referred to simply as "Digital Inoh's Map"), which is a georeferenced Inoh's large-scale map (scale 1: 36,000) [31]. In the Digital Inoh's Map, georeferencing was conducted to collate the curved sections and branching points of the existing roads and coasts on the GSI map. Control points were set up at shrines and temples, whose locations have been fixed until the present. At least 50 control points (at most 200 control points) were acquired for a single Inoh's large-scale map [31]. Figure 4 shows the non-surveyed coastlines in Inoh's map and the actual coastlines of those days (substituted by GSI map). The shape of each non-surveyed coastline is different from the shape of the actual coastline. This is because the non-surveyed coastline was drawn from the tip of the surveyed coastline by estimation (observation). Therefore, it is considered that the difference is larger further away from the surveyed coastline.

**Figure 4.** Example of a non-surveyed (estimated) coastline. The crimson line indicates the survey line corresponding to the road. Source: Digital Inoh's Map and GSI map. **Figure 4.** Example of a non-surveyed (estimated) coastline. The crimson line indicates the survey line corresponding to the road. Source: Digital Inoh's Map and GSI map.

It is preferable to use the topographic map of the same era, but there are no national maps showing land use or terrain in the late Edo period. Therefore, the GSI map was used in this study to supplement national scale topographic and land use information. We constructed the non-surveyed sections using the coastline data of the GSI map and combined the topographic and land use information with Inoh's map. The criterion for determining the topography and land use of the non-surveyed coast was based on the nearest adjacent map symbol in the GSI map (Figure 5). **Figure 4.** Example of a non-surveyed (estimated) coastline. The crimson line indicates the survey line corresponding to the road. Source: Digital Inoh's Map and GSI map.

seaward and landward of the current coastline. The surveyed coastline data include a small number of coastlines based on the shape of the georeferenced Inoh's large-scale map. In this study, the gaps generated from the surveyed coastline were not considered **Figure 5.** Examples of topographic and land use classification on the non-surveyed (estimated) coastlines: (**a**) cliffs, rocky beaches, and unknown, and (**b**) wetlands and mudflats. Non-surveyed coastlines were replaced by the current coastlines contained in the GSI map. Source: Digital Inoh's Map and GSI map. **Figure 5.** Examples of topographic and land use classification on the non-surveyed (estimated) coastlines: (**a**) cliffs, rocky beaches, and unknown, and (**b**) wetlands and mudflats. Non-surveyed coastlines were replaced by the current coastlines contained in the GSI map. Source: Digital Inoh's Map and GSI map.

indicates the example of derived gaps. The spatial difference (gap) between non-surveyed and actual coastlines was visualized as follows. First, the data for the coastlines were converted from line data to polygon data. Next, both polygon data were superimposed to derive the difference between the The spatial difference (gap) between non-surveyed and actual coastlines was visualized as follows. First, the data for the coastlines were converted from line data to polygon data. Next, both polygon data were superimposed to derive the difference between the

seaward and landward of the current coastline. The surveyed coastline data include a small number of coastlines based on the shape of the georeferenced Inoh's large-scale

because the purpose was to evaluate the accuracy of the non-surveyed coastline. Figure 6

because the purpose was to evaluate the accuracy of the non-surveyed coastline. Figure 6

indicates the example of derived gaps.

seaward and landward of the current coastline. The surveyed coastline data include a small number of coastlines based on the shape of the georeferenced Inoh's large-scale map. In this study, the gaps generated from the surveyed coastline were not considered because the purpose was to evaluate the accuracy of the non-surveyed coastline. Figure 6 indicates the example of derived gaps. *ISPRS Int. J. Geo-Inf.* **2021**, *10*, x FOR PEER REVIEW 6 of 13

**Figure 6.** Gap between the non-surveyed (estimated) coastline and the actual coastline. Source: Digital Inoh's Map and GSI map. **Figure 6.** Gap between the non-surveyed (estimated) coastline and the actual coastline. Source: Digital Inoh's Map and GSI map.

#### **3. Results 3. Results**

First, we identified the geographical features of the non-surveyed (estimated) coastlines. Then, we compared the difference (gap) between the non-surveyed (estimated) coastline and the actual coastline, focusing on the landscape features in the specified locations, given that the estimated coastline was explored from the land and sea. First, we identified the geographical features of the non-surveyed (estimated) coastlines. Then, we compared the difference (gap) between the non-surveyed (estimated) coastline and the actual coastline, focusing on the landscape features in the specified locations, given that the estimated coastline was explored from the land and sea.

#### *3.1. Relationship between Non-Surveyed Coastline Distribution and Topographical/Land 3.1. Relationship between Non-Surveyed Coastline Distribution and Topographical/Land Conditions*

*Conditions* Our spatial analysis reveals that 38.6% of the non-surveyed coastlines were cliffs, 25.7% were rocky beaches, and 6.2% were wetlands and tidal lands. The muddy coastline with tidal flats was also depicted on the basis of Inoh's estimation (including visual observation from the sea) without actual surveying, even on a flat coast. The results show that land reclamation modified 29.5% of the coastline over the last 200 years. However, Our spatial analysis reveals that 38.6% of the non-surveyed coastlines were cliffs, 25.7% were rocky beaches, and 6.2% were wetlands and tidal lands. The muddy coastline with tidal flats was also depicted on the basis of Inoh's estimation (including visual observation from the sea) without actual surveying, even on a flat coast. The results show that land reclamation modified 29.5% of the coastline over the last 200 years. However, the terrain of the coastlines of those days was unknown because there were no data.

the terrain of the coastlines of those days was unknown because there were no data. Figure 7 displays the distribution of non-surveyed coastlines by terrain and land classification. Cliffs and rocky beaches are prevalent in rias and peninsulas across Japan. Most of these coasts have kept the shapes unchanged during the last 200 years, and therefore the results of the analysis are valid. However, as the coastlines in wetlands and tidal lands were not surveyed in Inoh's map, the coastlines were not easy to identify. As is well known, various reclamation projects targeting wetlands and other regions had started in different parts of the country since the Meiji era (1868–1912), which had led to the rapid Figure 7 displays the distribution of non-surveyed coastlines by terrain and land classification. Cliffs and rocky beaches are prevalent in rias and peninsulas across Japan. Most of these coasts have kept the shapes unchanged during the last 200 years, and therefore the results of the analysis are valid. However, as the coastlines in wetlands and tidal lands were not surveyed in Inoh's map, the coastlines were not easy to identify. As is well known, various reclamation projects targeting wetlands and other regions had started in different parts of the country since the Meiji era (1868–1912), which had led to the rapid development of ports, industrial sites, and residential areas.

Cliffs are predominantly distributed along the cape's tip in the rias, and rocky beaches are adjacent to cliffs in the Tohoku region. The "unknown" coastlines are predominantly spread in the bay. "Unknown", which accounts for approximately 30% of all non-surveyed coastlines on the main islands, can be found scattered primarily near large cities and on intricate coasts. Specifically, "unknown" corresponds to harbors, factories,

idle land, and parks in the present.

development of ports, industrial sites, and residential areas.

**Figure 7.** Topographic and land use status of non-surveyed (estimated) coastline. Non-surveyed coastlines were replaced by the coastlines contained in the GSI map. Source: Digital Inoh's Map and GSI map. **Figure 7.** Topographic and land use status of non-surveyed (estimated) coastline. Non-surveyed coastlines were replaced by the coastlines contained in the GSI map. Source: Digital Inoh's Map and GSI map.

*3.2. Differences between Non-Surveyed Coastlines and Actual Coastlines* Figure 8 displays the distribution of the gaps in the non-surveyed coastlines. At first glance, it appears that the gaps are prominent on capes, especially in the rias. Bays tend to be less uncertain than capes. The non-surveyed coastlines with flat terrain have smaller gaps. As shown in ⑤ and ⑥ in Figure 8, the non-surveyed coastlines in western Japan Cliffs are predominantly distributed along the cape's tip in the rias, and rocky beaches are adjacent to cliffs in the Tohoku region. The "unknown" coastlines are predominantly spread in the bay. "Unknown", which accounts for approximately 30% of all non-surveyed coastlines on the main islands, can be found scattered primarily near large cities and on intricate coasts. Specifically, "unknown" corresponds to harbors, factories, idle land, and parks in the present.

#### have relatively high accuracy. In the next step, we examine the characteristics in the gaps, linking them with the *3.2. Differences between Non-Surveyed Coastlines and Actual Coastlines*

nized that the coastline facing the cliffs protruded more toward the sea.

terrain and land use conditions. First, terrain and land use gaps are calculated as the percentage of terrain and land use in the non-surveyed coastlines adjacent to the gaps. Table 1 shows the area of the gaps based on terrain and land use. As a result of comparing cliffs and rocky beaches, the total area of the gaps was 1.5 times larger for cliffs. In addition, the difference in the direction of the gaps (on the sea or landward of the current Figure 8 displays the distribution of the gaps in the non-surveyed coastlines. At first glance, it appears that the gaps are prominent on capes, especially in the rias. Bays tend to be less uncertain than capes. The non-surveyed coastlines with flat terrain have smaller gaps. As shown in <sup>5</sup> and <sup>6</sup> in Figure 8, the non-surveyed coastlines in western Japan have relatively high accuracy.

coastline) is more prominent for cliffs than for rocky beaches. The gap area on the seaward of the present coast is approximately 2 times larger than the landward. Inoh's team recog-

**Classification**

Rocky

Wetlands and tidal

**Figure 8.** Distribution of gaps in the non-surveyed (estimated) coastlines. The right maps show the enlarged views of each area. **Figure 8.** Distribution of gaps in the non-surveyed (estimated) coastlines. The right maps show the enlarged views of each area.

**Table 1.** Number of gaps and their areas in the non-surveyed (estimated) coastlines in Japan. **Total Seaward of the Present Landward of the Present Coastline Average Max**  In the next step, we examine the characteristics in the gaps, linking them with the terrain and land use conditions. First, terrain and land use gaps are calculated as the percentage of terrain and land use in the non-surveyed coastlines adjacent to the gaps.

**(km<sup>2</sup> ) Coastline (km<sup>2</sup> ) (km<sup>2</sup> ) Number Area (km<sup>2</sup> ) Number Area (km<sup>2</sup> )** Cliffs 192.5 1553 126.5 1550 66.0 0.1 6.7 beaches 128.8 <sup>1135</sup> 72.9 <sup>1090</sup> 55.9 0.1 7.8 468.0 60 451.3 79 16.7 1.5 78.4 Table 1 shows the area of the gaps based on terrain and land use. As a result of comparing cliffs and rocky beaches, the total area of the gaps was 1.5 times larger for cliffs. In addition, the difference in the direction of the gaps (on the sea or landward of the current coastline) is more prominent for cliffs than for rocky beaches. The gap area on the seaward of the present coast is approximately 2 times larger than the landward. Inoh's team recognized that the coastline facing the cliffs protruded more toward the sea.


#### total distance was 468.0 km, accounting for 10.9% of the whole non-surveyed coastline. The maximum and minimum length difference was greater than 40 km (max = 44.7 km, *3.3. Non-Surveyed Coastline Drawn by Sea-Based Measurement*

min = 0.1 km). The average length was 1.5 km per section. The average distance became shorter in the later years of the ten survey trips. This means that the measurement from the sea became more elaborate over time. The accuracy of the coastlines drawn by sea-based measurement (without actual surveying) was investigated, and the working section was visualized spatially using GIS (Figure 9). The exploration from the sea was often found in areas with complex coastal terrain. For example, in the Tohoku region, where rias were dominant, 61.4% was occupied by sea-based measurements.

**Figure 9.** Distribution of the map sections drawn by sea-based measurement. The right maps show the enlarged views of each area. Source: Digital Inoh's Map. **Figure 9.** Distribution of the map sections drawn by sea-based measurement. The right maps show the enlarged views of each area. Source: Digital Inoh's Map.

Table 2 shows the topographical characteristics of the non-surveyed coastlines drawn by sea-based measurement. Excluding "unknown", the rocky beaches (147.6 km, 31.5%) reveal the most significant distribution of the coastlines measured by boat (not actual surveying). Here, "unknown" includes small sandy beaches between cliffs, such as Ago Bay, and ports in bays without rocky beaches. On the other hand, sea-based measurement by boat was not generally carried out in wetlands and tidal lands. Measuring from the sea by boat was conducted 303 times in the whole of Japan. The total distance was 468.0 km, accounting for 10.9% of the whole non-surveyed coastline. The maximum and minimum length difference was greater than 40 km (max = 44.7 km, min = 0.1 km). The average length was 1.5 km per section. The average distance became shorter in the later years of the ten survey trips. This means that the measurement from the sea became more elaborate over time.

In Table 3, excluding "unknown," the total area of the gaps is the largest for rocky beaches (13.9 km<sup>2</sup> ). In addition, the direction of the gaps is more outstanding on the landward of the current coastline for both cliffs and rocky beaches. **Table 2.** Distance of non-surveyed (estimated) coastlines drawn by sea-based measurement in Japan. Table 2 shows the topographical characteristics of the non-surveyed coastlines drawn by sea-based measurement. Excluding "unknown", the rocky beaches (147.6 km, 31.5%) reveal the most significant distribution of the coastlines measured by boat (not actual surveying). Here, "unknown" includes small sandy beaches between cliffs, such as Ago Bay, and ports in bays without rocky beaches. On the other hand, sea-based measurement by boat was not generally carried out in wetlands and tidal lands.


**Classification Total (km) Rate of the Total (%) Table 2.** Distance of non-surveyed (estimated) coastlines drawn by sea-based measurement in Japan.

**Classification Total (km<sup>2</sup> ) Seaward of the Present Coastline Landward of the Present Coastline Average (km<sup>2</sup> ) Max (km<sup>2</sup> ) Number Area (km<sup>2</sup> ) Number Area (km<sup>2</sup> )** In Table 3, excluding "unknown," the total area of the gaps is the largest for rocky beaches (13.9 km<sup>2</sup> ). In addition, the direction of the gaps is more outstanding on the landward of the current coastline for both cliffs and rocky beaches.

Cliffs 6.4 165 1.9 165 4.5 0.1 0.6 Rocky beaches 13.9 184 3.0 167 10.9 0.1 0.8

tidal lands 1.3 9 0.4 8 0.9 0.1 0.6 Unknown 25.8 405 15.2 382 10.6 0.1 5.4

section of sea-based measurement).

Wetlands and


**Table 3.** Number of gaps and their areas in the non-surveyed (estimated) coastlines in Japan (the section of sea-based measurement). **4. Discussion** This section discusses the validity of our study and the contribution to historical GIS.

#### **4. Discussion** in the middle Meiji period (Ihoh estimated the coastline of Tokyo Bay along the coastal roads because the actual survey was not possible). Tokyo Bay was covered with tidal flats

rent GSI map.

tice.

This section discusses the validity of our study and the contribution to historical GIS. Inoh's map had the highest accuracy, but Inoh could not survey 25% of the total coastlines. The non-surveyed coastlines drawn by estimation were in places where actual surveying was impossible due to technical problems. in the early 19th century (which has been reclaimed into agricultural land, factory sites, housing estates, etc., since the modernization period). Using GIS, we superimposed the non-surveyed coastlines in Inoh's map, the old land use maps in 1881–1882, and the cur-

70.5% of these non-surveyed coastlines were in "cliffs, rocky beaches, wetlands and tidal lands". In this study, we classified the remaining 29.5% as "unknown". Unfortunately, we cannot obtain the exact coastlines of the "unknown". In the early 19th century, most of the unknown coastlines were in wetlands or tidal flats and are on reclaimed land today. A good example showing this situation is in Tokyo Bay (Figure 10). This figure shows the spatial relationship between the non-surveyed coastline and land use patterns in the middle Meiji period (Ihoh estimated the coastline of Tokyo Bay along the coastal roads because the actual survey was not possible). Tokyo Bay was covered with tidal flats in the early 19th century (which has been reclaimed into agricultural land, factory sites, housing estates, etc., since the modernization period). Using GIS, we superimposed the non-surveyed coastlines in Inoh's map, the old land use maps in 1881–1882, and the current GSI map. We found that the inland distributions of tidal flats and non-surveyed (estimated) coastlines were almost consistent. This approach would be applicable to the coastal areas of large cities such as Nagoya and Osaka and regions along the coast of the Sea of Japan, where land reclamation had been very active (fortunately, many land use maps are available for these regions). Therefore, the proportion of wetlands and tidal lands (6.2%) seems to be underestimated because the unknown (29.5%) includes swampy lands. In other words, non-surveyed (estimated) coastlines, which were in places where actual surveying was difficult, were expected to be more than 70.5%. In addition, since Japan belongs to the Circum-Pacific Belt, more than two-thirds of the country are mountainous [32]. Thus, steep topographical surfaces might have had a negative influence on the surveying prac-

**Figure 10.** Overlaying of the non-surveyed coastline and land use maps (1881–1882). The land use classification is based on the 1st Military District 1:20,000 Jinsoku-Sokuzu map. Source: Digital Inoh's Map and Geospatial Information Authority of Japan. **Figure 10.** Overlaying of the non-surveyed coastline and land use maps (1881–1882). The land use classification is based on the 1st Military District 1:20,000 Jinsoku-Sokuzu map. Source: Digital Inoh's Map and Geospatial Information Authority of Japan.

We found that the inland distributions of tidal flats and non-surveyed (estimated) coastlines were almost consistent. This approach would be applicable to the coastal areas of large cities such as Nagoya and Osaka and regions along the coast of the Sea of Japan, where land reclamation had been very active (fortunately, many land use maps are available for these regions). Therefore, the proportion of wetlands and tidal lands (6.2%) seems to be underestimated because the unknown (29.5%) includes swampy lands. In other words, non-surveyed (estimated) coastlines, which were in places where actual surveying was difficult, were expected to be more than 70.5%. In addition, since Japan belongs to the Circum-Pacific Belt, more than two-thirds of the country are mountainous [32]. Thus, steep topographical surfaces might have had a negative influence on the surveying practice.

As Figure 8 shows, overall, the non-surveyed (estimated) coastlines were more accurate in western Japan than in eastern Japan because of the solid financial support by the central government [3]. Focusing on the landscape features, the gap between the non-surveyed coastlines and the actual coastlines was the greatest in cliffs (Table 1). It means that the estimation of rocky beaches is more accurate in comparison with cliffs. Most rocky shores were in bays, and the distance from the estimation point was relatively short. Therefore, rocky shores were more accessible for Inoh's team.

On the other hand, the largest distributor of the non-surveyed coastlines in the seabased measurement was found on rocky beaches (Table 2). Comparison of the accuracy between the land-based and sea-based measurements shows that the error was minor for the sea-based measurement, regardless of the landscape features. Inoh's team often rented boats from villagers [33]; boat use might have been effective in reducing uncertainty because the coastline gaps estimated by sea-based measurements were generally small (as seen by comparing Tables 1 and 3).

According to the direction of the gaps, we can say that the seaward for the land-based estimates and the landward for the sea-based estimates were prominent, respectively (Tables 1 and 3). It suggests that Inoh perceived the coastline as further away than it was, regardless of forecasts from the land or sea. Therefore, we conclude that a specific spatial distortion existed in Inoh's cognition of the distance, more or less.

Our findings contribute to interpreting the uncertainty in the past coastlines and the use of GIS in historical maps. Firstly, as [34] states, geomorphological factors were not discussed in the previous historical GIS studies, although the land use pattern was wellconsidered. In contrast, we could examine the relationship between surveying situations and topography using Inoh's map and GIS techniques. Secondly, our study succeeded in delineating the surveying activity quantitatively on a nationwide scale. On the contrary, in previous studies, surveying was mainly analyzed using the survey logbook [33]. Thirdly, our study examined the coastline accuracy (i.e., Inoh's cognition) from the geographical conditions of the coasts. Thus, this study might be the first approach combining historical GIS and critical cartography [35].

#### **5. Conclusions**

For the discussion on the accuracy of Inoh's map, it is crucial to examine why these coastlines were not directly surveyed. This study focused on the non-surveyed coastline drawn by estimation in Inoh's map and investigated their distributions, physical conditions, and differences with actual coastlines at that time. Over 70% of all non-surveyed coastlines were situated in "cliffs, rocky beaches, wetlands, and tidal lands". We conclude that the non-surveyed coastlines were in places where actual surveying was not possible by using the surveying technology available at that time. The approach used in this study can be applied for examining the survey accuracy of historical maps in other countries.

Finally, we would like to note several remaining tasks and future research possibilities. First, this study assumed that geographical changes in coastal areas are small in cliffs and rocky beaches because man-made modifications to the cliffs and rocky shores are limited compared to wetlands and tidal lands. In this connection, we thought the present-day GSI map was substitutable as a base map showing steep topographical areas in the early

19th century. However, we must note that the coastline shape of cliffs and rocky beaches is continuously changing due to erosion by waves (at least 0.3m per year) [36]. Although our results show that the difference between the GSI and Inoh's maps tends to be larger than the effect of erosion (Table 1), a more sophisticated approach is necessary to check the influence of the erosion phenomena. Second, the Digital Inoh's Map (Professional Edition) we used was a georeferenced map, not the original. In addition, we regarded the surveyed coastlines to be consistent with the GSI map. However, we found there was distortion between them, and thus we must check the setting of the control points in the Digital Inoh's Map. Third, cliffs and rocky beaches in some coastal areas were directly surveyed by Inoh's team. Thus, it is an important task to investigate why actual surveying became possible in these coastal areas.

There is room to study Inoh's map from a historical GIS viewpoint. Inoh's survey logbook (diary) is also helpful in assessing Inoh's team's activities. As Inoh's map is valuable historical material among various maps prepared on a national scale, empirical scientific research on regional differences is required to develop quantitative historical geography.

**Author Contributions:** Conceptualization, Yuki Iwai and Yuji Murayama; writing—original draft preparation, Yuki Iwai and Yuji Murayama; writing—review and technical analysis for spatial analysis and visualization, Yuki Iwai; supervision, Yuji Murayama. All authors have read and agreed to the published version of the manuscript.

**Funding:** We used Grant-in-Aid for Scientific Research (C) 21K01027 to carry out this research.

**Acknowledgments:** The authors would like to express their gratitude to Kota Inohara (Tokyo Cartographic Co., Ltd.) for arranging the surveying data in "Digital Inoh's Map". We are indebted to the Inoh Tadataka Study Group members for the place name index and the translation of the survey diary into modern languages. We are also grateful for the insightful and constructive comments of the audience at the AAG Annual Meeting 2021.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


## *Article* **GIScience and Historical Cartography for Evaluating Land Use Changes and Resulting Effects on Carbon Balance**

**Canio Manniello <sup>1</sup> , Giuseppe Cillis 1,\* , Dina Statuto <sup>1</sup> , Andrea Di Pasquale <sup>2</sup> and Pietro Picuno <sup>1</sup>**


**Abstract:** Multi-chronological examination of territory using GIScience and historical cartography may reveal a strategic tool for investigating changes in land use and the surrounding landscape structure. In this framework, the soil plays a key role in ecosystem evolution, since it governs all the mechanisms at the basis of vegetal growth, as well as all components of the total environment contributing to the formation of a rural landscape, including the balance of carbon dioxide. The present study was developed using a GIS approach applied to historical maps and aims to assess the environmental impact of land-use change, with particular attention to its effects on agricultural soil and atmospheric carbon dioxide balance. Thanks to a comparison between historical cartographic maps of different periods, this geospatial approach has enabled the assessment of the evolution of the rural land of the study area in the municipality of Ruoti (Basilicata Region—Southern Italy). This area, indeed, has been affected by deep land-use transformations, mainly caused by agricultural activities, with a resulting impact on the atmospheric CO<sup>2</sup> balance. These transformations have been analyzed and quantified in order to contribute to the understanding on how the changes in land use for agricultural purposes have led to unforeseen changes in the rural landscape, ecosystems and the environment. The results showed that the greatest changes in land use were caused by the abandonment of large rural areas, resulting in the expansion of urban areas, a decrease in orchard and arable land (about less 25%), and an increase in woodland (more than 30%). These changes have resulted in a doubling in soil carbon fixation value. The final results have therefore confirmed that historical cartography within a GIS approach may decisively offer information useful for more sustainable agricultural activities, so as to reduce their negative contribution to climate change.

**Keywords:** historical maps; GIS; land use; carbon balance; rural landscape; total environment

#### **1. Introduction**

The natural dynamics that have followed one another slowly over millions of years on Earth have triggered processes and phenomena that have been the basis of the evolution of different habitats and living organisms. Now this natural balance has been strongly compromised by anthropogenic activity [1]. In fact, scientists speak of these last two centuries as the Anthropocene period, in which human beings have interfered with the slow natural dynamics, causing often irreversible changes to natural ecosystems and environment [1], alterations that can be experienced in all the different terrestrial environments. These phenomena, however, are only now being addressed with interest, but the lack of concrete activities of sustainable development and detailed knowledge of the impacts at the local level represent a further slowdown to the trigger processes of environmental and ecological reconstruction [2]. The resulting changes have been so severe that, in some cases, the entire structure of the terrestrial ecosystem itself was altered. It is necessary therefore to pursue

**Citation:** Manniello, C.; Cillis, G.; Statuto, D.; Di Pasquale, A.; Picuno, P. GIScience and Historical Cartography for Evaluating Land Use Changes and Resulting Effects on Carbon Balance. *ISPRS Int. J. Geo-Inf.* **2022**, *11*, 179. https:// doi.org/10.3390/ijgi11030179

Academic Editors: Wolfgang Kainz and Motti Zohar

Received: 2 February 2022 Accepted: 5 March 2022 Published: 8 March 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

a multidisciplinary approach that takes into account geographical, environmental and landscape factors as variables that interact with each other and with social and economic aspects [3]. In this contest, a detailed analysis of the changes carried out and the global monitoring of all ecosystems is necessary in order to propose appropriate environmental protection policies [4] aimed at reusing resources in the framework of a circular economy approach [5].

#### *1.1. New Spatial Planning Techniques*

Among the most important environmental consequences induced by land-use change from man is the loss of fertility [6], the modification of the surrounding rural landscape [7,8], and variation in atmospheric carbon dioxide [9]. The greatest loss in carbon storage per hectare results from the conversion of forests to cropland. Indeed, forests hold 20–50 times more carbon per hectare than cleared lands, and 100–200 MgC/ha may be lost to deforestation [10].

Since the beginnings of the industrialization period, CO<sup>2</sup> concentration levels have increased by 37%, from 260–280 ppm in 1880, to 360–380 ppm in 2001. Following the trend of the last 30 years, the current increase is almost 1.0–1.5 ppm per year [11,12]. More recent assessments by the Scripps Institutions of Oceanography (California) and the National Oceanic and Atmospheric Administration (Maryland) state that the world exceeded the 400 ppm threshold for the first time in 2013, reaching a record value of 420 ppm in May 2021. The heavy use of fossil fuels, on the one hand, and deforestation (and the resulting soil erosion) on the other are the main causes of the current strong qualitative and quantitative change in atmospheric gases, especially CO2, and the resulting phenomenon of global overheating.

This "greenhouse effect" of planetary dimension is the result of the sum of individual local conditions. Real and useful tools for the sensible planning of agricultural land by politicians and planners, which should also take environmental aspects into account, have only recently been developed. The relationship between agriculture, ecosystems, and environment has been proposed by some authors [13,14] as new contributions to turf and landscape design and land planning and management.

On the other hand, integrated GIS-based techniques are fundamental for the rational acquisition and analysis of forest and agricultural land data. GIS-based techniques, image processing, remote sensing and other new technologies for land development surveying, planning and management allow a more accurate analysis of rural landscape and environment today [15,16]. Indeed, besides vegetation, there are other elements related to the landscape, such as buildings, which should be properly considered in data processing [17]. If historical sources, such as cadastral registers, old maps, etc. [18,19], are available, it might also be possible to analyze the evolution of the rural landscape over time.

Geographical information systems (GIS) are excellent tools for landscape modelling, for knowing changes in vegetation and for performing three-dimensional analyses [20]. They allow the easy digitization of geographic information and land-cover structure, and they facilitate graphical representation [21,22].

#### *1.2. Historical Cartography for Land Transformation Analyses*

In order to evaluate territorial transformations and dynamics, it is possible to use different tools. One of these is represented by historical cartography, which can be integrated within the GIS environment [23,24]. The integrated approach can be used to obtain georeferenced information from old caldastral maps, topographic maps, military maps, aerial photos, and landscape and thematic maps, etc., and compare them with current digital geodata with high accuracy [25–28]. Indeed, all these historical documents are "speaking tools", able to inform us on life during previous centuries, including the modifications imposed by humankind on Earth. What is now crucial is our ability to detect, analyze, interpret and understand this information [29,30]. Therefore, these documents can now be

recovered and used to analyze the temporal evolution of non-urban land, environment and landscapes.

Pärtel et al. [31] described the landscape history of the largest calcareous semi-natural site in Estonia, showing a greater decrease of species in 20–40-year-old forests than in open grassland. Jordan et al. [32] investigated the influence of historical land-use changes on soil erosion and sediment transport in the Kali basin (Lake Balaton, Hungary) using historical maps from 1784 onwards. Bender et al. [33] investigated two study cases aimed at studying landscape changes since 1850 using cadastral maps and land registers. The results were presented at the plot level using a diachronic GIS, which provides valuable information for planning processes and conservation in changing cultural landscapes.

There have been many uses of historical topographic maps for the analysis of rural landscapes: Stauble et al. [34] for example, used them to reconstruct historical changes since the mid-19th century of the landscape of the canton of Valais in Switzerland, while Gimmi et al. [35] used them to analyze the change in wetlands in the canton of Zurich (Switzerland) over the last 150 years. The unique potential of these maps has also been exploited by the archeologists Williams [36] and Barclay [37]. Skaloš et al. [38] similarly studied long-term land-cover changes in the Czech Republic using old military survey maps [39]. Historical topographic maps, old vegetation maps and military survey maps have been used for the investigation of land-cover distribution and changes in landscape visual quality [29]. Historical cadastral maps were used in studies of Domaas [40] and Trpák [41].

More recently, Lievosky et al. [42], thanks to historical maps, have produced the first digital land use map for a 160-year period in the Carpathian Ecoregion, the Hungarian part of the Pannonian plains and the historical Moravia region in the Czech Republic as the study area. Finally, Valent et al. [43] analyzed land-use changes over a period of almost three centuries in the Myjava River basin, while Chen et al. [44] reconstructed land-cover changes in Taiwan between 1904 and 2015 from historical maps and satellite images.

This paper is intended to be a preliminary analysis of the assessment of the small-scale CO<sup>2</sup> budget based on the decrease/increase of natural resources due to land changes. Indeed, land-use dynamics over almost two centuries (from 1848 to 2017) were investigated by the comparison between historical maps, aerial photos and more recent orthophotos through the use of an open source software GIS in order to highlight the rural landscape changes, their connection with the natural cycles and anthropic activities, and the impacts on agricultural land and the physical environment. Finally CO<sup>2</sup> variations, by comparing chronologically different land-cover maps, has been evaluated. In the simple and expeditious survey presented in this paper, we have exploited the potential of spatialized historical data that allows for further detailed analysis of the CO<sup>2</sup> balance at the local level.

#### **2. Materials and Methods**

#### *2.1. Study Area*

The study area (about 25 km<sup>2</sup> ) is part of the municipality of Ruoti (Figure 1), located in the central–western part of the Basilicata region in southern Italy (40◦43005.4300 N, 15◦40032.1200 E), and the study area is bounded by the perimeter of the historical cartography of 1848 taken as a reference and is an area of strong interest and attention from regional and Italian planners. The Basilicata region represents much of the landscape variability of southern Italy. In addition to geological variability, the area of the region exhibits a high morphological variability, with the presence of surfaces of different periods and a great variability in the soils that have formed in these environments. This area is characterized by a hilly and mountainous landscape, and the altitudes range from 400 to 1000 m.

The municipality of Ruoti is situated on a hill dominating the course of the Avigliano river. In this territory, besides sheep and goat breeding, from whose milk excellent cheeses are obtained, cereal cultivation—in particular wheat, fodder, vegetables—is very widespread. Of particular value in this area are olive groves, orchards, and vineyards, from which the "Asprino" wine is produced. Industry is chiefly addressed to the food sector: the

*ISPRS Int. J. Geo-Inf.* **2022**, *11*, x FOR PEER REVIEW 4 of 14

cessing of straw and wicker [45,46].

chief product is dairy. In the field of craftsmanship, the area is renowned for the processing of straw and wicker [45,46]. nated by farmland, which is mainly present in the hills, with extensive pastures and vineyards in the northern area. The hills are covered with rich forests, consisting mainly of undergrowth, such as fir wood.

The study area mainly consists of agricultural land (57%), forested and semi-natural land (38%) and artificial surfaces (5%). The high hill landscape of the study area is domi-

are obtained, cereal cultivation—in particular wheat, fodder, vegetables—is very widespread. Of particular value in this area are olive groves, orchards, and vineyards, from which the "Asprino" wine is produced. Industry is chiefly addressed to the food sector: the chief product is dairy. In the field of craftsmanship, the area is renowned for the pro-

**Figure 1.** Geographical localization and altitude profile of the study area. The digital terrain model for the altitude map and the 2017 orthophotos have been downloaded from the geoportal of the Basilicata region under license IODL 2.0 [47]. **Figure 1.** Geographical localization and altitude profile of the study area. The digital terrain model for the altitude map and the 2017 orthophotos have been downloaded from the geoportal of the Basilicata region under license IODL 2.0 [47].

*2.2. The Cartography* This work is based on the integration and updating of the calculations proposed by Tortora et al. [20], since new methodologies for processing historical data have been developed, making the data more accurate. To understand the changes that the study area has undergone in the last two centuries, three different periods have been analyzed, i.e.,: The study area mainly consists of agricultural land (57%), forested and semi-natural land (38%) and artificial surfaces (5%). The high hill landscape of the study area is dominated by farmland, which is mainly present in the hills, with extensive pastures and vineyards in the northern area. The hills are covered with rich forests, consisting mainly of undergrowth, such as fir wood.

#### 1848, 1953, and 2017. *2.2. The Cartography*

2.2.1. 1848 Historical Map The most important historical document of the Ruoti municipality is represented by a historical map produced in 1848 (Figure 2). With reference to the study area and the surrounding territory, it represents a valuable source of information on land use during that period. This is an iconographic type This work is based on the integration and updating of the calculations proposed by Tortora et al. [20], since new methodologies for processing historical data have beendeveloped, making the data more accurate. To understand the changes that the study area has undergone in the last two centuries, three different periods have been analyzed, i.e.,: 1848, 1953, and 2017.

#### map realized in watercolor. The surveying techniques with which it has been realized are based on the topographic techniques and on the surveying instruments of the time. The 2.2.1. 1848 Historical Map

map, in paper format, is preserved at the State Archives of the Municipality of Potenza (Italy) and is freely available for consultation and digitalization. The most important historical document of the Ruoti municipality is represented by a historical map produced in 1848 (Figure 2).

With reference to the study area and the surrounding territory, it represents a valuable source of information on land use during that period. This is an iconographic type map realized in watercolor. The surveying techniques with which it has been realized are based on the topographic techniques and on the surveying instruments of the time. The map, in paper format, is preserved at the State Archives of the Municipality of Potenza (Italy) and is freely available for consultation and digitalization.

It shows the main rivers of the area (Fiumara di Ruoti and Fiumara di Avigliano), land classified as unirrigated farmland and irrigated territory (north of the area), alternation of olive groves and farmland in the central part of the area, and the distribution of vineyards in the eastern part, as well as part of the area classified as woodland. The legend located

in the western part of the map shows the territorial extent of the categories of vegetation represented here. The place names of different districts are also reported, but the farms and the roads within the territory are unfortunately not reported. *ISPRS Int. J. Geo-Inf.* **2022**, *11*, x FOR PEER REVIEW 6 of 14

**Figure 2.** Historical map of 1848 preserved by a protective film with some zoom on details. In the first detail (**A**), it is possible to observe the representation of olive groves and vineyards and part of the river zone; in the second detail (**B**), the representation of woodlands; in the third detail (**C**), instead the urban area of Ruoti and, in the fourth detail (**D**), the arable land and pastures. **Figure 2.** Historical map of 1848 preserved by a protective film with some zoom on details. In the first detail (**A**), it is possible to observe the representation of olive groves and vineyards and part of the river zone; in the second detail (**B**), the representation of woodlands; in the third detail (**C**), instead the urban area of Ruoti and, in the fourth detail (**D**), the arable land and pastures.

2.2.2. 1953 Aerial Photos The different land use categories for the year 1953 were obtained from an aerophotogrammetric survey conducted in the 1950s. The process of rectifying and georeferencing these frames was not easy, because the camera calibration certificate was not available. In These historical maps contain a lot of information, but in order to use them properly from a geographical point of view, their accuracy must be well assessed [48,49], so new approaches had to be applied to transform the maps from simple archival documents into geographic datasets [50].

this case, a series of georeferencing operations was carried out for each group of rural buildings, identifying, time by time, a sufficient number of ground control points (GCP) on more recent orthophotos (1988–2017) already georeferenced. A polynomial series correction was applied, which ensured the root mean square error remained even below 5 m in some cases [54]. This was possible because the study areas are not very extended in surface and do not have any particular features. 2.2.3. 2017 Orthophotos To obtain a land-use map for the year 2017, digital orthophotos were used, which match the image properties of a photograph with the geometric properties of a map. They allow the direct measurement of distances, areas, angles, and positions because the relief shift in the orthophotos had been removed, so that land features were shown in their actual ground position (Figure 3). First, they have been scanned and converted to a digital format; then they were georeferenced using a method already used for similar maps in a previous study [51]. The quality of geo-referenced map was assessed by overlaying the historical maps with the 2013 Regional Technical Map (CTR) of the Basilicata region at a scale of 1:5000. In some cases, due to the characteristics and realization techniques of historical maps, the changes in the landscape since the time of mapping, and taking into account the ground roughness of the study area, the quadratic mean (RMS) was over 100 m. Indeed, using a simple Helmert transformation (shift, rotation, change of scale) as usually used for topographic maps, the errors on the control points were very high. Therefore the historical map was treated as if it were an unorthorectified image. Increasing the number of control points in order to cover all the different areas of the study area, a second-order polynomial was used, so as to have an average RMS value of about 20 m. The average value was calculated on the RMS of each point, since QGIS allows the evaluation of RMS error for each point, both visually and in tabular format. Since the study area is small, it was not necessary to apply further transformations or assessments to improve the georeferencing result. In order to reduce it,

a second-order polynomial was used to obtain a RMS value of about 20 m [52], which also allowed a visual accuracy assessment, validating the georeferencing [53].

#### 2.2.2. 1953 Aerial Photos

The different land use categories for the year 1953 were obtained from an aerophotogrammetric survey conducted in the 1950s. The process of rectifying and georeferencing these frames was not easy, because the camera calibration certificate was not available. In this case, a series of georeferencing operations was carried out for each group of rural buildings, identifying, time by time, a sufficient number of ground control points (GCP) on more recent orthophotos (1988–2017) already georeferenced. A polynomial series correction was applied, which ensured the root mean square error remained even below 5 m in some cases [54]. This was possible because the study areas are not very extended in surface and do not have any particular features.

#### 2.2.3. 2017 Orthophotos

To obtain a land-use map for the year 2017, digital orthophotos were used, which match the image properties of a photograph with the geometric properties of a map. They allow the direct measurement of distances, areas, angles, and positions because the relief shift in the orthophotos had been removed, so that land features were shown in their actual ground position (Figure 3). *ISPRS Int. J. Geo-Inf.* **2022**, *11*, x FOR PEER REVIEW 7 of 14

**Figure 3.** The same portion of the study area in 1953 aerial photos and 2017 orthophotos. **Figure 3.** The same portion of the study area in 1953 aerial photos and 2017 orthophotos.

#### *2.3. Land Use Data Elaboration 2.3. Land Use Data Elaboration*

herbaceous plants).



that occupy irregular but significant areas.

With reference to each time period considered and the maps shown above, different categories of land use have been identified, taking into account the symbols present on the map and the level of detail of the cartographic base. Given the small size of the study area, a manual digitization of land use was chosen [51]. With reference to each time period considered and the maps shown above, different categories of land use have been identified, taking into account the symbols present on the map and the level of detail of the cartographic base. Given the small size of the study area, a manual digitization of land use was chosen [51].

In the 1848 historical map, eight land-use categories have been highlighted and identified. For the year 1953, on the other hand, 10 land-use categories have been identified, thanks to the analysis of georeferenced aerial photographs and with a higher level of detail. Finally, the colored orthophotos for the 2017 map resulted in 11 different land-use categories. Therefore, to standardize the data and allow for more direct comparison, the major land-use categories were grouped together and defined as "elements". Seven major "elements" were identified in the study area. For each of these elements, the total area in hectares (ha) and percentage (%) was calculated using the function GIS. Their change over the years was also calculated. The major elements are described below: In the 1848 historical map, eight land-use categories have been highlighted and identified. For the year 1953, on the other hand, 10 land-use categories have been identified, thanks to the analysis of georeferenced aerial photographs and with a higher level of detail. Finally, the colored orthophotos for the 2017 map resulted in 11 different land-use categories. Therefore, to standardize the data and allow for more direct comparison, the major land-use categories were grouped together and defined as "elements". Seven major"elements" were identified in the study area. For each of these elements, the total area in hectares (ha) and percentage (%) was calculated using the function GIS. Their change over the years was also calculated. The major elements are described below:




River zone: These were considered in their areal extent, i.e., the river bed and the vegetation present along the river, both in their linear extent. In this study area, the river has a seasonal flow, and it has a riverbed almost entirely covered by hygrophilous shrub vegetation. Streams and rivers cause erosion, with the consequent movement and displacement of sediments of the bottom of the bed or bank, in a watercourse, by the action of the water current. Most of the erosion process is caused by both rainwater and surface water flowing downstream. This process, which depends on hydraulic properties, sediment properties, and local characteristics, is able to change the morphological structure of

trees. - Arable land: Cultivated land with different crops (i.e., in this study: cereals); - Arable land: Cultivated land with different crops (i.e., in this study: cereals);

the area. Much of the present landscape is the result of a process of erosion.


River zone: These were considered in their areal extent, i.e., the river bed and the vegetation present along the river, both in their linear extent. In this study area, the river has a seasonal flow, and it has a riverbed almost entirely covered by hygrophilous shrub vegetation. Streams and rivers cause erosion, with the consequent movement and displacement of sediments of the bottom of the bed or bank, in a watercourse, by the action of the water current. Most of the erosion process is caused by both rainwater and surface water flowing downstream. This process, which depends on hydraulic properties, sediment properties, and local characteristics, is able to change the morphological structure of the area. Much of the present landscape is the result of a process of erosion.


#### *2.4. Carbon Dioxide Calculation*

With the aim to quantify the impact of land-use changes on the environment, with particular attention to air quality, we estimated the CO<sup>2</sup> time variations associated with the use of the crops (arable land, natural land, pastures, river zone, vineyards, olive groves, built-up areas) in the study area in the three different periods (1848, 1953, and 2017). CO<sup>2</sup> sequestration rates were calculated using the model in CO2FIX V.2 [55], a tool for the dynamic estimation of the carbon sequestration potential of forest management, agroforestry and afforestation projects. The C++ programming language was used to execute the model, and, as output, it returned data in tabular form. Thanks to it, we could evaluate the amount of carbon sequestered by forest stands and especially its evolution over time. This value is given by the sum of the carbon stored respectively in the living biomass, in the organic matter of the soil, and in wood products.

The carbon stored by living biomass is calculated using a forest cohort model or groups of individual stands considered homogeneous within the model. Carbon stored in the entire forest stand is expressed as the sum of the contributions of the individual existing stands. This approach considers several parameters such as competition, natural mortality, logging, and mortality due to logging damage. Soil carbon, on the other hand, considers five stock pools, three for litter and two for humus. Finally, carbon stored in wood products is determined by considering processing efficiency, by-product reuse (mainly for energy production), recycling, and disposal forms (wood by-products are left on the site with the potential to enter in the soil). This model has a wide applicability for both temperate and tropical conditions [55]. To initialize the model, several parameters and assumptions consistent with the software input characteristics [56] and the local area characteristics [15] were used. The main characteristics used for the forestry area were tree species, area, age, dominant height, standing volume, growth class, and the coordinates of the stand.

*Quercus cerris* is the most common specie in the study area. The rotation is 80 years long, and the maximum biomass in the stand is 2000 Mg/ha. Existing CO2FIX runs for comparable species were used to retrieve the allocation factor for foliage, branches, and root formation. The turnover (annual rate of mortality of the biomass component) was calculated to be 0.3 for foliage, 0.06 for branches, and 0.05 for roots.

The organic matter compartment in the soil is made up of dead wood, litter layers, and stable humus. A total carbon stock ranging from 32 to 134 Mg/ha and an average atmospheric carbon sequestration roughly 25 MgC/ha/yr were determined based on this investigation.

In the study area, the orchard areas are mostly vineyards with the occasional presence of an olive grove. For CO<sup>2</sup> calculation, the orchard area was compared to a forest of tall trees with a rotation of 20 years and the periodical removal of organic matter through agronomic procedures like pruning, with a turnover of 0.3 for foliage, 0.07 for branches, and 0.04 for roots. In an orchard, carbon balance is determined by the intrinsic structural and morphological characteristics of each species and is also influenced by population density, rearing system, and especially on the canopy and aboveground and underground woody organisms. Furthermore, in the case of a new plantation, the canopy must provide for a relatively small number of branches and roots and, consequently, primary production is net-positive, and the surplus of organic matter grows year after year until maturity, when dry matter increases and then approaches zero [56]. According to this theory, orchards sequester 7.25 MgC/ha/yr, shrubland 2.75 MgC/ha/yr, and arable land 3.6 MgC/ha/yr.

On the other hand, urban areas represent a source of CO<sup>2</sup> emissions from both municipal and industrial combustion; a yearly amount of 15.0 MgC/ha/yr of CO<sup>2</sup> release into the atmosphere was therefore estimated based on a report on the environmental state of Basilicata [57].

All the above-mentioned values of average atmospheric carbon sequestration were adopted for each one of the three time periods (1848, 1953, and 2017).

The data resulting from the implementation of the GIS gave the values reported in Table 1 expressed in terms of areas occupied by the different vegetation typologies and, applying their respective CO<sup>2</sup> sequestration rates, in terms of the absolute values of the annual sequestration of CO2. The balance of CO<sup>2</sup> does not include the effects of the agricoltural machinery, supplies, and transportation on CO2: in woodland, these factors are almost absent, while in the case of orchard and arable land, they depend strongly on crop techniques and in some cases are negligible.


**Table 1.** Land-use analysis from 1848 to 2017.

#### **3. Results and Discussions**

*3.1. Land Use*

Table 1 presents landscape use for each indication and for each different cartographic base; the comparison allowed the analysis of land-use changes from 1848 to 2017, covering a 169-year period and providing information on the historical persistence of land-use typologies together with their changes over time. The predominant land-use typologies of the site have been grouped, to better compare output data.

Visualizing the spatial and graph data expressed in Figure 4, we can observe a decrease during these 169 years in the arable land, vineyards and olive groves and an increase in natural land, built-up areas, pastures and river zones. The urbanized area has experienced considerable growth, especially after World War II. As evidenced by a widespread trend in the various areas of the region, the amount of land devoted to agriculture and arable farming has been significantly reduced, diminishing its role in the balance of natural ecosystems.

**Figure 4.** Land-use categories map and the evolution of main land use over the three different time **Figure 4.** Land-use categories map and the evolution of main land use over the three different time periods.

*3.2. Carbon Dioxide Balance* Information coming from the implementation of the GIS gave the values detailed in Tables 1 and 2 expressed in terms of areas involved by the different vegetation typologies and, applying their individual CO2 sequestration rates, in terms of absolute values of the Indeed, following the application of Commission Regulation (EEC) no 1272/88 of 29 April 1988 (detailed rules for the application of the aid scheme to encourage the setaside of agricultural land) [58] aimed at reducing prices, farmers have been encouraged to abandon the cultivation of cereals and other crops. This process has not only greatly increased but has accelerated in recent decades, due to two factors:

annual sequestration of CO2. The calculation has been performed by aggregating some classes belonging to similar categories: the orchard category, for example, consists of olive groves and vineyards, while the shrubland category includes both pastures and river 1. The depopulation of small villages, with the consequent abandonment of cultivated areas, especially in Basilicata, which made up of fragmented agricultural properties and very small farms [59]

zones. The balance of CO2 does not incorporate the impacts of the agricultural machinery,

periods.

2. The economic crisis of the agricultural sector at the national level, which has led to the abandonment of the less productive areas, which is disadvantageous from an economic point of view.

#### *3.2. Carbon Dioxide Balance*

Information coming from the implementation of the GIS gave the values detailed in Tables 1 and 2 expressed in terms of areas involved by the different vegetation typologies and, applying their individual CO<sup>2</sup> sequestration rates, in terms of absolute values of the annual sequestration of CO2. The calculation has been performed by aggregating some classes belonging to similar categories: the orchard category, for example, consists of olive groves and vineyards, while the shrubland category includes both pastures and river zones. The balance of CO<sup>2</sup> does not incorporate the impacts of the agricultural machinery, supplies, and transportation on CO2: in woodlands, these variables are essentially absent, while, in orchard and arable land, they are heavily influenced by crop practices and in some circumstances are minimal.


**Table 2.** Annual balance of CO<sup>2</sup> in the study area.

From the analysis of the results reported in Tables 1 and 2, it can be deduced that the greatest modifications in land use occurred after the abandonment of large areas that over time have become shrubby, going from 220.6 hectares in 1848 to 290.4 hectares in 2017, accompanied by a strong expansion of urban areas, going from 1.4% to 6.5% of the total surface. The orchard reduction was also significant, with reductions ranging from 22.3% (562.8 hectares) to 4.8% (120.6 hectares) of the overall surface.

There was also a decrease in arable land (from 1447 hectares in 1848 to 828.6 hectares in 2017 with a percentage decrease of 24.6%), while woodland increased from 254.5 hectares in 1848 to 1118.5 in 2017, a percentage increase of 34.27%. As a result of the different performance in terms of CO<sup>2</sup> fixation and relative to the investigated study area, all of these land changes caused a progressive increase in carbon dioxide sequestered by biotic agents embedded in the soil. We can argue that the sequestration of land carbon in 1848 was lower than in more recent periods and that, over time, the land carbon balance has improved, while the heavy emission of greenhouse-effect gases in the atmosphere by urban settlements were at the same time increasingly growing. In fact, the land carbon fixation value has increased over time from 15,718.9 MgC/ha/yr in 1848 to 19,810.7 MgC/ha/yr in 1953 and, finally, to 30,170.2 MgC/ha/yr in 2017. At the same time, massive emissions of greenhouse gases into the atmosphere by metropolitan areas were increasing.

This pattern might be considered a common occurrence also for many other areas of southern Italy or even elsewhere. This method represents a useful tool for rural landscape and environmental planning and management: the research example demonstrated that the careful analysis of land-use changes over the years could allow us to control CO<sup>2</sup> emissions into the atmosphere deriving from the diffusion of human activities. This information may be crucial for setting proper planning policies, in line with the Sustainable Development Goals (SDG), including Climate Action (SDG13), Life on Land (SDG15), etc.

#### **4. Conclusions**

Environmental sustainability must be a primary objective of extra-urban land planning. Sustainable rural development, at least in European countries, has been acknowledged by social awareness and sensibility and is regularly considered by new laws and regulations

aimed at protecting natural resources. In this contest, a precise study of performance differences and global ecosystems monitoring appears to be required in order to propose environmental protection regulations, which are critical elements for the sound planning of extra-urban area and for the civilized world's long-term growth. The alterations in land use—thanks to a comparison between historical cartographic maps of different periods and CO<sup>2</sup> fixation that have occurred in the municipality of Ruoti from 1848 to 2017—thanks to the CO2FIX model—have been shown in this analysis. It has shown how the results of the applied agronomic practices, in terms of CO<sup>2</sup> fixation, can contrast the heavy emissions of greenhouse effect gases into the atmosphere by urban settlements, demonstrating how proper rural-site management can efficiently balance environmental pollution caused by the human development.

This model has been a very useful tool for estimating the amount of carbon fixed by forests and soils, thanks to the relative ease of finding input data for the factors considered and the possibility of initializing the model without the need to have historical data. However, the correspondence between the data observed and that provided by the model is very important. CO2FIX V.2 provides a very accessible tool to estimate the carbon stored by forest stands. It is also very flexible because it can be adapted to different forestmanagement contexts. It presents, in general, some limitations: its input requires data and parameters that must always be updated and are especially peculiar to the study area. Consequently, in output, the values obtained are referable only and exclusively to this area, given the considerable variability of carbon content in forest ecosystems.

As far as the cartographic aspects are concerned, in consideration of the type of cartography and digital image used and of the difference in realization, the study assumes that, especially in the historical map, the precision of the surveys carried out could not guarantee a high accuracy that can be fully compared with modern cartography. But the historical GIS approach allows the limiting of these errors because the use of historical cartography is the only way to spatialize geographically information about land use and coverage before the realization of modern topographic cartography [27].

The application of this approach to the other environmental factors occurring in the total environment—including the atmosphere, lithosphere, hydrosphere, and biosphere, apart from the anthroposphere—would lead to a more comprehensive understanding of landscape development dynamics through its principal environmental components, assisting in the formulation of production-oriented policies capable of compensating for natural balance alterations and for the effective application of sustainable development.

In future studies, the model will be applied to other study areas, characterized by other conditions and management systems. Moreover, we will try to apply it at the landscape level to integrate the parameters within a GIS tool, or we will include other time steps and experiment with new silvicultural systems and forms of treatment that maximize the CO2-fixing capacity of forest and soils. In this perspective, forest owners and policy makers should equip themselves with appropriate planning and certification tools.

**Author Contributions:** The five authors equally shared their contribution to preparing this paper. Canio Manniello and Giuseppe Cillis proposed and developed the research design, methodology, manuscript writing, data analysis, and elaboration. Dina Statuto and Andrea Di Pasquale collected the data; Pietro Picuno supervised the work, provided additional comments on the results and interpretation, and reviewed and approved the final version. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


## *Article* **Historical Vltava River Valley–Various Historical Sources within Web Mapping Environment**

**Jiˇrí Krejˇcí \* and Jiˇrí Cajthaml**

Faculty of Civil Engineering, Czech Technical University in Prague, 166 29 Prague, Czech Republic; jiri.cajthaml@fsv.cvut.cz

**\*** Correspondence: jirikrejci@fsv.cvut.cz; Tel.: +420-224-354-646

**Abstract:** The article deals with a comprehensive information system of the historic Vltava River valley. This system contains a number of resources, which are described. For old maps, which are the basis of the whole system, their georeferencing and potential problems in creating seamless mosaics are described. Other sources of data include old photographs, which are localized and stored in the system, along with the definition point of the place from which they were probably taken. The vectorization of data is described, not only for area features used for the analysis of land-use changes, but also for the vectorization of contours. These were vectorized from old maps and are substantial for the creation of historic DEM. Vectorized footprints of buildings and vectors of other functional areas subsequently serve as a basis for the procedural modeling of the virtual 3D landscape. The creation of such a complex and broad information system cannot be described in one article. The aim of this text is to draw attention to a possible approach to the presentation and visualization of the historic landscape, along with links to important documents.

**Keywords:** information system; Vltava River; old maps; geolocation; photographs; HGIS

#### **1. Introduction**

Information about the historical landscape is a very important source for the analysis of land use and understanding of the issue of changes and development of territories over time. Traditionally, this information is obtained from old maps [1]. In addition, there are numerous of other materials on historical landscapes that can be used. These include old photographs, building plans, inventories of objects, thematic information (water management, social issues, transportation), or books (chronicles, documentary monographs) and historic film footage. All these materials form a memory of the landscape, which, as a whole, provides us with information about the area [2]. In the article, we deal with the creation of a comprehensive information system on the example of the Vltava River, which includes a number of sources. The aim is to connect all the materials with spatial information and to create a historical geographic information system (HGIS) [3]. The system of the former landscape around the river will be in the length of more than 350 km. In this context, the study of landscape or urban change is usually based on the HGIS method [4]. Maps always represent the basis of such an information system. In our case, these are old analogue maps of the river [5]. The oldest topographic maps depicting the river landscape date from the eighteenth century. Thus, our theme is limited to the period from the latter half of the eighteenth century; the greatest emphasis is aimed at the nineteenth and twentieth centuries, and for this time there already exist accurate maps based on geodetic bases.

#### **2. Materials and Methods**

#### *2.1. Theoretical Background*

Information systems focused on historical sites have become quite common. These are usually systems supporting tourism and the presentation of cultural heritage. As an

**Citation:** Krejˇcí, J.; Cajthaml, J. Historical Vltava River Valley–Various Historical Sources within Web Mapping Environment. *ISPRS Int. J. Geo-Inf.* **2022**, *11*, 35. https://doi.org/10.3390/ijgi11010035

Academic Editors: Motti Zohar and Wolfgang Kainz

Received: 31 October 2021 Accepted: 26 December 2021 Published: 4 January 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

example, we can cite the City of Jeddah for the promotion of its history through a webbased information system based on maps [6]. Very often we also encounter databases of historical objects that serve as inventories of cultural heritage [7]. In our project, we also focus on the inventory of interesting historical buildings around the river. These may comprise homesteads, inns, churches or chapels, but also bridges, ferries, weirs and dams. Hydrotechnic objects are thus a very interesting part of the project and represent an industrial heritage that should be included in the system. The database of hydrotechnic objects is also dealt with [8]. Among the infrequent sources focused on the study of historical rivers can be included [9,10]. Some other authors are interested in the river landscape as a specific type of terrain [11,12]. Important and famous objects of a particular period may be highlighted in a contemporary old map along with a popup information window, as in the Map of the City of Madrid from time of Cervantes [13]. A similar approach to ours is taken by the colleagues documenting the City of Granada [14]. They combine old maps with point layers of objects of interest and their photographs.

Old photographs are an integral part of our project. Due to the fact that the Vltava is the most famous Czech river, but also very important from the economic, touristic, transport, communication and geographic standpoint, there are a number of historical photographs documenting its development. For each photo, we create its reference point from where the shot was probably taken. The geolocation of photography can be found in some similar projects [15]. In addition to old photographs, it is advantageous to take current photographs, or rather use photographs of identical places from different time periods [16]. These comparative pairs are then preferably used within the information system.

A specific feature of old riverine maps of watercourses is mainly their linear character. This is also related to the more complex process of georeferencing, as the selected ground control points (GCPs) usually lie in the river line and global transformation equations fail. Then, it is necessary to use local transformation methods based on interpolation mechanisms, such as the Thin Plate Spline [17]. For photographs as well as for maps, it is advisable to georeference representatives from different time periods and allow their comparison within the application. In this case, it is appropriate to vectorize the course of the river so that the data can be easily overlapped [18]. It is advantageous to create vector layers where we deal with Land-Use Changes (LUC). In our project, we decided to create a complete vector data model, which will be used for analysis and visualization. Based on analytical calculations, it is possible to infer the stability of the landscape [19]. From the landscape-hydrology perspective, it is also possible to monitor the development of the pond network interconnected to the river, which has undergone a number of changes in history. In this case, there are not any specific hydrological maps, but usually historical military topographic mappings [20]. Driving forces for land-use and forest changes in the 19th century were studied as well [21,22]. It is also important to obtain information about the altitudinal component of the historic landscape, either from contour lines on the maps or in combination with current models acquired using LiDAR [23]. The DEM derived from them can be a suitable basis for visualizing the landscape in 3D. A way of reconstructing the historical 3D valley of the Vltava River was chosen, mainly because a large part of it was flooded after the construction of a cascade of a total of nine dams during the twentieth century. It is necessary to think about landscape concept and how the visual materials provide information for landscape reconstruction as well [24,25].

Creating virtual landscapes in 3D is currently a very popular visualization tool [26–28]. Extinct places, landscapes and cities can be revived in a virtual environment, as The Ancient City of Rome [29]. In addition to historic DEM, 3D objects (buildings, trees) could be added. The standard for modeling such objects is Trimble SketchUp [30]. By connecting DEM and 3D objects, a virtual landscape is created and can be rendered using a photorealistic visualization (for example the Lumion 3D software). Certainly, it is also possible to use historical aerial photographs and orthophotos. These can be used especially for endangered monuments that have been partly damaged in the last hundred years (for example, the Citadel in the City of Erbil) [31]. We also use the historical orthophoto, because the Vltava underwent

a significant transformation, especially in the twentieth century. Some already extinct monuments can be modeled in great detail using a variety of spatial data. The Roman Circus in Milan was reconstructed as an example of this approach [32]. If monuments are modeled in detail, including the DEM in the vicinity, it is possible to use visualization tools to present objects in different seasons and in different weather conditions, using photorealistic landscape rendering [33]. For a larger area, it is usually not possible to use the methods of classical 3D modeling for all objects in the scene (it may represent thousands of buildings). Then, procedural modeling might be used, which offers significant improvements over a uniform visualization of all elements. Using a set of rules, objects are modeled and textured, involving even an element of randomness. The result is a virtual landscape that can be detailed for only a few key monuments; the other objects are modeled procedurally [34]. Conversely, if we want to model some objects at the level of building elements, it is possible to use the BIM approach [35]. In recent years, a visualization technique using virtual reality (VR) has become a trend. It is used mainly in computer games that use one of the game engines (Unity, Unreal Engine). A typical application area of VR is archaeology. Here, VR can be used as a virtual museum of archaeological artefacts [36]. In the case of landscape modeling in VR [28], the situation is more exacting mainly due to the complicated data transfer between GIS and VR. In our project, the use of VR is planned, but we are only at the beginning of its utilization. However, there are examples of suitable applications, such as, for a landscape that has undergone an extensive transformation due to the mining industry [37].

#### *2.2. Data and Methods*

The Vltava River is an important river crossing most of the territory of Bohemia from south to north. It flows through the Capital City of Prague and is an important national symbol and motif in the works of painters, poets and musicians. The river, which originates in the Šumava forests, flowed through a morphologically diverse landscape, from the gradual Šumava plains, through narrow valleys in the Bohemian Forest and the flat landscape around the South Bohemian metropolis of Cesk ˇ é Budˇejovice to markedly cut valleys closer to the center of Bohemia.

It was also an important transport route, especially for the transport of salt and wood. Equally important was the use of water-power energy and there were many mills, sawmills and other objects on the river, later, for example, hydroelectric power plants. The river was dammed by weirs, including culverts allowing uninterrupted navigation. The traffic significance of the river was the reason for causeway works and caused the creation of a number of organizationally and technically demanding works (for example the canal lock/sluice in Županovice). A towpath led along the river bank, allowing transport boats to move upstream. Although the river itself was an important transport route and economic element, the river's surroundings, especially in hard-to-reach areas, were often very poor. Since the end of the nineteenth century, the charm of the river landscape has enchanted more and more people. Especially near Prague, the river has become a popular destination for tourists and paddlers, which triggered the development of tourism infrastructure.

The importance and ways of the economic use of the river have changed over time. Small weirs have been replaced by larger ones; mills were converted into power plants, and; from the end of the nineteenth century, there were projects for the comprehensive use of water-power and the construction of a dam system. This system of dams was built up during the second half of the twentieth century and it has significantly changed the Vltava riverine landscape and life in it, and interrupted the navigability of the river between Cesk ˇ é Budˇejovice and Prague. The goals of the construction and function of water works also developed over time, including the improvement of navigability, use of water energy, flood protection and providing water for the Temelín nuclear power plant. Works to make the river navigable are again, and currently, in progress.

The studied area represents a historical, cultural, organically developed landscape, the formation of which was significantly influenced by the predominant method of economic use. Due to the dominant form of use, we refer to it as the landscape of dam reservoirs [38]. It is characterized by the flooding of the original settlements connected to the riverbed and the construction of new settlements, usually with a dominant recreational function. The banks of the Vltava have suffered due to the construction of the system of reservoirs—Vltava Cascade. The swelling of the water surface in the monitored part of the river represents more than its two thirds, and it meets the above-described characteristics [39]. The level of the filled reservoirs avoided important historic towns, but rural settlements disappeared underwater. Both in the flat landscape where peat mining, flax growing, logging and floating wood were people's livelihood, and in the valleys and gorges lined with mills, raftsmen' inns, quarries, etc. Facilities for individual and mass recreation have sprung up on the banks. The transport of material on ships and rafts has been replaced by recreational transport [40,41].

Due to its importance, the Vltava River was the subject of interest of many engineers and cartographers. The first maps of all or a significant part of the river were created from the middle of the eighteenth century and mostly depict the most economically important part of the river from Cesk ˇ é Budˇejovice to Prague, or to the confluence with the Elbe River [42]. A number of these valuable and interesting maps was founded across various archival institutions. Most of them were forgotten and almost unknown. These maps vary in the scale, date and method of creation, but all of them represent important images of the Vltava river in the specific period. The overview of so-far researched maps is in Table 1. These manuscript river maps contain information important for navigation (sands, shoals, rocks, islands, weirs, rapids, and streams), crossing the river (fords, ferries, and bridges), the route of the tow-path and the economic background (mills, inns). Some maps also show measured transverse profiles or the velocity of the stream. The topographic content is limited to the river itself and narrow strips around the river, the elevation is shown by hatching (Figure 1). Due to the complexity and inaccessibility of the riverbed, older maps are not very accurate; however, since the first half of the nineteenth century, they have been based on accurate cadastral maps. Therefore, if we are writing about landscape reconstruction, we are usually focused on the first half of the nineteenth century.


**Table 1.** Manuscript river maps of the Vltava River, Source: author's elaboration.

For the description of the whole researched area, there are suitable map works, which cover the whole area uniformly. Due to the fact that the Vltava River and its surroundings belonged to the Habsburg Monarchy in the eighteenth and nineteenth century, it is possible to use the historical mappings of this empire [43,44]. For medium scales, it is possible to use military topographic surveys. An overview is presented in Table 2, including the time of origin of the used map sheets and suitable georeferencing methods.

Since 1918, it has been possible to work with topographic mappings of independent Czechoslovakia. Unfortunately, until 1957, no mapping was carried out to cover the entire territory of the state. For maps containing the surroundings of the Vltava River, military topographic mappings from 1953–1957 at a scale of 1:25,000 and from 1957–1971 at a scale of 1:10,000 are available. These mappings contain, in certain areas, already created dams and reservoirs built most often in the 1950s. For comparison with the current state of the

landscape, it is possible to use contemporary map works, especially the Basic Map of the Czech Republic created from the ZABAGED digital model [45]. *ISPRS Int. J. Geo-Inf.* **2022**, *11*, x FOR PEER REVIEW 5 of 21

> **Figure 1.** Depiction of a manuscript river map derived from the Ebert's map of the Vltava river from eighteenth century, Source: Archive of the Povodí Vltavy, State Enterprise (Prague, Czech Republic). **Figure 1.** Depiction of a manuscript river map derived from the Ebert's map of the Vltava river from eighteenth century, Source: Archive of the Povodí Vltavy, State Enterprise (Prague, Czech Republic).


For the description of the whole researched area, there are suitable map works, which **Table 2.** Topographic surveys of the Habsburg Monarchy in the research area, Source: author's elaboration.

oration. **Mapping Origin Date Scale Suitable Transformation Estimated Accuracy [m]**  First Military Survey 1764–1768 1:28,800 polynomial (control points) 200–300 Second Military Survey 1842–1852 1:28,800 projective (corner points) 20–30 Third Military Survey 1874–1880 1:25,000 projective (corner points) 20–30 Since 1918, it has been possible to work with topographic mappings of independent Czechoslovakia. Unfortunately, until 1957, no mapping was carried out to cover the entire territory of the state. For maps containing the surroundings of the Vltava River, military topographic mappings from 1953–1957 at a scale of 1:25,000 and from 1957–1971 at a scale of 1:10,000 are available. These mappings contain, in certain areas, already created dams and reservoirs built most often in the 1950s. For comparison with the current state of the If we are interested in a larger scale, it is possible to use cadastral maps. From this perspective, the most interesting is the Stable Cadastre, the mapping of which took place on the territory of Bohemia in the years 1826–1843. The scale of cadastral maps is 1:2880, the estimated accuracy is around 2–3 m. These maps were made on the basis of the Imperial Patent of 1817, which established the Stable Cadastre. Cassini-Soldner's projection was used as the map projection; for Bohemia, the Gusterberg datum of cadastre co-ordinates were used. Maps were created by using a method of the plane table in the fathom scale of 1:2880 directly in the field. Imperial Obligatory Imprints (original maps copies) were not used for plotting additional changes in the field, and thus depict the landscape in its original state from the mapping period (Bohemia 1826–1843). Maps were manually colored. Georeferencing these maps is relatively difficult due to mapping one cadastral unit by another on separate map sheets. Of course, it is also possible to use newer cadastral maps of the re-ambulated Stable Cadastre or the Czechoslovak cadaster of real estate (herein, already at a scale of 1:1000).

landscape, it is possible to use contemporary map works, especially the Basic Map of the Czech Republic created from the ZABAGED digital model [45]. If we are interested in a larger scale, it is possible to use cadastral maps. From this perspective, the most interesting is the Stable Cadastre, the mapping of which took place on the territory of Bohemia in the years 1826–1843. The scale of cadastral maps is 1:2880, In terms of elevation, military topographic surveys of the Habsburg Monarchy are inappropriate. First and Second Military Surveys contain only hachures; Third Military Survey contains contour lines, but with a relatively large interval. For the needs of modeling, the historic Vltava valley, only first edition of State Maps derived at a scale of 1:5000 are

the estimated accuracy is around 2–3 m. These maps were made on the basis of the Impe-

was used as the map projection; for Bohemia, the Gusterberg datum of cadastre co-

available. The planimetric map component in the first edition was derived from cadastral maps, while the altimetry component was derived from topographic maps in the S-1952 system, or, in some cases, also from topographic profiles of the 3rd military mapping. The altimetry, which was the subject of the successive processing, was marked in brown comprising contour lines, spot elevations and, potentially, technical or topographic hachures and descriptions. The elevations are in the Baltic Vertical Datum after Adjustment. The map sheets have dimensions of 2.5 × 2 km and are in the 90-degree angle sheet line system projected in the S-JTSK coordinate system. Their first edition is dated 1950–1953, which is a period just before the construction of most dams. Contour lines in these maps have a different quality and interval of contour lines. Therefore, it is not an ideal and homogeneous map work; however, no other source covers the whole area with contour lines.

An important work covering the entire area of interest is also an orthophoto from historical images of military aerial photography from the 1950s (mostly from 1952–1956) (Military Geographical and Hydrometeorological Office, Dobruška, Czech Republic). Aerial photography took place before the construction or filling of the Vltava Cascade (with the exception of Vrané and Štˇechovice Dams built in the 1930s). Therefore, it shows the original landscape and objects around the river. In the case of Slapy and Orlík Dams, it shows the landscape during the dam construction (dam construction, clearing of forests), but before filling the reservoirs.

In addition to these documents, a number of other maps and plans were created, mostly dealing with the regulation and possible use of the Vltava water-power. Most of these water structures were not built due to their complexity or were replaced by a newer project. Many regulatory works and small water structures disappeared or were hidden below the surface of large reservoirs created in the second half of the twentieth century. However, these are interesting constructions and technical solutions that are important to include within the information system. The origin and the owners of these plans and maps are different. They could be state regional archives, the archive of the Povodí Vltavy, the National Technical Museum, regional museums and private collections. Some of these documents are already digitized but, in most cases, the digitization is undertaken within the project.

For use within the river information system, it is important to properly georeference the maps. For newer maps of Czechoslovakia and the Czech Republic, the maps are created in the national JTSK system, and it is possible to transform them only with the use of the map sheets index and projective transformation using corner points. For older maps, it is usually necessary to derive a global transformation key for the transformation into one of the current reference coordinate systems. This was achieved for the Second and Third Military Surveys and the Stable Cadastre. Here, after the transformation of the map sheet index, it is possible to use a projective transformation and four corner points again. The situation is more complicated in the case of the Stable Cadastre, where it is necessary to resolve the contact of individual cadastral areas, as they were mapped separately on separate sheets.

The georeferencing of the First Military Survey maps or manuscript river maps is problematic. Here, it is necessary to use GCPs and a suitable type of transformation, which will ensure the best possible fitting for GCPs and, at the same time, the continuity of adjacent map sheets [46]. For river maps that contain only objects near the river, it is very difficult to find a suitable method. If the sections are drawn with larger errors, large image distortions occur. Many transformations also collapse due to the collinearity of GCPs. The identity of GCPs must be checked carefully. For some maps, it is necessary to accept the non-continuity of map sheets after the resulting transformation.

An excellent source of information about the historic landscape is also old photographs and postcards, or even paintings, drawings and views. In large numbers, they show already extinct or changed places and provide an easy idea of the appearance and changes of the landscape. For the needs of the information system and map applications, all photos are located and represented by one reference point (Figure 2). This allows us to view one site

in different documents, on old maps, current maps, historical and current photos from different places, and attain a comprehensive view of the development of the landscape. *ISPRS Int. J. Geo-Inf.* **2022**, *11*, x FOR PEER REVIEW 8 of 21

**Figure 2.** Visualized database of localized photographs and other places, Source: author's elaboration. **Figure 2.** Visualized database of localized photographs and other places, Source: author's elaboration.

In addition to the possibility of viewing raster maps in the information system, the conversion into a vector data model is suitable for selected layers and elements. At a minimum, it is suitable for the water surface of the river in order to be able to observe and analyze changes in the riverbed. The vectorization of a complete land use in a corridor about six km wide was chosen for cadastral maps. The area up to a distance of three km on each side of the river banks includes the area most affected by land-use changes in connection with the river. It is an area of cadastre units that are adjacent to or affected by the river (for example, flooded by a dam reservoir). Maps from the years 1843 (Stable Cadastre), 1953 (SMO5) and 2021 were chosen as time sections. The oldest period (1843) refers to the oldest possible high-quality, largescale maps. We assume that here we can see the characteristics of the pre-industrial landscape. The 1950s are the second period with an important situation before creating most of the dam reservoirs. Here, we can see the state of the landscape right before its flooding. The third period refers to the current state of the landscape. The current data use the RUIAN database, which contains a drawing of cadastral boundaries. Vector data can be used for the LUC analysis. They also serve as a basis for 3D visualizations, with the foot-There are a large number of photographs and postcards, which puts great demands on their selection and also on the identification of the photographed place. Above all, significant, high-quality and unique, and not commonly known photographs were selected. We went through many important archives, museums or personal collections to attain such photographs (National Archives of the Czech Republic, Prague, Czech Republic, State Regional Archives, Prague, Tˇreboˇn, both Czech Republic, regional museums in Horní Planá, Písek, Cesk ˇ ý Krumlov, Cesk ˇ é Budˇejovice, Pˇríbram, all Czech Republic). The disadvantage is the uneven distribution of photographs. Popular and easily accessible places were photographed often, hard-to-reach places not at all. A good knowledge of the Vltava landscape, as well as other historical and contemporary photographs, is essential for the correct location. The utilization of orthophotos from the 1950s, which capture the yet unflooded landscape and are also close to the time of the creation of the most interesting photographs (first half of the twentieth century), proved to be an essential tool for localization. The determination of the reference points of the photographs was carried out in an online map application using different map layers. Basic metadata are stored as reference point attributes available for each photo. The photographs were searched in the photographic collections of various archives, museums, institutions and private collectors.

prints of buildings being used for procedural modeling [34]. To vectorize the hypsography, we described the methodology of its conversion from SMO5 raster maps to the final 3D model of the whole valley [47]. It is a workflow in which the amount of colors is reduced first; then, automatic vectorization is performed in the ArcScan software. Subsequently, automatic vector adjustments are performed and information on contour heights is manually added. These contour lines, together with eleva-One of the goals of the information system is to show particular places, buildings and structures that have either disappeared or still exist. A database of important objectsconnected with the river and life around the river was established. These objects arewater management facilities (weirs, culverts, dams, limnigraphs), transportation (fords, ferries, bridges), economic structures (mills, sawmills, mills, power plants, paper mills, raft binding yards), sacral buildings (churches, chapels), tourist attractions (inns, restau-

The resulting data (georeferenced maps, vectorized layers, DEM) can be used to

tion spots, form the basis for the creation of the DEM of the historic Vltava valley.

rants and guesthouses), as well as major homesteads and hamlets. The structure of the database was designed with regard to the uniform storage of inconsistent information on different types of objects. The resources for selecting and identifying objects are diverse: literature, old maps, archival sources, existing databases (for example Available online: http://vodnimlyny.cz (accessed on 9 December 2021), and the knowledge of witnesses and researchers as well. The aim of the application is to offer basic information: location, description and purpose of the object, dating, and if there is also a photograph, reference to literature, selected plan documentation, link to archival and online resources for further research, with everything linked to one reference point.

In order to display the detailed situation in a specific locality (mill, paper mill, church), another interesting source are situational and construction floor plans of individual buildings or areas. These are very different in dating, processing, scale, and availability. They are used together with photographs and text descriptions to supplement information about selected important objects within the information system.

In addition to the possibility of viewing raster maps in the information system, the conversion into a vector data model is suitable for selected layers and elements. At a minimum, it is suitable for the water surface of the river in order to be able to observe and analyze changes in the riverbed. The vectorization of a complete land use in a corridor about six km wide was chosen for cadastral maps. The area up to a distance of three km on each side of the river banks includes the area most affected by land-use changes in connection with the river. It is an area of cadastre units that are adjacent to or affected by the river (for example, flooded by a dam reservoir).

Maps from the years 1843 (Stable Cadastre), 1953 (SMO5) and 2021 were chosen as time sections. The oldest period (1843) refers to the oldest possible high-quality, large-scale maps. We assume that here we can see the characteristics of the pre-industrial landscape. The 1950s are the second period with an important situation before creating most of the dam reservoirs. Here, we can see the state of the landscape right before its flooding. The third period refers to the current state of the landscape. The current data use the RUIAN database, which contains a drawing of cadastral boundaries. Vector data can be used for the LUC analysis. They also serve as a basis for 3D visualizations, with the footprints of buildings being used for procedural modeling [34].

To vectorize the hypsography, we described the methodology of its conversion from SMO5 raster maps to the final 3D model of the whole valley [47]. It is a workflow in which the amount of colors is reduced first; then, automatic vectorization is performed in the ArcScan software. Subsequently, automatic vector adjustments are performed and information on contour heights is manually added. These contour lines, together with elevation spots, form the basis for the creation of the DEM of the historic Vltava valley.

The resulting data (georeferenced maps, vectorized layers, DEM) can be used to model the virtual landscape in 3D. As mentioned, the footprints of buildings are used for procedural modeling (Figure 3). Various types of soil cover determine the texture of the surface (or it is possible to use the texture of the raster image of the map). The elevation used from DEM is usually exaggerated (2.5 times) to make the terrain more distinctive. The created virtual model is the basis for a potential creation of physical 3D models (3D printing technology or architectural modeling). The whole model can then be transferred to VR and presented with greater impact on users [37].

to VR and presented with greater impact on users [37].

**Figure 3.** Virtual landscape with procedurally modeled built-up areas, Source: author's elaboration. **Figure 3.** Virtual landscape with procedurally modeled built-up areas, Source: author's elaboration.

surface (or it is possible to use the texture of the raster image of the map). The elevation used from DEM is usually exaggerated (2.5 times) to make the terrain more distinctive. The created virtual model is the basis for a potential creation of physical 3D models (3D printing technology or architectural modeling). The whole model can then be transferred

The aim is to create a comprehensive system that allows us to display and analyze cartographic, photographic, textual and other materials in the context of their geographical location and their relationships. We use the Esri platform for the entire data processing. The georeferencing and vectorization is undertaken in ArcGIS Pro, along with the database creation and their filling as well. The City Engine product is used for procedural modeling. As part of VR testing, we opted for the Unreal Engine platform. ArcGIS Server is used for publishing map layers; the basic composition is created on ArcGIS Online and the application itself uses ArcGIS API for JavaScript 4.17. The aim is to create a comprehensive system that allows us to display and analyze cartographic, photographic, textual and other materials in the context of their geographical location and their relationships. We use the Esri platform for the entire data processing. The georeferencing and vectorization is undertaken in ArcGIS Pro, along with the database creation and their filling as well. The City Engine product is used for procedural modeling. As part of VR testing, we opted for the Unreal Engine platform. ArcGIS Server is usedfor publishing map layers; the basic composition is created on ArcGIS Online and the application itself uses ArcGIS API for JavaScript 4.17.

Within the information system, in addition to the introductory texts, emphasis is placed primarily on the presentation of the results through a combination of a 2D map application and 3D scenes. The 2D map application should enable basic orientation in space, search for localities, the display and comparison of each map layers, as well as the display and viewing of localized photographs and points of interest. These points contain not only basic descriptive information about the object, but, furthermore, they contain links to photographs, plan documentation and other data stored externally (website, archive fund, etc.). Within the information system, in addition to the introductory texts, emphasis isplaced primarily on the presentation of the results through a combination of a 2D map application and 3D scenes. The 2D map application should enable basic orientation in space, search for localities, the display and comparison of each map layers, as well as the display and viewing of localized photographs and points of interest. These points contain not only basic descriptive information about the object, but, furthermore, they contain links to photographs, plan documentation and other data stored externally (website, archive fund, etc.).

The 3D scene combines the created DEM of the historical landscape, vector models of land use in different time sections as possible surface layers, and a procedurally generated model of development (buildings). The procedural modeling of buildings based on their footprints and determined building rules allows a quick creation of a 3D model of thousands of individual buildings and, thus, completes the 3D landscape. We do not create historically accurate models of individual buildings via procedural modeling. The 3D scene combines the created DEM of the historical landscape, vector models of land use in different time sections as possible surface layers, and a procedurally generated model of development (buildings). The procedural modeling of buildings based on their footprints and determined building rules allows a quick creation of a 3D model of thousands of individual buildings and, thus, completes the 3D landscape. We do not create historically accurate models of individual buildings via procedural modeling. However, appropriately selected rules will enable the creation of a model of built-up areas that corresponds in the structure and style. Significant objects can be accurately modeled in CAD. The result is published as a 3WS web scene and made available to the public via ArcGIS Online. The aim of the web scene is to enable a smooth interactive movement and exploration of the now flooded Vltava River landscape in the entire area of interest.

Notably, 2D and 3D applications are complemented by other documents focused on the river landscape and its changes, aspects of life around the river, hydrotechnic objects and human use of the river during the past century. the river landscape and its changes, aspects of life around the river, hydrotechnic objects and human use of the river during the past century.

*ISPRS Int. J. Geo-Inf.* **2022**, *11*, x FOR PEER REVIEW 10 of 21

However, appropriately selected rules will enable the creation of a model of built-up areas that corresponds in the structure and style. Significant objects can be accurately modeled in CAD. The result is published as a 3WS web scene and made available to the public via ArcGIS Online. The aim of the web scene is to enable a smooth interactive movement and exploration of the now flooded Vltava River landscape in the entire area of interest.

Notably, 2D and 3D applications are complemented by other documents focused on

#### **3. Results**

**3. Results** 

The georeferencing of the maps used in the information system was relatively complicated. The maps of military topographic mappings were georeferenced for the whole territory of Bohemia within research on other projects. The mean errors achieved by individual mappings in the standard alignment can be found in the literature [48]. It can be seen that the First Military Survey achieves about ten times greater errors than the Second Military Survey. Therefore, it was necessary to use another proposed method [46]. A total of 250 map sheets and almost 7000 GCPs were used. When adjusted together by a polynomial transformation using the robust IRLS method, a mean error of 280 m was achieved. A preview of the resulting mosaic can be seen in Figure 4. The georeferencing of the maps used in the information system was relatively complicated. The maps of military topographic mappings were georeferenced for the whole territory of Bohemia within research on other projects. The mean errors achieved by individual mappings in the standard alignment can be found in the literature [48]. It can be seen that the First Military Survey achieves about ten times greater errors than the Second Military Survey. Therefore, it was necessary to use another proposed method [46]. A total of 250 map sheets and almost 7000 GCPs were used. When adjusted together by a polynomial transformation using the robust IRLS method, a mean error of 280 m was achieved. A preview of the resulting mosaic can be seen in Figure 4.

**Figure 4.** Complete mosaic of the First Military Survey, Source: author's elaboration. **Figure 4.** Complete mosaic of the First Military Survey, Source: author's elaboration.

Peripheral map sheets remain an unsolved problem, which could not be included in the compensation due to a small amount of GCPs. It will be necessary to add them to the mosaic manually, probably using a polynomial transformation using points on the edges of map sheets. Alternatively, the Grid Shift Binary method could be used for georeferenc-Peripheral map sheets remain an unsolved problem, which could not be included in the compensation due to a small amount of GCPs. It will be necessary to add them to the mosaic manually, probably using a polynomial transformation using points on the edges of map sheets. Alternatively, the Grid Shift Binary method could be used for georeferencing of the First Military Survey, which results in similar accuracy [49].

ing of the First Military Survey, which results in similar accuracy [49]. Newer map works (SMO5, Czechoslovak topographic maps) were georeferenced relatively easily using the map sheet layouts and corner point coordinates. For the SMO5 maps, the resulting mosaic has 334 map sheets. Although a projective transformation using the corners of map sheets was used, the affine transformation was also tested in order to eliminate potential gross errors (for example subtraction of the map corner). The mean errors of these transformations for all map sheets did not exceed five meters.

A total of more than 26,000 km of contours and more than 3000 elevation spots were vectorized using the above mentioned procedure of altimetry vectorization from SMO5 maps. The territory was analyzed to find the used contour interval. Almost 60% of the territory has a contour interval of ten meters, about 24% has an interval of only twenty meters. For ideal modeling, it would be better to have an interval of five meters or less. It is, therefore, necessary to consider that the model created by us is relatively rough. For DEM modeling, the points of the longitudinal profile of the Vltava River from 1940 were also used for refinement. The shoreline is, therefore, very well assigned a height. For broader relationships, it was necessary to link the created DEM with the current quality DMR5G (Digital Terrain Model of the Czech Republic of the 5th generation). It was necessary to smooth the data at the contact between the models so that the transition of both models would be smooth. We analyzed the quality of the DEM created by us using the current flood line of the Slapy, Orlík and Lipno reservoirs. A total of more than 26,000 km of contours and more than 3000 elevation spots were vectorized using the above mentioned procedure of altimetry vectorization from SMO5 maps. The territory was analyzed to find the used contour interval. Almost 60% of the territory has a contour interval of ten meters, about 24% has an interval of only twenty meters. For ideal modeling, it would be better to have an interval of five meters or less. It is, therefore, necessary to consider that the model created by us is relatively rough. For DEM modeling, the points of the longitudinal profile of the Vltava River from 1940 were also used for refinement. The shoreline is, therefore, very well assigned a height. For broader relationships, it was necessary to link the created DEM with the current quality DMR5G (Digital Terrain Model of the Czech Republic of the 5th generation). It was necessary to smooth the data at the contact between the models so that the transition of both models would be smooth. We analyzed the quality of the DEM created by us using the current flood line of the Slapy, Orlík and Lipno reservoirs.

errors of these transformations for all map sheets did not exceed five meters.

Newer map works (SMO5, Czechoslovak topographic maps) were georeferenced relatively easily using the map sheet layouts and corner point coordinates. For the SMO5 maps, the resulting mosaic has 334 map sheets. Although a projective transformation using the corners of map sheets was used, the affine transformation was also tested in order to eliminate potential gross errors (for example subtraction of the map corner). The mean

Figure 5 shows the mismatch of the line generated from the created DEM (blue line) and the current shoreline (red line). Figure 5 shows the mismatch of the line generated from the created DEM (blue line) and the current shoreline (red line).

*ISPRS Int. J. Geo-Inf.* **2022**, *11*, x FOR PEER REVIEW 11 of 21

**Figure 5.** Mismatch between modeled and right shoreline, Source: author's elaboration. **Figure 5.** Mismatch between modeled and right shoreline, Source: author's elaboration.

The georeferencing of the Stable Cadastre maps was the most difficult step in raster data processing. This is mainly because the resulting mosaic contains almost 2000 map sheets. These map sheets had to be georeferenced using both corner points and GCPs. As these are island maps, the connection of individual cadastral areas was a serious issue. These are masked on the map sheets according to cadastral boundaries so that there is no gap or overlap anywhere. Most map sheets were transformed using a projective or affine transformation. A polynomial transformation was used on fewer sheets. The mean errors of these transformations ranged in the order of units of meters, with a mean value of 2–3 The georeferencing of the Stable Cadastre maps was the most difficult step in raster data processing. This is mainly because the resulting mosaic contains almost 2000 map sheets. These map sheets had to be georeferenced using both corner points and GCPs. As these are island maps, the connection of individual cadastral areas was a serious issue. These are masked on the map sheets according to cadastral boundaries so that there is nogap or overlap anywhere. Most map sheets were transformed using a projective or affine transformation. A polynomial transformation was used on fewer sheets. The mean errorsof these transformations ranged in the order of units of meters, with a mean value of 2–3 m. Nevertheless, in some places it was not possible to eliminate completely the discrepancy (conflict) at the contact of the map sheets (Figure 6).

modeling in 3D.

ancy (conflict) at the contact of the map sheets (Figure 6).

**Figure 6.** Visual problems at the borders of cadastral areas, Source: author's elaboration. **Figure 6.** Visual problems at the borders of cadastral areas, Source: author's elaboration.

For our purposes, 3D modeling of buildings can be most easily divided into procedural and detailed in CAD (SketchUp). Accurate modeling of individual buildings brings a detailed graphical representation of an extinct building, but it is very time consuming. It is possible to create models of only carefully selected significant buildings by this method. An important basis for accurate modeling is the planning documentation, construction-historical survey of the building, drawings, and photographs. The model with a Both the mosaic of the Stable Cadastre and the SMO5 mosaic were used for the Land-Use vectorization. At the time of writing, the vectorization of the Stable Cadastre is complete, in the form of approximately 150,000 polygons. The vector data accurately follow the rasters in the form of a data model that uses the map key of the Imperial Obligatory Imprints of the Stable Cadastre. Of course, the vector data have passed through a complete topology check. The vector data will be used for standard LUC comparisons.

m. Nevertheless, in some places it was not possible to eliminate completely the discrep-

topology check. The vector data will be used for standard LUC comparisons.

Both the mosaic of the Stable Cadastre and the SMO5 mosaic were used for the Land-Use vectorization. At the time of writing, the vectorization of the Stable Cadastre is complete, in the form of approximately 150,000 polygons. The vector data accurately follow the rasters in the form of a data model that uses the map key of the Imperial Obligatory Imprints of the Stable Cadastre. Of course, the vector data have passed through a complete

All raster and vector layers will be available to the user within the 2D web map application. This application will be the basis for the information system. In addition to the 2D application, we are working on the development of a 3D application that will allow a better visual perception of the historic landscape. Vector data serve here as a basis for

surrounding terrain can then be modeled using the SketchUp or another CAD software. The Lumion 3D software, for example, can be used as a possible visualization tool, which enables the rendering of a 3D scene and animations in real time. In addition to photorealistic textures, a library of other elements (tree, grass) can be used, as well as the effects of sun exposure, time of day, weather, clouds, etc. The use and combination of effects was All raster and vector layers will be available to the user within the 2D web map application. This application will be the basis for the information system. In addition to the 2D application, we are working on the development of a 3D application that will allow a better visual perception of the historic landscape. Vector data serve here as a basis for modeling in 3D.

tested to achieve the most realistic visualization of an extinct site. The disadvantage is that in a detailed form, it is possible to present only a local model with its immediate surroundings. Due to the laboriousness, only the most valuable objects can be processed in this way. In Figure 7 is an example of the visualization of the parish complex in Červená nad Vltavou. The presentation in an online environment is problematic and the possibilities are limited to sharing a non-interactive video or images. For our purposes, 3D modeling of buildings can be most easily divided into procedural and detailed in CAD (SketchUp). Accurate modeling of individual buildings brings a detailed graphical representation of an extinct building, but it is very time consuming. It is possible to create models of only carefully selected significant buildings by this method. An important basis for accurate modeling is the planning documentation, construction-historical survey of the building, drawings, and photographs. The model with a surrounding terrain can then be modeled using the SketchUp or another CAD software. The Lumion 3D software, for example, can be used as a possible visualization tool, which enables the rendering of a 3D scene and animations in real time. In addition to photorealistic textures, a library of other elements (tree, grass) can be used, as well as the effects of sun exposure, time of day, weather, clouds, etc. The use and combination of effects was tested to achieve the most realistic visualization of an extinct site. The disadvantage is that in a detailed form, it is possible to present only a local model with its immediate surroundings. Due to the laboriousness, only the most valuable objects can be processed in this way. In Figure 7 is an example of the visualization of the parish complex in Cerven ˇ á nad Vltavou. The presentation in an online environment is problematic and the possibilities are limited to sharing a non-interactive video or images.

**Figure 7.** Parish complex in Červená visualized using Lumion 3D, Source: author's elaboration. **Figure 7.** Parish complex in Cerven ˇ á visualized using Lumion 3D, Source: author's elaboration.

The aim is a 3D visualization of the original landscape and buildings along the entire Vltava River, and thus it is necessary to choose a procedural modeling for the modeling of a huge number of buildings. This allows the generation of approximate models of individual buildings according to specified rules. Although the buildings themselves do not exactly correspond to the original, the overall structure style of the development does. Moreover, it is possible to create models of tens of thousands of buildings. It is essential to define the rules according to which individual buildings are drawn and, subsequently, select suitable textures for the visualization of partial parts of buildings. To determine the rules of common development, we use various sources, for example, the literature concerned with rural development in the nineteenth century [50,51]. A very valuable source is the information, drawings and floor plans of buildings that were created during the rescue survey (Figure 8.) carried out by the Institute of Ethnology of the Czech Academy of Sciences before the Orlík Dam had been dammed up. During the survey, all buildings in the area to be flooded were recorded, as well as buildings in the surrounding area. This created a set of data specifying the building style around the river. Due to the fact that a less developed area was involved, a number of buildings have been preserved in a form corresponding to the middle of the nineteenth century. Significant features of buildings were statistically evaluated and rules for building modeling were defined. The aim is a 3D visualization of the original landscape and buildings along the entire Vltava River, and thus it is necessary to choose a procedural modeling for the modeling of a huge number of buildings. This allows the generation of approximate models of individual buildings according to specified rules. Although the buildings themselves do not exactly correspond to the original, the overall structure style of the development does. Moreover, it is possible to create models of tens of thousands of buildings. It is essential to define the rules according to which individual buildings are drawn and, subsequently, select suitable textures for the visualization of partial parts of buildings. To determine the rules of commondevelopment, we use various sources, for example, the literature concerned with rural development in the nineteenth century [50,51]. A very valuable source is the information, drawings and floor plans of buildings that were created during the rescue survey (Figure 8)carried out by the Institute of Ethnology of the Czech Academy of Sciences before the Orlí<sup>k</sup> Dam had been dammed up. During the survey, all buildings in the area to be flooded were recorded, as well as buildings in the surrounding area. This created a set of data specifying the building style around the river. Due to the fact that a less developed area was involved, a number of buildings have been preserved in a form corresponding to the middle of the nineteenth century. Significant features of buildings were statistically evaluated and rules for building modeling were defined.

Existing buildings can also be modeled using photogrammetric or laser scanning and subsequent processing. Existing buildings can also be modeled using photogrammetric or laser scanning and subsequent processing.

Most of the available photographs show well-known and popular places; postcards and photographs are also described and, therefore, at least the approximate location (the nearest village, a significant landscape point) is known. The localization of some photographs, which were taken in little-known and less photographed places and are not described proved to be a problem. These photographs are often the only depiction of the landscape in an individual locality. A careful and laborious comparison of photography, landscape morphology, and old maps will allow at least an approximate localization of the photo. All photographs, as well as map plans from various sources, are provided with a watermark of its owner or administrator or the 'Vltava' project logo for licensing reasons. Obviously, there are tens of thousands of photographs that could be included in

the application. At the moment, we have processed over 1400 photographs, which will be supplemented by others in the coming years. *ISPRS Int. J. Geo-Inf.* **2022**, *11*, x FOR PEER REVIEW 14 of 21

> **Figure 8.** Depiction of rescue survey documentation, Source: Masaryk Institute and Archive of the CAS. **Figure 8.** Depiction of rescue survey documentation, Source: Masaryk Institute and Archive of the CAS.

> Most of the available photographs show well-known and popular places; postcards and photographs are also described and, therefore, at least the approximate location (the nearest village, a significant landscape point) is known. The localization of some photographs, which were taken in little-known and less photographed places and are not de-The database of objects of hydrotechnic, economic, transportation, sacral or touristic character contains over 1200 entries at the moment. As the archival research is still in progress, information and related documents to these objects will be improved. Additionally, the number of entries might slightly increase.

> scribed proved to be a problem. These photographs are often the only depiction of the landscape in an individual locality. A careful and laborious comparison of photography, landscape morphology, and old maps will allow at least an approximate localization of the photo. All photographs, as well as map plans from various sources, are provided with a watermark of its owner or administrator or the 'Vltava' project logo for licensing reasons. Obviously, there are tens of thousands of photographs that could be included in the application. At the moment, we have processed over 1400 photographs, which will be supplemented by others in the coming years. The 2D map application uses and displays most of the data mentioned above. The basic map layers consist of maps of the Stable Cadastre from 1843, SMO5 from 1953 and their vector models. For an easy land-use comparison, a layer from the present RUIAN data is available. An equally important layer is the historical orthophoto from the 1950s and the current orthophoto as well. For a broader view of the region, the layers of three military surveys are attached as well. Manuscript river maps also provide an interesting view of the river, the riverbed and its transformations. Situation and floor plans of important buildings are added as another layer, or they are tied to the important object points.

> The database of objects of hydrotechnic, economic, transportation, sacral or touristic character contains over 1200 entries at the moment. As the archival research is still in progress, information and related documents to these objects will be improved. Additionally, the number of entries might slightly increase. The 2D map application uses and displays most of the data mentioned above. The basic map layers consist of maps of the Stable Cadastre from 1843, SMO5 from 1953 and their vector models. For an easy land-use comparison, a layer from the present RUIAN data is available. An equally important layer is the historical orthophoto from the 1950s and the current orthophoto as well. For a broader view of the region, the layers of three In addition to common map tools, the map application itself is extended using ArcGIS API for JavaScript with more advanced functions. It allows filtering the content of point layers according to the type of object, comparing different map layers using user transparency settings, dividing the map window into two parts with different maps and using the swipe tool (Figure 9). The visualization of data and other image materials is solved with regard to user-friendly accessibility. This, as well as the entire map application, can be tested, evaluated and modified. Another goal is to connect the 2D and 3D content within one application, which will allow a more comprehensive view of the site and a closer understanding of all the context and changes in the transformed landscape.

> military surveys are attached as well. Manuscript river maps also provide an interesting view of the river, the riverbed and its transformations. Situation and floor plans of im-

points.

**Figure 9.** Map application comparing two map layers with the swipe tool, Source: author's elaboration. **Figure 9.** Map application comparing two map layers with the swipe tool, Source: author's elaboration.

In addition to common map tools, the map application itself is extended using ArcGIS API for JavaScript with more advanced functions. It allows filtering the content of point layers according to the type of object, comparing different map layers using user transparency settings, dividing the map window into two parts with different maps and using the swipe tool (Figure 9). The visualization of data and other image materials is solved with regard to user-friendly accessibility. This, as well as the entire map application, can be tested, evaluated and modified. Another goal is to connect the 2D and 3D content within one application, which will allow a more comprehensive view of the site and a closer understanding of all the context and changes in the transformed landscape.

The 3D view of the Vltava Valley is intended for presentation to the general public, and thus different ways of presenting the data were chosen: The 3D view of the Vltava Valley is intended for presentation to the general public, and thus different ways of presenting the data were chosen:


Detailed models of important objects rendered using Lumion 3D rendering software are visualized separately with the immediate surroundings using video animation or images. A great advantage here is the use of a library with a great number of natural and artificial surfaces and objects. Cloud movements, various daylighting, the nature of the weather and visibility, as well as movements of trees and grass when wind is blowing, are applied here. Thanks to this, there are impressive scenes that can bring inexperienced users a feeling of the modeled landscape. Detailed models of important objects rendered using Lumion 3D rendering software are visualized separately with the immediate surroundings using video animation or images. A great advantage here is the use of a library with a great number of natural and artificial surfaces and objects. Cloud movements, various daylighting, the nature of the weather and visibility, as well as movements of trees and grass when wind is blowing, are applied here. Thanks to this, there are impressive scenes that can bring inexperienced users a feeling of the modeled landscape.

Procedural modeling of landscape is used for 3D web map application. Here we use the Esri City Engine product, which allows us to present a full-fledged 3D model in the web environment, of course, with the limitations of the web environment. For example, an application must be loaded into memory at startup. This means a delay of a few seconds before viewing the landscape. Texture handling also has a big impact on the size and speed of the entire application. It is possible to use rather colored surfaces with simple effects (water ripples), but it is not possible to achieve photorealistic fidelity as in the case Procedural modeling of landscape is used for 3D web map application. Here we use the Esri City Engine product, which allows us to present a full-fledged 3D model in the web environment, of course, with the limitations of the web environment. For example,an application must be loaded into memory at startup. This means a delay of a few seconds before viewing the landscape. Texture handling also has a big impact on the size and speed of the entire application. It is possible to use rather colored surfaces with simple effects (water ripples), but it is not possible to achieve photorealistic fidelity as in the case of Lumion 3D software. As for the detail of individual objects, it is given by the rules for their generation. For buildings, for example, it is possible to set rules for the representation of roof types, building heights, textures of their facades, etc. When modeling, it is possible to set a degree of uncertainty with which the objects are modeled. Then, based on the rules, the buildings are created as 3D objects with the appropriate textures. There are more than 20,000 buildings generated in our scene. Other elements of the landscape are visualized using colorful textures of land use. Only some important objects (castles, chateaux, churches) are inserted in the form of a detailed model created in Trimble SketchUp.

The planned output also includes three large-format physical models of the landscape around large reservoirs, popular with tourists, aimed to present the extinct landscape in local places (regional museum, information center). The models have a scale of 1:6000 (or 1:8000, respectively), consisting of several parts and show an area of six by twenty four km (or ten by thirty-two km, respectively). Analogically to the virtual model, heights are 2.5 times exaggerated. Due to their dimensions, the models are divided into parts. The terrain is milled from permanent hardened polystyrene. In critical areas (river bank, urban area), it can be modeled by hand. The original method of the translation of the topographic content to the model by applying and adhering the surface texture printed on the elastic foil proved to be non-functional due to the high fragmentation and altitude differences of the terrain. For this reason, a more laborious method of manually redrawing the surface according to the image projected on the model was chosen. Detailed areas in the urban part of municipalities are solved by a sticker with a printed texture. The original river surface is shown in deep blue. The new water level of dam reservoirs is represented by a lightly tinted transparent modeling material, so that the flooded terrain underneath remains visible. The mass is self-leveling, and during pouring, it is necessary to proceed very carefully and within the limits of a truly flooded valley. Significant objects and localities can be easily identified on the large models using interactively controlled diodes. The models are currently in production.

#### **4. Discussion**

In the article, we presented our vision of a comprehensive information system of the historic landscape, specifically the Vltava River valley. Commonly presented information systems serve mainly to promote tourism [6]. In our work, we try to comprehensively process various materials so that they can be used not only by the public, but also by other experts (geographers, historians). It turns out that there are not too many information systems focused on historical river flows. It is possible to find more applications focused on databases of specific objects [8]. We focused not only on hydrotechnic objects, but also on other important objects of economic and social importance, for which we make resources and information available for further research. Our work is based on the concept of HGIS not only in classical 2D, but also in 3D [3]. We can state that the approach in which the position of the object is the basic identifier is still not a standard. In contrast to our approach, it is more often possible to find the classic approach of an information system based on a table-oriented database, where maps and locations are only additional data with links to external map applications. We consider our method more synoptic and approachable. Various materials, while geolocated, can be presented together in the map application giving a deep view into archival and current material for the particular location (area) or period. Metadata and attributes of these materials are still table-oriented. Therefore, the searching and filtering in data is easy to contrive. The processing and geolocation of archival material can be undertaken by cooperating experts in proper fields. Our approach is complex and thereby very time-consuming.

The georeferencing of maps and geolocation of photographs, our main workflow, is limited to raster data and 2D applications [14]. It is possible to observe problematic georeferencing of riverine maps, where only the close vicinity of the watercourse is depicted. For this reason, in our approach, we combine global transformation methods and local methods based on interpolation algorithms. Our approach to georeferencing military survey maps and Stable Cadastre, where we use a joint adjustment of map sheets with robust IRLS method, is quite unique as well. Of course, it is necessary to distinguish between the achieved accuracy of georeferencing and the accuracy of used old maps. The resulting mean transformation error gives us information about the estimated accuracy. This information should be kept in mind when further interpreting the information displayed on the maps. Therefore, the oldest maps with higher errors are used only in the raster form and are not suitable for vectorising data for further geographical analysis.

The processing and localization of old photographs is used in a number of projects [16]. Unlike them, we try to capture the probable place of taking a photograph on the map, rather than the place of interest in the photograph, and thus enable a more illustrative approach. The photographs are not available in the form of a classic web gallery, but they are available directly at the places in the GIS where the objects of interest are located. The link to external sources of information is then possible thanks to the unambiguous position of the object. Using old maps, their data model, important objects, building plans and photographs, we try to restore the valuable extinct landscape and present it to the public in a 2D map application and in a 3D environment.

When processing 3D scenes, we try not only for classical modeling of 3D objects [29]; we are focused rather on a procedural approach [34]. This makes it possible to create large-scale models in the order of tens of kilometers, even with buildings. Compared to single monument online visualization a web scene with thousands of buildings using photorealistic textures of objects and surfaces seems to be slightly problematic, so far [52]. It will be important to focus on the generalization and the use of different levels of detail (LOD) [53]. For visualization, we will try to work with photorealistic scenes [33]. The disadvantage of this procedure is the creation of only fly-through videos and the inability of the user to browse the entire model. Browsing in our model is enabled not only by 3D applications, but also by the intended VR application, as is the trend today [37].

The 2D web map application, 3D scene and thus the entire information system aims to summarize, synthesize, publish and attractively visualize a wide range of information, old maps, photographs, archival documents and research results in one place. The focus of the information system is very wide, and the main common denominator is the Vltava River.

The project itself is in its middle phase. Used methods seem to be working well for our needs. The majority of important data have already been processed. Databases of objects and photographs are sufficiently filled by hundreds of entries. Moreover, archival research is still in progress, and there will be more photographs or additional data to process; even any forgotten river map might be discovered. In the future, we will focus on the analysis of land use, the addition of more maps, the refilling of data entries of individual objects, the development and evaluation of map applications and, last but not least, on the methods of 3D visualization and use of VR.

There exist some other projects focused on the virtual reconstruction of the landscape [26,33]. Their approach is usually slightly different. They are creating 2D or 3D virtual environments for landscape reconstruction, especially for visualization purposes. We are using maps as a medium for other material connections. We would like to utilize maps as a center point of the information system. This is the uniqueness of the project.

#### **5. Conclusions**

As written above, the Vltava River is an important river for many reasons, outlined above. The motivation for our project was the missing complex information system about the history of the river. There are many visual materials such as maps, plans and photographs that are not processed nor used within combined information systems. The idea of our work was to digitize and process the most interesting items and bring them together into the digital environment. Providing easier access to the information in the future was our main driving force.

Our work documents the processing of various map, and photographic and descriptive materials, which are the basis for the future information system about the historic valley of the Vltava River. The aim was to present particular data components and the methodology of their processing so that they can be used for the information system. The system is based on the spatial determination of all available data, or their link to a spatial point object. Map as well as image data are processed on the Esri platform. In addition to the classical methods of georeferencing, more complex methods were used, especially the polynomial transformation with the conditions of the linkup of map sheets, supplemented by the robust IRLS statistical method. The proposed work procedure was used to vectorize the hypsography, which leads to a high-quality vector model of contour lines.

In addition to data processing methods, the text also provides an overview of the possibilities for visualizing the historic landscape. It is obvious that traditional 2D applications must be the basis of the whole system. It is appropriate to complement this application with a 3D visualization, as a virtual historical landscape. It can be viewed either in a web browser or as a VR application. If there is a suitable exhibition space, physical models created by 3D printing or architectural modeling methods can be used.

It is very important to supplement the information system with localized photographs of the historic landscape. The photos suitably complement the maps and map outputs, which are not intuitive for every user. Maps and photographs are suitably complemented by plans of important buildings, text documents and links to external sources.

The comprehensive information system revives the extinct historical Vltava landscape, presents valuable old maps and photographs and enables their mutual comparison. The database of important objects is a detailed list of cultural, hydrotechnic, technical and economic objects. Each record contains position information, basic descriptive information and other additional data.

The aim of the application is to bring closer the historical landscape, which has largely disappeared as a result of the construction of a cascade of a total of nine dams. Our approach is very interdisciplinary. We try to connect humanities (history, archiving) with a technical approach (georeferencing of maps, web map applications, 3D modeling). We think that our system will increase awareness of the entire area under study, and that it will serve as an entry point for research into the extinct surroundings of the Vltava River.

We see the novelty of our approaches mainly in complex methods of georeferencing old maps, in the method of using reference points of photographs, in the use of procedural modeling of the historical landscape for its visualization in 3D, and in plans to use VR as an access point for a virtual walk through the landscape. The scope of our work is very wide, so the text is a selection of descriptions of the methods and techniques used and an overview of some specific implementations.

**Author Contributions:** Conceptualization, Jiˇrí Krejˇcí and Jiˇrí Cajthaml; methodology, Jiˇrí Krejˇcí and Jiˇrí Cajthaml; software, Jiˇrí Krejˇcí; writing—original draft preparation, Jiˇrí Krejˇcí and Jiˇrí Cajthaml; writing—review and editing, Jiˇrí Krejˇcí and Jiˇrí Cajthaml; visualization, Jiˇrí Cajthaml; project administration, Jiˇrí Cajthaml. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by the Grant Agency of the Czech Technical University in Prague, grant SGS22 "Modern cartographic methods for processing, analysis, visualization and presentation of spatial data in 2D and 3D".

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


## *Article* **GisGCN: A Visual Graph-Based Framework to Match Geographical Areas through Time**

**Margarita Khokhlova 1,2 , Nathalie Abadie 1,\* , Valérie Gouet-Brunet <sup>1</sup> and Liming Chen <sup>2</sup>**


**\*** Correspondence: nathalie-f.abadie@ign.fr

**Abstract:** Historical visual sources are particularly useful for reconstructing the successive states of the territory in the past and for analysing its evolution. However, finding visual sources covering a given area within a large mass of archives can be very difficult if they are poorly documented. In the case of aerial photographs, most of the time, this task is carried out by solely relying on the visual content of the images. Convolutional Neural Networks are capable to capture the visual cues of the images and match them to each other given a sufficient amount of training data. However, over time and across seasons, the natural and man-made landscapes may evolve, making historical imagebased retrieval a challenging task. We want to approach this cross-time aerial indexing and retrieval problem from a different novel point of view: by using geometrical and topological properties of geographic entities of the researched zone encoded as graph representations which are more robust to appearance changes than the pure image-based ones. Geographic entities in the vertical aerial images are thought of as nodes in a graph, linked to each other by edges representing their spatial relationships. To build such graphs, we propose to use instances from topographic vector databases and state-of-the-art spatial analysis methods. We demonstrate how these geospatial graphs can be successfully matched across time by means of the learned graph embedding.

**Keywords:** historical visual sources; graph embeddings; geospatial descriptors; indexing and retrieval of historical data

#### **1. Introduction**

Historical visual sources, such as maps, engravings, drawings, or photographs, are the only visual representations of the past geographical realm still accessible. As such, they are extremely important sources of information for analyzing the past territory and its evolution. In recent years, many historical Geographic Information System (GIS) initiatives have thus been developed, essentially using maps or, more rarely, old photographs to reconstruct the past territory.

A first challenge faced by these works is the georeferencing of the historical visual sources. This operation aims at determining the geographic position of a given visual source. It can be performed using either a physical sensor model based on precise physical and geometrical information describing the sensor used for data capture or a correspondence model derived from a set of ground control points for which ground and image coordinates are known [1]. Whether performed manually [2] or automatically [3], this task requires providing enough ground control points to the system that estimates the georeferencing model. This thus implies identifying beforehand the area covered by each historical visual source, which is, in the case of old aerial photographs, sometimes very poorly documented, can be extremely difficult.

A second challenge lies in the extraction of useful information from the visual sources and their integration within a GIS data structure to make them usable. As an example, ancient place names can be extracted from old maps and stored with their respective

**Citation:** Khokhlova, M.; Abadie, N.; Gouet-Brunet, V.; Chen, L. GisGCN: A Visual Graph-Based Framework to Match Geographical Areas through Time. *ISPRS Int. J. Geo-Inf.* **2022**, *11*, 97. https://doi.org/10.3390/ ijgi11020097

Academic Editor: Motti Zohar

Received: 29 October 2021 Accepted: 24 January 2022 Published: 29 January 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

locations in some GIS data format. This operation can be performed manually, most often through collaborative approaches (see for example, https://geo.nls.uk/maps/gb1900/ (accessed on 28 October 2021) or https://geohistoricaldata.org/) (accessed on 28 October 2021), or semi-automatically, using collaborative correction or validation tools (such as http: //buildinginspector.nypl.org/ (accessed on 28 October 2021)), applied on data produced by automatic image segmentation and vectorization approaches [4,5].

The work presented in this article addresses the aforementioned challenge, that is to automatically identify and index the geographical area covered by some given historical visual source. We especially focus on old aerial photographs, but the approach could also be extended to other visual sources such as old maps. We see this work as the step towards the unified framework for analysis of various multidimensional cross-temporal sources.

Existing approaches designed to match an aerial photograph usually then compare it to a database of georeferenced photographs [6] using visual cues i.e., color, texture or luminance descriptors. It is, however, not a trivial task to match the aerial photographs available across years due to the structural and appearance changes occurring through time. New buildings and roads appear, old ones are demolished, rivers, and watercourses can change their beds, etc. Nevertheless, intuitively, these man-made and natural objects are more stable than pure visual cues, which strongly vary with the season and weather conditions on the acquisition day.

While our approach targets historical aerial image matching with georeferenced photographs, we propose a fundamentally different and novel way to encode the information captured by photographs: we propose to leverage the geometrical (normalized perimeter and eccentricity, shape, concativity, etc.) and topological properties of geographic entities represented by the photographs instead of pure visual cues. Indeed, we think that using the properties and spatial relations of geographic entities could outperform matching approaches based on visual cues, because they change less through time than the visual appearance of photographs. To use these properties of geographical entities, a first step is to extract them from photographs, using automatic vectorization tools for instance. However, to test our hypothesis, we propose to first dispense with this step and to use instances of vector geographic databases instead. Indeed, modern vector topographic databases are commonly captured from aerial photographs and are thus geometrically consistent with them. Hence, inspired by the recent progress in learning on graph data, we propose to model the scene depicting a geographic area as a graph, which is a very powerful and intuitive way to encode geoinformation. In this graph, geographic entities are represented by nodes, their properties are added to the nodes as labels, and their spatial relationships are represented by edges between the respective nodes. For example, an aerial photograph can be converted into an undirected graph, where rivers, buildings and roads are the nodes and edges represent their neighbourhood relationships. Graph representation is very promising because it can capture and standardize the semantics and relational information from various sources of data such as photographs, historical maps, or textual descriptions.

The contributions of this work are the following. First of all, we propose a novel approach for cross-time geographic area search, formulating the problem as a graph similarity matching task. Second, we create a multisource cross-time graph database which is the first benchmark of such a kind available for other researches. Third, we evaluate the classical key-points based search on our new dataset to establish a baseline. Fourth, we propose an approach to learn the graph similarity using a Graph Neural Network based Siamese-like architecture.The tests we performed show that the proposed method can be successfully used for cross-time data matching and that it is a promising avenue as there is still a lot of possible improvement for the more challenging cross-year cross-source retrieval task.

The remainder of this article is structured as follows. First, we present the problem definition and basic notions that motivates our proposal in Section 2. We give the background on the relevant research areas in Section 7. This section completes the research paper and gives the overview of the methods, on which this work is based; however, it is not necessary for understanding of the proposed method. We outline the proposed framework which

uses a Graph Convolutional Network (GCN) as backbone and a siamese structure to learn the graph embeddings in Section 3. We introduce a new across-time graph dataset retrieved from the historical photographs and the intuition behind the conversion of visual data to structured one in Section 4. In the next section we describe the experiments performed to evaluate the proposed approach. The last section, Section 8, outlines the future works.

#### **2. Problem Statement**

In our project, we aim at developing a unified framework for the analysis and indexing of various types of visual sources, seeking the optimal common representation of different features coming from archive photographs, historical maps, modern vectorized maps etc. This article focuses on a subtask in the scope of the project: the cross-time matching of geographical areas, based on the similarity of the geographic entities that are located in them. This subtask is particularly difficult as the landscape may have changed over time. To evaluate our approach, we present in this article a use case scenario of cross-time matching of aerial vertical photographs, based on the similarity of the geographic entities they represent.

Aerial image retrieval is particularly difficult, as almost all images represent the same type of geographic entities: buildings, road, watercourses, fields, forest, etc. At that resolution, the main characteristic that distinguishes one geographic area from the others is the spatial distribution of the geographic entities it contains. Moreover, the majority of rural areas contain no distinguishable objects at all, being occupied by the forests or agricultural fields. The three closest directions to the problem at hand are visual landmark cross-view research, retrieval analysis of visual information in various historical-based studies, and more generally Remote Sensing (RS).

The first group of methods commonly requires the objects present in the image to be captured with a higher resolution and contain specific landmarks. Their features can then be learned to be matched across views by means of a Convolutional Neural Network (CNN) as it was done in [7,8]. However, archive aerial images often do not have a significant level of details and high resolution and cannot be matched using these methods. Some recent works proposed to learn the relationship between ground level appearance and overhead appearance and land cover attributes from sparsely available geotagged ground-level images [9–11]. The methods have demonstrated their potential for the problem of determining the real-world geographic location across visual sources. However, they do not particularly target across-time scenarios and perform poorly on unseen locations. Geospatial zone visual research and retrieval in the context of aerial images was analyzed to a much lesser extent, with some examples of existing software such as [6] where aerial image databases are extensively searched and matched to geolocate a parcel of land data in the USA. The second group of methods covers conceptually different methods, but it is always assumed that the images are at least partially georeferenced as in the photogrammetry-based analysis of historical aerial photography [12,13]. The analysis of the changes is commonly based on visual characteristic dynamics. For example, ref. [12] quantify volumetric changes along sandy beaches from archival imagery using photogrammetry-based analysis. In this scenario, the temporal changes are the ones which should be tracked and quantitatively or qualitatively analyzed using historical data. Remote-Sensing gathered lots of attention of researchers recently, thanks to the availability of accurate data of Earth observations. While the approaches and applications vary significantly, the common characteristic of RS is the large volume, large variety, large velocity, large veracity and large value, which raises awareness about the importance of large-scale image processing, fusion, mining and indexing and retrieval. Hence, the indexation of the images plays a important role to solve the RS problems [14,15]. Historical data possess similar processing challenges and will benefit from a set of appropriate indexing and retrieval tools.

Figure 1 shows the unified framework for analysis of various multidimensional crosstemporal visual sources proposed in this work. As shown this figure, our global approach starts with vectorization step of the historical visual sources. This step can be performed

either manually, or with an automatic semantic segmentation method. It aims at extracting the main geographic entities represented in the visual sources, their characteristics—such as shape, size, orientation, nature, etc., and their spatial relations. Inspired by the recent advances of machine learning graph-based approaches and the achievements of CNNbased visual matching algorithms, we propose to use a graph representation derived from this semantic segmentation process results to match aerial images across time. Similarly, some visual-based methods such as [9] rely on a prior segmentation, which is made by the means of CNNs in a separate step. The segmentation usually allows to detect important and significant objects present on an image. While CNNs are showing state-of-the-art results in segmentation, this step is still error-prone. Thus it remains a question whether the possible errors in this step are impacting a lot results of the following cross view matching task. In this work, we exclude the segmentation task to concentrate on the data representation and matching steps, hence we choose to use manually annotated dataset. This allows us to exclude the potential source of errors due to the segmentation and evaluate the feasibility of the matching part of the algorithm in the controlled setup. Therefore, in the remainder of this article we will focus on the steps 1 to 4 of the GisGCN approach.

**Figure 1.** The proposed approach GisGCN consists of several steps: (1) selection of the POI and surrounding geographical zones for cross-time data; (2) for each geographic zone, representation of its geographic entities as a connected labelled graph; (3) for each graph, GCN-based embedding learning with shared weights; (4) similarity-based retrieval of geographic areas across time.

In step 1, Points Of Interest (POI) are sampled using the characteristic of the geographic entities present on the historical visual sources. In step 2, the geographic entities in the area of a predefined size around a POI are used to form the nodes of a graph *G*. The edges are built based on the spatial relationships between the geographic entities and nodes are labelled with semantic and geometric attributes of their respective geographic entity. Given the set of graphs G={*G*0, *G*1, . . . , *GT*} representing the territory around a POI through years, our purpose is to retrieve nearly identical geographic areas across time for the query graph by employing a similarity measure. Step 3 thus aims at learning an embedding vector *D* of the geographical graph structure, which can take into account the node attributes and topological relations between the nodes, be robust to the noise and changes in the data coming from different sources or different dates, and be compact. In step 4, these embeddings are used for fast search and retrieval in a large database of thousands of graphs. Figure 2 shows the training and evaluation scheme used in our article. We train on the dataset of known geographic area correspondences using contrastive learning approach. At the evaluation step, the pre-trained network is used to obtain a descriptor from a query graph, which is then matched to a database by the means of K nearest neighbors (*K*NN) search.

**Figure 2.** The main steps of the proposed approach GisGCN on graphs: training and evaluation.

We called the proposed approach GisGCN, where Gis stands for geographic information science which deals, among other things, with geographic information representation and analysis, and GCN is simply the term used for a neural network that operates on graphs.

#### **3. GisGCN: A GCN-Based Siamese Model for Geographic Area Embedding Learning**

We suggest that the topological information encoded in a graph structure can be essential when we need to distinguish between two different geographic areas with a similar set of objects and similar geometric attributes. Therefore, we propose a novel learning pipeline using a GCN [16] network to learn to match geographic areas represented as graphs across time. In the original work, GCN model was proposed to perform a node classification task for big sparse graphs. Our model aims to learn an embedding space for variable size geographic graphs by exploring the notion of deep graph matching.

#### *3.1. Building Graphs to Represent Geographic Entities and Their Spatial Distribution*

For each geographic zone of interest, the graph representing its spatial configuration can thus be summed up by the following equation:

$$\mathbf{GC}\_{\mathcal{C}} = \left(\mathbb{R}\_{\mathcal{C}}, \boldsymbol{X}\_{\prime}\boldsymbol{A}\right) \tag{1}$$

where:


Such a graph can be built from the geographic vector data captured in pre-requirement step of the pipeline proposed in Figure 1. POI can be chosen according to user's criteria, based on geographic entities properties. The size of the geographic areas built around the POI may vary depending on the scale of analysis of each use case. Node features are attributes describing geographic entities properties. They may be provided by the input vector dataset, if available, or computed from the geometric properties of the geographic entities by means of classical spatial analysis methods [17]. Edges are derived from spatial relationships between geographic entities: they may be simple neighbourhood relationships, computed by means of a Delaunay Triangulation [18] for example, or topological relationships implemented in the graph as labelled edges.

#### *3.2. Model Architecture*

Given two graphs *G*<sup>1</sup> = (*X*1, *A*1) and *G*<sup>2</sup> = (*X*2, *A*2), where *A* is a connectivity matrix and *X* is node parameters, we want a model that produces the learning function *f* : *G* → *D* through the GCN with learnable parameters *w*, to compare them next with the similarity score *s*(*D*1, *D*2) between them in a new vector space. The encoding function *f* takes the *A* and *X* values of a current POI and all geographic entities within the reference area *R* as input and outputs the embedded geospatial contextual information. Our model allows to convert graphs to descriptors, which enables efficient retrieval with fast nearest neighbor search data structures.

In this work, we propose to adapt Siamese networks to handle graphs to learn their embeddings. Indeed, a Siamese network is particularly effective in training a model with few examples for each class, which is most often the case for historical visual sources. The Siamese network consists of two identical networks (with shareable weight parameters). In our case, each of the networks is essentially a GCN with maxpooling, as depicted in Figure 3. This graph matching embedding model is inspired by GNN [16], and comprises four main parts:


The architecture is schematically shown in Figure 3. The aggregation layers follow the formulation of the GCN by [16] and are defined as:

$$X^{L+1} = \sigma(AX^L \mathcal{W}) \tag{2}$$

where *A* is the normalized and modified as in [16] adjacency matrix and *W* are weights to be trained and *σ* is a ReLu activation function. As in the original work by [16] we use two layers of propagation, so that the representation for each node will accumulate information in its local 2-hop neighborhood.

**Figure 3.** Schematic architecture of the model proposed to train the graph embeddings.

After we obtain the final node representations, we aggregate across them to get graphlevel representations. This is implemented by a simple maxpooling followed by a MLP operation that reduces the node representation into a single vector and then transforms it:

$$D = MLP\_G(\operatorname{max} pooling\_{x \in n}(X\_i^l))\tag{3}$$

where *X* are the learned graph nodes *n* features.

The proposed architecture mainly differs from [16] in the point (3) where we do not calculate the node-level features, but compute a graph-level representation instead by performing a maxpooling operation over the nodes in a graph to obtain the whole graph descriptors *D<sup>G</sup>* similar to [19]. The pooling layer maps the input graph of any structure and size to a fixed size-structured output.

During training, the embedding model will jointly reason about the graph structure as well as the graph features to come up with an embedding that reflects the notion of similarity described by the training examples. The proposed Siamese GCN model is endowed with the contrastive loss to train on the data with the ground truth correspondence. The NT- Xent [20] loss function for a positive pair of examples of matching geographic areas through time (*i*, *j*) is defined as:

$$l = -\log \frac{\exp(\text{sim}(D\_1, D\_2)/\tau)}{\sum\_{k=1}^{2N} \mathbf{1}\_{k \neq i} \exp(\text{sim}(D\_1, D\_2)/\tau)}\tag{4}$$

where: *τ* is the temperature, sim(*D<sup>i</sup>* ,*D<sup>j</sup>* )—cosine similarity, *i*, *j*—two graphs in batch of size *N*. The final loss is computed as an arithmetic mean across all positive pairs, both (*i*, *j*) and (*j*, *i*), in a mini-batch.

$$L = \frac{1}{2N} [l(2k - 1, 2k) + l(2k, 2k - 1)]\tag{5}$$

Following the idea of [20], we propose to create batches of random graphs to train the model. However, instead of altering them to use as an input for the second branch of the Siamese GCN, we take the graphs representing the same geographic area but from a different time frame to form positive samples. Then the loss encourages the embeddings for the same geographic area to be closer in the embedding space in terms of the cosine distance; and the embeddings of different areas to be farther apart.

#### *3.3. Similarity-Based Retrieval of Geographic Areas*

Once the model has been trained, it can be used to produce further graph embeddings for areas of unknown location for which there is a historical visual representation at a certain date. The resulting embeddings can then be compared to those pre-computed during model training using state-of-the-art measures such as cosine similarity or L2 distance. If the searched area is represented in our training database, then its embedding should obtain a higher similarity score with those of the graphs representing the same area at different dates, thus allowing the location of the searched area. Ideally, it should obtain a higher score with the temporally closest embedding, provided that the available historical representations have the same level of detail, thus allowing to estimate the valid time of the visual historical source from which it is extracted.

#### **4. Dataset Preparation**

Our target use case deals with ancient aerial photographs picturing French landscape evolution through time: we therefore need annotations for each set of photographs representing the same area at different dates. Besides, aerial image retrieval is particularly difficult, as almost all images represent the same type of geographic entities: buildings, road, watercourses, fields, forest, etc. At that resolution, the main characteristic that distinguishes one geographic area from the others is the spatial distribution of the geographic entities it contains. To maximize our chances, we thus focus on the geographic areas where points of interest are located.

#### *4.1. Input Data Selection*

The input dataset in our experiments is made of temporal snapshots of a topographic vector database which precisely correspond to vertical aerial images taken from three French departments (the so-called "departments" are administrative divisions of the French territory): Moselle, Bas-Rhin, and Meurthe-et-Moselle in four different years (2004, 2010, 2014 and 2019). The data used to create our dataset originate from the French Mapping Agency (IGN) [21]. They have been annotated manually, following the same data capture rules and in an incremental way through time. They are thus very homogeneous, and especially very geometrically consistent through time (See Figure 4). We have selected the following types of geographical entities for their visual salience in the landscape and their durability: roads, railways, waterways, buildings, airports, sports facilities, and cemeteries. For each department and each temporal snapshot, we leverage the geometric representations of geographic entities as well as their associated attributes information to build a graph describing their overall spatial distribution on the territory.

**Figure 4.** Temporal snapshots overlapped with an orthophotograph from 2020. Buildings captured at different dates are shown with green and white strokes for 2004 and rose and beige solid color for 2019. Unless the terrain changes, they are perfectly identical from one database version to the next.

To test the potential for generalization of our method, we need more heterogeneity between geometries, as if they had been produced by some automatic segmentation algorithm. Therefore, we also added data for the same geographic entity types, taken from Open Street Map [22] (OSM) for two of the departments (Moselle and Meurthe-et-Moselle). OSM data are produced by crowd-sourcing, without precise data capture rules, and therefore, they are more likely to have heterogeneous geometries. They are only available in the most recent edition dated 2020. In addition to these cross-source data, we also experiment with randomly generated noise added to geometric attributes across the single source dataset.

#### *4.2. Defining Geographic Zones of Interest*

To ensure that the areas represented by our graphs contain enough geographical entities to identify them, we build our graphs around POIs: building of religious nature, historical objects and monuments, castles or forts, local governmental buildings, buildings with the sport functions, railway stations, airports, etc. We take the exact geometric centers of POIs, randomly adding a small shift to its easting and northing coordinates up to 10 m. An example of the geographic zones selected for the database is shown in Figure 5, where the area bounding boxes of dimensions 200 × 200 m are superimposed on the aerial footage images.

For each vector database edition (namely, 2004, 2010, 2014 or 2019), we randomly shift the center of the bounding boxes up to 20% so that the represented area is not the same one for a more realistic retrieval scenario. Some of the selected POIs happened to be close to each other, so the bounding boxes around them overlap. If the overlapping area is more than 50%, we consider the resulting graphs as representing the same geographic zone and assign this zone the same unique identifier. This results in a small number of zones which have nonunique correspondences in the database (see the Table 1 for the details). We reassign the labels after the graph creation using the Union Find algorithm [23]. An example of nonunique zone is shown in Figure 5: see the two overlapping bottom-left bounding boxes.

For each of these zones of interest, the semantic information available in the topographic vector data and their geometries are used to create a relational graph *GCe*(*R<sup>e</sup>* , *X*, *A*) with attributes *X*, as described below.

#### *4.3. Building the Graph*

For each POI *e*, we build a graph describing its geographical context *GC<sup>e</sup>* , i.e., the geographic entities located in its surrounding geographic area and theirs spatial relationships.

#### 4.3.1. Creating Nodes and Edges

Each node *V* of the graph represents a single geographical entity. In our experiments, we choose the following geographic entity types as nodes for the graph: buildings, roads

topological segments, watercourses topological segments, railroad topological segments, airports, sport facilities, cemeteries. We choose them as they are quite perennial, clearly visible, and structurally important landmarks in vertical aerial images.

**Figure 5.** An example of the bounding boxes of an equal surface around POIs. Note that if two bounding boxes overlap more than 50%—we consider them representing the same geographic zone.

We consider the geographic coordinates of the nodes as a finite planar set to compute the neighbourhood relationship between them and build the graph edges. Toussaint et al. [24] compared three methods to connect and cluster a set of points. Based on their analysis, we selected Delaunay triangulation [18] between the centroids of nodes geometry, as the method provides the biggest amount of edges between locally neighboring points. We also experimented with the Range Neighboring Graphs (RNG) and provide the code for RNG graph creation along with the Delaunay one. Both methods are selected to guarantee the formation of connected graphs.

#### 4.3.2. Adding Semantic Labels to Nodes

Graph nodes have discrete labels representing their nature: river, road, railroad, religious buildings, castle, fort, tower, arc, monument, cemetery, sports ground, normal building, public buildings, and airport. In the case where some precise label is not available, which is sometimes the case for buildings, the label 'normal building' is applied.

#### 4.3.3. Adding Geometric Attributes to Nodes and Edges

Each node also has geometric attributes *X* (normalized perimeter and eccentricity). Many other geometrical attributes are commonly used in spatial analysis, such as general orientation, mean axes of geometric forms, surface descriptors, and various shape descriptors [17]. However, we limited this research to the simplest ones which do not require any orientation information nor high level of details to compensate for the different level of details in the annotations across time and data sources and to avoid making assumptions on the orientation of the scene. Eccentricity in our case is simply *E* = *<sup>L</sup> <sup>W</sup>* , where *L* and *W* are the length and width of the geometry's minimal bounding box. Normalized perimeter is simply *P<sup>n</sup>* = *<sup>P</sup> <sup>H</sup>*×*<sup>W</sup>* .

We also store the edge attributes in the form of the normalized distance between the nodes. The other option is to use the angles but that will make this feature dependent on the orientation and we seek our graphs to be rotation and scale-invariant: hence, the angles between objects cannot be used as edge weights in our scenario. Note that we do not use the edges later on in our tests but leave them for future research, limiting this work to the most challenging case scenario. The graphs in our dataset are therefore undirected and unweighted.

The code designed for graph creation from the vector shapefile data is available on an open repository available here: https://www.alegoria-project.fr/en/GENR\_dataset (accessed on 28 October 2021).

#### *4.4. Discussion on the Resulting Graphs Dataset*

An example of the resulting graphs representing the same area across time is shown in Figure 6. Note the difference in the corresponding landscape and graph structure within the two corresponding dates with a 15 years gap. We observe that in some cases the changes are very significant. It is even more the case for cross-source graphs, even when the time gap is small. Two examples of the resulting cross-source graphs are shown in Figures 7 and 8. Visually, the biggest changes are in the categories of objects such as roads and watercourses, and there are less changes in the buildings, although there are nodes of the category "building" which appear and disappear across time.

Apart from using matching geographical areas from the three departments mentioned before, we also add some clutter data without correspondences from a fourth French department, namely, Côtes-d'Armor, to make the indexing and retrieval scenario for the same source data more challenging. Table 1 presents the data in the graph database we developed for geographic area retrieval across time. Note that the maximum number of nodes in a single graph is 150. This is done intentionally, we just removed several geographic zones with a bigger number of vertices to limit the final graph size. Similarly, we removed all graphs with less than 3 nodes. This choice is explained by the selection of the GCN model later in this work.

**Figure 6.** An example of the resulting graphs representing the same geographical area at two different time points. (**Left**): 2004. (**Right**): 2019. The detailed geometries are shown for the reference, their categories are color-coded: red = building, orange = road segment, yellow = building from a special category (such as church, monument, castle, etc).

**Figure 7.** Resulting graphs for graph #17. (**Left**): OSM 2020, (**Middle**): IGN 2019. Node labels are color coded: red = buildings, orange: roads, violet: rail. (**Right**): superimposed IGN (orange) and OSM (blue) graphs.

#### *4.5. Graph Dataset Statistics*

An example of the final graph data characteristics and distribution for Moselle department is summarized in Table 2. It is interesting to see the difference between the data distribution across the years. Note also that there is a particular change in the number of nodes and edges in 2014—this is probably due to some change in the manual vector data capture process since we use the same code to convert the vector data to graph representations for all years and both sources. The graph characteristic distributions for different databases with a small gap seem to be rather similar, with less nodes for OSM graphs than for IGN graphs, which goes along with our observation that IGN data have a higher level of detail on the geographic extent of our dataset.

**Figure 8.** (**Left**): Superimposed graphs OSM 2020 (blue) and IGN 2004 (orange). (**Right**): superimposed graphs OSM 2020 and IGN 2019. Same geographical area across years and sources.

**Table 1.** The characteristics of the proposed dataset for cross-time geographic area retrieval. Note that we mostly deal with single geographic area correspondence across years.


**Table 2.** Single department statistics example.


As explained above, our dataset does not always contain the same geographic areas per each year, i.e., the bounding boxes are shifted for a more realistic scenario. We thus provide the statistical analysis of the similarity of the attributes between matching graph data. To that end, Intersection Over Union (IOU) between two graphs representing the same geographic area across time is used:

$$IOU\_{G\_1, G\_2} = \frac{2 \ast \sumXe\_{G\_1} == Xe\_{G\_2}}{\sum e\_{G\_1} + e\_{G\_2}} \tag{6}$$

where *X*(*X*1, *X*2) are geometric attributes of the nodes representing geographical entities *e*. The resulting distributions of the IOU between matching graphs coming from the same data source are shown in Figure 9. The obtained distributions are not always normal, however; analyzing the resulting histograms, we can see that the number of graphs with smaller values of IOU is bigger when the difference in time is bigger. Since the landscape and thus the geographical entities it contains are more likely to have evolved over a long period of time than over a short interval, it seems logical that the differences between graphs of very different valid times are greater than those between graphs close in time.

**Figure 9.** IOU histograms for matching graph geographic areas obtained for Meurthe-en-Moselle department. If all the nodes attributes match exactly (injective matching) across graphs, the IOU value will be 1.

Figure 10 demonstrates the cross-source statistical analysis for IOU over node attributes. While calculating IOU for cross-source graphs, we introduced a new parameter *d*, which corresponds to the number of digits after the decimal point in the continuous attribute values. We observed that varying the *d* value, we obtain rather different IOU histograms, as shown in Figure 11a. Using the histograms and based on the experimental results from the next section, we choose *d* = 3 for the cross-source data experiments. Moreover, with a visual examination and further histogram analysis, we discovered that the buildings have the most persistent geometric features across the two databases, and the rivers the least. This is probably due to the fact that the same cadastral plan was used to annotate buildings for OSM and IGN data sources. The histograms in Figure 11b summarise the IOU for the attributes of different categories. Still, across time, the average IOU score for cross-source data stays much lower than the same source one.

**Figure 10.** IOU histograms for matching graph geographic areas obtained for Meurthe-et-Moselle department for two graph sources: IGN 2019 and OSM 2020. The precision level of attributes is fixed at 3.

**Figure 11.** IOU (**a**) histograms between IGN and OSM data 2019–2020, with varying *d*. (**b**) IOU histograms for matching graph geographic areas obtained for Meurthe-en-Moselle department for two graph sources (IGN 2019-OSM 2020) and 4 different categories: red—ordinary buildings, yellow—all buildings (including monuments), blue—rivers and black—roads.

We have thus created a cross-source cross-time dataset containing graph representations of matching geographic areas. For the generalisability, the data are taken from two different sources: Open Street Map [22] (OSM) and French Mapping Agency (IGN) [21]. The selection of these particular dates for our experiments was dictated by the availability of the aerial photographs and their corresponding vector data. We confirmed statistically the assumption that the bigger is the time gap, the more different the resulting graphs for the same territory are. We used as much intermediate dates as possible to evaluate how the accuracy of our approach is dependent on the number of human-made landscape changes through years. These final cross-time cross-database graphs representing nearly identical geographical zones are then used to assess the performance of the proposed GisGCN method. Overall, we expect our approach to work for the case when there is a significant number of objects preserved on given geographic zone across years, and we are aware that the performance can degrade when time frame increases.

#### **5. Experiments and Model Evaluation**

In this section, we present two ways of comparing the data produced from different editions of the BD TOPO<sup>r</sup> database to find those representing the same area of the geographical realm. The first, unsupervised, simply compares the attribute values computed as presented in Section 3. The second uses graph embeddings produced using the network presented in the previous section.

#### *5.1. Attributes Based Nonsupervised Similarity Search*

We use a nonsupervised baseline to determine whether the geometric attributes of the scene objects alone are enough for the graph matching task, without using any topological information (i.e., graph representation) and any learning. We employ the *K* Nearest Neighbors (*K*NN) to retrieve top 5 matching results and report map average precision value (map@5):

$$map\circledast \mathbf{K} = \frac{\sum\_{n=1}^{N} P\_{av} \circledast \mathbf{K}}{N} \tag{7}$$

where *N* is number of queries, *Pav* is the average precision for a single query, *K* = 5.

We use Facebook AI Similarity Search (Faiss) library [25] to retrieve the geographic areas across time. Faiss is designed to search for multimedia documents that are similar to each other using the *K*NN algorithm. We use the *L*2 distance measure to retrieve the most similar geographic areas across years based on the local geometric features and semantics of each object present in the scene. Other similarity metrics available in Faiss proved to be less efficient experimentally.

The results of the Faiss-based similarity search are summarised in Table 3. The obtained map@5 scores are quite high, which means that the semantic and geometrical attributes are representative enough to describe the geographical areas. Nevertheless, there is still room for improvement that can potentially be gained by using the topological information between the nodes, i.e the graph representation. It is interesting to see that when the returned data are actually wrong. Note that even if the actual correspondence geographic area from 2004 contains many entities present in 2019, the other areas were returned instead. This example shows the limitation of the attributes-only search, when no relational information about the scene is used. In the same moment, we can see that there are also significant graph structure modifications across years, which we ideally want to be robust to in the graph matching scenario. An example of the correctly returned data is shown in Figure 12.

Overall, all our results show that the smaller the time gap between the query and the database is, the better are the map@5 retrieval results. This seems logical since the time frame difference should be connected to the number of landscape changes occurred. There is a sudden drop in the retrieval accuracy between the cross-sources data with a 10 years gap.

#### *5.2. Node Attributes Robustness in the Presence of the Noise*

In our so-called *same source* database, we have very precise node attributes with the precision of six decimals that are similar thanks to access to the manually and incrementally captured vector data. It is especially the case for the tests with data coming from the same source, and the corresponding map@5 results are higher than the cross-source map@5. If the same information should have been extracted from the images automatically, the node attribute similarity would have dropped due to the errors in the segmentation and vectorization stages. To simulate this realistic scenario, we made a group of tests by decreasing the number of decimals in the node attributes up to one, two, and three decimals, and by adding normally distributed noise with *µ* = 0 and *σ* = [0.1, 0.01, 0.001] to the queries attributes.


**Table 3.** map@5 Faiss similarity search results, same source database with clutter graphs.

**Figure 12.** An example of the correct top 5 similar graphs in 2010 returned for the query from 2020, cross source database. The node colors represent the semantics of geographic entities.

Faiss method results with added Gaussian noise and less precise query features are summarised in Table 4.


**Table 4.** No learning retrieval baseline using Faiss similarity search and modified node attributes.

#### *5.3. GisGCN Model*

We consider the two following scenarios presented in Table 5 for our graph similarity learning problem:



**Table 5.** Two scenarios used in our experiments.

Both scenarios are deployed for the *same source* and *cross source* data. In addition, we try to directly transfer the learned model to the new data source to further establish its performance on this data source. We compute the similarity between the final descriptors using a simple L2 similarity metric in the vector space from Faiss library, and evaluate the results using map@K metric. Note that the similarity in the loss function we use is cosine similarity; however, we found that L2 distance worked better at the inference time. We also report the average retrieval time for a single query.

#### *5.4. Model Parameters*

Throughout the experiments, we fixed the dimensionality of graph embeddings to 512, trying the following commonly used values: 128, 256, 512. Our tests on "same source" data have shown that when the final descriptor size is lower than 512, the learning capacities of the model stay the same, but the generalization capacities for the cross-time matching on new time frames are much lower. We can obtain the same map@5 values on the training set with all descriptor sizes. However, the validation map@5 reaches the plateau and stops increasing for lower map@5 values in case of the smaller descriptor size. We observed that the loss function follows a similar trend, and the smaller the descriptor size, the faster the model starts to overfit the training data. The first FC layer has a constant size and maps the features to 128 dimension space.

We use two modified GCN layers in the backbone of our model. The weights in the GCN layers are equal to 512 as well, although we experimented with different parameters from 128 and up to 1024. We observed that the smaller number of weights leads to underfitting on our data. The dimensionality of FC layer added to GCN layers is also equal to 512.

#### *5.5. Training Parameters*

In the graph preprocessing step, the discrete labels of the nodes were one-hot encoded, and the continuous attributes were left intact and concatenated with the one-hot encoded ones. Since the graphs have different number of nodes, we use padding to create graphs of the same size as input to the network. The maximum number of nodes in a single graph was thus set to 150. An example of the graphs for the same zone across years is shown in Figure 13. The first fully connected layer with learning weights of the model maps the features to 128 dimensions.

The training graph pairs are selected during the run time and are shuffled randomly at the end of each epoch. The network is trained for approximately 200 epochs until the moment the validation map@5 accuracy reaches the plateau and starts to decrease. In the case of a single department training, 100 epochs were sufficient to rich the maximum validation accuracy. Further training leads to overfitting to the training set, so we take the model which shows the best validation score and report its results. We used the batch size of 64 graphs through all experiments.

We do not use any data augmentation techniques and NT-Xent loss allows to avoid hard sample mining commonly used for Siamese network training. The temperature parameter in the contrastive loss is equal to 0.5. Adam [26] optimizer is chosen for the optimization of the learning weights. The learning rate is equal to 0.1<sup>3</sup> , with a decay 0.1<sup>5</sup> and a multistep learning rate scheduler, decreasing after 80 and 120 epochs with gamma 0.05. All training is performed on a CPU and does not require excessive calculating power, mainly thanks to the small size of our graphs. It takes about 12 h to train the model end-to-end on our dataset with 5732 graphs across two distinct years from the database, along with calculating map@5 precision for training and validation data after each epoch. A single branch of the Siamese model is used in the inference step.

**Figure 13.** An example of the graph representation of the nearly identical geographical zone across 4 years.

#### **6. Discussion**

Table 6 summarises the results for the cross-time learning scenario for our GCN-based descriptors. We report the performance of the global and local descriptors, taking the resulting node embeddings before the maxpooling layer. We observed that the latter takes much longer to compute and perform less well than the global ones, which corresponds to our learning objective. On average, the retrieval time for the query represented as a single descriptor is twice smaller than the one obtained with Faiss local features earlier on the same computer. The obtained map@5 values are lower than the ones obtained with a purely feature-based similarity search from our baseline.

Next, we evaluated if the *same source* trained network can be used for *cross source* retrieval scenario by querying the OSM data with the IGN data pretrained model. The obtained map@5 values are quite low, although they correspond to the baseline results in the case of a big time gap between the two sources. Therefore, we can conclude that the training procedure is under specified and the model does not generalize well to the new data. This result leaves a lot of margins for further improvement, even in case of data with a small time gap between them.

Table 7 shows how the noise added to the query graph attributes affects the results at the inference stage. Here we can see that, in contrast to the baseline method, GCN-based descriptor is relatively robust to noise, with map@5 values decreasing up to 10% in contrast with the 70% decrease of Faiss noisy local feature search. We only added noise to a single source data, assuming that the cross-source data already contain a significant amount of noise due to their capture process.

**Global Local Query Year db Year Role map@5 t per q map@5 t per q** same source 2019 2010 training 0.660 0.017 0.303 0.880 2019 2014 validation 0.554 0.017 0.317 0.412 2019 2004 testing 0.5522 0.017 0.371 1.14 2010 2004 testing 0.684 0.019 0.638 0.97 cross source 2020 2004 testing 0.056 0.01 n/a n/a 2020 2010 testing 0.062 0.01 n/a n/a 2020 2014 testing 0.044 0.01 n/a n/a 2020 2019 testing 0.085 0.01 n/a n/a

**Table 6.** The map@5 results for the global and local descriptors. Note that training and validation same-source training data do not contain clutter, but the testing data does. Cross-source test is performed without clutter data. Time is indicated in s.

**Table 7.** The map@5 results for the single source database global and local descriptors with noise *σ* = 0.01.


Table 8 shows the potential for generalization of our GCN-based descriptors for the case when data come from the same source and are trained and tested on different geographical areas. Note that the map@5 is reported for the similarity search in two (training) and one (testing) departments correspondingly. We kept the hyperparameters tuned for the cross-time learning scenario previously, and can see that the network is capable to create a meaningful descriptor for the new region it did not see during training, so it generalizes rather well when the data come from the same source.

**Table 8.** The map@5 results for the global and local descriptors, cross-department learning.


Finally, Table 9 shows the *cross source* cross-time learning results. The model was trained from scratch, but we left all parameters intact from the previous tests. The results are significantly lower than those obtained by the *same source* scenario, which demonstrates the difficulty of the aimed task, especially for data with a significant time gap. It is hard to say what is the reason of the decreased performance in this case, since both the geometric attributes and graph structure can be affected a lot in the case of the *cross source* data. Nevertheless, taking into account the fact that our dataset ground truth consists of mostly

single correct matches, the results are interesting and confirm the intuition that the graph representations can be used for this challenging task.


**Table 9.** The map@5 results for the global and local descriptors, cross-year correspondence learning, cross-source data.

Table 10 shows the generalizability of our GCN-based descriptors for the case when data are coming from two different sources. Note that the map@5 is for the similarity search in one (training) and one (testing) departments correspondingly. We got a lower map@5 score than the one for the same source cross-year data but a higher result.

**Table 10.** The map@5 results for the global and local descriptors, cross-department learning and testing, cross-source data.


The obtained results show that the proposed method works and outperforms the baseline approach for scenarios dealing with less precise attributes and *cross source* data: hence it has more potential for the real-case geographical graph matching task we aim at. Table 7 proves that the resulting descriptors are robust to noise and retrieve the correct geographic area in 50% of the queries when dealing with same source data (i.e., trained and tested on the same geographic area). Note that with the unsupervised baselines, noise is dramatically affecting the map@5 results while GCN model proves robust. In a real-world scenario, where it is yet impossible to fully automatically produce perfect segmentation results, and even human annotation might vary from person to person and from database to database, this property is crucial. However, the *cross source* tests have shown quite low map@5 values, but demonstrated that our approach is not generalizable enough at the current state. We explored different ways to represent geographical data as graphs (Delaunay graphs, RNG graphs), and experimented with different parameters of our Siamese GCN model and discovered that the graph architecture and network parameters are very dependent on each other but we did not discover a particular trend characteristic for all tests performed. Saying that, we assume that although the first results obtained are not yet on the point with classical methods, it is a path to a new research direction.

#### **7. Background**

The idea of using distinct geographic entities to automatically index, match—and optionally geolocalize—aerial images are not novel. Early research in Computer Vision envisaged it long time ago [27], and there are recent works for RS image-based indexation and retrieval [28,29]. However, the majority of the works dedicated to image matching and retrieval do not explicitly use the graph structure to represent the spatial relations between geographic entities and do not aim for cross-time matching as we do in this research inspired by the recent advances in structural learning. In the following we give the background of the various research areas, on which the current work is based.

#### *7.1. Representing Geographic Information as Graph*

Graph-based representations are traditionally used for road, train, or watercourse network representation, and their associated computations such as routing applications. They are also commonly used for spatio-temporal geographic data [30]. But graph-based representation of places and landscapes can also reveal important insights in scenarios such as scene geolocalization or geographic information retrieval.

Initially coined by Google to describe an enhancement of its search engine with semantics [31], the term "knowledge graph" is widely used today to refer to any graph-based representation of general purpose knowledge, such as the big knowledge bases of the Linked Open Data (LOD) cloud, namely DBPedia [32], Yago [33], etc. From the early days of the Web of data, geographic data have played a central role in the LOD cloud, as they provide an intuitive way to link datasets from different fields such as life sciences, humanities, heritage, media, social networks, etc. [34]. Following the Ordnance Survey example [35], many national mapping agencies have published their geographic data in compliance to Web of data good practices and standards [36–39]. In order to go a step further than directly translating geographic vector data into a RDF graph, [40] proposes a set of metrically refined approximate topological relations to enrich a geographic knowledge graph and improve its question answering capabilities. Geographic knowledge graph summarization is largely covered in [41]. The work covers topics such as understanding, representing, and reasoning about POI but does not propose ways to learn geographic area descriptors. Finally, Trisedya et al. [42] propose an entity alignment model for knowledge graphs based on the popular earlier Trans-E approach [43] which models relationships by interpreting them as translations operating on low-dimensional embeddings of the entities. However, the method is only applicable for 1 to 1 linking between two graph databases and uses triples instead of graphs to learn the embeddings, so it cannot be directly used in our scenario. Finally, Ling Cai et al. [44] propose a unified GCN encoder framework, named TransGCN, which can learn entity embeddings and relation embeddings simultaneously. The work is similar to ours in the idea to use the graph embedding to encode the neighboring structure; however, the authors target node embedding and not a whole graph embeddings as we do.

#### *7.2. Graph Kernels and Graph Distances*

Similarity search based on graph kernels is a well-known research subject with a great number of various kernels proposed for various specific cases and data types, many are available in Grakel library [45]. The similarity is typically defined by either exact matches (full-graph or subgraph isomorphism) [46], random walks or paths on graphs [47], propagation of the information in the graph structure [48] or others. A recent survey on graph kernels can be found in [49]. It should be noticed that the kernels themselves are hand-designed and motivated by graph theory, and only some of them are designed to handle continuous attributes on edges and nodes of a graph. Graph kernels can be formulated as first computing the feature vectors for each graph, and then taking the inner product between these vectors to compute the kernel value, no learning involved with the exception of [50]. Graph kernels have shown themselves as very efficient tools for graph comparison, but often take a significant time to compute.

Hand-engineered or learned graph distances are very similar to graph kernels. Common choices include spectral distances and distances based on node affinity. [51] compares commonly used graph metrics and distance measures, and demonstrates their ability to discern between common topological features found in both random graph models and real world networks. Many of the classical graph kernels are also based on the graph distances [45,52]. Recently, graph distances caught researchers attention, with recent works exploiting attention mechanisms to make learnable metrics [53].

#### *7.3. Graph Embeddings*

Recently different Machine Learning and Deep Learning algorithms were proposed for graph data. The data mining community has a strong interest in (knowledge) graph summarization because graph structure is ubiquitous: all kinds of data from social networks and up to research publications can be represented as graphs. A popular idea is to learn the embeddings for nodes [54] or even the whole graph [55,56] based on their features and structure. There is also an example of the use of such a node embedding method in a geographical context. Yan et al. [57] propose to estimate the similarity and relatedness of place types with their surroundings using place types embeddings. However, all these algorithms are based on the models coming from text processing, so they were designed to produce embeddings such that nodes with similar network neighbourhoods are embedded close together: the nodes are handled as words taken from a vocabulary. Moreover, these methods can handle structure or label information but not both at the same time, which limits their application for our scenario.

#### *7.4. Convolutional Graph Networks*

In the past few years, graph neural networks (GNNs) have emerged as an effective class of models for learning representations of structured data and for solving various supervised prediction problems on graphs. Such models are invariant by design to permutations of graph elements and compute graph node representations through a propagation process which iteratively aggregates local structural information. Nodes on isomorphic graphs (with the same node and edge features) will have the same representations regardless of the ordering. GNN networks have different architectures and can be roughly classified into several categories: spectral methods [58] perform graph convolution by employing the eigen vectors of the graph Laplacian as the transformation matrix, methods that work in the spatial domain [16] and methods complementary to GNNs and agnostic to the choice of a GNN itself (i.e., pooling, attention) [59]. GNNs have been successfully used in many domains from drug discovery [60] to social network classifications [61]. Independently of the network nature, the common task accomplished by them is the supervised learning of the node embeddings. These node representations are then either used directly for node classification, or pooled into a graph vector for graph classification. Problems beyond supervised classification or regression are relatively less well studied for GNNs. Xu et al. [62] proved that with Convolutional Neural Networks, we can measure the graph's similarity such as Weisfeiler-Lehman similarity test. However, the graph isomorphism problem is not very relevant for our use case. More general graph similarity learning approaches were recently proposed by [63,64]. These learned models can adapt to the desired metric and are potentially interesting for our target scenario. However, Li et al. [63] demonstrate the performance of the method on graphs with only minor changes and no node attributes. In this paper, we focus on representation and similarity metric learning for attributed graphs representing the same geographic area across time.

#### *7.5. Siamese Networks*

Siamese network architectures aim to construct embeddings, where two extracted features corresponding to the same real world entity are more likely to be similar than features representing different real world entities [65]. They are a popular choice for scenarios dealing with so-called one-shot learning problems, when a single training sample is available for each class. The efficiency of Siamese networks was previously demonstrated for visual object tracking [66], person re-identification [67], cross-view image matching [68] and other tasks. Siamese networks can also be used for graph similarity learning as demonstrated in [63]. The closest to ours is the recent work of [69] where the authors successfully use a Siamese Graph Convolutional architecture for indexing and retrieving remote sensing images represented as region adjacency graphs.

We follow the similar idea to use the descriptive power of graph representation along with a Siamese-based GCN; however, the graph creation process differs from region adjacency graph (RAG) approach [69], and the architecture we propose is conceived for our type of data and corresponding features. The final scenario also differs: we want to retrieve the exact location and not the similar classes, hence we deal with a more challenging problem with many classes and a few examples per class (mostly, a single correspondence). Moreover, our end goal is to make an image-to vector data correspondence.

#### **8. Conclusions**

With the growing availability of large volumes of historical and modern geographic data of different modalities (image, structured, or textual data, etc.), the development of a unified framework for their joint analysis is of great interest. In this work, we started to move towards such a framework by exploring models for geographic area similarity learning, which is itself a relevant research direction.

This article proposes an approach to the problem of vertical cross-time image indexing and retrieval from a new point of view: we interpret the geographic entities represented in photographs and their geometric and semantic properties as a connected graph. We then proposed a novel deep learning-based method to learn graph representations of geographic entities and their spatial relationships and compare them across time. To test this approach, we created two original *same source* and *cross source* datasets. The proposed GCNbased model is currently outperformed by the unsupervised method based on attributes similarity we used as a baseline. But in contrast to this baseline, our model is robust to the presence of noise in the attributes, which makes it credible in a real-world scenario, such as indexing and retrieval of automatically segmented and vectorized aerial images or even an heterogeneous vector databases matching task or a spatio-temporal pattern retrieval application. The integration of GIS data and Machine Learning allowed us to successfully match the geographical areas across time, obtaining a correct match in more than 50% of cases, up until across 15 years. Besides, the proposed approach can be directly used to learn graph embeddings in any attributed graph similarity problem.

We see the proposed approach to be used by historians and archivists working with large amounts of historical visual sources and looking for ways to retrieve the images photographs, maps, engravings, etc.—representing the same geographical zones, or to locate these images when poor location indication is given. However, using our approach is not straightforward as aerial photographs or ancient maps are rarely segmented, and the proposed GisGCN method requires the segmentation as the prior step. It can be done manually for one image at the query time provided that the database is coming from geographical vector data as it is done in this work. Exploiting image and vectorial data together is more complicated in the case of matching across oblique aerial images sets for example. We see two possibilities here. The automatic segmentation CNN-based methods are gaining more and more accuracy, and they can be used to annotate the data as a prior step. Another approach would be using the proposed graph-based learning methods, but getting the graph nodes by the means of object detection. Further, some additional nodes can be introduced capturing the textual annotations, partly present on many historical data. However, it will bring a new difficulty of encoding textual features such as toponyms, a challenge we leave for future research. Another interesting direct application of our work is the search for particular spatial configurations (buildings wedged between a road and a river, for example) regardless of location, useful for professions that study land use, space organization and urbanization.

Lastly, there are still a number of interesting challenges to resolve: to improve on the efficiency of the matching models, to study different matching architectures, to adapt the GCN so it can use graphs of different sizes, and to apply our model to new application domains. We think that the graph models with attention can work efficiently in the aimed application, and we plan to adopt an attention mechanism in the future. Another possible direction is to improve the graph representation, which could lead to better retrieval results. We leave these directions for future research. We also provide the new *cross source* and *crosstime* dataset along with dataloaders (including an adaption of the dataloader for Pytorch

Geometric library [70]). We hope that this work can spur further research in geographical graph matching and provide the first benchmark for learning on graphs for cross-time area matching.

**Author Contributions:** Conceptualization, Nathalie Abadie, Valérie Gouet-Brunet and Margarita Khokhlova; methodology Nathalie Abadie and Margarita Khokhlova; software Margarita Khokhlova; validation Nathalie Abadie and Valérie Gouet-Brunet; formal analysis Valérie Gouet-Brunet, Nathalie Abadie and Margarita Khokhlova; investigation, Nathalie Abadie and Margarita Khokhlova; resources, Nathalie Abadie; data curation, Nathalie Abadie; writing—original draft preparation, Margarita Khokhlova; writing—review and editing, Nathalie Abadie and Valérie Gouet-Brunet; visualization, Margarita Khokhlova; supervision, Nathalie Abadie, Valérie Gouet-Brunet and Liming Chen; project administration, Valérie Gouet-Brunet; funding acquisition, Valérie Gouet-Brunet. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work is supported by ANR, the French National Research Agency, within the ALEGORIA project (https://www.alegoria-project.fr (accessed on 28 October 2021)), under Grant ANR-17-CE38-0014-01.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The dataset used in this work can be found at [71].

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


## *Article* **Impact of the Cartographer's Position and Topographic Accessibility on the Accuracy of Historical Land Use Information: Case of the Second Military Survey Maps of the Habsburg Empire**

**Krzysztof Ostafin \* , Małgorzata Pietrzak and Dominik Kaim**

Institute of Geography and Spatial Management, Faculty of Geography and Geology, Jagiellonian University, Gronostajowa 7, 30-387 Kraków, Poland; malgorzata.pietrzak@uj.edu.pl (M.P.); dominik.kaim@uj.edu.pl (D.K.) **\*** Correspondence: krzysztof.ostafin@uj.edu.pl; Tel.: +48-12-664-6498

**Abstract:** Historical maps are critical for long-term land use reconstructions; however, quantifying the uncertainty involved in comparing historical maps with recent data remains a considerable challenge. To date, many works have focused on the technical aspects of comparing historical and contemporary materials, but the potential sources of uncertainty inherent in historical data remain poorly understood. In this paper, we analyze the impacts of the topographic accessibility and cartographer's field position on the content quality of historical Austrian second military survey maps by referring to independent census data. Our results show that the topographic accessibility and visibility from the cartographer's surveying table points had very little impact on the map content quality and that the surveying table point locations were uniformly distributed throughout the area, regardless of the landscape conditions. These findings demonstrate that the second military survey maps can be seen as valuable and consistent historical data sources, making them especially useful for long-term land use research in Central Europe.

**Keywords:** historical maps; uncertainty; land use; visibility; topographic accessibility; Central Europe

#### **1. Introduction**

To date, dynamic land use changes have had substantial impacts on biodiversity [1], climate [2], and human well-being [3]. To assess the dynamics of such changes, long-term land change studies must be performed to facilitate a comparison with recent analyses. However, due to the limited availability of fully comparable remote sensing or statistical data [4], researchers often use regional archival maps to obtain historical land use information [5–7]. Since historical maps differ in scale, have inconsistent legends, or were produced with different aims, there is always uncertainty involved in using them for long-term comparisons [8–10]. The conceptual framework designed by Leyk et al. [11] to assess the level of uncertainty involved in such analyses distinguishes three independent domains of uncertainty for subsequent land use modelling—those inherent in historical data, those caused by data processing, and those dependent on the application. The assessment conducted by Leyk et al. [11] to evaluate the importance of these domains showed that the first domain—the uncertainty inherent in historical data (production-oriented uncertainty)—seems to be the most relevant.

In addition, automatic feature extraction techniques have been employed recently to acquire spatial information from historical maps referring mainly to the transformationoriented uncertainty [12–14]. Nevertheless, the quality of the original historical cartographic work, instructions, equipment, or procedures and their impacts on the map content and subsequent analyses remain poorly understood. For instance, Jenny et al. [15] created a software package capable of illustrating local map distortions that helped to understand the spatial distribution of potential errors in a broad range of map series; however, the

**Citation:** Ostafin, K.; Pietrzak, M.; Kaim, D. Impact of the Cartographer's Position and Topographic Accessibility on the Accuracy of Historical Land Use Information: Case of the Second Military Survey Maps of the Habsburg Empire. *ISPRS Int. J. Geo-Inf.* **2021**, *10*, 820. https:// doi.org/10.3390/ijgi10120820

Academic Editors: Motti Zohar and Wolfgang Kainz

Received: 22 September 2021 Accepted: 1 December 2021 Published: 4 December 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

sources of distortion were not identified. Recently, Ostafin et al. [16] indicated that changes in the content of neighboring map sheets of the same type are connected to the map sheets' authorship changes. Gimmi et al. [17] analyzed the contents of archival topographic maps against the true situations presented on orthophotos and terrestrial imagery serving as a proxy of the cartographer's perspective and indicated greater discrepancies between the maps and reference data in less topographic accessible and less visible locations. However, analyzing the historical topographic accessibility or visibility of a location and their impact on the map content quality is usually difficult due to the unavailability of historical road networks available in geographic information system (GIS) formats and a lack of knowledge regarding the precise position of the cartographer' in the field. Although accessibility is very often used as an important variable explaining environmental and social phenomena [18–20], usually proxy data are used in the form of current road networks, Euclidean distances to the settlements, or proximity to the waterways, instead of original, historical data [6,21].

In this paper, we were able to overcome both of the abovementioned obstacles and verify the impacts of these spatial determinants on the accuracy of Austrian second military survey maps covering a large portion of 19th century Europe [22,23]. The maps are a very popular source of historical information in the region and have thus far been used in many landscape-change analyses in different parts of Central Europe. The maps have been found useful for research on land use change [23–25], riverbed evolution [26], pond reconstruction [27], and landslide analyses [28]. Although the Austrian second military survey maps were produced as a form of generalization of the detailed cadastral maps (1:2880), little is known about the extent of changes, updates, and improvements done by the cartographers between the times of cadastral mapping (1844–1854 for Galicia) and military mapping (1861–1864 for Galicia). One can assume that the changes could be reflected in the independent census data collected in a similar period, however, due to a relatively short temporal framework of the map preparation, the potential updates could be also spatially determined and conducted mainly in the most accessible areas. Although reliability of the archival cartographic sources covering the Habsburg Empire are carried out [10,29,30], their authors indicate the necessity of further research, especially those focused on the comparability to other sources. It would be particularly useful in reconstruction of the land use over long time [31]. Similarly, like for the landscape painters, optimal location characterized by the high visibility enabled proper landscape reconstruction [32]; we also assumed that for the cartographers verifying and updating detailed cadastral information in order to create military maps, visibility played crucial work. Although general rules of the measurements conducting [33,34], or organizational structure [35–37] of the Austrian military cartography is well recognized, little is known about the impact of the local landscape conditions, or human factor on the quality of their work [16]. This might be especially important taking into account the size, weight, and quality of the former measuring equipment [38]. Specifically, we want to verify the role of the historical topographic accessibility and cartographer's field position on the accuracy of second military survey maps, based on the information from exactly the same cartographic sources, offering no time delay, which gives high potential in assessing their impact. In doing so, we sought to answer the following questions:


#### **2. Study Area**

Galicia was one of the largest provinces of the former Habsburg Monarchy (Austrian Empire). The region belonged to the Kingdom of Poland until 1772 and is currently divided

between Poland and Ukraine. To assess the impact of the cartographer's position and topographic accessibility on the accuracy of historical map information, we chose three test areas characterized by diverse landscape conditions within the region. The first test area (3 × 3 map sheets) was located in the western part of Galicia in the Carpathian Foothills and the Beskidy Mountains. The second test area (4 × 3 map sheets) was mainly located in the lowlands of the Sandomierz Basin. The third test area (3 × 4 map sheets) represented the eastern part of the region, covering the Outer Eastern Carpathians and Dniester Plain [39] (Figure 1). In total, these three test regions covered an area of 7813 km<sup>2</sup> , including the complete envelopment of 500 cadastral communes according to the historical boundaries during the 1860s (Table 1). Empire). The region belonged to the Kingdom of Poland until 1772 and is currently divided between Poland and Ukraine. To assess the impact of the cartographer's position and topographic accessibility on the accuracy of historical map information, we chose three test areas characterized by diverse landscape conditions within the region. The first test area (3 × 3 map sheets) was located in the western part of Galicia in the Carpathian Foothills and the Beskidy Mountains. The second test area (4 × 3 map sheets) was mainly located in the lowlands of the Sandomierz Basin. The third test area (3 × 4 map sheets) represented the eastern part of the region, covering the Outer Eastern Carpathians and Dniester Plain [39] (Figure 1). In total, these three test regions covered an area of 7813 km2, including the complete envelopment of 500 cadastral communes according to the historical boundaries during the 1860s (Table 1).

Galicia was one of the largest provinces of the former Habsburg Monarchy (Austrian

*ISPRS Int. J. Geo-Inf.* **2021**, *10*, x FOR PEER REVIEW 3 of 17

**2. Study Area** 

**Figure 1.** Study areas within the boundaries of historical Galicia. **Figure 1.** Study areas within the boundaries of historical Galicia.



#### **3. Materials and Methods 3. Materials and Methods**

The main materials we used in the study were the Austrian second military survey maps of 1861–1864 (1:28,800) (Figure 2) in addition to census data (from 1869) and land use data (from 1857) published at the commune level that were employed as the reference data (Figure 3). For the three test areas, we were able to compare the numbers of houses recorded in both data sources, while for test area 1, we also performed a comparison of the forest cover (based on information from 218 cadastral communes). The data on exact house locations from the maps were acquired from the building database of the second military survey, available for the area already in shapefile format. The publication related to this database discusses the geometric accuracy and generalization of the representation of buildings in relation to accurate cadastral maps and other factors [40]. Forest cover data The main materials we used in the study were the Austrian second military survey maps of 1861–1864 (1:28,800) (Figure 2) in addition to census data (from 1869) and land use data (from 1857) published at the commune level that were employed as the reference data (Figure 3). For the three test areas, we were able to compare the numbers of houses recorded in both data sources, while for test area 1, we also performed a comparison of the forest cover (based on information from 218 cadastral communes). The data on exact house locations from the maps were acquired from the building database of the second military survey, available for the area already in shapefile format. The publication related to this database discusses the geometric accuracy and generalization of the representation of buildings in relation to accurate cadastral maps and other factors [40]. Forest cover data were acquired as a result of manual vectorization within the FORECOM project [41]. The geometric accuracy of forest areas on the second military survey maps

has been verified in publications analyzing forest changes and fragmentation in the Polish Carpathians [42,43]. The corresponding census data on housing units were acquired from the 1869 census [44,45]. Housing units in the census are presented separately for the villages and small towns (ger. *Orts-Gemeinden*) and large estates (landed estates or manor farms) (ger. *Gutsgebiete*). Summarizing the abovementioned categories yielded the results for cadastral communes (ger. *Katastralgemeinden*). Census data on forest cover at the village level were acquired from other sources [46], where the data were collected as a result of cadastral measurements conducted in Galicia, mainly in 1844–1854. All the statistics are presented for cadastral communes separately for small and great ownership. Adding these two ownership categories together and changing the original planar units into hectares (1 Vienna morg [ger. *Joch*] = 0.57546 ha) [47] adequately modified the data for a comparison with the spatial map-based data. The commune boundaries were acquired by manual vectorization directly from the georeferenced sections of the second military survey maps (Figure 3). Some cadastral communes for test area 1 were aggregated due to slight changes in the commune list for 1857 compared to the cadastral communes from the maps; e.g., the lower and upper parts of the commune were merged into one. The differences in the number of buildings between the map and the census were calculated for each cadastral communes by comparing the number of buildings (points) on the map and from the census data. The differences in forest areas were calculated in a similar way. *ISPRS Int. J. Geo-Inf.* **2021**, *10*, x FOR PEER REVIEW 5 of 17

**Figure 2.** Austrian second military survey maps: **1**—cartographer positions (surveying table points), **2**—forest areas, **3**—buildings (houses in red), **4**—first-class roads, and **5**—maintained roads*.* tively; [52]). **Figure 2.** Austrian second military survey maps: **<sup>1</sup>**—cartographer positions (surveying table points), **2**—forest areas, **3**—buildings (houses in red), **4**—first-class roads, and **5**—maintained roads.

All the points were then used to create a visibility map based on the ALOS Global Digital Surface Model (Version 3.1) [51] with an original spatial resolution of 1 arc second; the data were resampled to 30 m and then transferred to the Lambert azimuthal equal

To assess the relations between the quality of the materials and either the topographic accessibility or the cartographer's visibility, we conducted a set of simple regressions; the difference between the map content and the census records was taken as the dependent variable, while topographic accessibility (mean distance to main roads) and visibility (invisible areas as a percentage of a particular area) were taken as the independent variables. The analyses were conducted at the map sheet level (n = 33) and cadastral commune level (n = 500), and the relations between topographic accessibility and visibility at the cadastral commune and map sheet levels were quantified by using Pearson's correlation (*r*), since the datasets followed the normal distribution (skewness *μ*3 = 0.96 and *μ*3 = 1.15, respec-

surface model resolution on the analysis by testing a light detection and ranging (LIDAR) based 1 m resolution digital elevation model (DEM) available only for contemporary Poland and compared it with the ALOS-based analysis. In the analysis, we used a variable reflecting invisible areas as a percentage of the area under study (i.e., the test area, map sheet, or cadastral commune). In addition to the simple division of visible/nonvisible areas, we also created a frequency mask where the output records indicate the number of times each cell location in the input surface raster can be seen by the input observation

locations. The visibility analysis was performed in ArcMAP 10.8.

**Figure 3.** Schematic of the proposed workflow. **Figure 3.** Schematic of the proposed workflow.

**4. Results**  *4.1. Initial Comparison of Digital Elevation Models*  Our visibility analysis was based on a globally available DEM model since the study area is located in the territory of Poland and Ukraine and because the more detailed LI-DAR-based DEM did not cover the whole study area. In accordance with the ongoing The historical topographic accessibility was analyzed based on the road network presented on the second military survey maps. The network we used consisted of four main road categories that were passable in the mid-19th century all year round [48]. The data were acquired in the form of GIS files [49] and used to create a 'distance to main roads' variable. For each cadastral commune, we then indicated a mean distance to main roads as a proxy of the topographic accessibility used in the subsequent steps of the analysis.

discussion, to ascertain whether a more detailed model might have an impact on the results [53,54], we tested the ALOS and LIDAR-based DEMs for two study areas located in Poland: one in a mountainous region and one in the lowlands (test areas 1 and 2, respectively, each with an area of 210 km2). In the mountainous terrain (test area 1), the difference in invisible areas between the two models was only 1.22% (Figure 4), while in the lowlands, the difference was higher at 9.21% (in both cases, more invisible areas were indicated in the 30 m resolution model than in the LIDAR-based model). Nevertheless, since only the ALOS model was available in test area 3 (most of which is mountainous terrain located in Ukraine), we decided to employ the ALOS model for the latter part of the analysis. *4.2. Visibility Analysis*  To assess the cartographer's visibility in the field, we first manually acquired all the map symbols indicating the cartographer's position (surveying table points) (ger. *Tischstand*, *Fixpunkte* [48]), i.e., the places where the measurements were acquired and where the cadastral maps used as a basis for creating the second military survey maps were updated. The second military survey maps were created by reducing the selected cadastral map sheet content: land cover, settlements, and communication lines. The basic materials were prepared by the Cadastre General Directorate (ger. *Generaldirection des Grundsteuer-Katasters*) [35,50]. The sections of the maps prepared in this way were then verified in the field, where topographic details and elevation were added to the maps; this process also verified the triangulation network and helped update the map content changes. In the countries where cadastral maps were available, each of the military cartographers had one field assistant to help him or her, and in the areas where cadastral maps were not available, the number of field assistants increased to three people [35].

According to the analysis, the density of surveying table points for each of the test areas was uniform at 0.27 points/km2. A visibility analysis based on the surveying table points indicated that in test area 1, 9.63% (22,021 ha) of the terrain was not visible to the cartographer, while in test areas 2 and 3, the corresponding percentages were 37.25% (102,970 ha) and 13.21% (36,503 ha), respectively. When analyzing invisible areas in smaller reference units, on average, 21% of each map sheet and 18% of the cadastral com-Contrary to the road-based visibility or other proxy-based variables used in many studies, the surveying table points considered herein indicate the exact position that the cartographer reached and used when working on the data. The map symbols indicating the cartographer's position were acquired from all 33 map sheets analyzed in the study and the 2-km buffer around each of the three test areas. We discovered 717, 866, and 861 cartographer positions for test areas 1, 2, and 3, respectively.

mune area were not visible. In test areas 1 and 3, the largest invisible regions were mainly steep slopes obscured by the ranges of the mountains and foothills. In contrast, due to the different landscape characteristics, the invisible areas in test area 2 were located mainly in

All the points were then used to create a visibility map based on the ALOS Global Digital Surface Model (Version 3.1) [51] with an original spatial resolution of 1 arc second; the data were resampled to 30 m and then transferred to the Lambert azimuthal equal area (LAEA) (EPSG 3035) reference system. Additionally, we verified the impact of the surface model resolution on the analysis by testing a light detection and ranging (LIDAR)-based 1 m resolution digital elevation model (DEM) available only for contemporary Poland and compared it with the ALOS-based analysis. In the analysis, we used a variable reflecting invisible areas as a percentage of the area under study (i.e., the test area, map sheet, or cadastral commune). In addition to the simple division of visible/nonvisible areas, we also created a frequency mask where the output records indicate the number of times each cell location in the input surface raster can be seen by the input observation locations. The visibility analysis was performed in ArcMAP 10.8.

To assess the relations between the quality of the materials and either the topographic accessibility or the cartographer's visibility, we conducted a set of simple regressions; the difference between the map content and the census records was taken as the dependent variable, while topographic accessibility (mean distance to main roads) and visibility (invisible areas as a percentage of a particular area) were taken as the independent variables. The analyses were conducted at the map sheet level (n = 33) and cadastral commune level (n = 500), and the relations between topographic accessibility and visibility at the cadastral commune and map sheet levels were quantified by using Pearson's correlation (*r*), since the datasets followed the normal distribution (skewness *µ*<sup>3</sup> = 0.96 and *µ*<sup>3</sup> = 1.15, respectively; [52]).

#### **4. Results**

#### *4.1. Initial Comparison of Digital Elevation Models*

Our visibility analysis was based on a globally available DEM model since the study area is located in the territory of Poland and Ukraine and because the more detailed LIDARbased DEM did not cover the whole study area. In accordance with the ongoing discussion, to ascertain whether a more detailed model might have an impact on the results [53,54], we tested the ALOS and LIDAR-based DEMs for two study areas located in Poland: one in a mountainous region and one in the lowlands (test areas 1 and 2, respectively, each with an area of 210 km<sup>2</sup> ). In the mountainous terrain (test area 1), the difference in invisible areas between the two models was only 1.22% (Figure 4), while in the lowlands, the difference was higher at 9.21% (in both cases, more invisible areas were indicated in the 30 m resolution model than in the LIDAR-based model). Nevertheless, since only the ALOS model was available in test area 3 (most of which is mountainous terrain located in Ukraine), we decided to employ the ALOS model for the latter part of the analysis.

#### *4.2. Visibility Analysis*

According to the analysis, the density of surveying table points for each of the test areas was uniform at 0.27 points/km<sup>2</sup> . A visibility analysis based on the surveying table points indicated that in test area 1, 9.63% (22,021 ha) of the terrain was not visible to the cartographer, while in test areas 2 and 3, the corresponding percentages were 37.25% (102,970 ha) and 13.21% (36,503 ha), respectively. When analyzing invisible areas in smaller reference units, on average, 21% of each map sheet and 18% of the cadastral commune area were not visible. In test areas 1 and 3, the largest invisible regions were mainly steep slopes obscured by the ranges of the mountains and foothills. In contrast, due to the different landscape characteristics, the invisible areas in test area 2 were located mainly in large river valleys. In the mountainous and foothill areas featuring diverse landforms and high relative elevation differences, the invisible areas were usually compact, while in the lowlands, the pattern of invisible areas was scattered as relatively small patches (Figure 5). 5).

**Figure 4.** Invisible areas defined from the surveying table points (red) for the 30 m resolution ALOS model (**upper** map) and 1 m resolution LIDAR model (**lower** map). **Figure 4.** Invisible areas defined from the surveying table points (red) for the 30 m resolution ALOS model (**upper** map) and 1 m resolution LIDAR model (**lower** map).

high relative elevation differences, the invisible areas were usually compact, while in the lowlands, the pattern of invisible areas was scattered as relatively small patches (Figure

*ISPRS Int. J. Geo-Inf.* **2021**, *10*, x FOR PEER REVIEW 8 of 17

**Figure 5.** Invisible areas (black) and cartographer positions (red) in the three test areas (1, 2, 3 as in Figure 1).

Figure 1).

According to the analysis, 7.77% of the forests were not visible to the cartographers, while 6.44% of their area was visible from one surveying table point, and 85.79% of the forested area was visible from more than one surveying table point (Figure 6). A similar analysis of houses indicated that 17.7% (11,037 houses) were invisible to the cartographers, 11.6% were visible from one surveying table point (7264 houses), and 70.7% (44,201 houses) of all houses were visible from more than one surveying table point (Figure 7).

forested area was visible from more than one surveying table point (Figure 6).

**Figure 5.** Invisible areas (black) and cartographer positions (red) in the three test areas (1, 2, 3 as in

According to the analysis, 7.77% of the forests were not visible to the cartographers, while 6.44% of their area was visible from one surveying table point, and 85.79% of the

*ISPRS Int. J. Geo-Inf.* **2021**, *10*, x FOR PEER REVIEW 9 of 17

**Figure 6.** Forested areas in test area 1 that were invisible (black), visible from one surveying table point (yellow), and visible from two or more surveying table points (green). **Figure 6.** Forested areas in test area 1 that were invisible (black), visible from one surveying table point (yellow), and visible from two or more surveying table points (green).

A similar analysis of houses indicated that 17.7% (11,037 houses) were invisible to the cartographers, 11.6% were visible from one surveying table point (7264 houses), and 70.7% (44,201 houses) of all houses were visible from more than one surveying table point (Figure 7).

#### *4.3. Impacts of Topographic Accessibility and Visibility on the Map Content*

On average, the surveying table points (n = 2444) were located 2974 m from a main road. There was almost no relation between the mean distance to main roads and the number of surveying table points located in the cadastral communes (*r* = 0.12, *p* < 0.05). Furthermore, our analysis indicated a relatively weak correlation between topographic accessibility and visibility at both the map sheet level (*r* = 0.36, *p* < 0.05) and the commune level (*r* = 0.39, *p* < 0.05).

The quality of the maps reflected by the difference in the number of houses observed on the map sheets and expressed as a percentage of the corresponding census data was not explained by either the visibility (R<sup>2</sup> = 0.02) or the topographic accessibility (R<sup>2</sup> = 0.01) at the commune level. A similar analysis for forests indicated that the difference in the forested area between the maps and census data at the commune level also could not be explained by either visibility (R<sup>2</sup> = 0.003) or topographic accessibility (R<sup>2</sup> = 0.02).

*ISPRS Int. J. Geo-Inf.* **2021**, *10*, x FOR PEER REVIEW 10 of 17

**Figure 7.** Houses in the three test areas (1, 2, 3 as in Figure 1) that were invisible (black), visible from one surveying table point (yellow), and visible from two or more points surveying table points (green). **Figure 7.** Houses in the three test areas (1, 2, 3 as in Figure 1) that were invisible (black), visible from one surveying table point (yellow), and visible from two or more points surveying table points (green).

#### **5. Discussion**

#### *5.1. Difference between Census and Map Data*

When analyzing the differences between the number of houses vectorized from the archival maps and the number of houses recorded in the census data, we must bear in mind that the difference in publication time between the maps and census is 5–8 years. Moreover, the target user and the purpose of creation also diverge between the two sources: the maps were created for military purposes to enable the efficient movement and stationing of military troops, while the census served as a tool for the development of sociodemographic policy. Nevertheless, it is difficult to assess the impact of the temporal difference between these sources over large areas. A cross-validation conducted on the local level and based on independent sources including, for example, parish documents would be helpful in answering this question. The overall difference in the number of houses based on the census (71,755) and that from the map-based vectorization (62,502) was substantial at 9253 houses. In the first test area, there were fewer houses according to the census than there were according to the maps (23,076 vs. 23,396, respectively). For the second and third test areas, however, the numbers based on the census were higher (test area 2, 31,298 vs. 27,580; test area 3, 17,381 vs. 11,526). These differences were especially high in the small towns. For instance, in Drohobych (ukr. oó, pol. Drohobycz), the difference was approximately 60% (1816 houses in the census vs. 1085 houses on the maps); in Boryslav (ukr. oa, pol. Borysław), the difference reached an astounding 479% (537 houses according to census vs. 112 on the maps) (Figure 8). The differences mentioned above might be partly due to cartographic generalization, where a large number of houses could not be presented on the maps in detail. However, in the 19th century, part of test area 3 (especially Boryslav) became a globally important centre of the oil industry, which triggered highly dynamic development [55] that was not fully captured by the materials. Since the maps are, in general, older than the census data with regard to housing units, a higher number of houses recorded in the statistical data may reflect the general sociodemographic trend resulting in increases in the overall population [56] and number of housing structures, especially in the Carpathian part of the region [57]. In 1857, the census for the districts of Galicia indicated 4,597,470 inhabitants and 760,181 houses [58], while 12 years later, there were 5,418,016 inhabitants and 857,702 houses [44], thereby reflecting the dynamics of the economic development and demographics.

The temporal difference between the census data and maps in terms of the forest cover (10–20 years, depending on the area) are even higher than that in the case of housing units. Regarding the forested area, the census data [46] had to be slightly aggregated, as the administrative commune coverage differed at the time of publication when compared to the maps. We analyzed 218 spatial units based on the 1857 census (presenting data on forest cover) vs. 222 units in 1869 (containing data on the number of houses). The statistical data present areas up to 1 morg (0.57546 ha), while the minimal area of the map-based forest cover patch is 0.02 ha. The overall difference in the forest area based on the census (50,423 ha) and the one from the map-based vectorization (48,518 ha) was 1905 hectares (Figure 9). Being aware that the differences between the census and the map-based information are substantial, we argue that the role of visibility and topographic accessibility in terms of second military survey maps might be lower than commonly expected in the research on other historical maps quality. Our research indicated uniformly distributed surveying table points, regardless the landscape conditions, which indicates limited role of the above-mentioned variables on the map content quality.

*ISPRS Int. J. Geo-Inf.* **2021**, *10*, x FOR PEER REVIEW 12 of 17

**Figure 8.** Differences in the number of buildings between census and maps in cadastral communes in the three test areas (1, 2, 3 as in Figure 1). **Figure 8.** Differences in the number of buildings between census and maps in cadastral communes in the three test areas (1, 2, 3 as in Figure 1).

#### The temporal difference between the census data and maps in terms of the forest *5.2. Impact of the Materials Used in this Study on the Final Results*

cover (10–20 years, depending on the area) are even higher than that in the case of housing units. Regarding the forested area, the census data [46] had to be slightly aggregated, as the administrative commune coverage differed at the time of publication when compared to the maps. We analyzed 218 spatial units based on the 1857 census (presenting data on forest cover) vs. 222 units in 1869 (containing data on the number of houses). The statistical data present areas up to 1 morg (0.57546 ha), while the minimal area of the map-based forest cover patch is 0.02 ha. The overall difference in the forest area based on the census (50,423 ha) and the one from the map-based vectorization (48,518 ha) was 1905 hectares (Figure 9). Being aware that the differences between the census and the map-based information are substantial, we argue that the role of visibility and topographic accessibility in terms of second military survey maps might be lower than commonly expected in the research on other historical maps quality. Our research indicated uniformly distributed surveying table points, regardless the landscape conditions, which indicates limited role of the above-mentioned variables on the map content quality. The main source of the information used in this study was second military survey maps, which served as the source of the road networks, houses locations, forest cover, and surveying table point locations. The maps were pre-processed before acquiring the data, and one of the steps was georeferencing. Geometric correction and georeferencing to the Poland CS92 (EPSG 2180) or UTM zone 34N (EPSG 32,634) coordinate system were obtained using 2nd-order polynomial transformation and then transferred to the Lambert azimuthal equal area (LAEA) reference system. This step was accomplished by using at least 20 control points per map sheet, and the average value of the root mean square error (RMSE) in most cases ranged between 10 and 30 m, only occasionally exceeding the upper threshold. On the one hand, we are aware that changing the surveying table point locations within the RMSE values might have an impact on the resulting visibility maps; however, the upper RMSE value is close to 1 mm at the original scale of the map (1:28,800). Hence, the comparison of the visibility maps, being a result of potential changes in the cartographer's position, is questionable, as the symbol representing the cartographer's position on the maps was more dependent on the drawing precision of the original map than on the RMSE value itself. This example supports the findings of Leyk et al. [11]

regarding the role of production-oriented uncertainty in overall uncertainty assessments of historical land use maps. *ISPRS Int. J. Geo-Inf.* **2021**, *10*, x FOR PEER REVIEW 13 of 17

**Figure 9.** Differences in the forest area between census and maps in cadastral communes. **Figure 9.** Differences in the forest area between census and maps in cadastral communes.

*5.2. Impact of the Materials Used in this Study on the Final Results*  The main source of the information used in this study was second military survey maps, which served as the source of the road networks, houses locations, forest cover, and surveying table point locations. The maps were pre-processed before acquiring the data, and one of the steps was georeferencing. Geometric correction and georeferencing to the Poland CS92 (EPSG 2180) or UTM zone 34N (EPSG 32,634) coordinate system were obtained using 2nd-order polynomial transformation and then transferred to the Lambert azimuthal equal area (LAEA) reference system. This step was accomplished by using at least 20 control points per map sheet, and the average value of the root mean square error (RMSE) in most cases ranged between 10 and 30 m, only occasionally exceeding the upper Our analysis also showed that using different DEMs might have an impact on the visibility analysis. However, as our test areas were located in two countries, we had to use a model available in both of them. That is why we employed the ALOS DEM. Moreover, most of the test areas were located in mountainous terrain, and the differences between the models, according to our tests, were relatively low. Being aware that a more detailed DEM could yield different visibility maps, we must remember that the ALOS model showed, on average, more invisible areas than the LIDAR-based model, so using the former model over the whole study area would potentially result in even less of an impact of visibility on the quality of the maps. Such a result would confirm our findings that the quality of the second military survey maps did not depend on the surveying table point locations.

threshold. On the one hand, we are aware that changing the surveying table point locations within the RMSE values might have an impact on the resulting visibility maps; however, the upper RMSE value is close to 1 mm at the original scale of the map (1:28,800). Hence, the comparison of the visibility maps, being a result of potential changes in the cartographer's position, is questionable, as the symbol representing the cartographer's position on the maps was more dependent on the drawing precision of the original map than on the RMSE value itself. This example supports the findings of Leyk et al. [11] regarding the role of production-oriented uncertainty in overall uncertainty assessments of historical land use maps. Our analysis also showed that using different DEMs might have an impact on the visibility analysis. However, as our test areas were located in two countries, we had to use a model available in both of them. That is why we employed the ALOS DEM. Moreover, most of the test areas were located in mountainous terrain, and the differences between the models, according to our tests, were relatively low. Being aware that a more detailed DEM could yield different visibility maps, we must remember that the ALOS model showed, on average, more invisible areas than the LIDAR-based model, so using the former model over the whole study area would potentially result in even less of an impact The potential drawback of our work might be that the impact of land cover on the cartographer's visibility is omitted. However, we hypothesize that the cartographer's position was chosen carefully such that the impact of the land cover was assessed before performing the measurements. The height of vegetation, although roughly discernible from cadastral maps, cannot be determined from the second military survey maps used herein. As the forest height was diverse at the time, assigning one value for the whole forest would introduce additional uncertainty, which we decided to avoid. Usually, archival map quality assessments are conducted when the old map content can be directly compared with other cartographic sources from different periods [23], maps from similar periods but at different scales [10,29,59], or other contemporary data sources such as LIDAR point clouds to reconstruct land use changes [31]. In such situations, quality assessments of the historical maps are not the main aim of the research [24,60], and, thus, the quality is not questioned as long as the map content confirms the general trends, e.g., of the land use trajectories over time. In our analysis, we wanted to transcend such assessments to focus on the uncertainties that are dependent on the acquisition processes of historical data. Although these problems are very important for the conscious use of original sources, these issues are rarely discussed [61–63].

#### of visibility on the quality of the maps. Such a result would confirm our findings that the **6. Conclusions**

quality of the second military survey maps did not depend on the surveying table point locations. The potential drawback of our work might be that the impact of land cover on the In this paper we aimed to explain the role of accessibility and visibility on the potential updates and improvements of the second military survey maps, as a generalization of

cartographer's visibility is omitted. However, we hypothesize that the cartographer's po-

performing the measurements. The height of vegetation, although roughly discernible

the detailed cadastral mapping. Our analysis showed that topographic accessibility and visibility did not impact the quality of the land cover information from 19th century second military survey maps. Most likely, the high quality of the second military survey maps was achieved mainly by using very detailed, large-scale cadastral maps as their basis, which had never been achieved theretofore, not only for Galicia, but also for other parts of the Habsburg Empire. However, contrary to the cadastral mapping, the second military survey maps also offer analysis of the map content in relations to other topographic features. Historical maps offer the exact location of the objects in space, which cannot be acquired from detailed census data. The scale of the second military survey maps is appropriate for conducting data analyses at the level of cadastral communes, allowing them to be compared with census values on a large scale [64,65]. These 19th century cadastral maps are still commonly used for historical and geographical research, e.g., when studying land tenure and their positional accuracy has been confirmed, e.g., when compared in GIS software with new large-scale datasets [66]. Our results might be of importance for research conducted on large portions of the former Habsburg Empire (including, for example, contemporary northern Italy, Austria, Czechia, Slovakia, Hungary, Slovenia, Croatia, Romania, southern Poland, or western Ukraine), where second military survey maps are available, although their production-oriented uncertainty is still insufficiently recognized.

Overall, it is important to analyze the uncertainty inherent in historical data, as an increasing number of historical maps are becoming accessible electronically and the methods of data acquisition are developing rapidly. As a consequence, historical, mapbased information is becoming more accessible and is being implemented in a greater number of analyses of different environmental studies.

**Author Contributions:** Krzysztof Ostafin and Dominik Kaim, conceptualization, methodology, formal analysis, and writing—original draft preparation; Małgorzata Pietrzak, writing—review and editing. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Ministry of Science and Higher Education, Republic of Poland under the frame of "National Programme for the Development of Humanities" 2015–2021, as a part of the GASID project (Galicia and Austrian Silesia Interactive Database 1857–1910, 1aH 15 0324 83).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


## *Article* **Application of GIS Tools in the Measurement Analysis of Urban Spatial Layouts Using the Square Grid Method**

**Łukasz Musiaka and Marta Nalej \***

Faculty of Geographical Sciences, University of Lodz, Kopci ´nskiego 31, 90-142 Łód´z, Poland; lukasz.musiaka@geo.uni.lodz.pl

**\*** Correspondence: marta.nalej@geo.uni.lodz.pl; Tel.: +48-426-35-45-68

**Abstract:** The principal aim of this paper is to present the capabilities of newly developed GIS tools for measurement analysis of urban spatial layouts, using the square grid method. The study of urban morphology and metrology is a multistage process, which involves the metrological analysis of town plans. The main research step is the determination of measurement modules of town layouts, using the square grid. By using GIS software, the authors developed a new tool, named HGIS Tools, which allow to create any number of modular grids with the selected cell size that corresponds to urban units of distance and surface area. When investigating a large number of towns and cities, this offers a significant improvement of the research procedure. The paper presents a test of the tool's potential on the example of regular medieval towns from the area of the former Teutonic Order state (currently the territory of Poland), diversified in terms of size, genesis and morphometrics. The obtained results confirmed that HGIS Tools allowed to determine the hypothetical measurement module of the layout of the cities studied. The results were consistent with the analyses of other authors carried out with the traditional grid-square methods. The test of the HGIS Tools showed their significant potential in conducting morphometric analyses of spatial arrangements of urban spatial layout on a larger scale.

**Keywords:** historical GIS; HGIS; GIS tools; fishnet; grid; urban morphology

#### **1. Introduction**

The purpose of this paper is to present the possible use of GIS tools for the measurement analysis of urban spatial layouts, using the square grid method. The functionality of tools available in the ArcMap application from ESRI is insufficient to conduct morphological and morphogenetic research. The authors had to create appropriate tools, using the Python programming language. The operational capabilities of the new tools were presented on examples of medieval towns of the Teutonic Order in Prussia (now in Poland): Nowy Staw, Elbl ˛ag and Grudzi ˛adz. The main criterion for the selection of the research group was the regularity of spatial planning and similar period of origin and territorial vicinity, resulting from being located by the same investor—the State of the Teutonic Order. The article is an introduction to a study of the measurement of a larger group of cities in the area of former Prussia, and its main purpose is to present and test the tool.

Scientific studies of the origins and processes of the formation of the spatial layouts in towns and villages date back to the 19th century [1]. Today, the morphological and morphogenetic research of settlement unit is mostly conducted by geographers, architects, historians and art historians [2]. Morphological studies of the settlement are also the domain of archaeology [3,4]. Each of these disciplines uses its own research apparatus, which emphasizes the different methodological aspects and sources. Regardless of the methodology adopted, retrogressive investigation of spatial arrangements faces many difficulties. The main barriers include the shortcomings of the source base, both written and cartographic, which often forces the hypothetical nature of many analyses [5,6]. The common denominator of morphological and morphogenetic settlement research is

**Citation:** Musiaka, Ł.; Nalej, M. Application of GIS Tools in the Measurement Analysis of Urban Spatial Layouts Using the Square Grid Method. *ISPRS Int. J. Geo-Inf.* **2021**, *10*, 558. https://doi.org/ 10.3390/ijgi10080558

Academic Editors: Motti Zohar and Wolfgang Kainz

Received: 13 July 2021 Accepted: 15 August 2021 Published: 17 August 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

the analysis of archival cartographic material [1,7–14]. However, we can be faced with numerous challenges when attempting to use archival maps and plans to obtain data concerning the spatial development of a given settlement unit in the past [15]. There can be also listed bad technical conditions and decreased legibility of archival documents, lack of carthometry or faithful mapping of the surface on a map plane, and the inaccuracy of archival cartography [5,16–19]. Metrological studies of modern base maps alone, not backed up by archaeological and architectural research, may lead to erroneous conclusions, due to the changing ownership divisions and development over the ages [20]. The process of 'town-plan analysis', which uses post-medieval cartography as a source from which to reconstruct medieval urban topographies, in contrast to the Conzenian school, is often criticized [15,21]. For these reasons, results based on cartographic sources have to be confronted with other data sources [2]. Analyses of non-cartographic sources conducted by historians and art historians are especially important here. For dozens of years now, the significance of archaeology and architecture in settlement research has been growing, as they serve to verify the assumptions and theories of both geographers and historians [22,23]. Therefore, it seems that the best results can be achieved by combining different research approaches with a critical analysis of the results obtained. Unfortunately, in practice, such studies are rare.

In urban morphogenetic and metrological research, the main aim is usually to find the roots of the form, to explain its genesis, and then to recreate its past development [24]. In morphological studies, the primary issue in the spatial research of urban layouts (measurement analysis) is the determination of modular rules for its measurement, and the discovery of the units of distance and area measurement used when plotting them in the field. The next stage in the metrological analysis is the search for the modular structure of field parcels, plots in urban settlements, urban blocks, building sites and other elements [25].

One of the methods used for metrological analyses of urban measurements is the square grid [25–31]. It is used as an auxiliary method when searching for modular divisions of a settlement, overlaying it on the layout under analysis [26,32]. The grid may be wellmatched or not. The grid has to be drawn over a copy of the town plan, using trial and error. To examine a larger group of towns, a whole set of grids of different gauges and scales could be drawn by hand. Slightly translucent paper or transparent film for back-lit projectors was used for this purpose. The results of the analysis, the chartered spatial layout within a grid of a certain gauge, was presented on a separate sheet. Figure 1 shows the use of a typical manual graph paper for that purpose.

**Figure 1.** Hand-drawn sketch of the most important elements of the city layout on the background of the squares grid, example of Wieliczka [33].

By using a modular grid for a town plan and conducting detailed metrological measurements, the shapes, size and proportions of markets, as well as plots and urban blocks can be analyzed. Based on empirical research, the initial dimensions of urban plots used to determine the most widely used measuring modules in the city centers can be recreated. By focusing on the issues of sizes and proportions of cities, modular layout schemes can be created for groups of cities with common geographical and chronological scopes [34,35].

The medieval surveyors may have plotted a certain grid of squares in the field when initially measuring for the city. Alternatively, they could have used the tradition of ancient geometry [25,34]. Such a method is an excellent way to facilitate the spatial planning of a regular medieval city. Very rarely can we find almost perfectly geometric examples of chartered cities because the initial plan is deformed by terrain limitations, such as a bend in the river, a steep hill, previous buildings and roads, or ownership divisions [36]. New towns, founded on the 'raw root' (in cruda radice), or in previously uninhabited locations, had a better chance of having a regular grid of streets and blocks [1]. This method of simple medieval parceling allows us to overlay rod or rope grids and other dimensions onto the plans of different cities.

We can distinguish two types of grids: primary—actual measurements and a grid of essential elements of the urban layout—which facilitates the discovery of the composition idea behind the urban settlement, an analysis of the elements of layout of individual systems, such as the location of the market, the ratio between its size and the overall area of the city, as well as an analysis of the repeatability of market dimensions in other elements of the city. Additionally, it makes it easier to compare different planning concepts [25].

The progress in the computerization of research processes may facilitate the morphological and morphometrical studies of urban centers. Over the last twenty years, new spatial analysis methods have appeared, utilizing GIS tools, as well as spatial and historical tools—Historical Geographic Information Systems (Historical GIS or HGIS) [37–44]. New technologies allow for improvements of traditional research methods, introducing brand

new data processing methods, and even expanding the current research perspectives. The application of GIS in studying the morphology and morphometry of settlements will surely not solve all methodological problems but may significantly facilitate analyses and create an opportunity to combine different approaches on a common plane.

The paper is organized as follows. The first section "Introduction" provides a brief overview of the traditional grid method used in urban morphology research. The next section, "Materials and Methods", describes the study area and presents the new HGIS Tools set. The "Results" section contains analysis of usefulness of HGIS Tools in morphological studies, which is presented based on studies of selected cities. The article closes with "Discussion" and "Conclusions", which summarize the possibility of using GIS tools for the analysis of historical urban layouts with the use of square grids.

#### **2. Materials and Methods**

#### *2.1. Study Area*

The examined towns of Elbl ˛ag, Grudzi ˛adz and Nowy Staw are located in the northern part of present Poland [45–47]. The urban centers selected to test the tool are characterized by a regular spatial layout, typical for cities of the late Middle Ages. They were granted city rights in the years 1246 (old town Elbl ˛ag), 1291 (Grudzi ˛adz) and 1343 (Nowy Staw). In Elbl ˛ag, instead of a rectangular or square market, there is a functioned widened street, parallel to the waterfront. Nowy Staw, founded under Chełmno law, also has an elongated market square, but the reason for this was the important function of serving the trade of the Vistula delta area. In Grudzi ˛adz, we have a layout of a classical Gothic town with an orthogonal spatial layout with a market in the center, in place of one of the quarters. From the 13th and 14th centuries until 1466, the towns and cities were located within the territory of the State of the Teutonic Order in Prussia. They were located under the rule of the Teutonic Order, so they are related not only in terms of chronology (time of foundation), but also territorially.

Figure 2 illustrates the geographical location of the towns analyzed with the HGIS tool. The studied centers are placed both within the present-day and medieval political boundaries.

#### *2.2. Data Sources*

The research mainly used digitalized city plans and orthophotomaps in various scales (Table 1). The time range of selected cartographic materials is 1710–2020. All cartographic materials were available in the form of reproductions in monographic publications (Nowy Staw) or as digitized source materials collected in the Historical Atlas of Polish Towns (Elbl ˛ag, Grudzi ˛adz), and in the POLONA digital repository (Elbl ˛ag). A contemporary orthophotomap from the portal mapy.geoportal.gov.pl (Grudzi ˛adz), which was also used in this study. The cartographic materials represent the historic center of Elbl ˛ag and Grudzi ˛adz. In the case of Elbl ˛ag, only the spatial layout of the Old Town was taken into consideration. The New Town area was excluded from the analysis. In Grudziadz, plans covering the area of the medieval town surrounded by walls were selected for analysis. In the case of Nowy Staw, due to its small spatial scale, the archival plans included the area of the entire town, together with its suburbs, and the central area, without the suburbs and the adjacent rural development of Nytyska Wie´s.

**Figure 2.** Location of the study area [48].



The historical cartographic materials were spatially georeferenced. Archived plans, not available in digital form, were scanned at a resolution of 400 dpi and saved in \*.tif format. Then, using the georeferencing tool available in the ArcMap software, the analyzed material was calibrated. The materials used in the study did not include information about the spatial reference system assigned to them, which would make it possible to perform mathematical operations converting the coordinates of points from one coordinate system to another, calibrating the plans. Therefore, a common orthophotomap was used as a

reference map, by determining the Ground Control Points (GCP), they were georeferenced in the PL-1992 [6,53] coordinate system currently in use in Poland. Root Mean Square Error (RMS) was used as a measure of the spatial fit quality. RMS values and the number of ground control points are summarized in Table 2. Based on literature findings, it was considered that the achieved accuracies were sufficient for further research [54–56].

**Table 2.** Characteristics of the quality of spatial matching of source materials, summary of Ground Control Points (GCP) and Root Mean Square Errors (RMS).


Additionally, in the research, publications and sources concerning the use of square grid methods and the application of the GIS method in historical studies, as well as those describing aspects of the history and metrology of the studied cities, were used.

#### *2.3. GIS Tools*

Conducting metrological analyses using computer techniques is similar to traditional 'analogue' methods. The primary stage of research involves the determination of measurement modules using a square grid. In the case of GIS tools, this stage has to be preceded by the spatial alignment (georeferencing) of archival cartographic materials. Calibration of the historical map is a vast subject, which has already been discussed at length in the literature [6,18,57–68].

Initially, the square grid for metrological analyses was created using ArcMap 10.7, especially the geoprocessing tools in the ArcToolBox module. Primarily, the Create Fishnet tool available in the Data Management Toolbox, Sampling set, was used. The Create Fishnet tool is available in all versions of ArcMap and in ArcGIS Pro, whereas in the open source software QGIS, there is a tool called "Create grid", which has similar functionality but additionally allows to create diamond and hexagon cell shapes [69]. It was discovered that the functionality of these tools is insufficient to conduct morphometrical research, as it does not allow full customization of the overlaid grid to the required shape and orientation of the polygon study area, which limits their usability, especially in the case of a larger number of analyzed towns plans. The authors used Python programming language to create their own set of tools and added it to the ArcToolbox in ArcMap 10.7. The ArcPy site package was used, as it allows for creating scripts in Python to enable access to geoprocessing tools, as well as additional functions, classes and modules. This allows for significant improvements, even to the more complex data processing. For ease of access, the scripts (tools) were placed in the HGIS Tools box (Supplementary Materials). The set consists of two scripts: HGIS Fishnet and HGIS Fishnet Rhombus.

The first and primary analytical tool used to create a square grid is the Fishnet script. It is based on the CreateFishnet\_management (...) function of ArcPy site package, which requires several parameters. The primary innovative function of the HGIS Fishnet tool is the ability to precisely align the grid with the area of study.

The functionality of the script allows the creation of any number of grids with different cell sizes. What is especially important for metrological research is that the grids can be applied to any location of the city plan that serves as a base layer for analyses. The location

must be defined by the user as rectangular polygon saved in the shape file. This allows for acceleration and significant automation of work, with minimum user requirements. The user has to specify only four parameters: working space—location, where the output layer will be stored, the name of the output layer, the location of the shape file containing the extent of the study area and the cell size (Figure 3). The functionality of this tool is enhanced by adding the ability to choose the cell size from a list that contains medieval distance measures used in the area of study. In addition, the HGIS Fishnet tool allows for differentiating the width and height of the cell size by filling in the 'Multiplicity' field.


**Figure 3.** Description of the HGIS Fishnet tool interface [70].

In the next stages of the square grid creation, invisible to the user, a buffer is created around the area of study (Figure 4). The buffer is sized according to the cell gauge, which allows for its later precise alignment to the city plan (base layer) to determine its measurement module. The coordinates of the four end points for the maximum spatial buffer range are then downloaded and used to apply rotation to the grid, in line with the rotation of the area of study. The result is a grid with the chosen cell size and proper rotation, aligned to the extent of study area, expanded by the initially created buffer.

Not all medieval chartered towns were identically measured in space. This could have been caused by an uneven terrain, the course of rivers, an earlier communication system, or the existence of earlier ownership divisions. Therefore, another tool called HGIS Fishnet Rhombus was created to allow for the formation of a grid based on a study area that is not rectangular. HGIS Fishnet Rhombus allows the user to create a rhomboid grid. The principle operation of this tool is different than that of HGIS Fishnet. Grid geometry (angles α and β) (Figure 5) is based on the drawn polygon (tetragon) that includes the study area. The tool uses mainly the UpdateCursor (...) function of the ArcPy site package. The script creates construction lines, parallel to the boundaries of the study area, at intervals defined by the user. The user can choose the cell size from a list or define it directly by entering a value in meters. Then, the construction lines are changed into a grid. A certain limitation of this tool is that the areas of study are of an irregular shape, in some specific cases, the grid may not completely cover the area. However, this problem can be solved by slightly enlarging the polygon that marks the boundaries of the area under investigation.

**Figure 4.** HGIS Fishnet tool workflow stages: (**a**) buffering study area, (**b**) collecting of the four end points of the buffer maximum range coordinates, (**c**) rotation of the grid in accordance with the study area rotation, (**d**) trimming of the grid to the extent of the study area.

Although the operation of both tools from the HGIS Tool box is based on the assumptions of the Create Fishnet tool from the Data Management Tools Sampling set, available in ArcToolbox (ArcMap 10.7), significant changes were made. The HGIS Tools enable better matching of the grid to the study area, while reducing the need to specify many parameters by the user. Table 3 presents a comparison of the parameters and functionality of the standard Create Fishnet tool and HGIS Tools.

**Figure 5.** HGIS Fishnet Rhombus tool workflow stages: (**a**) creating lines referring to the area of study at a distance corresponding to the adopted measure unit and determining the size of the angles of the study area range, (**b**) copying and shifting horizontal and vertical lines within the range, (**c**) grid line coverage of the study area, (**d**) converting the line grid to a polygon grid and trimming of the grid to the extent of the study area.

**Table 3.** Comparison of the parameters and functionality of the standard Create Fishnet tool (standard ArcToolbox—ArcMap 10.7) and HGIS Tools.



**Table 3.** *Cont.*

Geospatial approaches for the analysis of historical cartographic material requires their proper preparation and the use of geoprocessing tools. The procedure of preparation of historical cartographic materials for analysis with the use of HGIS tools was summarized in a schematic workflow given in Figure 6.

**Figure 6.** Schematic workflow of preparation of historical cartographic materials for analysis with the use of HGIS tools.

#### **3. Results**

The functionality of the HGIS Tool set was presented by overlaying a number of square grids on archival and modern maps of the chosen towns. The examples of towns studied were selected to represent both different morphogenesis, size and importance during the medieval period and today, and to have slightly different features of spatial planning of their old towns. Nowy Staw is characterized by a very elongated town square, whereas Elbl ˛ag is characterized by a grid layout whose actual layout module, due to the adaptation to local terrain features, is neither square nor rectangular. Finally, Grudzi ˛adz represents the typical, fully developed form of a Gothic town with a regular market square. Such a selection of the settlement units allowed for the analysis of various functionalities of the tested tool.

The analysis started with Nowy Staw, a small town in Gda ´nsk Pomerania, chartered over an earlier village, which was built near a pre-charter settlement. The most likely date of charter made by the Teutonic Knights under Chełmno law is 1343. In the Middle Ages, Zuławy was already a very productive agricultural area, due to the fertile alluvial soil ˙ forming the Vistula delta. Its considerable distance from Gda ´nsk and Elbl ˛ag did not allow a hike from one's residence to the cities to sell crops and back home in a single day. Hence, there was a need to create a new commercial center to serve local agriculture in the middle of the Vistula delta.

The commercial function of the town is emphasized by its large marketplace, measuring 250 by 45 m (initially probably 60 × 10 rods), which occupied a large portion of the chartered town. The first map, on which a square grid was overlaid, shows Nowy Staw (Neuteich) in 1802 after the fire. Figure 7 shows the process of creating a square grid using the HGIS Fishnet tool. The first step involved the determination of the study area whose shape would resemble the marketplace and include the whole town presented by the map (Figure 7a). Next, the tool was launched, with 1 Chełmno rod, used in Silesia and Pomerania in the 14th century, as the cell size [25]. A buffer equal to the cell size was created around the study area to ensure the inclusion of the whole area in the grid and allow small corrections of its location. Then, a square grid of adopted cell size was created within the buffer (Figure 7b). Additionally, the same tool was used to overlay a second grid, with a cell size 10 times larger than 1 rod, over the study area (Figure 7c). The square (cell size) of 10 rods was called 'w˛ezysko' (snakepit). Overlaying two grids over ˙ the map allowed for the notation of significant elements of spatial planning, such as the marketplace, an urban block, and the streets (Figure 7d).

This process was repeated for the oldest preserved city plan from 1711. The aim of this operation was to compare the layout and measurement of Nowy Staw before the great fire of 1802 (Figure 8a) and after the regulation of the town following this event (Figure 8b).

An analysis of the oldest preserved plan did not require the recreation of grids and tedious drawing. The digitized and calibrated maps had only to be added to the ArcMap project. By comparing the plans of Nowy Staw from the years 1711 and 1802, it was determined that the earlier plan showed that the layout was not perfectly plotted, so the intended courses of streets, the shape of the marketplace and urban plots was not preserved throughout the existence of the town. This was probably influenced by the local topography, the valley of the Swi˛eta river, the existence of Nytyska Wie´s that preceded the ´ town, and the pre-charter communication system. As a result of the regulation, the layout of Nowy Staw was cleaned and made more geometric, but the primary elements did not vanish (Figure 8a,b).

**Figure 7.** The use of HGIS tools. Example of Plan von der Königlichen Immediat Stadt Neuteich (Town plan) 1802: delimitation of the study area (**a**), creation of grid (cell size 1 rod) (**b**), imposition of a modular grid (cell size 10 rods) (**c**), zoom to the study area and searching for the principles of the modular layout of the city (**d**).

**Figure 8.** Comparison of the grid created using HGIS Tools for the Gründ Riss der Stadt Neuteich (Town plan) 1710, before the fire of year 1711 (**a**) and Plan von der Königlichen Immediat Stadt Neuteich (Town plan) 1802, after the regulation of the city (**b**).

Elbl ˛ag is an interesting example of an old chartered city. It received the Lübeck law in 1246, which was often given to ports. The spatial layout of such cities, in which the street network was perpendicular or almost perpendicular to the main commercial street or the market square from the port wharf is called a grate layout [24]. Elbl ˛ag is an interesting example in which the axis of the streets and urban plots is inclined in relation to the primary route at an angle other than 90 degrees. Such a layout of the streets and plots makes using a regular square grid practically useless. The HGIS Fishnet Rhombus tool developed by the authors allowed for plotting a differently shaped grid, namely, a rhombus grid. A comparison of the traditional and new grids is shown in Figure 9. Despite graphical differences between the maps, the analysis of the module measurements has confirmed that a study conducted using GIS tools is similar to the one conducted using traditional methods.

The HGIS Tools square grid method can be successfully applied not only to analyzing archival town plans but also to modern cartographic materials, such as orthophotomaps. This is shown using a series of plans of Grudzi ˛adz, founded on the riverside slope of the right bank of the Vistula river (Figure 10). The city was most likely founded in the mid eleventh century, and chartered in 1391 under Chełmno law. To this day, it has preserved much of its initial, regular spatial layout.

**Figure 9.** Comparison of the research results conducted using traditional method (**a**) [71] and HGIS tools, example of Plan de la Ville d'Elbing (Town plan) (**b**).

**Figure 10.** An example of the possibility of applying HGIS tools. Grid imposed on archival Grundriss der Stadt Graudenz Anno 1772 (Town plan) (**a**), Town plan (Grudzi ˛adz) from the year 1897 (**b**), and a contemporary orthophotomap (**c**).

#### **4. Discussion**

The aim of the study was to empirically verify the usefulness of the created tools for the measurement analysis of the historical cartographic materials. The analysis was carried out using the HGIS Tools toolbox. The conducted research confirmed the usefulness of HGIS Tools for generating the grid for morphological and morphometric studies for the example of selected urban centers. Grids of squares are widely used in the analysis of urban areas [27,35]. The use of grids of excavation squares is among the basic research tools in archaeology [31]. Various types of square grids are also used in other analyses, such as studies of the spatial development of cities [72–78], cities' landscape change and land use [79–81]. They are also used in studies on urban climate, build-up areas [76,82–85] and on issues related to the quality of life of inhabitants [86,87] or population density [88]. In most of the cases cited, grids were created using GIS tools, including those available in ArcMap software. Not surprisingly, these tools have come into widespread use also on

HGIS grounds, including urban studies [89]. It should be noted that the progress in the computerization of research processes may facilitate the morphological and morphogenetic studies of urban centers. In recent years, Historical GIS appeared as a new spatial analysis methods, utilizing GIS tools, as well as spatial and historical tools [90]. New technologies allow for improvements of the traditional research methods, introduce brand new data processing methods, and even expand current research perspectives. According to Gregory et al. [91], there are three main benefits of GIS for historical research. Firstly, spatial locations of objects can be used to organize the data and different datasets in various ways. Secondly, data can be visualized in a variety of ways, using conventional maps or different technics. Thirdly, spatial analysis tools can be used to investigate the data. On a larger scale, GIS tools have been used in the study of metrology, urban spatial layouts and archaeology since the early 1990s [11,31,35,37,38,92,93]. The application of GIS in studying the morphology and morphometry of settlements will surely not solve all methodological problems [5,6,94] but may significantly facilitate analyses and create an opportunity to combine different approaches on a common plane.

The biggest advantage of the proposed HGIS Fishnet tool is the ability to easily and quickly adjust the grid cell size to the calibrated town plan, whereas traditional methods force the necessity of drawing separate grids for plans of different scales. The tool also provides custom alignment of the grid to the area of study and the ability to automatically convert historical measurement units to modern ones and vice versa, which is a considerable improvement, compared to the standard ArcToolBox 'Create Fishnet' tool. Additionally, the tool provides ease of comparison and analysis of results for many cities, and the possibility to attach and overlay several historical and modern maps (including Digital Terrain Model—DTM and orthophotomaps) on the map under analysis.

A certain limitation in the use of the HGIS Fishnet tool may be the requirement for proficient use of GIS software. However, the popularity of the use of Geographic Information Systems as research tools allows us to assume that the majority of users will not have problems with its application. A more serious problem seems to be the appearance of errors in the spatial alignment of the historical maps. Maps should be calibrated using control points and the type of geometric correction that ensures the lowest possible Root Mean Square Error (RMS) and deformation. Map calibration distortions can have a significant impact on the performance of the tools because they are unevenly distributed over the map surface [18]. This can make it difficult to size the grid cells proportionally to the size of the distortion. The values of the distortions and their impact are an individual feature of historical cartographic materials resulting from the adopted parameters of calibration. Proper preparation of source materials and taking into account the occurring calibration errors allows for the use of tools and an efficient, time-saving analysis. Another limitation in the use of tools is their incompatibility with software other than ArcMap. The tools do not work with ArcGIS Pro and opensource software, such as QGIS. However, it is possible to develop the tool and adapt it for use in other types of software.

Selected town plans were used to demonstrate the functionality of the new HGIS Fishnet tool. Thanks to the tool, it was possible to confirm information on medieval measurements used in the planning of Elbl ˛ag and Grudzi ˛adz, and to determine probable measurements for Nowy Staw. These results may indicate the utility and usefulness of the developed tool in metrological studies of towns. According to the authors, the HGIS Fishnet tool can be used for the analysis of various regular cities, dating from ancient times to the present. These can be cities from different geographical regions as well as periods of urban development, which have a regular spatial layout [27,35].

#### **5. Conclusions**

This paper presents the possible use of the HGIS Tools set developed by the authors for measurement analysis of urban spatial layouts, using the square grid method. The 'Create Fishnet' tool offered in ArcMap 10.7 did not yield satisfactory results, as it did not allow the custom alignment of the grid to the required shape and orientation of the study area, nor the ability to create a slanted grid in which the cell is not regular, but rhombus shaped. For this reason, a new HGIS Tools set was developed. The tools were tested using plans of regular chartered towns from the former monastic State of the Teutonic Order.

An analysis of comparison of features offered by the two methods of square grid plotting showed that the HGIS Tools set allows better functionality and speed. This can be particularly useful when study of a larger number of towns or cities is needed. Any number of grids with the desired cell size may be applied to city plans. A comparison of results of layout and measurement studies using GIS methods should also be considered to be easier and faster. The primary drawback of the tool proposed by the authors is the need to familiarize oneself with the GIS software, even to a minimal extent. Especially important is the knowledge of the editing tools used to create a detailed area of study to be filled with a grid of basic fields.

Quick creation of multiple grids was also tested, as well as their application to any calibrated map in different configurations. The usefulness of these tools is proven by tests conducted on town plans that were already studied using traditional square grids.

In conclusion, it should be noted that the use of traditional methods, as well as the new version of the tool proposed here to create square grids is just one of the elements that support the process of analyzing the layout, measurements and morphology of medieval or any other cities. As mentioned at the beginning of the paper, the methodology used in historical geography, based on analyses of archival town plans, should be confronted with historical sources and verified through archaeological research.

The authors hope that the proposed HGIS Tools will be useful in modern research on urban layout, measurements and morphology. The presented paper is an introduction to further research dedicated to the metrological analysis of a larger group of cities, using the proposed tool.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/article/ 10.3390/ijgi10080558/s1. HGIS Tools box form ArcMap 10.7. The set consists of two scripts: HGIS Fishnet (Fishnet regular.py) and HGIS Fishnet Rhombus (Fishnet rhombous.py).

**Author Contributions:** Conceptualization, Łukasz Musiaka and Marta Nalej; methodology, Łukasz Musiaka and Marta Nalej; software, Marta Nalej; validation, Łukasz Musiaka and Marta Nalej; formal analysis, Łukasz Musiaka and Marta Nalej; investigation, Łukasz Musiaka and Marta Nalej; resources, Łukasz Musiaka and Marta Nalej; data curation, Marta Nalej; writing—original draft preparation, Łukasz Musiaka and Marta Nalej; writing—review and editing, Łukasz Musiaka and Marta Nalej; visualization, Łukasz Musiaka and Marta Nalej; supervision, Marta Nalej; project administration, Łukasz Musiaka; funding acquisition, Łukasz Musiaka and Marta Nalej. Both authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Publicly available datasets were analyzed in this study. These data can be found here: https://polona.pl/ (accessed on 15 October 2020), http://atlasmiast.umk.pl/ atlasy/grudziadz/ (accessed on 16 October 2020), https://mapy.geoportal.gov.pl (accessed on 10 January 2021).

**Acknowledgments:** This work was supported by Faculty of Geographical Sciences at University of Lodz.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


## *Article* **Coupling Historical Maps and LiDAR Data to Identify Man-Made Landforms in Urban Areas**

**Martino Terrone <sup>1</sup> , Pietro Piana <sup>2</sup> , Guido Paliaga 3,\* , Marco D'Orazi <sup>4</sup> and Francesco Faccini 1,3**


**Abstract:** In recent years, there has been growing interest in urban geomorphology both for its applications in terms of landscape planning, and its historical, cultural, and scientific interest. Due to recent urban growth, the identification of landforms in cities is difficult, particularly in Mediterranean and central European cities, characterized by more than 1000 years of urban stratification. By comparing and overlapping 19th-century cartography and modern topography from remote sensing data, this research aims to assess the morphological evolution of the city of Genoa (Liguria, NW Italy). The analysis focuses on a highly detailed 1:2'000 scale map produced by Eng. Ignazio Porro in the mid-19th century. The methodology, developed in QGIS, was applied on five case studies of both hillside and valley floor areas of the city of Genoa. Through map overlay and digitalization of elevation data and contour lines, it was possible to identify with great accuracy the most significant morphological transformations that have occurred in the city since the mid-19th century. In addition, the results were validated by direct observation and by drills data of the regional database. The results allowed the identification and quantification of the main anthropic landforms. The paper suggests that the same methodology can be applied to other historical urban contexts characterized by urban and architectural stratification.

**Keywords:** urban geomorphology; anthropogenic landforms; old maps; contour lines; Genoa

#### **1. Introduction**

Detailed knowledge of landforms, both natural and human, contributes to the understanding and awareness of geo-hydrological risk [1–7]. However, these landforms are often no longer immediately identifiable, particularly in densely urbanized contexts, where they have been progressively obliterated by urban growth [8,9]. This is the case of Mediterranean cities, characterized by ancient urban history, where the identification of the various phases of landscape modification is crucial for urban planning. This process helps us to identify hidden or neglected risks due to diverted or covered streams and superficial landforms (slope cuttings, retaining walls) [10–17].

Urban geomorphology [8,14,18,19] is a new discipline aimed at investigating and quantifying urban landscape changes under a historical perspective, for instance, by comparing multi-temporal cartographies in order to identify landforms modifications and volumes' shape changes through centuries.

However, while precision and cartographical accuracy of a modern cartographical product are defined by standardized and known validation processes, the same concept does not apply for historical cartography, since map making techniques have changed through time and often lack precision and accuracy [19–21].

**Citation:** Terrone, M.; Piana, P.; Paliaga, G.; D'Orazi, M.; Faccini, F. Coupling Historical Maps and LiDAR Data to Identify Man-Made Landforms in Urban Areas. *ISPRS Int. J. Geo-Inf.* **2021**, *10*, 349. https:// doi.org/10.3390/ijgi10050349

Academic Editors: Motti Zohar and Wolfgang Kainz

Received: 31 March 2021 Accepted: 14 May 2021 Published: 18 May 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

This makes it necessary to study their evolution and production techniques [21–24]. More particularly, maps produced between the late 18th and the early 19th century progressively lost their links with landscape art. Increasingly, they were produced using scientific surveying techniques, which entailed astronomic studies as well as geodetic and topographic field observations [25]. At the same time, these maps depict landscape settings of the pre-industrial age, before urban sprawl and substantial landforms modifications took place. In this period, the map is conceived as a mathematical representation of the landscape, the adoption of contour lines being a particularly meaningful example [26].

Traditionally associated with the military field, maps went through substantial technological improvements in the 18th century, when a new cartography based on modern instruments and geometric rules was increasingly functional to national states. The first systematic mapping of France on a scale of 1:86,400 using a measuring apparatus was carried out during the 18th century by members of the Cassini family [27]. The British Ordinance Survey undertook national mapping from 1791 [28], and topographical drawing and mapping were an important part of army and navy officer's education in 18th-century military schools. Cartography was also functional to the establishment of overseas colonial empires in the 18th and 19th century, and to the modernization of newly conquered areas, as was the case of Italy in the Napoleonic period [29]. The methodology of production was used and improved until the 1980s when significant adoption of digital techniques in data acquisition has emerged [30,31]. Particularly valuable examples of precise geographical representation in Italy are detailed cadastral plans of some pre-unitary states, including the Grand Duchy of Tuscany, the Kingdom of Lombardy-Venetia, the Sardinian Kingdom, and the Kingdom of the Two Sicilies [32,33], whose astronomic observatories allowed optimal geodetic measurements, similarly to those in northern and central Europe states.

In Liguria (north-western Italy), the first example of fairly accurate cartographical representations are the maps by Matteo Vinzoni, appointed official surveyor of the Commissari della Sanità of the Republic of Genoa in 1722. By using the alidade and the plane table, Vinzoni depicted the whole Ligurian coast in his "*Pianta delle due Riviere della Serenissima Repubblica di Genova divise ne' Commissari della Sanità*". In addition to representing one of the first systematic surveys of an entire state in the Italian context, Vinzoni's work is used as one of the sources to monitor coastal dynamics in Liguria [13,32]. The main cartographical product of early 19th-century north-western Italy was the *Carta degli Stati Sardi in Terraferma di S. M. il Re di Sardegna*, divided into 91 sheets on a scale of 1:50.000, of which 86 were published between 1851 and 1859 [33,34]. The map is the product of accurate field surveys carried out in the previous decades by the Piedmontese cartographers, who produced a series of 112 highly accurate manuscript maps (1:9'450, 1:10'000, and 1:20'000) commonly called *Minute di Campagna*.

With the unification of the Kingdom of Italy, the newly founded Istituto Geografico Militare (IGM) (1861) carried out a complete mapping of the entire Italian territory at various scales (1:10'000; 1:50'000; 1:100'000) using the technique of contour lines for the depiction of elevation [35]. In order to cover the entire Italian territory, the Bessel ellipsoid was adopted as datum and the reference median was that of Monte Mario 12◦27008,40000 East of Greenwich.

This work aims to assess volumetric changes of anthropic landforms in an urban environment in the last two centuries in the city of Genoa (NW Italy) through computational analysis and cartographical comparison. This time frame is considered appropriate, as maps produced in the early 19th century are the oldest available documents for a scientific reconstruction of pre-urban sprawl morphology. This means that the earliest mapped anthropic landforms represent the total amount of pre-existing elements as a result of the stratification of several urban phases, particularly of the Middle Ages and Modern Age.

In particular, this paper could be divided into two targets:


#### **2. Materials and Methods**

#### *2.1. Study Area*

Genoa is a port city of the Mediterranean and capital of Liguria (north-western Italy). The municipality, whose current size dates back to the 1926 "Grande Genova" project, is 240.29 km<sup>2</sup> wide and the coastline 35 km long. The most ancient core of the city, the yellow area in Figure 1, was smaller (c. 11.45 km<sup>2</sup> ) and was roughly located between the city's two main streams: the Bisagno in the eastern part, and the Polcevera in the West. From a geological point of view, the study area is largely constituted by marly limestones, marls, and clays in Flysch sequences of Upper Cretaceous, while the western part, in the Polcevera Valley, is dominated by shales of Lower Cretaceous [36–38]. The two formations are found both as outcrops and bedrock. Another important lithotype, today largely hidden by the urban structure, is constituted by stiff fissured Pliocene clays, exclusively found in the historical center [37–42]. Morphology is characterized by many small and steep catchments, mostly oriented orthogonally to the coastline, originating from the neotectonic activity [42].

**Figure 1.** (**A**) Study area and areas of interest. (**B**) The *geological sketch map* box: (1) waterfront embankments; (2) alluvial deposits; (3) slope deposits; (4) stiff fissured clays; (5) shales; (6) marly limestone with clayey shales interlayers.

The first permanent settlement dates back to the pre-Roman period. In the Middle Ages, Genoa developed around the old city core with an East–West direction along the coastline [43]. Urban sprawl over the hills, with a South–North direction, and along the coast developed in more recent times [44] (Figure 2). In the early 19th century, the two main valleys, the Polcevera to the West and the Bisagno to the East, were characterized by an agricultural landscape with a dense pattern of small fields across the floodplain and terraces over the slopes with a prevalence of vineyards and olive orchards [45,46]. Terraces are only partly visible in the current landscape and have been largely obliterated by more recent urban structures, decametric retaining walls, and high blocks of flats of the second half of the 20th century [47]. Urban expansion brought hydraulic works in the two valleys, with progressive narrowing of the riverbeds and the establishment of new

urban areas in the floodplains in the early 20th century up to the 1930s [12]. It is estimated that loss of agricultural surfaces between 1820 and 2014 in the Bisagno Valley was 57%, while woodland and new urban areas increased respectively by 248% and 750% [45]. In the second post-war period, the Promontory of San Benigno, exploited for the extraction of marly limestone as building material from the Middle Ages, was completely excavated. The complete removal of San Benigno Promontory allowed territorial continuity between Genoa center and its western suburbs, with intense urban saturation of many areas, including the former rock quarry [34].

**Figure 2.** Urban sprawl in western Genoa (**a**) Vinzoni (1773) "Il dominio della Serenissima Repubblica di Genova in terraferma (Riviera di Ponente)"; (**b**) Carta degli Stati Sardi in Terraferma di S. M. il Re di Sardegna "Minute di campagna scale", scale 1:9,450 period 1815–1823; 1:25,000 scale topographic map of the Italian Military Geographic Institute of 1878 (**c**) and 1934 (**d**).

#### *2.2. Materials*

Analyzed material can be divided into two subgroups: historical cartography and current cartography.

The most important document is the "*Carta Generale di Difesa di Genova*" by engineer Ignazio Porro (1848). For the purpose of this research, we analyzed the digitalized version of black and white photographs of the 74 original plates (scale 1:2'000) of the "*Carta Generale di Difesa di Genova*" depicting a large part of the current urban area. Porro, who was engaged by the Kingdom of Sardinia for the production of the map, can be considered as a pioneer of topographical surveying and is the inventor of the Porro prism for binoculars.

The map represented with high detail the system of military fortifications and the morphology of the Genoa area, providing topographical information for urban planning aimed at improving the military defenses of the city. Starting from data acquired by Napoleonic topographers and unpublished documents by astronomer Baron De Zach (1816–1817) and with the help of engineer Giulio D'Andreis, Ignazio Porro carried out a new triangulation using mountains around Genoa as reference points and applying for the first time in Italy the tacheometry technique, which he appositely improved for the purpose of this project [48].

The final work consists of 74 plates on a scale of 1:2'000 with a reference system related to the Lanterna (Lighthouse of Genoa Harbour), with 10-m-interval contour lines.

This cartographical product is one of the first Italian examples of topographical representation with contour lines. Overall, it is a metric planimetry based on relative coordinates and with the North direction coaxial to the plates' vertical squares

The area analyzed in this research is in the central part of the city and the coast nearby and is covered by 7 of the 74 plates produced by Porro. These are plates n◦ 42, 52, 53, 62, 63, 72, and 73, according to Porro's numbering system.

Notably, the analysis focuses on particularly interesting sections from a historical and urban geomorphological point of view (see Figure 1): (I) Sant'Agata Old Bridge (included within sheets n◦42 and 52); (II) Circonvallazione a Monte (included within sheets n◦52, 53, 63); (III) San Benigno Promontory and via Digione area (both included within sheets n ◦ 62 and 72; and (IV) Polcevera Viaduct area (i.e., Morandi Bridge area), entirely included within sheet n◦ 73.

Cases I (Sant'Agata Old Bridge) and IV (Polcevera Viaduct) are located in alluvial plains respectively of the Bisagno and Polcevera Streams, while II (Circonvallazione a Monte) and III (San Benigno Promontory and via Digione area) are hillside case studies, very significant for the substantial urban and morphological modifications they went through during the years, particularly due to rock excavation.

Modern cartography consists of sheets of the Regional Technical Map of Liguria on a scale of 1:5'000 ("*Carta Tecnica Regionale"* ed. 2007), digital terrain models (DTMs) with a 1-m resolution and equivalent scale of 1:1'000, and orthophotos with a 10-cm resolution. The last two constitute the LiDAR aerial survey that Genoa Municipality commissioned in September/October 2018 across the municipal area (Figure 3).

**Figure 3.** (**a**) Porro Plate n◦52. (**b**) The same area today, overlay of Regional Technical Map and DTM LiDAR.

In addition, stratigraphic information from the database of Regione Liguria [49] and boreholes promoted by the *Special Commissioner for the reconstruction of Polcevera viaduct*, after the Morandi Bridge disaster [50], were used. These are fundamental data sources in order to reconstruct the deepness of the contact interface between landfill material and natural rock/sediment. These point data were integrated with the digitalized contour lines and are the input for interpolation calculation.

#### *2.3. Methods*

The methodology is summarized in Figure 4. Two different GIS software were used, QGIS and GRASS, while the accuracy of the old cartography was defined with MapAnalyst.

QGIS was used for the georeferencing and digitalization of contour lines, while GRASS was used to generate the digital elevation model (DEM) and for statistical and algebraic analysis of the difference of DEMs (DoDs) system [51–54].

The analysis initially seeks to define the intrinsic quality of maps, considering that the accuracy error of old maps on the x-y plane (*RMSEH*) is given by the sum of various factors summarized in the following relation (1):

$$RMSE\_H = \left( RMSE\_s^2 + M \cdot RMSE\_m^2 + M \cdot RMSE\_d^2 \right)^{0.5} \tag{1}$$

where *RMSE<sup>S</sup>* is the error introduced by the ground survey, *RMSE<sup>m</sup>* is the map error, *RMSE<sup>d</sup>* is the digitization error, and *M* is the map scale factor [55,56].

As Porro himself described in his relation [57], topographical surveying was affected by an azimuthal instrumental error (RMSE) of 0◦ 0' 70'', which meant a real deviation of 0.33 m on the x-y plane at a distance of 1 km.

Moreover, digital images are the product of a process of paper maps photography that has very likely produced perspective distortions (*RMSE<sup>d</sup>* ). In addition, it has been observed that original images show deformations also due to paper deterioration processes (*RMSEm*).

#### - Accuracy analysis with MapAnalyst

The first passage entailed a monitor comparison between the old map and the current orthophoto/topographical map aimed at identifying the perimeters of medieval, Renais-

sance, and 18th/19th-century palaces still found in Genoa's urban fabric. Once reference points of both maps were found, ground point control (GCP) positioning was carried out. These are common points between plates of the Porro Map and the current cartography, using firstly orthophotos and secondly OpenStreetMap only as auxiliary images for references of place names.

MapAnalysts simulates a georeferencing process, which consists in imposing a geographical reference system known to the target image. This dislocation causes variable image distortions depending on the desired reference system. However, the output provided by MapAnalyst is not a ground control point (GCP) but a statistical basis on the old map accuracy compared to the modern one.

For this analysis, three geo-referencing algorithms, available in MapAnalyst, were applied: *Linearised Helmert*, the *Affine transformation with 5 parameters*, and the *Affine transformation with 6 parameters*.


Once the statistical basis was obtained, the same procedure was applied using QGIS, identifying for each image the same configuration and the same number of GCPs allocated in MapAnalyst.

Images in QGIS were georeferenced with the linearized Helmert algorithm. This algorithm does not distort an old image; rather, it roto/translates it, minimizing the deviation between the old image and the geo-referenced one. The adopted reference system is Gauss-Boaga Roma 40 (EPSG 3003).

In general, it is possible that some points might be higher than RMSE. In this case, a threshold value was arbitrarily decided, slightly higher than the error estimate under which the variance between Porro and modern topography can be considered acceptable.

In cases where the distortion was higher than the threshold, additional GCP were assigned by iteration until the desired effect was reached.

This was done for each single map section, which falls within a single or two or more plates. The process allowed us to minimize the distortion error as each plate is separate and can be geo-referenced singularly.

Subsequently, contour lines were digitalized through a linear geo-vector by assigning the relevant elevation above sea level. Once this operation was completed, the vertexes of these isolines were extracted.

In the hillside area (Circonvallazione a Monte, San Benigno Promontory, and Via Digione), this operation was sufficient and complete.

In floodplain areas (Sant'Agata Old Bridge and Polcevera Viaduct) where contour lines are not frequent or not visible due to the presence of buildings, another data source was represented by borehole stratigraphy (Figure 5) [49,50].

In particular, the ancient interface surface between landfill material and natural rock/sediment was considered, but only in cases where the borehole stratigraphy indicated post-industrial landfills that included traces of cement, industrial bricks, etc.

Moreover, since riverbed areas in the Porro Map lack spot elevations, it was decided to conceptually match the old topographical surface with the current situation. This allowed us to bind the interpolation to the thalweg in order to maintain the morphology of typical cross-sections of Ligurian riverbeds.

Elevations of the current riverbed were extracted from LiDAR and assigned to the old riverbed as spot elevations. In this way, the same vector file gathered elevations extracted from contour lines: at the interface between landfill material and natural sediment/rock substrate and in the two riverbeds, bearing in mind that no significant natural process occurred, such as severe erosion or intense sediment accumulation due to flooding, that changed the floodplain terrace elevations between Porro's survey epoch and the late 1800s, when intense urban sprawls began, unlike the phenomena described in [58–60].

**Figure 5.** An example of a borehole with translated stratigraphy used for this research.

#### - Interpolation and DoDs with GRASS GIS

Since hand-drawn graphical elements of the Porro map have an equivalent resolution of 2 m (1 mm to the naked eye on a scale of 1:2'000), previously obtained elevation data were interpolated on a raster map with a 1-m resolution, also to reduce potential vertical errors due to horizontal offset between the old and current map, according to a methodology explained by Stoker et al. [61].

This was carried out in GRASS GIS through three different steps: (1) the *v.surf.rst* command [62], an algorithm based on spline with smooth and tension, used for preliminary interpolation that allowed us to convert elevations into a first raster dummy surface; (2) from this surface, a vector points cloud with casual location was extracted with the *r.random* command; (3) these points were interpolated a second time with the *v.surf.rst* command, obtaining a smooth surface, not affected by classic interpolation errors, such as quad-three, gradient gaps, etc. [63].

Moreover, once output raster maps of the five case studies were obtained, an assessment on vertical unreliability was carried out according to the following relation [54]:

$$
\varepsilon\_{\mathbb{Z}} = \varepsilon\_H \cdot \tan \llcorner \tag{2}
$$

where *ε<sup>Z</sup>* is the elevation error caused by horizontal offsets on a sloping surface, *ε<sup>H</sup>* is the horizontal error, and α is the local slope.

In slopes characterized by α < 45◦ gradient, potential vertical offsets between the old and current cartography are less than the horizontal one. In contrast, if angle α > 45◦ , the quantification of the vertical offset is crucial as this is higher than the horizontal one.

The maximum vertical offset RMSEZij MAX [64] was calculated, corresponding to the cell with the highest slope degree, according to relation (3), in order to have full knowledge of the limits of the extracted output:

$$\text{RMSE}\_{\text{Zij}\,\text{MAX}} = \text{RMSE}\_{H} \tan \text{ S}\_{\text{ij}\,\text{MAX}} \tag{3}$$

where *RMSE<sup>H</sup>* is the horizontal error for the map, and Sij MAX is the maximum slope in a i-row and j-column grid cell of the output raster map.

Furthermore, sea level rise in the period 1884–2009 [65] was considered in order to quantify the heights' vertical offset between modern and 1830s sea level measurements.

The resulting old DEMs were subtracted from the 1-m-resolution 2018 LiDAR of Genoa Municipality, obtaining the various DoDs, where positive areas correspond to fillings and negative ones to excavations [15].

#### **3. Results**

The methodology adopted produced the following results:


From a first modelling process through *MapAnalyst*, it was possible to quantify potential errors by applying the linearized Helmert transformation, the affine transformation with five parameters, and the affine transformation with six parameters.

It was observed that the accuracy in old cartography increases in consolidated urban areas; this is due to the higher number of reference points that Porro and his team could use, compared to hilly areas. This also explains why the number of GCP couples is higher in plates of urban areas. This makes them more reliable than plates of rural areas characterized by less GCP couples. The GCP minimum number for each couple of images is 16 as shown in Table 1, alongside the error estimations for each plate.

**Table 1.** Georeferenced plate numbers with standard deviation, RMSE computed with various algorithms, and the number of used GCP.


An acceptable deviation threshold of 10 m on a x-y plane between old and current cartography was arbitrarily established, a slightly higher value than the RMSE (Affine 5) of sheet n◦73 (the least accurate and located in areas with no urban settlements at the time of Porro).

Interpolation obtained through *spline with smoothing and tension* at a 1-m resolution increases the smoothness of old topographical surfaces, generally characterized by a lower gradient than current ones. In addition, by interpolating contour lines it was possible to follow riverbeds and/or erosion landforms of ancient, often ephemeral, streams, now completely obliterated by anthropic landforms.

The total surface of the calculated raster maps is 1.79 km<sup>2</sup> , of which only 0.2% have a gradient > 45◦ , and the cell with the highest gradient has a value of α = 84◦ so that RMSEZijMAX = 28.02 m, obtained through relation (3) with RMSE<sup>H</sup> = 4.09 (Sheet 72, Helmert in Table 1).

In the remaining part (99.8%), the vertical offset value is less than the horizontal one.

The case studies, all less than 2.5 km from the coastline, lie within 150 m asl. In particular, the Sant'Agata Old Bridge and Polcevera viaduct areas are respectively included between 5.2–60.3 m asl and 4.92–82.89 m asl, while via Digione and San Benigno Promontory are located respectively between 18.02–111.69 m asl and 10.02–92.41 m asl. Circonvallazione a Monte lies between 57.47 and 130.74 m asl.

The heights' vertical offset between modern and 1830s sea level measurements is about 25 mm with an average rate of 1.1 mm/year. Considering an average slope of 16.94◦ for the old DEM hilly areas, as obtained only by the Porro's isolines interpolation, an offset

of 0.25 m leads to an average x-y shifting of 0.9 m, which may be considered as a negligible value.

For what concerns the DoDs, outputs can be divided into three different macro-groups. The first group is characterized by a high frequency of positive differences. This is the case of floodplain areas, including Sant'Agata Old Bridge (Figure 6) and Polcevera Viaduct (Figure 7). As already discussed, the preliminary assumption that old and current surfaces of Bisagno and Polcevera riverbeds precisely overlap means that the elevation difference in the DoDs equals to zero. Instead, there is a significant difference for what concerns current artificial embankments, which is maximum in their immediate proximity and tends to decrease with distance, reaching almost zero values near the foothill. The average thickness of the alluvial plain is around 2.45 m.

**Figure 6.** Sant'Agata Old Bridge: (**a**) Old surface; (**b**) DoD with the medieval and modern trace.

**Figure 7.** Polcevera Viaduct: (**a**) Old surface; (**b**) DoD.

The case of Polcevera Viaduct is even more complex due to the numerous anthropic landforms (*sensu* Rosenbaum) (i.e., railway, industrial sites, etc.).

Along the slopes, local fillings with differences between 10–15 m are found.

The second group is characterized by a greater occurrence of negative differences. This is the case of via Digione and San Benigno Promontory (Figure 8). Except for slight positive differences of c. 2–3 m in the slope's upper side, significant negative differences, around 60–75 m, are noticed.

**Figure 8.** San Benigno Promontory and Via Digione: (**a**) Old surface; (**b**) DoD; (**c**) Via Digione cross-section; (**d**) San Benigno Promontory cross-section.

The third group is represented by Circonvallazione a Monte (Figure 9). In this case, it was not possible to identify areas characterized by a high frequency of positive and negative values. Amongst the negative differences there are values below 20 m, while the positive ones can reach 12–15 m. Compared to previous cases, areas with values around zero in Circonvallazione are rare.

**Figure 9.** Circonvallazione a Monte: (**a**) Old surface; (**b**) DoD.

#### **4. Discussion**

The analysis shows how an average offset < 10 m between the current cartography and the "*Carta Generale di Difesa di Genova*" can be reached by repeating the geo-referencing process. This value can be considered acceptable to distinguish between natural and anthropic landforms at large-scale maps (1.10'000–1.5'000).

The difference between anthropic and natural landforms, however, is only visible for the period between present day and the second half of the 19th century. Anthropic landforms before this period can be easily confused with natural morphologies. For this reason, the Porro cartography and derived products can be integrated by using older cartographies, although these can lack geometrical accuracy. In addition, paintings, topographical views, and historical photographs can offer insights into past landscapes [66].

For instance, the raster of Figure 10a shows a gently sloping area along the left side of the stream with a regular pattern and constant NE-SW direction. From the analysis of the Porro cartography, it is not possible to establish whether this is a natural morphology or a filling landform of the period before 1848 [15].

However, the painting of Figure 10b, made in 1853 from previous *en plein air* sketches, shows that this can be a filling of the floodplain area of the Polcevera Stream confined by a consistent stone retaining wall. Here, starting from 1850, various railways lines connecting Genoa to northern Italy were established. The section of Figure 7b, near the viaduct, shows that the thickness of the two railway embankments is around 5 m.

The area of Sant'Agata Old Bridge is also characterized by dense urbanization, although the evolution of anthropic and natural landforms is less complex. Here, a consistent filling dating back to the period between late 19th and early 20th century is found, related to the Bisagno Torrent narrowing. The bridge of Sant'Agata (Figures 6b and 11a) was made in the medieval period along the track of the old Aurelia road and consisted of 28 arcades for a total length of 360 m. Currently, the new bridges are only 50 m long, providing evidence of the significant narrowing of the Bisagno Stream riverbed [12]. Alongside Polcevera Stream, the Bisagno is currently one of the mostly exposed to geo-hydrological risk Italian rivers. Their narrowing is one of the reasons of the frequent flood events of the two streams. In the 1950–2015 period, 5 floods with discharge between 700 and 1000 m3/s occurred [12,45,67].

**Figure 10.** (**a**) The output of Polcevera viaduct, in purple dotted line the Campasso area; (**b**) Left bank, elevated area with embankment sustaining wall in the Campasso area (view by Carlo Bossoli, 1853).

**Figure 11.** (**a**) Sant'Agata Bridge in the background (painting of Giuseppe Bisi, c. 1825); (**b**) Ancient Genoa: outside the city walls, over the hills, the site where Circonvallazione a Monte currently lies was characterized by a rural landscape (painting, unknown author, c. XV-XVII century); (**c**) San Benigno Promontory not yet excavated and the Lighthouse (*La Lanterna*) in the background; (**d**) Via Digione landslide (1968).

During intense rainfall events (>50 mm/h, [44,45]), minor sink-holes due to collapses of culvert-arches of underground streams took place [42], something that is also very likely to occur in the Circonvallazione a Monte area.

Circonvallazione a Monte crosses an area characterized by 11 minor catchments less than 3 km<sup>2</sup> in size [42]. As shown by Figure 9, filling and excavation anthropic landforms analyzed with DoD are quantitatively commensurate to current topography, constituted by excavations and significant decametric retaining walls (Figure 11b). These walls, built on limestone bedrock, are not monitored nor is their deterioration state known. Local building collapse episodes due to the subsidence of significant landfill materials thickness is indicated in small valleys that have been filled since the second half of the 19th century [42].

However, the most significant excavation case is San Benigno Promontory, where an entire slope was removed. Only considering the southernmost part of the slope, from the excavation face edge, more than 10.79·10<sup>6</sup> <sup>m</sup><sup>3</sup> of rock were removed between the late 19th century and 1960 (Figure 11c). This value was obtained by DoD in Figure 8, excluding the area of via Digione [68]. Here, another limestone quarry, completely excavated and characterized by the establishment of new buildings in the post-war period, was affected by a landslide in 1968, with 19 fatalities (Figure 11d).

Such a product, in addition to having multidisciplinary scientific interest, is also a useful tool for stratigraphy analysis of the buildings' foundation levels and for the production of civil protection maps [69] or multimedia products aimed at raising public awareness on geo-hydrological risk in urban contexts [70–72].

#### **5. Conclusions**

Digital image acquisition of analogic photographs of the "*Carta Generale di Difesa di Genova*" map; its error analysis, georeferencing, digitalization, and contour lines interpolation; and DEMs difference between the derived product and the current DTM (2018) allowed us to develop a specific methodology useful for the study of urban geomorphology. The analysis demonstrated that this specific old cartography can be compared with modern large-scale topographical maps. It provides evidence of the potential use of such historical material to identify and quantify anthropic landforms in urban geomorphology survey and mapping. In particular, it is worth mentioning cadastral maps with no elevations or contour lines and old views depicting past landscapes.

This analysis involved five sample areas that were compared with the topography of Ignazio Porro. The result can be divided into three different categories:


Potential further research can involve broader geographical areas and time frames; from a size point of view, the whole area of Genoa can be covered using the 77 plates made of the Porro map. Assessments of the city's morphometric changes can be carried out by analyzing more recent maps, for instance, geo-referencing and interpolating the IGM (Istituto Geografico Militare) 1:25'000 scale maps [14,71] produced between the late 19th and early 20th century in order to gain a precise city evolution with a 30-year interval.

This approach can be used for other cities that, like Genoa, are characterized by urban and architectural stratification worthy of being discovered and studied. This process would necessarily entail the use of a sufficiently old map or set of maps, which would depict with accuracy the landscape before urbanization took place.

**Author Contributions:** Conceptualization, Francesco Faccini and Martino Terrone; methodology, Martino Terrone and Guido Paliaga; software, Martino Terrone and Guido Paliaga; validation, Pietro Piana, Francesco Faccini and Marco D'Orazi; formal analysis, Martino Terrone; investigation, Martino Terrone and Francesco Faccini; resources, Martino Terrone and Francesco Faccini; data curation, Pietro Piana, writing—original draft preparation, Martino Terrone and Pietro Piana; writing—review and editing, Pietro Piana and Guido Paliaga; visualization, Marco D'Orazi; supervision, Francesco Faccini; project administration, Martino Terrone and Marco D'Orazi; funding acquisition, Martino Terrone and Francesco Faccini All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was financed by Programma Operativo Regione Liguria 2014-2020 - Fondo Sociale Europeo–Asse 3 "Istruzione e Formazione". RLOF18ASSRIC/60/1 "Tecnologie per il monitoraggio integrato e la mitigazione del rischio geomorfologico da frana nella gestione delle reti di sottoservizi" (Scientific coordinator: Francesco Faccini).

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **References**


MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*ISPRS International Journal of Geo-Information* Editorial Office E-mail: ijgi@mdpi.com www.mdpi.com/journal/ijgi

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34 Fax: +41 61 302 89 18

www.mdpi.com ISBN 978-3-0365-4123-5