CSIS Discussion Paper No. 34
A Computational Procedure for Making Seamless Map Sheets
Takashi Sato*, Yukio Sadahiro**, and Atsuyuki Okabe**
JANUARY, 2001

*Department of Urban Engineering, University of Tokyo**Center for Spatial Information Science, University of Tokyo

Center for Spatial Information Science
University of Tokyo
7-3-1, Bunkyo-ku, Tokyo 113-8656, Japan

January 11, 2001
A Computational Procedure for Making Seamless Map Sheets(seamless map sheets)

Takashi Sato*, Yukio Sadahiro**, and Atsuyuki Okabe***Department of Urban Engineering, University of Tokyo**Center for Spatial Information Science, University of Tokyo
Contact:

Phone:
Fax:
E-mail:Yukio Sadahiro
Center for Spatial Information Science andDepartment of Urban Engineering, University of Tokyo
7-3-1, Hongo, Bunkyo-ku, Tokyo 113-8656, Japan+81-3-5841-6273+81-3-5841-8521
sada@okabe.t.u-tokyo.ac.jp

Keywords: seamless map sheets, GIS, spatial database, transformation, computational
procedure
Acknowledgment: This paper is an extended version of Sato et al. (2000) which has
been published in Japanese. Authors express thanks to those who gave us valuablecomments on the earlier version of the paper.

A Computational Procedure for Making Seamless Map Sheets
Abstract
Map sheets are often used as basic spatial units for managing spatial data produced
from paper maps. This often results in incompatibility between adjacent map sheets spatial objects do not cross the boundaries smoothly and even the boundaries themselves
do not match their neighbors exactly. To solve the problem this paper proposes a
computational procedure for making seamless map sheets. Line objects digitizedseparately in different map sheets are considered, which are frequently used to representroad networks, gas pipelines, and boundaries of polygon objects. The procedure consistsof three steps: 1) extraction of end nodes, 2) detection of matching nodes, and 3)
transformation of the map sheet. Each step goes interactively so that unexpected errors
can be avoided by human observation. To test the validity of the procedure, map sheets
are combined containing the road network data of Tokyo 23-ku area, Japan.

1 Introduction
Paper maps are usually compiled as separate map sheets. In general, rectangularlattices are used as basic units of map sheets; in Japan, for example, 1/25,000
topographic maps are provided as 1.5 * 2.5 km rectangular maps. Landsat TM images
use 180 * 180 km square lattices.
This tradition still exists in spatial database - the digital form of paper maps. Map
sheets are often used as basic spatial units for managing spatial data produced from paper
maps. Maps are digitized separately and compiled as a set of map sheets in GIS.
However, it is obviously desirable for GIS users to use sheetless and seamless spatial
data because they wish to study regions larger than one map sheet. This complaint applies
to not only local analysis of geographical phenomena such as urban landuse and traffic
patterns but also sciences on a global scale, say, earth environment science and global
climatology.
Unfortunately, seamless spatial data cannot be easily obtained when paper maps are
separately digitized (Beard and Chrisman 1988). In theory, spatial data can be seamless if
individual map sheets are correctly digitized - if their reference to the earth's surface andthe location of spatial objects are correctly recorded in GIS. In practice, however, various
types of errors creep in map digitizing process such as displacement of control points,
distortion of papers, and human errors in manual operations. These errors cause
disagreement among map sheets: spatial objects do not cross the boundaries smoothly andeven the boundaries themselves do not match their neighbors exactly.
There are two approaches to obtain seamless data from separate map sheets: global
and local transformations. In either case one or several map sheets are transformed tomeet the adjacent ones. The basic difference between them lies in the scale oftransformation. Global methods transform the whole area of maps so that all the spatial
objects in map sheets may be affected. Local methods, on the other hand, work onlyaround map boundaries, preserving spatial objects near the center of map sheets.
In global methods, transformation is usually limited to the affine transformation to
avoid unnatural deformation of maps. Because of this global methods inevitably fail tomatch all the end nodes with those in adjacent map sheets (Petersohn and Venderohe
1982; White and Griffin 1985). This implies that sliver polygons still remain after thetransformation if applied to polygon data. Local expansion and shrinkage of paper maps
cannot be retrieved appropriately.
Among local methods, the nearest matching method is most frequently used in GIS.
This method moves end nodes to their nearest ones so that lines are exactly connected
between map sheets. A fuzzy tolerance is used to avoid excessive displacement of nodes.
Since the nearest matching method is easy to understand and its implementation is1straightforward, it is popular in geography and cartography; some systems can
automatically merge map sheets by the nearest matching method.
An extension of the nearest matching method is an algorithm called ZIPPER
proposed by Beard and Chrisman (1988). It uses two parameters to detect matchingnodes and to control the distortion area of maps. Consequently, it smoothly transforms
the original data to produce seamless maps containing no sliver polygons.
Unlike global methods, local methods cannot appropriately retrieve global distortionof maps which can be caused by systematic errors in digitizing process and numericalcalculation. Therefore, as Beard and Chrisman (1988) noted, a combination of global andlocal methods seems a desirable option because it can resolve both global and local errors.
Following the line of this hybrid approach, this paper proposes a new computational
procedure for making seamless map sheets. The procedure consists of three steps: 1)
extraction of end nodes, 2) detection of matching nodes, and 3) transformation of the map
sheet. Each step goes interactively so that unexpected errors can be avoided by humanobservation.
In the following we consider line objects without attribute information including
identification numbers; this type of line data are frequently used to represent road
networks, gas pipelines and boundaries of polygon objects. We assume that spatial data
are managed in rectangular lattice, and that one of the map sheets is distorted so that it has
to be transformed to match its surrounding sheets (Figure 1). The procedure, however, isalso applicable to matching of two adjacent sheets. In the following sections, each step of
the procedure is described successively. In Section 5, a set of road network maps are
combined to test the validity of the procedure. The last section summarizes theconclusions with discussion.
Figure 1 A distorted map sheet (thick lines) that has to be transformed to match its
surrounding sheets.
2 Extraction of end nodes
Suppose a rectangular lattice used for managing a set of map sheets containing
line objects. Let R0 be a rectangular region of the lattice whose corresponding map sheet
is distorted. We call the region R0 the standard unit. Our objective is to transform themap sheet to make it compatible with its surrounding sheets. Note, however, that thisdoes not mean that the surrounding sheets are perfectly correct but that they are relatively
correct compared with the target sheet.
The region covered by the distorted map sheet is called the distorted unit denoted
by R'. The distorted unit may not be rectangular because map distortion does not always2preserve the original shape of maps. The region that originally has to be covered by the
distorted map is called the true unit which is denoted by R. The true unit is not always
identical with R0 because errors usually exist in the surrounding sheets. Figure 2 shows
an example of the three units R0, R, and R’. The process of merging map sheets can beinterpreted as to find a one-to-one mapping from R' to R taking R0 into consideration.
Figure 2 An example of the standard unit R0 (gray-shaded area), true unit R (solid lines),and distorted unit R’ (broken lines).
Boundaries of R 0 , R and R ' are called the standard, true, and distorted
boundaries, and denoted by ∂R 0, ∂R, and ∂R’, respectively. In general, the standard
boundary ∂R0 is explicitly recorded in spatial database or at least computable from known
coordinates. The true and distorted boundaries, on the other hand, are not recorded sothat we have to estimate them from the data.
Line objects are distributed in the distorted unit R'. The number of lines linked to anode is called the order of the node. Nodes are classified by their order, and the set ofnodes of order i in R' is denoted by Ωi. The set Ω1 represents the set of end nodes. Let#(Ω i) be the number of nodes in Ω i, and Ω be the set of all the nodes in R', that is,
U Ωi .
i

The first step to merge map sheets is to extract end nodes, candidates for nodes tobe matched with those in adjacent map sheets. They are classified into two types:boundary nodes and interior nodes. The former is a set of nodes that are expected to be
exactly on the distorted boundary ∂R’. They usually have their corresponding nodes in
the adjacent map sheets and consequently quite important in merging map sheets. The
latter is also a set of nodes that seem to match nodes in the adjacent sheets though locatedinside ∂R’.
Extraction of boundary and interior nodes is a computational iterative process: wefirst set the initial candidates and update them iteratively until a predetermined condition issatisfied. The initial candidates are extracted from Ω1 as follows. Since boundary andinterior nodes are usually concentrated in the "neighborhood" of the standard boundary∂R0, the ratio of #(Ω1) to #(Ω) is extremely high near ∂R0. Mathematically, the ratio isrepresented by# (Ω1 ∈ B( x )) # (Ω1 ∈( R' \ B( x )))r( x ) =−,(1)# (Ω ∈ B( x )) # (Ω ∈( R' \ B( x )))where B(x) is the inner buffer region of distance x from ∂R0. The value of x maximizing
Equation (1) is denoted by xmax, and the inner buffer region B(xmax) is considered to be3

the "neighborhood" of ∂R 0 . End nodes in B(x max ) and outside of R 0 are the initial
candidates for boundary nodes, while those in B(2xmax) \ B(xmax) are candidates forinterior nodes.
Figure 3 Initial candidates for boundary (white circles) and interior nodes (double
circles). Gray-shaded area represents the inner buffer region B(xmax).
Initial candidates for boundary and interior nodes usually contain inappropriate
nodes, that is, nodes that do not have corresponding nodes in the adjacent map sheets. To
eliminate them we iteratively estimate the distorted boundary ∂R’ from candidate nodes
and update them by dropping the farthest one from the estimated boundary ∂Rˆ' until all
the candidates for boundary nodes are within the distance xmax/2 from ∂Rˆ' . The distorted
boundary is estimated by fitting the line that minimizes the sum of the square distance tothe candidates for boundary nodes. Figure 4 illustrates an example of this process.
Figure 4 Update of candidates for boundary nodes. Broken lines represent the estimatedboundary ∂Rˆ' . Black circles indicate boundary nodes rejected during the update process.
3 Detection of matching nodes
Having extracted end nodes in R’, we have to detect their corresponding nodes in
the adjacent map sheets. In theory, all of order 1 nodes in the adjacent sheets are
candidates for corresponding nodes. However, it is practically sufficient to examine only
the nodes in the "neighborhood" of ∂R0 which is given by the buffer operation.
Among the end nodes in R’ we first examine only the boundary nodes. We assign
them to the nearest nodes in the adjacent map sheets, and adopt node pairs that are thenearest neighbors with each other (Figure 5a). To improve the accuracy of matching, wethen drop the p percent (say, five percent) of the pairs whose node-to-node distance isrelatively long.
We next investigate the interior nodes and the boundary nodes rejected in theprevious step together. One can also include the nodes rejected in the update process ofcandidates for boundary nodes, if necessary. The procedure is quite similar; we assign
the nodes to the nearest nodes in the adjacent map sheets, and adopt node pairs that are
the nearest neighbors with each other (Figure 5b). We again drop the p percent of thepairs whose node-to-node distance is relatively long. The node pairs detected in these two
steps are considered to be "matching" (see Sato 2001 for details).
Detection of matching nodes plays a key role in merging map sheets. Consequently,
it is desirable for GIS users to examine the result of detection visually and revise it if

4necessary. The system developed by Sato (2001) provides the function of visualizing
matching nodes and their displacement vectors.
Figure 5 Detection of matching nodes. Solid lines represent line objects in the distorted
map sheet while broken lines in its adjacent sheet. White and black circles indicate theboundary and interior nodes, respectively. (a) Investigation of the boundary nodes. Nodepairs encircled by broken ellipses are the nearest neighbors with each other. (b)Investigation of the interior nodes and the boundary nodes rejected at the previous step.
Using the matching nodes detected, we can estimate the true and distorted units R
and R’. Let Ni' be the ith node in the distorted map sheet detected in the above procedure(i=1, ..., n; Nn+1'=N1'). We assume that the nodes are numbered counterclockwise withrespect to the gravity center of the standard unit R 0 denoted by G. The node in theadjacent sheets corresponding to Ni' is denoted by Ni. The polygons N1N2N3...Nn andN1'N2'N3'...Nn' are good estimators of R and R’, respectively.
4 Transformation of the map sheet
Having detected matching nodes, we finally transform the distorted map sheet toobtain seamless spatial data. Transformation of the map sheet proceeds by two steps: aglobal transformation and a set of local transformations.
In general, local transformations are not desirable because they often yield unnatural
distortion at the boundary of regions where different transformations are applied, even if
the transformation functions are smoothly continuous across the regions. On the other
hand, seamless map sheets cannot be obtained only by the global transformation. We thus
adopt a hybrid approach of these two methods in order to obtain seamless map sheets
without unnatural local distortions.
In the global transformation, nodes in the distorted map sheet are moved by theaffine transformation that minimizes the sum of the square distance between matchingnodes. The affine transformation includes not only the similarity transformation but also
the mirroring and shearing transformations. If the mirroring and shearing transformations
are to be avoided, the similarity transformation should be chosen which permits only thetranslation, rotation, and scaling of spatial objects. The transformation can also be limited
only to the translation and rotation to keep the map scale.
The global transformation is followed by a set of local transformations. The
distorted unit R’ is divided into a set of triangular regions each of which is given by G,
Ni', and Ni+1' (Figure 6, broken lines). Nodes are added to the line objects that cross theboundaries of triangles. The local transformation function is determined for each5

triangular region - transformation from the triangle ∆GN i'N i+1 ' into the trilateral∆GNiNi+1 (Figure 6, solid lines). Note that the resulted regions may not be triangles but
may contain curved edges; it depends on the transformation functions.
Figure 6 Tessellations of R (solid lines) and R’ (broken lines).
Let us consider the transformation from ∆GNi'Ni+1' into ∆GNiNi+1. Set the originof the coordinate system at G. The locational vector of a point in ∆GN i'N i+1 ' isrepresented by two parameters s and t:
x' = s{tx i ' +(1 − t )x i +1' } (0 ≤ s, t ≤ 1) ,(2)
where xi’ is the locational vector of Ni'. Let Ti,i+1(x') be the transformation function atx', that is, the point at x' in ∆GN i'N i+1 ' is moved to T i,i+1 (x'). The displacement
function ∆Ti,i+1(x') is then given by∆Ti , i +1 ( x' ) = Ti , i +1 ( x' ) − x' .
(3)
We adopt the transformation whose displacement function is∆Ti , i +1 ( x' ) = γ (s){t∆Ti , i +1 ( x i ' ) + (1 − t )∆Ti , i +1 ( x i +1' )},(4)

where s and t are given by Equation (2). The function γ(s) is a monotonically increasing
function of s that satisfies
γ (0) = 0 ,(5)andγ (1) = 1 .
Substituting Equations (3) into (4), we obtain the transformation function
Ti , i +1 ( x' ) = s{tx i ' +(1 − t )x i +1' } + γ (s){t ( x i − x i ' ) + (1 − t )( x i +1 − x i +1' )}= t{γ (s)x i + (s − γ (s))x i ' } + (1 − t ){γ (s)x i +1 + (s − γ (s))x i +1' }

(6)

.

(7)

The transformation given by Equation (7) keeps the topology of line objects in R’; thegravity center G is not moved while the vertices Ni' and Ni+1' are exactly moved to Ni
and N i+1 , respectively; line segments parallel to N i 'N i+1 ' are not curved after
transformation. To keep the topology of spatial objects is critical in GIS because it
represents the structure of spatial objects. The monotonicity of γ(s) assures that thedisplacement of spatial objects increases monotonously with s; points near the gravity
center G move slightly while those around the edge N i’N i+1’ move greatly. Figure 7
shows an example of the transformation.
Figure 7 Transformation of ∆GN1'N2' (broken lines) into ∆GN1N2 (solid lines).
There are a wide variety of functions satisfying the conditions mentioned above. In
this paper we use the power function6γ ( s ) = sα .

(8)

The parameter α determines the range of transformation; a small value of α indicates thatthe transformation works globally; as the value of α increases the transformation becomes
local, that is, only the spatial objects near the boundaries are moved. The trilateral∆GNiNi+1 becomes a triangle when α=1.
The parameter α is so determined as to maximize the correlation coefficient of thebidimensional regression where the points in N1'N2'N3'...Nn' and N1N2N3...Nn are thedependent and independent variables, respectively (for details of the bidimensional
regression, see Tobler 1994 and Nakaya 1997). The correlation coefficient ρ is given byn∑∫ρ = 1 − i =1n2x' ∈GN i ' N i +1 '

∑∫
i =1

x' − Ti , i +1 ( x' ) dx'
2

x' ∈GN i ' N i +1 ',(9)x' − x i , i +1' dx'

where x i , i +1' is defined by
x i , i +1' =∫ x' dx'∫ dx'

x' ∈GN i ' N i +1 '

.

(10)

x' ∈GN i ' N i +1 'The correlation coefficient ρ represents how close the transformation fromN 1 'N 2 'N 3 '...N n ' to N 1 N 2 N 3 ...N n is to the similarity transformation. Consequently, alarge value of ρ is desirable because the similarity transformation does not yield
unnnatural distortion of spatial objects. Integral terms in Equations (9) and (10) arepractically calculated by numerical integration.
Local transformation functions are calculated for individual triangular regions, andapplied to all the nodes of line objects including those added in the triangulation of R’.
This, however, often changes straight lines running across triangles into polygonal lines
because the transformation function is different among the regions. We thus finally
remove the redundant nodes on the boundaries of triangles to restore straight lines. The
topology of line objects is examined at each step of node removal to preserve the original
topology. Figure 8 shows an example of the transformation process.
Figure 8 Transformation of line objects. (a) Original (distorted) lines and triangles, (b)triangulation (nodes are added), (c) transformation, (d) removal of redundant nodes.
5 Empirical study
To test the validity of the procedure proposed above, we combine map sheetscontaining the road network data of Tokyo 23-ku area, Japan. About 250 map sheets
cover the Tokyo 23-ku area, each of which is a rectangle of 1.5 * 2.5 km.
We have transformed three map sheets that do not match their surrounding sheets,7

but show only the result of one of the sheets S0 due to the limitation of space. Figures 9
and 10 shows the data in the original sheet and some mismatches between S0 and its
adjacent map sheet.
Figure 9 Original data of the map sheet. Solid and broken lines indicate the road networkin the S0 and its adjacent sheets, respectively.
Figure 10 Mismatches between S0 and its adjacent map sheet.
To transform S 0 we first extracted its end nodes, candidates for nodes to bematched with those in adjacent map sheets, using the method proposed in Section 2. The
distance xmax which determines the "neighborhood" of ∂R0 was set to 0.240m, whichyielded 110 initial candidates for boundary nodes and 26 interior nodes. The initial
candidates for boundary nodes were then updated by the iterative estimation of the
distorted boundary, and as a result, 83 boundary nodes were extracted. Figure 11 shows
the boundary and interior nodes in S0.
Figure 11 Boundary and interior nodes extracted in S0. White and black circles indicate
the boundary and interior nodes, respectively.
Though the distance xmax was considerably small, extraction of end nodes seemed
successful by visual observation. We thus accepted the results and proceeded to detect
their corresponding nodes in the adjacent map sheets. At the first step of detection, only
the boundary nodes in S0 were examined. Among 83 boundary nodes, 67 nodes wereassigned to their corresponding nodes in the adjacent sheets. We then dropped the five
percent of the pairs whose node-to-node distance is relatively long. At the second step,
we investigated the interior nodes and the boundary nodes rejected together, and foundcorresponding nodes for 33 interior and boundary nodes among 42. We again dropped
the five percent of the pairs whose node-to-node distance is relatively long, and obtained
100 matching node pairs. Figure 12 shows several examples of the nodes that failed tofind their corresponding nodes.
Figure 12 Nodes that failed to find their corresponding nodes (black circles). Solidand broken lines indicate the road network in the S0 and its adjacent sheets, respectively.
Using the matching nodes detected we finally transformed the map sheet S0 by
global and local transformations. Global transformation is so determined as to minimize

8the sum of the square distance between matching nodes. After the global transformation,
the map sheet S0 is divided into 100 triangular regions as shown in Figure 13 in each ofwhich a local transformation is applied.
Figure 13 Tessellations of the map sheet S0.
The parameter α in the function γ(s) which determines the range of transformationis set to 0.280 that maximizes the correlation coefficient of the bidimensional regression
mentioned earlier. Since the value is relatively small, the transformation is performed
locally around the boundaries. Local transformations were followed by removal ofredundant nodes on the boundaries of triangles to restore straight lines. Figure 14 shows
several examples of the removal of redundant nodes, which effectively demonstrates theusefulness of this operation. Figure 15 shows the final result of transformation in which
the road network is exactly compatible between adjacent map sheets.
Figure 14 Examples of the removal of redundant nodes. Road networks before(broken lines) and after (solid lines) the removal.
Figure 15 Final production of the transformation.
6 Concluding discussion
In this paper we have proposed a computational procedure for making seamless
map sheets. The procedure consists of three steps: 1) extraction of end nodes, 2)
detection of matching nodes, and 3) transformation of the map sheet. It is applicable tospatial data representing line objects managed in the rectangular lattice such as road
networks and gas pipelines. The procedure was empirically applied to the map sheetscontaining the road network data of Tokyo 23-ku area, Japan, and successfully combined
the map sheets into seamless spatial data.
Advantages of our approach over existing methods are summarized as follows.
1) Our method uses both global and local transformations to merge map sheets.
Consequently, it works for maps distorted in either scales; unlike global
methods, it produces exactly seamless maps; unlike local methods, it keeps
unnatural deformation of maps to a minimum.
2) The rectangular boundaries of map sheets effectively work as a key information
to extract end nodes. They greatly contribute toward extracting correct nodes,
and consequently, improve the efficiency of the whole process.

93) Our method does not use attribute information of line objects including
identification numbers, which provides its wide applicability. This property
makes it applicable to polygon data whose boundaries often do not haveattribute information; we can merge map sheets containing polygon objects bytreating their boundaries as line objects.
4) Unlike existing local methods, the transformation function works not only near
the map boundaries but also over the whole map sheet. This provides smoothtransformation of spatial objects.
5) The map transformation keeps the topology of original lines. This property is
indispensable especially to line objects such as road networks where thetopology is of great importance.
One of the characteristics of our method is, as mentioned above, that it uses therectangular boundaries of map sheets as a key information to extract end nodes. This isbecause spatial data are usually managed in rectangular lattice. However, the method canbe applied to other systems of data management such as the longitude-latitude system; we
can replace the linear function used in estimation of map boundaries by curvilinear
functions to treat various systems. In addition, our method works even if the map sheet
system is not known. Though the accuracy of end node extraction decreases, we can still
merge map sheets by interactively checking end nodes extracted.
Our method can be extended further, because it does not make full use of
information provided by spatial objects. For example, if lines have attribute information,
it will drastically improve the accuracy of the procedure. In road network data, names ofmajor streets are often recorded as an attribute, which are quite useful for matching node
extraction. Even the data of road width would be helpful in finding node pairs. The shapeof line objects also gives a clue to find matching nodes, because lines have to be
combined smoothly at the map boundary. Though it was sufficient in the empirical study
to use only the locational information of end nodes, taking additional information into
account is an important subject for future research.

10

References
Beard, M. K., Chrisman, N. R., 1988: Zipper: a localized approach to edgematching.
American Cartographer 15: 163-172.
Nakaya, T., 1997: Statistical inferences in bidimensional regression models.
Geographical Analysis 29: 169-186.
Petersohn, C., Venderohe, A., 1982: Site-specific accuracy of digitized property maps.
Auto-Carto 5: 607-619.
Sato, T., Sadahiro, Y., and Okabe, A., 2000: An interactive method of detecting andvisualizing the property of spatial data for joining spatial data units. Papers andProceedings of Geographic Information Systems Association 9: 399-402.
Sato, T., 2001: An interactive method for joining spatial data units. Master’s Thesis,
Department of Urban Engineering, University of Tokyo.
Tobler, W. R., 1994: Bidimensional regression. Geographical Analysis 26: 187-212.
White, M., Griffin, P., 1985: Piecewise linear rubber sheet map transformation.
American Cartographer 12: 123-131.

1

Figure 1

R0

Figure 2

distorted boundary
standard boundary

xmax

Figure 3xmax

standard boundary(a)xmax/2(b)

xmax/2(c)Figure 4

standard boundary

(a)

(b)

Figure 5

N13 N13'

N1

N12'

N11'

N12

N1'

N10

N11

N10'

N9'

N2
N2'

N9
G

N8
N4

N3'
N3

N4'

N5

N8 '

N6'
N5'

Figure 6

N6

N7'

N7

N 1'
N1

G
N2'
N2

Figure 7

(a)

(b)

(c)

(d)

Figure 8

Figure 9

Figure 10

Figure 11

Figure 12

Figure 13

Figure 14

Figure 15