A Raster-Based Methodology to Detect Cross-Scale Changes in Water Body Representations Caused by Map Generalization

In traditional change detection methods, remote sensing images are the primary raster data for change detection, and the changes produced from cartography generalization in multi-scale maps are not considered. The aim of this research was to use a new kind of raster data, named map tile data, to detect the change information of a multi-scale water system. From the perspective of spatial cognition, a hierarchical system is proposed to detect water area changes in multi-scale tile maps. The detection level of multi-scale water changes is divided into three layers: the body layer, the piece layer, and the slice layer. We also classify the water area changes and establish a set of indicators and rules used for the change detection of water areas in multi-scale raster maps. In addition, we determine the key technology in the process of water extraction from tile maps. For evaluation purposes, the proposed method is applied in several test areas using a map of Tiandi. After evaluating the accuracy of change detection, our experimental results confirm the efficiency and high accuracy of the proposed methodology.


Introduction
The multi-scale representation of spatial data is a classic topic in cartography that requires many generalization operations, such as simplification, aggregation, typification, and displacement [1][2][3][4][5][6], to be realized. Water is necessary for economic, social, and environmentally sustainable development. For the multi-scale representation of water areas, generalization operations such as the selection of objects, the simplification of boundaries, and the aggregation of multiple objects are usually applied, which leads to changes in water areas in terms of the quantity, shape, and distribution. To solve the quality problem of multi-scale maps, a change detection method can usually be applied.
Change detection refers to the process of identifying differences through a comparative analysis of the state of the same geographic entity acquired from images or databases in different periods [7], and it is an important method to maintain spatial data and update spatial databases [8,9]. Consistent multiresolution spatial data are required in many applications, such as rapidly transmitting spatial data over the internet [10], querying multi-resolution spatial data [11], and extracting and integrating spatial data information with different levels of detail [12].
On one hand, the traditional spatial data change detection methods ignore the influence of cartographic generalization, which cannot be applied for the change detection of cross-scale water body representations. For example, many scholars have detected changes in water systems at a single or adjacent scale, such as lake change detection [13][14][15], river change detection [16,17], coastline or coastal zone change detection [18][19][20], and water system [21] change detection. In addition, some mature This kind of technique for solving the problem by starting with the general scope and then focusing on the details, or by starting from the whole and moving to the partial, is the so-called hierarchical method. Change is a kind of difference information in spatial data. This paper argues that people should consider hierarchical features in the process of acquiring regularity; then, they can be aware of the essence of change information from simple to complex scales and better master the change rules of the spatial data from the macro to micro scales. Using hierarchical change detection, we can obtain the right amount of change information at different scales and avoid cognitive difficulties.

Division of Detection Hierarchy
To see details and focus their attention on the main problem without interfering with minor issues, people always hierarchically abstract subjects in an orderly way when thinking. As the level of detail exceeds a certain degree, the more things a person can see, the less things they can express [35].
The detection of detailed change information shows that we detect spatial change information between two uniform or adjacent scales. However, to better grasp the global change situation of spatial data, we need to observe and analyze changes in multi-scale spatial data from the macro view. Hence, this article attempts to study multi-scale spatial data changes at different spatial scales.

Body Layer
Spatial entities within the same scope have different detailed expressions at different map scales. Multi-scale spatial data can be seen as a collection of data composed of a group of spatial objects on the order of different map scales. This collection of data can be called a spatial series. Each object of the spatial series corresponds to an expression of the spatial entity at a certain map scale. A spatial series of the same geographical object can be expressed mathematically: This expression can be abbreviated as {X n , n ∈ N}. Similar to the relationship between the samples and the sample observations, we can use Formula (2) to express n ordered sets of observations of the spatial series in Formula (1). Formula (2) can be called the observation series of length n: x 1 , x 2 , · · · , x n (2) In the body layer, we analyze and research the change information of the spatial data from macroscopic and global perspectives, and the rules or information included in the orderly observations x 1 , x 2 , · · · , x n of the geographical objects are discovered.
For example, when the observed value x represents the area of a lake, namely, x n = Area n , through the variance calculation of this set of observations {Area 1 , Area 2 , · · · , Area n }, we can obtain the dispersion degree of area change in a lake in a multi-scale map. If the variance D(Area) is relatively large, the area of the lake would have a high dispersion degree and noted fluctuations in n scales. If the variance D(Area) is relatively small, the area of the lake would have a low dispersion degree and small fluctuations in n scales.

Piece Layer
In the map expression space, the scale range of every spatial entity is limited [36]. Taking a building as an example, buildings are usually expressed in the form of a polygon in large-scale maps, a simplified rectangle in medium-scale maps, and a point on small-scale maps. On maps with a scale of 1:250,000 or smaller, buildings can no longer be expressed unless they are special. This kind of change, with a substantial difference from a polygon to a rectangle or from a rectangle to a point, is called a mutation in this paper. In the piece layer, the mutation information is detected and expressed.
For a set of spatial series X 1 , X 2 , · · · , X n , the purpose of detection in the piece layer is to group this set of series, which makes the differences obvious among different groups, while the differences are less obvious among objects in one group. For example, for a set of spatial series {X 1 , X 2 , X 3 , X 4 }, if the difference ∆ 12 between object X 1 and object X 2 is relatively small, and the difference ∆ 34 between object X 3 and object X 4 is also relatively small, while the difference ∆ 23 between object X 2 and object X 3 is relatively large, then X 1 and X 2 can be classified as one group, and X 3 and X 4 can be classified as one group. Therefore, in the process of piece layer detection, we should focus on the larger difference ∆ 23 between the groups instead of the smaller difference ∆ 12 or ∆ 34 in one group.
As shown in Figure 1, the change between S 1 and S 2 is due to the process of boundary simplification of the double river. This kind of change is a slow change. The change after S 2 is a process of transforming the geometrical dimension, which can be seen as a mutation. The change between S 3 and S 4 is due to the process of shape simplification of a single river. This kind of change is a slow change. After S 4 , the river disappears, which can be seen as a mutation again. In the piece layer, we focus on how to detect the mutation change between S 2 and S 3 instead of the slow change between S 1 and S 2 . The reason is that the mutation change produces more differences than slow changes.
Sensors 2020, 12, x FOR PEER REVIEW 4 of 19 difference ∆ between object and object is relatively small, and the difference ∆ between object and object is also relatively small, while the difference ∆ between object and object is relatively large, then and can be classified as one group, and and can be classified as one group. Therefore, in the process of piece layer detection, we should focus on the larger difference ∆ between the groups instead of the smaller difference ∆ or ∆ in one group.
As shown in Figure 1, the change between and is due to the process of boundary simplification of the double river. This kind of change is a slow change. The change after is a process of transforming the geometrical dimension, which can be seen as a mutation. The change between and is due to the process of shape simplification of a single river. This kind of change is a slow change. After , the river disappears, which can be seen as a mutation again. In the piece layer, we focus on how to detect the mutation change between and instead of the slow change between and . The reason is that the mutation change produces more differences than slow changes.

Slice Layer
After grouping a set of spatial series , , ⋯ , in the piece layer, the difference between adjacent scales in each group is smaller. Relative to the piece layer, in the slice layer, this kind of change information is detected in every "piece" with a smaller difference. For example, in the spatial series { , , , } in Section 2.3, and are in one group, and the difference ∆ between and is smaller. In the slice layer, this kind of change with a smaller difference, such as ∆ , is handled.
In the process of map generalization [37][38][39][40][41][42], boundary simplification [43][44][45] usually produces changes in the slice layer. As shown in Figure 2, a famous simplification algorithm [46] is used for the smoothing curve. Figure 2 shows the smoothing results (red lines) of originalboundaries (blue lines) with different tolerances. This kind of detailed change with a smaller difference in Figure 2 is a slice layer change.

Change Classification and Detection of Water Area
It is necessary to consider the change type in change detection time. In the past few years, scholars have classified the change type from different perspectives depending on the type of object. Claramunt subdivided the evolution of a single entity into eight types: appearance and disappearance, movement, rotation, expansion, shrinkage, and deformation [47]. Raza and Kainz treated land blocks as the analysis object and subdivided the change types of areal objects into eight types: appearance, disappearance, size change, shape change, movement, transformation, merging, S1 S2 S3 S4 S5

Slice Layer
After grouping a set of spatial series X 1 , X 2 , · · · , X n in the piece layer, the difference between adjacent scales in each group is smaller. Relative to the piece layer, in the slice layer, this kind of change information is detected in every "piece" with a smaller difference. For example, in the spatial series {X 1 , X 2 , X 3 , X 4 } in Section 2.3, X 1 and X 2 are in one group, and the difference ∆ 12 between X 1 and X 2 is smaller. In the slice layer, this kind of change with a smaller difference, such as ∆ 12 , is handled.
In the process of map generalization [37][38][39][40][41][42], boundary simplification [43][44][45] usually produces changes in the slice layer. As shown in Figure 2, a famous simplification algorithm [46] is used for the smoothing curve. Figure 2 shows the smoothing results (red lines) of originalboundaries (blue lines) with different tolerances. This kind of detailed change with a smaller difference in Figure 2 is a slice layer change.
Sensors 2020, 12, x FOR PEER REVIEW 4 of 19 difference ∆ between object and object is relatively small, and the difference ∆ between object and object is also relatively small, while the difference ∆ between object and object is relatively large, then and can be classified as one group, and and can be classified as one group. Therefore, in the process of piece layer detection, we should focus on the larger difference ∆ between the groups instead of the smaller difference ∆ or ∆ in one group.
As shown in Figure 1, the change between and is due to the process of boundary simplification of the double river. This kind of change is a slow change. The change after is a process of transforming the geometrical dimension, which can be seen as a mutation. The change between and is due to the process of shape simplification of a single river. This kind of change is a slow change. After , the river disappears, which can be seen as a mutation again. In the piece layer, we focus on how to detect the mutation change between and instead of the slow change between and . The reason is that the mutation change produces more differences than slow changes.

Slice Layer
After grouping a set of spatial series , , ⋯ , in the piece layer, the difference between adjacent scales in each group is smaller. Relative to the piece layer, in the slice layer, this kind of change information is detected in every "piece" with a smaller difference. For example, in the spatial series { , , , } in Section 2.3, and are in one group, and the difference ∆ between and is smaller. In the slice layer, this kind of change with a smaller difference, such as ∆ , is handled.
In the process of map generalization [37][38][39][40][41][42], boundary simplification [43][44][45] usually produces changes in the slice layer. As shown in Figure 2, a famous simplification algorithm [46] is used for the smoothing curve. Figure 2 shows the smoothing results (red lines) of originalboundaries (blue lines) with different tolerances. This kind of detailed change with a smaller difference in Figure 2 is a slice layer change.

Change Classification and Detection of Water Area
It is necessary to consider the change type in change detection time. In the past few years, scholars have classified the change type from different perspectives depending on the type of object. Claramunt subdivided the evolution of a single entity into eight types: appearance and disappearance, movement, rotation, expansion, shrinkage, and deformation [47]. Raza and Kainz treated land blocks as the analysis object and subdivided the change types of areal objects into eight types: appearance, disappearance, size change, shape change, movement, transformation, merging,

Change Classification and Detection of Water Area
It is necessary to consider the change type in change detection time. In the past few years, scholars have classified the change type from different perspectives depending on the type of object. Claramunt subdivided the evolution of a single entity into eight types: appearance and disappearance, movement, rotation, expansion, shrinkage, and deformation [47]. Raza and Kainz treated land blocks as the analysis object and subdivided the change types of areal objects into eight types: appearance, disappearance, size change, shape change, movement, transformation, merging, and splitting [48]. The simplest outcome for change detection is straightforward; that is, the target data are either changed or not [49]. However, it is necessary to take the influence of map generalization into account when we detect changes in water objects at multiple scales. The existing classification methods do not consider the problem of scales, and they are not appropriate for the change detection of water objects at multiple scales.
There are two kinds of reasons for water changes at multiple scales. First, map generalization can cause differences between different scales. Water systems can be divided into two types without considering the type of water: linear water and areal water. Linear water includes a single river system and a double river system. The generalized operation of the former includes the river selection number and chooses what and how to simplify the system. The generalized operation of the latter includes the extraction of the double river central axis and the simplification of a single river after extraction. Areal water includes lakes and reservoirs. The generalized operation of areal water includes selection and simplification. Second, the differences caused by the actual measurements may lead to water changes at different scales.
By analyzing the causes of changes, the change types of water areas at multiple scales are divided into 9 categories in this paper, for a total of 14 subtypes, on the basis of existing change types. For all of the change types, we establish a set of indicators and rules to detect them and use hierarchical thought to classify and express the change information in three layers: the body, piece, and slice layers.

Definition of Detection Indicators
The types of water area change are diverse. It is difficult to identify and differentiate the various change types when we only use some basic indicators of geometric objects, such as length, area, shape, and center. For the purpose of detecting and inferring the change types of water areas, we define a series of new indicators for identifying them, and those different indicators are combined for use in this paper. These inference indicators include the following nine types: matching type, length, area, shape index, geometric center, overlap degree, buffer overlap degree, geometric type, and hole state. The calculation methods of these indicators are as follows: Matching type MT: The matching types of an object area mainly include five kinds: MT = 1 : 0, MT = 0 : 1, MT = 1 : 1, MT = 1 : m, and MT = n : 1, in which m, n, 0, and 1 represent the number of objects.
Length L: The length is used to describe the boundary of the water area. Before the length of the object is calculated, it is necessary to extract the edge pixels of the object area. Then, the length of the water object is calculated by counting the number of pixels Num, namely, L = Num. The length difference ∆ L before and after the change can be represented as ∆ L = L a f ter − L be f ore .
Area A: The precise area of the region is calculated by counting the number of pixels Num, namely, A = Num. The area difference ∆ L before and after the change can be represented as: Shape index SI: The shape index is used to describe the shape of a linear object or the boundary of an areal object. The common methods to describe shape and calculate similarity include Fourier descriptors [50,51] and the angle difference integration method [52]. In this paper, the included angle chain method [53] is used to describe the shape and calculate the similarity of the linear object or the boundary of the areal object. In this method, as shown in Figure 3, a curve is modelled by a number of linked equi-length line segments, for example, a length of 50 pixels at levels 11 and 12 in the experimental section and a sequence of codes using the included angles between a pair of neighboring line segments are used to represent the curve. The representation is invariant to rotation, scaling, and translation [53]. The Euclidean distance of the included angle chain code in n-dimensional space can be used to measure the similarity between matching objects.  Suppose that the set of equi-length line segments of curve Q is , , ⋯ , ; then, the included angle chain code of curve Q can be denoted by A, and A is defined as shown below: In this equation, denotes the angle from line segment to line segment in the counter-clockwise direction. The shape difference ∆ before and after the change can be represented as In this equation, Because the angle number of two linear objects may be different, for example, linear object A with n angles and object B with angles, in the actual calculation, by combining other detection indicators, we only need to calculate ∆ using the first n angles with the same starting point and counter-clockwise direction to meet the identification requirement.
Geometric center GC: The geometric center of the irregular polygon is usually equal to its center of gravity when its density is uniform, and the location of the center of gravity is only related to the shape of the object. Assuming that f(i, j) is a matrix that represents an image, in which i and j are the row and column values of this matrix, respectively, the coordinates of geometric center ( , ) can be calculated as follows: The difference in geometric center ∆ before and after the change can be represented as ∆ = , , where dis is the Euclidean distance between two points. Overlap degree OD: For each group of matching water area objects, the overlap degree OD is the ratio between the overlapping pixel number and the total pixel number after overlaying, namely, = ⁄ . As shown in Figure 4, curve A is composed of 11 pixels, curve B is composed of 11 pixels, and the overlapping portions between curves A and B appear red. It is easy to see that the overlap degree of A is = = 2 ⁄ 11 = ⁄ 18.2%, and the overlap degree of B is = ⁄ = 2 11 = 18.2% ⁄ . As shown in Figure 4, after overlapping, when the pixels of one object A are in the buffer of pixel p, p is one pixel of object B, and then pixel p is also seen as an overlapping pixel, even if there is no pixel overlapping with p in object A. Suppose that the total overlapping pixel number of matching objects is ; the total pixel number of matching objects is , and then the buffer overlap degree of the matching objects can be expressed as = ⁄ . In Figure  4, when the buffer radius r is equal to 3 pixels, the buffer overlap degree of curve A and curve B is Suppose that the set of equi-length line segments of curve Q is {c 1 , c 2 , · · · , c n }; then, the included angle chain code of curve Q can be denoted by A, and A is defined as shown below: in this equation, α i denotes the angle from line segment c i to line segment c i+1 in the counter-clockwise direction. The shape difference ∆ SI before and after the change can be represented as Because the angle number of two linear objects may be different, for example, linear object A with n angles and object B with n + m angles, in the actual calculation, by combining other detection indicators, we only need to calculate ∆ SI using the first n angles with the same starting point and counter-clockwise direction to meet the identification requirement.
Geometric center GC: The geometric center of the irregular polygon is usually equal to its center of gravity when its density is uniform, and the location of the center of gravity is only related to the shape of the object. Assuming that f (i, j) is a matrix that represents an image, in which i and j are the row and column values of this matrix, respectively, the coordinates of geometric center GC(i c , j c ) can be calculated as follows: i c =m 10 /m 00 , j c = m 01 /m 00 , The difference in geometric center ∆ GC before and after the change can be represented as ∆ GC = dis GC a f ter , GC be f ore , where dis is the Euclidean distance between two points.
Overlap degree OD: For each group of matching water area objects, the overlap degree OD is the ratio between the overlapping pixel number Num lap and the total pixel number Num all after overlaying, namely, OD = Num lap /Num all . As shown in Figure 4, curve A is composed of 11 pixels, curve B is composed of 11 pixels, and the overlapping portions between curves A and B appear red. It is easy to see that the overlap degree of A is OD A = Num lap /Num all = 2/11 = 18.2%, and the overlap degree of B is OD B = Num lap /Num all = 2/11 = 18.2%.
As shown in Figure 4, after overlapping, when the pixels of one object A are in the buffer of pixel p, p is one pixel of object B, and then pixel p is also seen as an overlapping pixel, even if there is no pixel overlapping with p in object A. Suppose that the total overlapping pixel number of matching objects is Num lap ; the total pixel number of matching objects is Num all , and then the buffer overlap degree of the matching objects can be expressed as BOD = Num lap /Num all . In Figure 4, when the buffer radius r is equal to 3 pixels, the buffer overlap degree of curve A and curve B is  Geometric Type GT: The geometric type of a water object is classified into linear and area objects. When the object l is linear water, = 0, and when the object l is a water area, = 1.
Hole State HS: When there are no holes in the water area, = 0; otherwise, = 1.
Before the identification of change types, some simple operations, such as moving, scaling, and overlaying, may be applied to the water objects to distinguish all of the change types. In addition, the parameter values used in the actual experiment are not fixed. The selection of these parameter values depends on the image resolution and different application scenarios.

Change Types
In this paper, the change types of water areas are divided into 10 categories, for a total of 15 subtypes, as shown in Table 1. The 10 categories are no change, appearance, disappearance, movement, rotation, scaling, boundary generalization, morphology change, derivative change, and hole change. Among them, the morphology change includes appending, distribution, and shape changes; the derivative change includes merging and splitting; and the hole change includes appearance, disappearance, and structure change of the hole.  Geometric Type GT: The geometric type of a water object is classified into linear and area objects. When the object l is linear water, GT l = 0, and when the object l is a water area, GT l = 1.
Hole State HS: When there are no holes in the water area, HS = 0; otherwise, HS = 1. Before the identification of change types, some simple operations, such as moving, scaling, and overlaying, may be applied to the water objects to distinguish all of the change types. In addition, the parameter values used in the actual experiment are not fixed. The selection of these parameter values depends on the image resolution and different application scenarios.

Change Types
In this paper, the change types of water areas are divided into 10 categories, for a total of 15 subtypes, as shown in Table 1. The 10 categories are no change, appearance, disappearance, movement, rotation, scaling, boundary generalization, morphology change, derivative change, and hole change. Among them, the morphology change includes appending, distribution, and shape changes; the derivative change includes merging and splitting; and the hole change includes appearance, disappearance, and structure change of the hole.

Detection Rules
We define a series of detection rules to detect the change types of The detection rules are made up of operations and indicators. Fo operations, such as moving, scaling, or overlaying, may be used to de indicators are selected and combined. Each change type correspond Table 2

Detection Rules
We define a series of detection rules to detect the change types of water areas at multiple scales. The detection rules are made up of operations and indicators. For each type of change, some operations, such as moving, scaling, or overlaying, may be used to detect it. Then, several of the nine indicators are selected and combined. Each change type corresponds to an indicator combination. Table 2 shows the detection rules of all 15 water area change types presented in this paper.

Detection Rules
We define a series of detection rules to detect the change types of water areas at multiple scales. The detection rules are made up of operations and indicators. For each type of change, some operations, such as moving, scaling, or overlaying, may be used to detect it. Then, several of the nine indicators are selected and combined. Each change type corresponds to an indicator combination. Table 2 shows the detection rules of all 15 water area change types presented in this paper.

Detection Rules
We define a series of detection rules to detect the change types of water areas at multiple scales. The detection rules are made up of operations and indicators. For each type of change, some operations, such as moving, scaling, or overlaying, may be used to detect it. Then, several of the nine indicators are selected and combined. Each change type corresponds to an indicator combination. Table 2 shows the detection rules of all 15 water area change types presented in this paper.   Table 2.

Hierarchical Expression of Water Area Changes
For the changes in every water object area at multiple scales, we establish a form of hierarchical expression for displaying the detection results. As shown in Figure 5, in the body layer, by analyzing the discrete degree of the area, the perimeter, or the center of the water area, the global fluctuations in the changes in the water area at multiple scales can be expressed. The metrics of global fluctuation can be the maximum deviation, the mean deviation, the standard deviation, and so on. In the piece layer, the mutation information of the water area at multiple scales is expressed. The corresponding mutation types include appearance, disappearance, morphological change, derivative change, and the appearance or disappearance of holes. It is difficult to use quantitative indicators to measure the differences in the mutations. Therefore, in this layer, the mutations can be seen as qualitative changes.
Sensors 2020, 12, x FOR PEER REVIEW 10 of 19 layer, the mutation information of the water area at multiple scales is expressed. The corresponding mutation types include appearance, disappearance, morphological change, derivative change, and the appearance or disappearance of holes. It is difficult to use quantitative indicators to measure the differences in the mutations. Therefore, in this layer, the mutations can be seen as qualitative changes.
In the slice layer, the detailed change information of the water area at multiple scales is expressed. The corresponding change types include no change, movement, rotation, scaling, boundary generalization, and structural change of the hole. It is easy to use quantitative indicators to measure the differences in the detailed changes. Therefore, in this layer, the changes can be seen as quantitative changes. Among them, the difference in the movement change can be measured by distance, the difference in the rotation change can be measured by the angle, the difference in the scaling change can be measured by the ratio, the difference in boundary generalization can be measured by shape similarity, and the difference in the structure change of the hole can be measured by length or area.
In other words, from a cognitive rule perspective, the hierarchical methods proposed in this paper are used to detect and express the changes in water area at multiple scales (from global to local, from macro to micro, and from qualitative to quantitative), which enables us to grasp the essence of change information from rudimentary to profound and better master the laws of changes at multiple scales. .

Water Area Extraction and Reconstruction
How to extract water from tile elements is a problem that cannot be ignored. Therefore, before change detection, we focus on how to extract independent water areas from the tile map. Compared with the colors of the other elements in a map, water areas are usually expressed with a uniform blue color. Therefore, it is effective to separate water areas from other geographical elements on one map when using a given color threshold. However, due to the overlay of other elements on the map, such as roads and annotations, the extracted water areas contain obvious fractures or hollow phenomena. In addition, the boundary colors of the water areas are usually not the same color but have gradual change effects. As a result, the boundary of the water areas needs to be processed specifically. Considering these two factors, this paper uses the following steps to extract water areas from tile maps [9].
(1) Extract water areas with color segmentation In the slice layer, the detailed change information of the water area at multiple scales is expressed. The corresponding change types include no change, movement, rotation, scaling, boundary generalization, and structural change of the hole. It is easy to use quantitative indicators to measure the differences in the detailed changes. Therefore, in this layer, the changes can be seen as quantitative changes. Among them, the difference in the movement change can be measured by distance, the difference in the rotation change can be measured by the angle, the difference in the scaling change can be measured by the ratio, the difference in boundary generalization can be measured by shape similarity, and the difference in the structure change of the hole can be measured by length or area.
In other words, from a cognitive rule perspective, the hierarchical methods proposed in this paper are used to detect and express the changes in water area at multiple scales (from global to local, from macro to micro, and from qualitative to quantitative), which enables us to grasp the essence of change information from rudimentary to profound and better master the laws of changes at multiple scales.

Water Area Extraction and Reconstruction
How to extract water from tile elements is a problem that cannot be ignored. Therefore, before change detection, we focus on how to extract independent water areas from the tile map. Compared with the colors of the other elements in a map, water areas are usually expressed with a uniform blue color. Therefore, it is effective to separate water areas from other geographical elements on one map when using a given color threshold. However, due to the overlay of other elements on the map, such as roads and annotations, the extracted water areas contain obvious fractures or hollow phenomena. In addition, the boundary colors of the water areas are usually not the same color but have gradual change effects. As a result, the boundary of the water areas needs to be processed specifically. Considering these two factors, this paper uses the following steps to extract water areas from tile maps [9].
(1) Extract water areas with color segmentation The water colors of different tile maps are different. However, the colors of water in one map are uniform. Therefore, it is effective to separate water from other geographical elements in one map when using a given color threshold. As shown in Figure 6, the left figure is a tile map of Tiandi, and the right figure is the result after extraction by using the color threshold. However, this extraction method is not suitable for scanned paper maps with various color values. Due to the overlay of the roads, the extraction area is not continuous, but there is a certain fracture. The water colors of different tile maps are different. However, the colors of water in one map are uniform. Therefore, it is effective to separate water from other geographical elements in one map when using a given color threshold. As shown in Figure 6, the left figure is a tile map of Tiandi, and the right figure is the result after extraction by using the color threshold. However, this extraction method is not suitable for scanned paper maps with various color values. Due to the overlay of the roads, the extraction area is not continuous, but there is a certain fracture.
(2) Fracture connection of water areas Due to the overlay between water and other geographical elements on the map, such as roads, water areas may contain fracture phenomena. As shown on the left side of Figure 6, the roads overlay the river, which leads to the fracture phenomenon in the extracted river. To connect the fracture, this paper uses an algorithm based on a convex hull to connect the fracture. The basic principle of this algorithm is shown below. First, the convex points of each connected domain should be calculated. Then, when the distances between convex points from two connected domains are less than a given threshold, these convex points should be connected. Finally, the connection lines should be converted into pixels to generate new rivers. The code is shown on the left side of Figure 7. The algorithm is based on the MATLAB environment, and the meanings of the corresponding functions and variables can be found in the MATLAB official help documentation. The right side of Figure 7 is the result after using the algorithm to connect the fracture. It was found that the algorithm maintains the original width and tendency information of the river well.  (3) Hole filling of the water areas (2) Fracture connection of water areas Due to the overlay between water and other geographical elements on the map, such as roads, water areas may contain fracture phenomena. As shown on the left side of Figure 6, the roads overlay the river, which leads to the fracture phenomenon in the extracted river. To connect the fracture, this paper uses an algorithm based on a convex hull to connect the fracture. The basic principle of this algorithm is shown below. First, the convex points of each connected domain should be calculated. Then, when the distances between convex points from two connected domains are less than a given threshold, these convex points should be connected. Finally, the connection lines should be converted into pixels to generate new rivers. The code is shown on the left side of Figure 7. The algorithm is based on the MATLAB environment, and the meanings of the corresponding functions and variables can be found in the MATLAB official help documentation. The right side of Figure 7 is the result after using the algorithm to connect the fracture. It was found that the algorithm maintains the original width and tendency information of the river well.   (3) Hole filling of the water areas For some larger planar water areas, because of the overlay of features, the extraction areas may contain holes, as shown in the middle of Figure 8. To solve this problem, it is necessary to determine if the neighborhoods of the feature pixels, such as the dark blue annotation pixels in Figure 8a, contain water pixels. Pixels with other inconsistent colors will be considered as non-feature pixels and finally removed, such as the red road pixels in Figure 8b. Assuming that the width of a neighborhood is 2 pixels, and it belongs to eight connected types, if the neighborhoods of the feature pixels contain water pixels, then the feature pixels are marked as water pixels. Otherwise, the feature pixels are not marked. In an actual experiment, a relatively high accuracy can be realized when the width is set to 2 or 3 pixels. The Figure 8c is the result of using this method to fill holes. Although this method of hole filling of the water areas is simple and limited, it can satisfy the basic experimental requirements because the size and color of features on the same map, even in different areas, are usually consistent. For some larger planar water areas, because of the overlay of features, the extraction areas may contain holes, as shown in the middle of Figure 8. To solve this problem, it is necessary to determine if the neighborhoods of the feature pixels, such as the dark blue annotation pixels in Figure 8a, contain water pixels. Pixels with other inconsistent colors will be considered as non-feature pixels and finally removed, such as the red road pixels in Figure 8b. Assuming that the width of a neighborhood is 2 pixels, and it belongs to eight connected types, if the neighborhoods of the feature pixels contain water pixels, then the feature pixels are marked as water pixels. Otherwise, the feature pixels are not marked. In an actual experiment, a relatively high accuracy can be realized when the width is set to 2 or 3 pixels. The Figure 8c is the result of using this method to fill holes. Although this method of hole filling of the water areas is simple and limited, it can satisfy the basic experimental requirements because the size and color of features on the same map, even in different areas, are usually consistent.
In addition, for a general raster map, the matching of water areas can be realized according to the pixel coordinates of the same features. For a tile map in this experiment, the water areas at different levels of detail are matched by overlaying operations. Because a tile from a map at a high level is always divided into four tiles with the same size at an adjacent low level, overlapping water areas can be treated as matched objects by twofold magnification.

Study Area
In this paper, the experimental data used to detect changes in water areas at multiple scales come from the map of Tiandi [54], China. Because the lakes basically have no changes before the 7th level and after the 13th level in the Tiandi map, we chose the data with the same map extent from the 7th (1:5,000,000) to the 13th (1:100,000) level of the Tiandi map; the detailed scales of the levels between the 7th and 13th levels can be found on the official website [54]. The area is located in western Qinghai Province, China. There are many lakes and few rivers in this area. The range of the area is between 34.34° and 36.62° N latitude and between 89.97° and 92.87° E longitude. The data at the 7th level consists of a piece of tile, and the size and coordinates of this tile are 256 × 256 pixels and = 96, = 50 , respectively. With the increase in level, the number of tiles also increases in accordance with pyramid rules. The data in the highest level, the 13th level, consist of 4096 tiles. The coordinate range of those tiles is 6144 ≤ ≤ 6207 and 3200 ≤ ≤ 3263. The size of the image consisting of those tiles is 16,384 × 16,384 pixels. Figure 9 shows the tiles in the 11th level. The left panel shows the original tile data, from which In addition, for a general raster map, the matching of water areas can be realized according to the pixel coordinates of the same features. For a tile map in this experiment, the water areas at different levels of detail are matched by overlaying operations. Because a tile from a map at a high level is always divided into four tiles with the same size at an adjacent low level, overlapping water areas can be treated as matched objects by twofold magnification.

Study Area
In this paper, the experimental data used to detect changes in water areas at multiple scales come from the map of Tiandi [54], China. Because the lakes basically have no changes before the 7th level and after the 13th level in the Tiandi map, we chose the data with the same map extent from the 7th (1:5,000,000) to the 13th (1:100,000) level of the Tiandi map; the detailed scales of the levels between the 7th and 13th levels can be found on the official website [54]. The area is located in western Qinghai Province, China. There are many lakes and few rivers in this area. The range of the area is between 34.34 • and 36.62 • N latitude and between 89.97 • and 92.87 • E longitude. The data at the 7th level consists of a piece of tile, and the size and coordinates of this tile are 256 × 256 pixels and x = 96, y = 50, respectively. With the increase in level, the number of tiles also increases in accordance with pyramid rules. The data in the highest level, the 13th level, consist of 4096 tiles. The coordinate range of those tiles is 6144 ≤ x ≤ 6207 and 3200 ≤ y ≤ 3263. The size of the image consisting of those tiles is 16,384 × 16,384 pixels. Figure 9 shows the tiles in the 11th level. The left panel shows the original tile data, from which we can see the overlay phenomenon between the lake and the roads or features. The image on the right shows the situation after extraction. We find that the results after fracture connection and hole filling are good.

Result and Analysis
(1) Detection of change types According to the proposed hierarchical detection method, after multi-scale target matching, we detect the changes in the study area by calculating the corresponding indicators and using certain operations. The calculation method of each indicator is shown in Section 3.1, and the detection rule of each change type is shown in Section 3.2. Table 3 shows the results of detection.

Result and Analysis
(1) Detection of change types According to the proposed hierarchical detection method, after multi-scale target matching, we detect the changes in the study area by calculating the corresponding indicators and using certain operations. The calculation method of each indicator is shown in Section 3.1, and the detection rule of each change type is shown in Section 3.2. Table 3 shows the results of detection. Table 3 shows that from level 7 to level 10, the change types and number of changes are lower. The main change types are boundary generalization, morphology change, and hole change. There are no other change types. From level 10 to level 12, the change types and number of changes increase significantly. The main change types are appearance, disappearance, morphological change, derivative change, and hole change. Among them, the proportion of appearance changes is the highest, accounting for 56.84% of all changes. The proportions of derivative change and hole change are second, accounting for just 5.6% and 3.1%, respectively. The proportions of disappearance and morphological changes are the lowest, accounting for only 1.1% of all changes. However, the proportion of disappearance depends on the direction of scale change. From level 12 to level 13, the change types and number of changes decrease significantly. Among them, the proportion of appearance change is still the highest. The proportion of hole change is the second highest. In general, the number of changes is smaller from level 7 to level 10. The number of changes increases significantly from level 10 to level 12 and decreases significantly from level 12 to level 13. Of all the change types, the appearance change is the greatest, accounting for 84.8% of all changes. The proportion of hole changes is the second highest, accounting for 8.1% of all changes. The numbers of morphology change and derivative change are lower; together, they only account for 5.8%. The remaining detected change types disappear and are identified as boundary generalization changes. In addition, some other change types, such as moving, rotation, scaling, derivative change 1, and hole change 2, are not detected. To verify the validity and reliability of our method, we checked all 15 change types by overlaying and visual judgement methods and found that the accuracy of change type detection is 91.1%. Mainly, the wrong types come from the morphology and hole changes. In addition, as the extraction of water areas is simple and limited, many detection errors occur to a certain degree.   Figure 10a shows the results of change detection from level 11 to level 12 using the proposed method, where red represents the change type of appearance, black represents the change type of disappearance, green represents the type of morphology change 2, brown represents the type of derivative change 2, yellow represents the type of hole change 1, and orange red represents the type of hole change 3. The thresholds used for change type detection are as follows: ∆ SI = threshold ≤ 20 • , and ∆ GC = threshold ≤ 100. We also applied the map overlay [55] method for the change detection of this region. The corresponding results of change detection can be found in Figure 10b, where red represents the positive change of water areas, and black represents the negative change of water areas. Compared with the traditional map overlay method, two advantages of the proposed method should be addressed. (1) As shown in Figure 10, the proposed method can effectively detect many different types of changes, such as morphology change, hole change, and derivative change, while the traditional map overlay method can only divide the changes into two types, namely, positive and negative changes, as shown in the black box. (2) As shown in Figure 11 in the next subsection, the changes of water areas can be hierarchically expressed, which enables the reader to progressively acquire change information at different levels of detail, while the traditional map overlay method cannot achieve this. (2) Hierarchical expression of change In turn, a water area can be expressed in the body, piece, and slice layers by using the proposed hierarchical ideas. For example, consider the case of lake A marked in Figure 11.
In the body layer, the situation of global change is measured. As shown on the left side of Figure 11, the changes in the representative area value from level 7 to level 13 are displayed by gradient colors. The deeper the red color of the area is, the larger the fluctuation in the area change is. Green represents the lakes without increases and decreases in their area within a tolerance range. The figure shows that the larger the area of a lake is, the more powerful the fluctuation in the area change is. When the area of a lake is smaller, it is not easy to observe larger changes in the area. The areal variance in lake A in the box is larger. Therefore, the deeper the red color is, the larger the area fluctuation in lake A in the body layer is.   (2) Hierarchical expression of change In turn, a water area can be expressed in the body, piece, and slice layers by using the proposed hierarchical ideas. For example, consider the case of lake A marked in Figure 11.
In the body layer, the situation of global change is measured. As shown on the left side of Figure 11, the changes in the representative area value from level 7 to level 13 are displayed by gradient colors. The deeper the red color of the area is, the larger the fluctuation in the area change is. Green represents the lakes without increases and decreases in their area within a tolerance range. The figure shows that the larger the area of a lake is, the more powerful the fluctuation in the area change is. When the area of a lake is smaller, it is not easy to observe larger changes in the area. The areal variance in lake A in the box is larger. Therefore, the deeper the red color is, the larger the area fluctuation in lake A in the body layer is. Figure 11. Hierarchical expression of change information. Lake A is a typical sample. A (a) (b) Figure 11. Hierarchical expression of change information. Lake A is a typical sample.
(2) Hierarchical expression of change In turn, a water area can be expressed in the body, piece, and slice layers by using the proposed hierarchical ideas. For example, consider the case of lake A marked in Figure 11.
In the body layer, the situation of global change is measured. As shown on the left side of Figure 11, the changes in the representative area value from level 7 to level 13 are displayed by gradient colors. The deeper the red color of the area is, the larger the fluctuation in the area change is. Green represents the lakes without increases and decreases in their area within a tolerance range. The figure shows that the larger the area of a lake is, the more powerful the fluctuation in the area change is. When the area of a lake is smaller, it is not easy to observe larger changes in the area. The areal variance in lake A in the box is larger. Therefore, the deeper the red color is, the larger the area fluctuation in lake A in the body layer is.
In the piece layer, the mutation information is expressed. As shown in Table 4, according to the calculation method proposed in this paper, the mutations in lake A are found at levels 8-9, 10-11, and 12-13, where the mutation types are morphology change 1, hole change 1, and morphology change 2, respectively. Therefore, the lakes at seven levels are divided into four groups. The change between groups is larger, but the change within one group is smaller. In the slice layer, after the expression of lake A in the piece layer, the detailed changes in the objects in one group are expressed. As shown in the image on the right side of Figure 11, boundary generalization occurred in levels 7-8. There is no change in levels 9-10. Hole change 3 occurs in levels 11-12. The difference in each kind of change before and after the overall change is shown in Table 5. Table 5. Changes in the slice layer. Change type boundary generalization no change hole change 3 no change difference ∆ SI = 9 • 0 ∆ A = 230pel 0

Conclusions
During mental processes, humans usually solve problems by starting with the general scope and then focusing on the details hierarchically [34]. In this paper, we proposed the hierarchical change detection principle at multiple scales from the perspective of spatial cognition. We also classified water area changes, established a set of indicators and rules for the change detection of water areas in multi-scale raster maps, and determined the key technology for water extraction from tile maps. Compared with the traditional water change detection methods at the same or adjacent scales, the following conclusions can be drawn: (1) The proposed change detection method of water areas contains more detailed change types and can even examine the changes generated from map generalization, while the traditional method cannot.
(2) The representation of water area changes is hierarchical, which enables readers to progressively acquire change information at different levels of detail. Thus, the cognitive burden caused by displaying too much information can be well avoided.
As the proposed method for the change detection of water areas is suitable for any binary image, it can be applied for other map sources or other locations when the corresponding map elements are extracted. For some open source tile maps that have data quality problems, such as OpenStreetMap, the proposed change detection method is very helpful for data quality improvement. For example, when volunteers upload map data to the OpenStreetMap database, the proposed method can be used to filter out inconsistent data with great changes and ensure the correctness of input data. However, indexes that are used to describe the change patterns in the body layer, such as the standard deviation in area, are relatively simple, and finding more effective and reasonable measure indexes is worthy of further consideration. Furthermore, because the collapse operation is usually applied to rivers with long and narrow shapes, the detection indicators and methods are not designed for the collapse of areal lakes. In addition, we only applied the tile map with regular color values as the experimental data, and complex data sources such as remote sensing images should be explored in the future. As time series change is also an important research topic, future research should focus on the seasonal changes of water areas using tile map data.

Author Contributions:
The individual contributions of the authors are as follows: the paper was drafted by Y.S. and subsequently revised and approved by T.A. All authors have read and agreed to the published version of the manuscript.