Georeferencing Historical Maps at Scale

Historical maps are a cartographic record of where a place was and perhaps still is. They offer a window into the past that can support indigenous communities to reconnect with their histories, language and places. Specifically, historical cadastral maps provide a record of the evolution of land interests during the colonisation of Aotearoa New Zealand in the 19th and 20th centuries. Given that georeferencing is typically a manual and time-intensive task, these rich resources are often inaccessible to most except for those who are familiar with cartography or the local histories of the places that were mapped. By labelling points in the map and calculating their real world coordinates, georeferencing enables a broader level of accessibility to the map image and the histories embedded therein. We propose an algorithm Koki Tauriterite (translated simply as determining angle equality) that leverages the ease of detecting intersecting lines in the map image and digital cadaster to perform robust record linkage between the two sources for the purposes of large scale automatic georeferencing of historical maps in Aotearoa New Zealand.

For New Zealand Māori, historical maps allow one to locate traditional places of food gathering, villages and burial grounds that are no longer visible in today’s landscape due to successive generations of te muru me te raupatu (land confiscation and dispossession), where historical and contextual knowledge of said places generally exists in the minds of a few [4]. Therefore historical maps enable one to do what most cannot; locate these cultural sights of significance in space. Attempts to revitalise and make accessible this knowledge are visible around the country and the utilisation of computational tools are aiding in that process of revival, for example Ngāi Tahu’s community generated atlas tool Kā Huru Manu [11].

The processes that lead to the creation of the map in the Aotearoa New Zealand context represent a traumatic colonial history for Māori communities [4]. Power structures in the colony which sought to benefit the colonising power were perpetuated while simultaneously assimilating indigenous relationality to place via the creation of the map [12]. The role that maps played in the dispossession of Māori and their resources through the Native Land Court and other various mechanisms of the state [6] highlights a contradiction evident with the elevation of the map as an important archival record of Māori language and geographic history [5], particularly (for one of us) as a descendant of Parihaka ploughmen who were imprisoned without trial and forced into hard labour for ploughing their own confiscated land [14].

However, those same symbols of paternalistic control can also play an important role in the retrieval, preservation and eventual dissemination of Māori relationality and re-connection with place, and subsequently aid in helping to make visible that rich history to the communities from which those places belong and vise versa. Therefore, the rematriation [10] of historical maps as records of place, via the provision of wider, open access to them for indigenous communities represents a form of spatial data justice [15] that inspires and motivates this work.

The task amounts to record linkage between these two very different sources of data. Source 1 is a digitised image of a map – those of interest are not currently georeferenced. The state of the map and quality of digitisation is variable and, in addition, neither overall orientation nor scale can be assumed for the map since there are often scarce metadata records to accompany the map. This necessitates an approach that is therefore scale and rotation invariant. Source 2 is essentially a large list of polygons and associated geospatial information, henceforth “the cadaster”. The task of georeferencing a map amounts to finding plausible corresponding locations for an individual map, within the (large) cadaster. However the two sources use completely different representations in their raw data (images, and polygons respectively). Contrast this with astrometry for example, where both image and stellar catalogue naturally result in lists of 2d vectors which are able to be compared [13]. An effective solution will involve making a defensible correspondence between features identifiable in both map images and the cadaster.

Record linkage of cadastral data and historical map images presents a significant challenge due to inherent discrepancies between analogue and digital-borne data [16]. While both sources purportedly depict the same geographic area, they capture it at different points in time, resulting in variations in geometry present in both the survey record and the map image. Temporal changes in land use as well as cartographic conventions complicate the process of establishing accurate correspondences between the two sources of data [17]. In Aotearoa New Zealand, the cadaster is a chronological record of land title administration that has been recorded, digitised and made public by the government agency Toitū Te Whenua Land Information New Zealand (LINZ), a process which began in the mid-1980s [19]. This is the primary dataset that contains millions of polygons that represent all current parcels in Aotearoa New Zealand [18].

Herein lies another major challenge in conducting historical record linkage. While the current survey record reflects the present-day cadastral landscape, the historical map offers a snapshot of a past configuration of land ownership. Consequently, the LINZ Parcels Dataset contains numerous parcel records that are absent from the historical map due to subsequent partitions of new parcels. Therefore, the historical map may depict parcels that no longer exist in the contemporary LINZ Parcels Dataset. These disparities necessitate a robust record linkage approach accounting for the inevitable attrition of records across the temporal divide. Record linkage in this case being the process of matching and combining records from multiple sources into a single space and then subsequently assessing whether or not they are in fact describing the same object.

Cadastral features such as parcel boundary intersections are persistent to changes over time. Parcels of land are more often partitioned than they are amalgamated as a result of accelerated urban development in Aotearoa New Zealand during the 20th century [1]. Although the polygonal shape that is created by the parcel is generally orthogonally segmented over time, the original shape of the parcel remains identifiable in the cadaster. The original parcel may simply now have a number of lines drawn through it, for example parcel 22 in 1 is partitioned in the cadaster. These straight lines that depict parcel boundaries and the adjacent parcels they intersect with will be the basis for our analysis since there is no complete ground truth yet developed that could represent all land parcels from a particular point in time when a given map was created. Therefore a common representation space (defined below) is required to be developed upon which comparison can be performed between the two sources.

Our proposed algorithm makes use of the fact that straight lines are particularly salient in both sources (respectively via the Hough Transform [7, 3] for images and simple geometric processing for the cadaster). This supports ready identification of line intersections in either source, and characterisation by their internal angles. Agreement between sets of internal angles is useful but not diagnostic: it does not provide sufficient specificity on its own to match maps to locations. However any pair of intersections (which we call a “dyad” in what follows) furnishes another pair of angles defining their relative orientation, which provides a much more finely grained basis on which to argue for a match.

A “distance” that is low for good matches is simply the sum of squared differences in these two sets of angles, $\mathrm{\Theta }$ from the map and $\mathrm{\Phi }$ from a potential site in the cadaster:

Although intuitively appealing, Equation 1 requires some justification. In Section 3, we derive it from considering a Bayes factor for the probability of a match being correct, given the respective angles $\mathrm{\Theta }$ and $\mathrm{\Phi }$, and taking careful account of their various dependencies.

In our pursuit of an effective method for georeferencing historical map images using data from the cadastral record, we explored several approaches before developing our current line intersection-based algorithm. These preliminary investigations, while ultimately not preferred, did highlight the challenges inherent in this domain. Use of image segmentation and shape recognition techniques to match full parcel polygons is one such option, but our experience has been that these methods are not robust enough to the variability inherent to historical maps, as well as the fundamental differences in spatial object representation between the cartographer’s approach and modern GIS practices that form the current cadastral record. The discrepancies in boundary delineation and feature abstraction between these two distinct technologies rendered direct shape-matching approaches unreliable. Subsequently, we explored point-based methods, employing point-finding algorithms coupled with triangulation techniques such as the Delaunay triangulation [8]. This approach, while theoretically appealing, was also found to be unsuitable, largely due to the noise present in the map images. The arbitrary addition and removal of points caused by image artifacts and inconsistencies in feature representation led to unstable and unreliable triangulations.

We also considered leveraging toponyms extracted through Optical Character Recognition (OCR) techniques in order to conduct direct georeferencing of the historical map. As is evident in Figure 2, this approach faced limitations due to the temporal disconnect between historical and contemporary place names. Many toponyms present on historical maps have since fallen into disuse or been replaced in the process of the colonisation of Aotearoa New Zealand [2]. Text on historical maps also tends to be heavily occluded, and/or follow the geography of the feature that it is describing, introducing even more complexity into the text extraction process [9]. Nevertheless, we found value in using the more persistent toponyms, such as major town names, street names or key geographic features as a means to estimate the general location of the map and thereby reduce the search space for our dyad matching pipeline (described below). These exploratory efforts underscore the complexity of the task at hand and the need for robust, noise-resistant methods in historical map analysis and record linkage.

In either source, we refer to a single intersection of two straight lines as a “monad”, to contrast with “dyads” (which are pairs of intersections) that provide the main discriminating features on which our matching algorithm is based – see Figure 3 for an example. Where the intersecting line segment does not extend past another but instead generates a “T” junction, we can still treat this as an “X” shape in effect. Each monad has a characteristic angle, defined to be the smaller of the two internal angles. This is denoted ${\theta }_i$ if the source is the map, and ${\varphi }_i$ if it is from the cadaster. Clearly ${\theta }_i\le \pi /2$ and similarly for ${\varphi }_i$.

Internal angles generated by intersecting parcel boundaries are biased toward $90\mathrm{°}$, as is evident from Figure 4. Similarly parcel boundary intersections in the cadaster result in internal angles close to $0\mathrm{°}$, due to small line segments generated by parcels that, for instance, follow the course of a winding river or are located in urbanised areas. Furthermore, the very definition of an intersecting parcel boundary in the cadaster lends itself to complexities with respect to feature extraction since it is difficult to determine whether or not a given point in the coordinate sequence of a polygon doesn’t merely sit along a straight line. For instance in most general contexts a square shaped polygon would likely produce four unique coordinates, but in the cadaster this cannot be assumed and so two parcels could intersect along an almost straight line, producing an angle very close to $0\mathrm{°}$. The very different nature of the two datasets (the historical map image being cartographic and the cadaster being generated computationally via survey measurements) leads to different complexities with respect to feature extraction since different methods for feature are required for the different datasets. Figure 4 confirms this very skewed distribution across the entire LINZ Parcels Dataset. The higher likelihood of seeing orthogonal or close to parallel lines reduces their discriminating power for the proposed algorithm, which relies on “suspicious coincidences” of angular relationships to attribute matches between dyad pairs.

A dyad is a set of two unique monads from the same source-this induces two further angles that are relative bearings to the line joining the monads, as illustrated in Figure 3. The range for the relative angle so defined is ${\theta }_r\le \pi -{\theta }_i$. Collecting these angles then, a given dyad $D$ extracted from the map image source, generated by the two monads $A$ and $B$ consists of $\mathrm{\Theta }=({\theta }_i^A,{\theta }_r^A,{\theta }_i^B,{\theta }_r^B)$, and similarly for a dyad extracted from the cadaster (which is a possible match for $\mathrm{\Theta }$) we have $\mathrm{\Phi }=({\varphi }_i^{A^{\prime }},{\varphi }_r^{A^{\prime }},{\varphi }_i^{B^{\prime }},{\varphi }_r^{B^{\prime }})$. Monads and dyads are the primary features extracted from both sources. This extraction allows for comparison to be made between the two sources thereby acting as a common representational space for matching across the sources.

We first pre-filter candidate matches on the basis of monad evidence alone, in order to reduce the computational burden of matching each component part of the dyad. This reduces the search space for subsequent, more computationally intensive comparisons. The pipeline can be summarised as follows:

1. 1.

   Monad matching: For each monad in the map and cadaster, identify putative matches where internal angles differ by no more than $\sigma$ degrees.

2. 2.

   Relative angle calculation: For each plausible pairing of monads, compute the relative bearing angles ${\theta }_r$ and ${\varphi }_r$.

3. 3.

   Correct for edge cases: Apply Algorithm 1 to ${\theta }_r$ and ${\varphi }_r$.

4. 4.

   Dyad comparison: For each candidate dyad pair, compute a similarity score $d(\mathrm{\Theta },\mathrm{\Phi })$ based on all internal and relative angles.

5. 5.

   Best match selection: Select the dyad match with the lowest distance score as the most likely correspondence.

The first step filters on the ${\theta }_i$ angles of each monad from both sources to generate a set of putative monad-monad matches where pairs meeting the filtering condition $|{\theta }_i-{\varphi }_i|\le \sigma$ are considered, with $\sigma$ set to $6\mathrm{°}$. This exclusion seeks to increase computational efficiency, although we note some potential unlikely matches are overlooked. The algorithm then calculates the relative bearing ${\theta }_r$ or ${\varphi }_r$ between pairs of monads leveraging the identified putative ${\theta }_i$ or ${\varphi }_i$ angle matches.

Algorithm 1 is applied at this point to relative angles, in order to correct for dramatic shifts in ${\theta }_r$ that can occur, should the noise present in the map image (or the cadaster) completely change the ${\theta }_r$ angle (or ${\varphi }_r$ respectively). As depicted in Figure 3, this needs to be dealt with because the ${\theta }_r$ angle from the connecting line to the next line segment, proceeding clockwise, can either become very small or very large given only a small amount of noise. To ameliorate this fragility we check if ${\theta }_r$ is very close (in either direction) to the line segments that constitute the monad, or 0. If so, we generate an edge case alternative value for ${\theta }_r$, referred to as ${\theta }_{r_e}$, which is then later computed during the distance step alongside ${\theta }_r$. The best score resulting from the comparison of these two is treated as the true ${\theta }_r$. This algorithm effectively handles the edge cases for ${\theta }_r$ (and similarly ${\varphi }_r$), adjusting ${\theta }_{r_e}$ based on its proximity to the edge of the monad’s line segments, while accounting for a noise factor.

Having generated putative monad matches between the two sources, the dyad selection stage can begin. This step creates the list of relative bearing angle combinations between filtered pairs of monads that each have an inter-source putative match. The algorithm then finds potential inter-source dyad matches and can therefore start to build out potential valid inter-source tuples of dyads where each constituent monad has a putative monad match that is contained in the alternative putative dyad.