Using High-Performance Computing to Address the Challenge of Land Use/Land Cover Change Analysis on Spatial Big Data

Abstract

Land use/land cover change (LUCC) analysis is a fundamental issue in regional and global geography that can accurately reflect the diversity of landscapes and detect the differences or changes on the earth’s surface. However, a very heavy computational load is often unavoidable, especially when processing multi-temporal land cover data with fine spatial resolution using more complicated procedures, which often takes a long time when performing the LUCC analysis over large areas. This paper employs a graph-based spatial decomposition that represents the computational loads as graph vertices and edges and then uses a balanced graph partitioning to decompose the LUCC analysis on spatial big data. For the decomposing tasks, a stream scheduling method is developed to exploit the parallelism in data moving, clipping, overlay analysis, area calculation and transition matrix building. Finally, a change analysis is performed on the land cover data from 2015 to 2016 in China, with each piece of temporal data containing approximately 260 million complex polygons. It took less than 6 h in a cluster with 15 workstations, which was an indispensable task that may surpass two weeks without any optimization.

1. Introduction

Land surface changes result from numerous driving forces, including tropical deforestation, rangeland modification, agricultural intensification, urbanization, social-political factors and the worldwide interconnectedness of places and people [1,2,3,4,5]. Telecoupling provides a new avenue of research that enables natural and social scientists to understand and generate information for managing how humans and nature sustainably coexist [6,7,8,9,10]. With global issues such as climate change, surging population growth and continued energy shortages, more focus is being directed toward global environmental changes [10,11,12]. As landscapes are enormously complex systems, most methods to study the driving forces decompose the system into different functional components and then describe the state of the landscape, the processes within the landscape, and the reactions of the landscape over time [3,13,14]. In addition, monitoring the locations and distributions of land cover change is important for establishing links between policy decisions, regulatory actions and the resulting land-use activities [15,16]. Land use/land cover change (LUCC) analysis has become a fundamental component of environmental change and sustainability research [17,18]. As stated by Li [19], except for the analysis of the driving mechanism of LUCC, other existing studies are currently dedicated to designing reasonable methods to detect the changes more accurately [11,20,21,22] or predicting the change trends of the spatial-temporal progress of LUCC [23,24,25,26], which can be performed at a fine scale [27,28,29,30,31,32] or a global scale [33,34,35,36,37].

Over the last few decades, dramatic changes have occurred throughout the world. LUCC is fundamentally a spatial process, resulting from the complex interaction of social, ecological and geophysical processes [38]. Land use maps indicating the status are more detailed for agricultural and urban land use classes than for others. Satellite data has always been the primary source for obtaining information about the earth’s land surface, and a wide variety of techniques ranging from unsupervised algorithms to parametric supervised algorithms to machine learning algorithms (MLA) have been used to classify the land cover over large areas from satellite [39,40]. Generally, post-classification comparison and multi-date composite image change detection are the two most commonly used methods in change detection [41,42,43], and the techniques of satellite remote sensing and GIS can also be used to examine the spatial and temporal patterns of the LUCC [44]. Moreover, after classification, overlaying two layers in different time series, in raster or vector format, can directly reflect the changes in pixels or polygons [45,46,47]. In addition, a hybrid classification by combining different types of methods or incorporating different data sources may also provide a better choice [48,49]. In recent years, the classification accuracy of these methods has been greatly enhanced with the increase of spatial resolution, such as very high-resolution (VHR) remote sensing images [50], and the optimization of knowledge-based approaches which play a role in use-specific photointerpretation and machine learning-based algorithms [51,52]. In the near future, it is likely that more and more robust and highly adjusted approaches will emerge for recognizing objects from images.

When multi-temporal images are available, how to detect the changes among these images has become a disparate, highly variable and ever-expanding area of research [53]. After classification, the land surface can be easily partitioned into different classes of land use/land cover parcels. On this basis, some research on change analysis addresses post classification methods, termed as a post-classification comparison [54,55]. For such comparisons, pixel-by-pixel [56,57] or object-by-object comparisons [58,59,60] are usually conducted. As suggested by Zhou et al. [59], object-based approaches provide a better means for change detection than pixel-based methods because they can provide a more effective way to incorporate spatial information and expert knowledge into the change detection process. With the resultant objects in polygon layers from the segmentation, land parcels are further split through GIS overlay functions, and then the classification can be well preserved at the object level [1,61]. Generally, two types of spatial models are widely used when representing the earth’s phenomena, vector and raster, which are the data formats for performing spatial operations [62]. Therefore, to detect the changes from multi-temporal data, GIS operations can be used on raster images, vector polygon layers, or a combination of both representing the land classifications.

With the advances in sensor and computer technology, the approach for collecting, managing and analyzing the spatial data is going through a revolution. Monitoring geographical conditions can provide consistent measurements of landscape conditions, allowing detection of both abrupt changes and slow trends over time [54]. National geographical conditions are viewed as an important part of the basic national conditions, including spatial distribution and relationships among natural elements, cultural elements, and socio-economic information [63]. Routine monitoring provides a persistent data source for the LUCC analysis, a very large amount of data and the complicated computing processes require that high-performance computing (HPC) technologies should be used to conquer these challenges. In geographical information science, geocomputation refers to the adoption of a large-scale computationally intensive approach to the problems of physical and human geography in particular, and the geosciences in general [64]. In essence, HPC is used as an important component of geospatial cyberinfrastructure (GCI) and is critical to large-scale geospatial data processing and problem-solving [65,66]. In the era of spatial big data, cluster computing, cloud computing, and GPU computing play complementary, overlapping and even crucial roles in the processing of geospatial data [67]. In fact, the utilization of HPC systems in spatial big data processing has become more and more widespread in recent years [68], such as image classification [69], land use change analysis [70], urban growth simulation [71], etc. To perform land use change analysis over massive spatial data at a national scale, a very heavy computational load is often unavoidable, especially when processing the multi-temporal land cover data with finer resolution using more complicated procedures. Particularly, to solve this geospatial problem, a theoretical approach towards capturing the computational loads can be used to direct the two necessary steps: spatial decomposition and task scheduling [72,73].

This paper presents a new approach to this problem: a cluster computing paradigm for conducting the data moving, clipping, overlay operation, area calculation and the matrix building. As a validation case, we used it to quickly detect the land use changes over China from multi-temporal layers. It should be noted that the overlay operation is performed on digital maps in a vector format for different temporal land cover layers. However, we guess that the overlay operation could also be applied to raster maps with small changes.

The remainder of the paper is organized as follows. Section 2 introduces the detailed method, including the LUCC procedures, the manner of spatial decomposition and the stream scheduling mechanism. Next, two individual experiments based on the estimating factors are given in Section 3. Section 4 provides a detailed discussion of the obtained results and finally the main conclusions of the work are presented in Section 5.

2. Materials and Methods

2.1. Procedures of LUCC Analysis Based on Multi-Temporal Vector Objects

In China, the Geographical Conditions Monitoring (GeoCM) project aims to monitor all kinds of indexes for the land surfaces in a dynamic and quantitative way, thereby achieving geographical descriptions of the spatial distributions and spatiotemporal changes of natural, economic and social factors. Since 2013, China has initiated the first national survey of geographical conditions. By the end of 2015, the complete land cover/land use information had been collected through interpretation from high-resolution imagery. Since 2016, it has entered a routine monitoring stage, and the land cover data for 2016 and 2017 are now both available.

The produced land cover data are all in vector format, which can therefore clearly depict the vegetation coverage, water and wetlands, desert and bare land, traffic networks, residential sites and facilities, etc. In 2015, for example, over 260 million complex polygons were collected, and each of the produced polygons belonged to a specific class of land cover. At a rough level, ten classes of land parcels are involved, i.e., farm-lands, garden lands, wood-lands, grass-lands, buildings, roads, constructed bodies, digging lands, desert lands, and waters, which form a complete coverage of the land area. Meanwhile, at the finest level, there are more than 80 classes involved. The detailed classification is provided in Supplementary Materials.

With the high volume of land cover data, a natural workflow is described in Figure 1, which involves five procedures, i.e., data moving, clipping, overlay of multi-temporal layers, area calculation and transition matrix building.

The data moving is the first procedure (P₁), which downloads a partition of the data from the distributed file sharing system to the local storage, according to some spatial boundaries. Particularly, the volume of the partial data to be moved should be predetermined with reference to the computing power of the local CPU resources. In fact, the following procedures should also be considered for assessing the computational load. As the original data is produced and stored in physical chunks, some redundant chunks should also be moved to a local computer, which ensures data integrity when performing the LUCC analysis according to semantic partitions, such as administrative divisions. As data moving is typically an input/output (I/O)-intensive operation, the memory consumption and time usage are both closely related to the size of the dataset.

With the available local dataset, the clipping procedure (P₂) partitions them into some semantic chunks for the subsequent overlay operation (P₃) and area calculation (P₄). Herein, the semantic chunks are usually bounded by administrative divisions. As the data clipping may cause heavy pressures both on the CPU and hard disks, its time usage is often highly relevant to the number of polygon vertices. According to our experiences, when managing a large dataset, the main proportion of the consumed time in this procedure is used in computation, instead of I/O.

The overlay operation creates a new layer according to two temporal layers produced in the last procedure (P₃). More importantly, the attributes of both original layers are assigned to the new layer by the intersection of the polygons. In other words, the new layer represents the LUCC information from different years. This operation is typically computationally intensive, and the time usage is determined by the number of polygon vertices to a great degree.

Area calculation is conducted to calculate the area for each of the intersected land parcels. Specifically, the area for each polygon is stored in a table row that contains the relevant information, as shown in

Table 1. Structure of the results for area calculation.

No.	Field Name	Notes
1	Region Code	Spatial identifier of each land parcel after intersection
2	Polygon ID	Unique identifier of each land parcel after intersection
3	Land Cover Class (from)	Land cover class before change
4	Land Cover Class (to)	Land cover class after change
5	Area	Spherical area of the land parcel after overlay operation
6	Other	A spare field

. Naturally, area calculation is closely related to the count of the vertices in each polygon. After the above four procedures, a transition matrix is built in an MS-SQL database (P₅), which can produce the matrixes for different levels of administrative districts and different levels of land cover classes through GROUP BY and SUM operators. In comparison with the first four procedures, the time used in this procedure is negligible.

2.2. Graph-Based Spatial Decomposition

2.2.1. Administrative Level Decomposition

To perform LUCC analysis over large areas, administrative districts provide a direct and effective way for partitioning the large-scale landform, and a hierarchical structure can be drawn from the study area. Shown in Figure 2 is a 4-level hierarchy of administrative districts in China.

From the top down, the country can be decomposed into some lower-level areas hierarchically. Meanwhile, from the bottom up, the transition matrix can be built level by level. Herein, each administrative district at the same level can be viewed as an independent partition of the country, and moreover, they can be processed simultaneously. However, different divisions usually have different areas, different numbers of land parcels, different numbers of polygon vertices, etc., which can be used to evaluate the time complexity of an algorithm or procedure and thus contribute to measuring the computational load for each of the of involved procedures. As the computational load may vary significantly in different procedures, a better strategy for load measuring requires referencing both the data and relevant operations in a combined way.

In fact, the administrative division is for the purpose of administration, and most countries in the world have three or more levels. Therefore, this initial decomposition can also be used in managing the data for other countries.

2.2.2. Computational Load Evaluation and Representation

Under the administrative level decomposition, each district in the hierarchical level involves a certain amount of computational load. How to quantitatively measure the load is a primary step for grouping the districts and scheduling the procedures. Among the five procedures in the LUCC analysis in Figure 1, the first four procedures are the most time consuming due to heavy loads in computations or I/O operations, and are thus performed in distributed workstations. Theoretically, the time usage in data moving and area calculation should present a linear correlation with the data volume to be processed, i.e., the number of polygons or the constituting vertices. Both the clipping and overlay are essentially intersection operations. According to Becker, et al. [74], the optimal time complexity of the overlay operation on two sets of polygons is O (N × log N + K), where N is the number of line segments, and K is the number of line intersections. When K is significantly smaller than N, a normal selection for measuring the computational loads is the number of vertices composing the polygons. In addition, I/O consumption is another factor that cannot be ignored, which is an impediment to the scalable speedups when parallelizing the overlay operations in a shared file system [75,76,77]. It is noteworthy that the jitter of the I/O process is unavoidable, which is often caused by the resource contention within multicore architectures, the communication, the kernel process scheduling and the cross-application contention [78].

Wang and Armstrong [73] proposed a theoretical approach which transforms the field-based and object-based spatial dataset to a spatial computational domain and then uses a composite function to represent the I/O consumption and computing time as compu-bands.

f = df + of

(1)

where df is a data-centric function, which transforms knowledge about spatial data characteristics in the context of spatial operations into either memory or I/O band. of is an operation-centric function, which takes into consideration the spatial operations that are directly or indirectly related to spatial data characteristics and transforms knowledge about the characteristics of spatial operations into a computing time band.

For the LUCC analysis, multiple procedures can be identified. Herein, we concentrate on how to finish the entire workflow as soon as possible, and the first four procedures are represented as four transforming functions. f₁(p) is a typical I/O operation, and its time usage is mainly determined by the data volume of a data partition p, which can be measured by the number of polygons, the number of geometry vertices, or the areas. f₂(p) and f₃(p) measure the computing time consumed in clipping and overlaying operations, and they have a similar measurement for the time complexity [79]. f₄(p) measures the area calculation, and the number of polygon vertices is naturally used as the main influencing factor.

$$f_{1,2,3,4}\left(p\right)\propto \left\{L_O\left(p\right),L_V\left(p\right),L_A\left(p\right),\dots \right\}$$

(2)

where L_O(p) is the number of polygons of p, L_V(p) is the number of polygon vertices of p, and L_A(p) is the area of p. Moreover, the Pearson correlation coefficient (PCC) is used to measure the linear correlation between the possible variables and the time usage, and then determine which factor should be selected to measure the loads. Finally, a compressive measurement is defined based on these functions.

$$T=T_{P1}\left(f_1\right)+T_{P2}\left(f_2\right)+T_{P3}\left(f_3\right)+T_{P4}\left(f_4\right)$$

(3)

where T_P1, T_P2, T_P3 and T_P4 are fitting functions used to convert the values calculated from (2) to computing time units (e.g., seconds).

During the entire workflow, the computational loads result from both the inner part of a district and the neighboring districts. Specifically, the inner loads are closely related to the inner part of the data and the operations performed on it, while the outer loads are closely related to the neighboring divisions due to the clipping operations. According to Kang and Lin [80], each district can be represented as a graph vertex with weights (i.e., the load produced from P₁ to P₄), and the relations among the districts can be represented as graph edges with weights (i.e., the load in clipping neighboring districts). After this representation, a graph that completely illustrates the data distribution and the intensity of the computational loads in the LUCC analysis can be obtained. This graph can be formally modeled by the following formula:

$$\mathrm{Graph}\left(LUCC\right)=\mathrm{G}\left(\mathrm{V}\left[i,f_v,w_v\right],\mathrm{E}\left[j,f_w,w_e\right]\right)$$

(4)

where V is a graph vertex set, and E is a graph edge set. Each vertex in V contains three elements: i (the identifier for the vertex), f_v (the vertex function) and w_v (the vertex weight). f_v denotes the specific operation upon a particular object, and w_v denotes computational load of the function on that object. Similarly, each edge in E contains three elements: j (the identifier for the edge), f_e (the edge function), and w_e (the edge weight). f_e is used to merge these results from f_v, and w_e denotes the computational load of the edge function on the spatial relation between two objects. In fact, f_e is optional and can be empty if output dependence does not exist. w_v and w_e may vary greatly for different vertices and edges due to the spatial heterogeneity. For more details about the graph representation, please refer to the definition of data dependence in spatial operations [80], which also provides an operable means for generating the graph vertices and edges.

Figure 3 illustrates a sample instance by simulating seven districts, which are all represented as graph vertices. The neighboring relations are represented as graph edges. The background grids are the originally produced chunks, each of which is filled by many land class parcels. In addition, the value in each grid is the corresponding load value, and the load of a district can be obtained by summing load values from these topologically contained or intersected grids. As a matter of fact, the shapes of these grids are often irregular, and they provide complete coverage of the entire area.

In Figure 3, the weights both exist on the graph vertices and the edges. As shown in

Table 2. Graph weight evaluation.

Type	Vertex ID	Weights
Vertex	A	WA = 2 + 5 + 2 + 6 + 3 + 5 + 6 + 5 + 4 + 3 + 5 + 6 + 6 + 4 + 2 + 4 + 5 + 5 + 5 + 5 + 3 + 4 + 4 = 99
	B	WB = 3 + 3 + 5 + 6 + 3 + 6 + 5 + 3 + 2 + 4 + 4 + 4 + 5 + 1 + 6 + 6 + 4 + 5 + 6 + 6 = 87
	C	WC = 6 + 5 + 5 + 2 + 5 + 3 + 5 + 5 + 6 + 5 + 6 + 5 + 6 + 5 + 6 = 75
	D	WD = 3 + 5 + 4 + 6 + 5 + 6 + 5 + 5 + 6 + 1 + 6 = 52
	E	WE = 5 + 5 + 5 + 3 + 4 + 4 + 5 + 4 + 5 + 6 + 6 + 2 + 4 + 2 + 2 + 3 + 6 + 5 + 5 + 6 + 6 + 2 + 6 + 6 + 5 + 2 + 6 + 5 = 125
	F	WF = 2 + 4 + 2 + 1 + 2 + 6 + 6 + 2 + 5 + 3 + 6 + 2 + 5 + 2 = 48
	G	WG = 5 + 6 + 4 + 6 + 1 + 6 + 5 + 2 + 5 + 4 + 1 + 4 + 3 + 5 + 3 + 2 + 3 + 5 + 5 + 2 + 2 + 5 = 84
Edges	A	WAB = WBA = 3 + 3 + 2 + 5= 13; WAE = WEA = 3 + 4 + 4 + 5 + 5 + 5 = 26
	B	WBC = WCB = 6 + 3 = 9; WBD = WDB = 3 + 4 + 6 + 6 = 19; WBE = WEB = 5 + 4 + 5 + 6 + 6 = 26
	C	WCD = WDC = 3 + 6 + 5 + 5 = 19; WCG = WGC = 5 + 6 = 11
	D	WDE = WED = 6; WDG = WGD = 6 + 1 + 6 + 5 = 18
	E	WEF = WFE = 2 + 4 + 2 + 6 + 6 + 2 + 2 = 24; WEG = WGE = 6 + 5 + 5 + 5 = 21

, the weight on each vertex is obtained by summing the loads of the related grids, while the weight on the edges is obtained by summing the loads of the crossing grids that intersect with the neighboring boundaries.

2.2.3. Graph Partitioning of the LUCC Analysis Graph

When the calculations are described as graphs in which the vertices and edges represent the computation and data dependence, they can be divided among the CPU cores of a parallel computer by partitioning the vertices of a graph [81]. In mathematics, the graph partitioning problem is defined in the form of a graph G = (V, E), with V vertices and E edges, such that it is possible to partition G into smaller components with specific properties [82]. Graph partitioning arises as a pre-processing step to divide and conquer algorithms, where it is often viewed as an effective way to divide things into roughly equal-sized pieces [83]. Partitioning a graph into equal-sized components while minimizing the number of edges between the different components can contribute to a balanced task assigning strategy in parallel computing [84,85]. Thus, the computational load can be divided among CPU cores. Moreover, some typical spatial operations have been optimized using the graph-based divide and conquer method by representing the geometry-level computation as vertices and dependence as edges, and then partitioning and encapsulating the sub-graphs into sub-tasks [80]. Specifically, each processing task over a spatial partition is viewed as a sub-task that may have a dependence on the other tasks over the neighboring spatial partitions.

After transforming the task of the LUCC analysis into a graph, it can thus be partitioned into several parts that can be simultaneously processed on different CPU cores. Two requirements should be met to complete the task efficiently, i.e., inter-independence and balanced load allocation. The inter-independence requires that each partitioned part can be processed independently without referencing the neighboring partitions in the computation, while balanced load allocation requires different partitions be processed in approximately equal time [73,86]. As geospatial science problems often have constraints in spatial and temporal dimensions (e.g., topological relations and order relations), relevant principles should be systematically considered in exploiting the parallelism toward optimizing the problems [86,87,88,89]. Two general rules can be concluded: (1) the independence among different parts from the physical phenomena are more intensive for those that are closer (localization); (2) the data access, computation and communication all contribute to the spatial distribution of the computational loads (heterogeneity). The localization enables the geographical domains to be divided and conquered, while the heterogeneity in the distribution of the computational loads requires that a reasonable load balancing strategy is adopted to partition the problem.

Many algorithms have been developed to find a reasonably good partition. Typical methods include spectral portioning methods [90], geometric partition algorithms [91,92,93], and the multilevel graph partitioning algorithms [94,95,96,97]. The multilevel graph partitioning algorithms work by applying one or more stages. Each stage reduces the size of the graph by collapsing the vertices and edges, partitioning the smaller graph, then mapping back and refining this partition of the original graph [97]. In many cases, this approach can provide both fast execution speed and very high-quality results. One widely used example of such an approach is called METIS [84], which was developed in the Karypis lab (http://glaros.dtc.umn.edu/gkhome/). METIS provides a set of programs for partitioning graphs, partitioning finite element meshes, and producing fill-reducing orderings for sparse matrices. The algorithms implemented in METIS are based on the multilevel recursive bisection, multilevel k-way, and multi-constraint partitioning schemes. Herein, the multilevel k-way partitioning scheme is selected as the partitioner in the proposed method. METIS was chosen because METIS’s multilevel k-way partitioning algorithm provides the advantageous capabilities of the minimized resulting sub-domain connectivity graph, enforced contiguous partitions and minimized alternative objectives.

2.3. Architecture and Scheduling

2.3.1. Master-Slave Architecture for Distributed Computers

We implemented a system employing a dynamic load balancing strategy, as shown in Figure 4. The system enables a workflow to support load measuring on the master side, spatial decomposition on the master side, data moving from the master side to multiple slave sides, data clipping on the slave sides, overlay operation on the slave sides, area calculation on the slave sides and the final matrix building on the master side.

Measuring the computational loads in the original chunks is performed on the server side. After this step, the LUCC analysis graph is constructed on the master side, and then partitioned to a certain number equal to that of the slaves. Next, the data partitions are moved from the master to the slave sides. Thus, data clipping, overlay operation and area calculation can be performed concurrently on each of the slave sides. Finally, the resultant data will be transferred to the master side.

2.3.2. Stream Scheduling of Tasks

With the above system, communications between any pair of workstations are supported. Furthermore, we adopted a stream scheduling method to exploit the parallelism among the slave workstations. Herein, each administrative division was viewed as an independent compute unit for scheduling, as illustrated in Figure 5.

Different administrative divisions often have different computational loads. The different computational loads mean the time usage for different compute units may vary greatly due to load differences. A process pool, whose size is set to the number of CPU cores on each of the slave workstations, is first initialized to support the computation. Moreover, in the pool, different units are scheduled asynchronously. As illustrated in Figure 5, four compute units were initiated at startup. When a unit finished, a new one was scheduled to push into the pool. When all the units finished, the whole task was terminated.

2.3.3. Implementation

On the side of each slave workstation, the stream scheduling was realized in the Python language and used to fully exploit the parallel performance in the process pool, which was supported in the package named Multiprocessing. A shared table was established on the server side to distribute the jobs to different workstations to work smoothly. At the same time, the working status of the workstations is continuously transferred to the server side.

Moreover, to support all the spatial operations from P₁ to P₃, we used the ArcGIS API for Python, which is a widely used Python library for working with maps and geospatial data. It provides simple and efficient tools for sophisticated vector and raster analysis, geocoding, map making, routing and directions. Furthermore, a series of error checks were developed to guarantee that all the procedures had correct inputs and outputs. Error checking enabled the workflow to work smoothly. For P₄, an area calculation was developed in the C++ language, and the File Geodatabase API provided by ESRI Company was used to access the spatial data. In comparison with the Python language, a very optimal performance could be obtained.

In addition, an MS-SQL database was used as the message database for exchanging the command directives, and the log and progress information between the master and the slave workstations.

3. Performance Evaluation

3.1. Testing Environment

A computer cluster was used to support the LUCC analysis, which included 15 workstations running Microsoft Windows 7 (×64) with a six-core Intel Xeon E5-1650 CPU (3.50 GHz) and 16 GB DDR4-2133. The total number of CPU cores was 90. In addition, Hyper-Threading Technology makes a single physical core appear as two logical cores, and the maximum number of threads for a workstation was 12. These workstations were connected through a ten-thousand-megabyte bandwidth switch.

3.2. Determination of the Factors in Estimation

To select the proper factors in measuring the computational loads, we used the PCC coefficient to measure the linear correlation between the possible factors and the time usage in each of the procedures. Values near 1 indicated that the curve shapes have high comparability. This is defined as follows.

$$\rho =\frac{\mathrm{cov}\left(X+Y\right)}{{\sigma }_X{\sigma }_Y}$$

(5)

where cov is the covariance, ${\sigma }_X$ is the standard deviation of $X$, ${\sigma }_Y$ is the standard deviation of $Y$. Herein, we randomly selected 89 counties in China, and the testing result is provided in the Supplementary Materials. Shown in

Table 3. Comparison of PCC values from different factors (The green items indicated higher correlations).

Factors\Time	TP1	TP2	TP3	TP4	T
F	0.969167	0.530038	0.568205	0.740030	0.582032
N	0.718055	0.973142	0.994402	0.964224	0.990629
A	0.513562	0.790036	0.825275	0.779288	0.809465
N × log(N)	0.706571	0.976781	0.994454	0.960296	0.991878

is the correlation analysis result obtained through calculating the PCC between the time usage in each of procedures and the potential influencing factors.

Herein, F is the count of polygons in the layers, N is the number of vertices that constitute the polygons. According to Table 3, F was the most influencing factor that determined T_P1 in data movement, and thus it was selected for measuring the loads in this procedure. As the land cover data were stored in vector objects, the movement of the data layers was realized through applying the query condition to each of the vector features. Even N should be a more effective parameter in theory for assessing the data volume, F showed a closer linear relationship with time T_P1. Thus, f₁ was set to F. Herein, a linear curve-fitting method was used to capture the trend of time with the changes in F. As illustrated in Figure 6, the time usage in this procedure was modeled using the following formula.

$$T_{P1}=0.0001743\times F+14.395$$

(6)

As the data clipping (P₂) and overlay (P₃) are essentially intersection operations, the time usage should be measured by the number of vertices. In P₂, the number of the vertices constituting the clipping boundaries was much less than that of the land cover polygons. Thus, factor N could be used in the time complexity format to measure the computational load, i.e., f₂ = N × log(N). In addition, clipping the two polygon layers involved more data I/O operations, and thus a linear fitting curve may not well match. In Figure 7, a second-order curve-fitting method was used to capture the trend of time with the change in the number of vertices.

Then, we have

$$T_{P2}={9.0}^{-17}\times {\left(N\times log\left(N\right)\right)}^2+0.0000007064\times N\times log\left(N\right)+2.5858$$

(7)

Meanwhile, in P₃, the overlay of the different layers involved two polygon layers from different years. Therefore, the number of vertices from the two participating layers was used, i.e., f₃ = N × log(N). As the time usage was mainly determined by the computation, we found a linear fitting curve was more suitable than other curves, as illustrated in Figure 8.

Then, we have

$$T_{P3}=0.0000007214\times N\times log\left(N\right)$$

(8)

In P₄, the computational load in the area calculation was positively proportional to the N, and a linear relation was most suitable, i.e., of₄ = N. The fitting curve was shown in Figure 9.

Then, we have

$$T_{P4}=0.0000008585\times N+16.631$$

(9)

Finally, the total time used in the four procedures can be measured as follows.

$$T=T_{P1}+T_{P2}+T_{P3}+T_{P4}$$

(10)

It is important to note that these time functions (T_P1, T_P2, T_P3 and T_P4) were all obtained under a certain testing hardware environment, including the CPU, the physical memory, the disk, etc. Different hardware conditions, including the CPU cores, memory bandwidth, disk transfer rate, etc., may cause certain differences in these functions. After transforming the related factors into an estimable function of the time, the computational load can be measured according to the data size and the related procedures.

3.3. Experimental Results

3.3.1. Testing data

The testing data covered all 31 provinces and autonomous regions in mainland China (Hong Kong, Macao and Taiwan regions were not covered), and 2858 counties were involved. Specifically, the data were organized into 2858 chunks, each of which was stored in an individual File Geodatabase. Moreover, these data chunks were located in a shared-disk file system, and each computer could freely access any part of the data independently. According to the time functions proposed in Section 3.2, we estimated the computational loads for all these chunks. Shown in Figure 10 is the estimated time usage for all the counties. Clearly, there were obvious differences among these counties due to the differences in the number of geometries and number of vertices. It could be seen that the computation on some large counties cost more than two hours, while computation on some small counties cost less than one minute.

In general, larger counties often had higher loads, while for some areas with complicated land cover conditions, the time usage per square kilometer showed a large difference. In other words, the computational loads in different locations were different. As shown in Figure 11, it can be concluded that in eastern and southern China, the time density was relatively higher. In fact, these areas were usually fragmented, and the representative polygons were often narrow and twisted. Thus, more points were essential for constraining the boundaries of the land cover parcels.

3.3.2. Results

Two individual experiments were conducted to verify the performance of the proposed method. The first one was to test the performance of the LUCC analysis over 7 provinces on a single workstation, which could validate the load measurement and the scalability of the program under a multi-core environment. The second one was performed in a cluster environment, which aimed at testing the validity of the spatial decomposition and the entire efficiency.

In the first experiment, Shanxi, Heilongjiang, Shanghai, Zhejiang, Shandong, Xizang and Qinghai were selected, and detailed testing results are provided in Table S3 in the Supplementary Materials. First, seven provinces were tested in a serial manner towards testing the time loads. Shown in

Table 4. Comparison of the observed time and estimated time.

Order	Name	Number of Counties	Time Usage (s)		Difference
			Observed	Estimated	Absolute	Ratio
1	Shanxi	117	54,158.173	51,664.532	2493.640	4.6%
2	Heilongjiang	128	39,550.255	37,914.468	1635.787	4.1%
3	Shanghai	16	2743.127	2529.132	213.995	7.8%
4	Zhejiang	90	28,432.039	21,489.041	6942.998	24.4%
5	Shandong	137	35,075.643	27,635.252	7440.391	21.2%
6	Xizang	74	74,659.131	77,058.161	2399.029	3.2%
7	Qinghai	46	32,129.256	34,644.451	2515.195	7.8%
-	Total/Avg	-	266,747.624	252,935.037	23,641.035	8.9%

are the observed and estimated time usage.

The maximum difference ratio was approximately 24.4%, and the minimum difference ratio was 3.2%. On the whole, the average difference ratio was approximately 8.9%. Considering some unstable factors, such as the disk transfer rate and geometry complexity, a further improvement in accuracy for measuring the loads should be difficult to achieve. Actually, the absolute differences could be further reduced when performing the test in a multicore CPU. Furthermore, we tested the performance improvement when giving different computing processes (CPU cores), as illustrated in

Table 5. Performance improvement in a multicore CPU environment (P is the number of processes).

Order	Name	Time (s)					Speedup
		P = 1	P = 2	P = 4	P = 6	P = 12
1	Shanxi	54,158.173	29,733.984	18,042.815	12,955.973	8994.652	6.021
2	Heilongjiang	39,550.255	27,368.083	19,710.915	13,951.313	10,218.074	3.871
3	Shanghai	2743.127	1658.247	1083.171	877.641	681.561	4.025
4	Zhejiang	28,432.039	15,515.316	8185.116	6340.825	4595.972	6.186
5	Shandong	35,075.643	19,257.271	10,199.049	7279.027	5601.673	6.262
6	Xizang	74,659.131	43,597.761	26,888.392	18,922.136	15,557.339	4.799
7	Qinghai	32,129.256	18,202.262	12,600.972	9032.178	6538.687	4.914
-	Total/Avg	266,747.624	155,332.924	96,710.430	69,359.093	52,187.958	5.111

, and the performance curves are shown in Figure 12.

In the second experiment, the proposed method was tested on a cluster. As all the 2858 counties would be processed, spatial decomposition should be conducted first. The target of partitioning the task of the LUCC analysis was to distribute it to different slave workstations, and then coordinate them to complete the entire task. Therefore, a partitioning method that can both maintain the balancing loads and the smallest possible cost of cutting should be pursued. Before partitioning the data, the LUCC analysis graph was first constructed by estimating the loads in the focused counties and the neighboring counties. The loads in the counties focused upon contributed to the major computational loads from P₁ to P₄, which were represented as graph vertices, while the loads in neighboring counties mainly contributed to the data moving load in P₁, which were represented as graph edges.

Thus, we obtained the LUCC graph that completely described the task. Before partitioning the graph, we stored the vertices and edges as a graph file by referencing the METIS soft manual (http://glaros.dtc.umn.edu/gkhome/views/metis). Moreover, to minimize the edge cut costs and the total communication volume, multilevel k-way partitioning was used. The LUCC graph was partitioned into 15 parts, as shown in Figure 13, to be in accordance with the testing environment.

After partitioning the graph, 15 sub-graphs were formed, and the related cutting edges were also recorded. As shown in

Table 6. Computing time loads for the partitioning domains and total loads.

Computer ID	Graph Loads (s)	Cutting Loads (s)	Costs Ratio	Task Loads (s)	Balance
192.168.3.1	78,737.800	7834.684	9.95%	86,572.484	6.52%
192.168.3.2	78,945.698	2513.631	3.18%	81,459.329	6.13%
192.168.3.3	78,444.883	6984.443	8.90%	85,429.326	6.43%
192.168.3.4	83,135.226	13,525.892	16.27%	96,661.118	7.28%
192.168.3.5	82,848.954	8524.572	10.29%	91,373.526	6.88%
192.168.3.6	80,668.012	7148.928	8.86%	87,816.940	6.61%
192.168.3.7	83,146.467	6905.600	8.31%	90,052.067	6.78%
192.168.3.8	83,042.980	6721.402	8.09%	89,764.382	6.76%
192.168.3.9	82,401.041	11,897.392	14.44%	94,298.433	7.10%
192.168.3.10	80,064.289	11,565.392	14.45%	91,629.681	6.90%
192.168.3.11	79,672.611	12,012.495	15.08%	91,685.106	6.90%
192.168.3.12	78,525.364	2017.754	2.57%	80,543.118	6.06%
192.168.3.13	78,805.905	3550.061	4.50%	82,355.966	6.20%
192.168.3.14	80,616.409	5116.888	6.35%	85,733.297	6.45%
192.168.3.15	82,511.703	10,519.652	12.75%	93,031.355	7.00%

, 15 sub-tasks were built. Then, we estimated the time loads in the sub-graphs and the cutting edges, respectively. Moreover, the total loads for each sub-task were obtained by summing the sub-graph loads and the cutting-edge loads.

According to Table 6, the LUCC analysis task was decomposed into 15 sub-tasks. These tasks were distributed in a relatively even manner. In addition, extra costs caused by decomposition were approximately 10%. The proportions of the computational loads ranged from 6.06% to 7.28%.

Table 7. Performance improvement in a cluster (All the CPU cores were used).

Slave ID	Counties	Task Loads (s)	Observed Time (s)	Speedup
192.168.3.1	176	86,572.484	16,713.703	5.180
192.168.3.2	212	81,459.329	21,245.131	3.834
192.168.3.3	189	85,429.326	18,003.370	4.745
192.168.3.4	327	96,661.118	18,944.093	5.102
192.168.3.5	294	91,373.526	17,064.811	5.354
192.168.3.6	333	87,816.940	18,688.152	4.699
192.168.3.7	207	90,052.067	17,794.495	5.061
192.168.3.8	268	89,764.382	18,774.252	4.781
192.168.3.9	197	94,298.433	19,458.336	4.846
192.168.3.10	106	91,629.681	16,454.875	5.569
192.168.3.11	88	91,685.106	18,048.192	5.080
192.168.3.12	137	80,543.118	17,363.678	4.639
192.168.3.13	71	82,355.966	18,177.326	4.531
192.168.3.14	111	85,733.297	19,216.351	4.461
192.168.3.15	142	93,031.355	16,692.415	5.573
Avg	190.5	88,560.409	18,175.945	4.872

shows the results of the entire testing data. Obviously, a large proportion of the time was reduced on each of the slave workstations, and the average time load was approximately 5.049 h. The performance improvement ranged from 3.834 to 5.573, and the maximum gap in time was approximately 4790 s, with reference to the maximum task load on workstation 192.168.3.4 and the minimum task load on workstation 192.168.3.10.

When the total task was performed on a single workstation without adopting any parallel optimization method, the estimated time was more than 14 days by summing all the estimated graph loads (i.e., the total task was performed on a single workstation, and no cutting loads were incurred). Whereas, in a cluster with the proposed method, the total time usage was approximately 5.901 h considering the maximum task load. On the whole, the speedup was approximately 57.028 in the cluster in comparison with the estimated time on a single workstation.

4. Discussion

Generally, time complexity is simply a measure of the time it takes for a function or expression to complete its task, as well as the name of the process to measure that time. However, the hardware environment should also be considered. This is because different CPU frequencies, memory bandwidth and disk transfer rates will no doubt exert an influence on the time usage. Therefore, there exists a gap between time complexity and actual time usage. Through a polynomial regression analysis, we built the relevance between the actual time usage and the data volume factors, such as the number of geometries and vertices. Moreover, the estimated time usage from different procedures could be transformed into a unified function. According to the testing results in the Supplementary Materials, the PCC coefficient for the observed time and the estimated time was more than 0.99. Thus, the estimated time function can exactly reflect the actual time usage by referencing the number of geometry vertices or the number of polygons. According to the analysis of the testing data, spatial heterogeneity is a common phenomenon that is determined both by the data and the corresponding procedures. Herein, two individual experiments were provided. The first one was performed on a single workstation, while the second one was performed on a cluster.

The results from the first experiment indicated that an average speed improvement of 5.1 times could be achieved over the sequential version with 12 processes. When the number of allocating processes was fewer than 6, a stable improvement could be achieved. However, there was no continuous performance improvement when the allocating processes were more than 6. The main reason was due to the I/O bottlenecks. When multiple processes read the data from disk or dump data to the disk with conflicts, the overall performance can easily decline due to the I/O speed limit of the storage drive. In addition, multiple procedures that were working on different CPU cores may also have caused some bad effects on the computation speed when assessing the shared memory space. Therefore, to alleviate the competition from slow disks, a flexible method is to adopt a faster drive, such as a solid-state drive (SSD). If conditions allowed, loading all the data into the main memory could definitely provide an ideal way of conquering large-scale data processing. According to our early tests [79], using a streaming framework that can coordinate the I/O events with the computation may also provide an effective way of alleviating this bottleneck. Thus, an individual procedure should be allocated to continuously reading the data and producing a continuous data block stream; then, other procedures can concurrently process these available data blocks acquired from the stream.

When performing the task on a single workstation, we focused on how to fully exploit the parallelism from multicore CPUs. In a cluster environment, the difficulty lies in how to allocate the computational loads to different workstations. Therefore, a quantitative representation of the task should first be accomplished. According to the results from the second experiment, the average difference ratio of the partitioned loads maintained a relatively low level, which demonstrated that the estimation and the partitioning method could be used to decompose the computational loads. Furthermore, two temporal land cover layers covering mainland China were decomposed into multiple sub-datasets with even loads. Thus, this extremely time-consuming task was solved in a very short period of time.

In addition, two improvements should also be made in the near future. One regards the lack of assessment on how the geometry complexity influences the time consumption. Though time functions have given an approximate estimation of the observed time, some deep rules might have been covered up. Another improvement is in terms of applicability. All the tests were performed in a homogeneous environment, and all the workstations had an identical configuration. Actually, heterogeneous environments are very common, and some adjustments should be made.

As a matter of fact, the earth’s surface is changing at an amazing speed, and thus, more and more large-scale, short-period and high-resolution data are becoming available with the development of sensor technology. The proposed method was only used to achieve a two-temporal-layer LUCC analysis, and no geomorphic layers were involved for a deep-level analysis. When dealing with more temporal layers and thematic layers, the time costs will definitely increase sharply. Therefore, methods for how to perform the change analysis for a large-scale region, such as a country, a continent or even the globe, especially using high-performance computing facilities, is still a challenge that should be addressed.

5. Conclusions

This work presented a novel graph-based spatial decomposition to conquer the very heavy computational load when overlaying multi-temporal land cover data. This approach quantitatively represented the computational loads as graph vertices and edges, and balanced subtasks were formed through graph partitioning. Basically, PCC analysis was first conducted towards selecting the determining factors that influence the time usage. Then, stream scheduling was used to achieve an efficient parallel LUCC method. Together with a master-slave architecture connecting the distributed workstations and a stream scheduling method, we performed two experiments for testing the performance improvements on a single workstation and a cluster.

Clearly, selecting the proper factors and regression methods can contribute to a unified representation that quantitatively transforms each of the procedures in a complex workflow to a time function with the data size as measurable variables. On this basis, graph-based spatial decomposition gives a balanced load allocation method, which not only considers the volume of the data to be processed, but also the spatial operations performed on the data. As a main research object, the LUCC analysis is thus optimized efficiently. It is worth mentioning that the consumed time was reduced from a two-week period to less than 6 h, for detecting the two-temporal changes over China.

In the near future, research should also be done to validate the feasibility of approach to raster maps and explore the potential influence of the geometric complexity that may also governs the parallel efficiency to some degree. In addition, a heterogeneous environment should also be considered.