GIScience Theory Based Assessment of Spatial Disparity of Geodetic Control Points Location

Abstract

Geodetic networks provide a spatial reference framework for the positioning of any geographical feature in a common and consistent way. An even spatial distribution of geodetic control points assures good quality for subordinate surveys in mapping, cadaster, engineering activities, and many other land administration-oriented applications. We investigate the spatial pattern of geodetic control points based on GIScience theory, especially Tobler’s Laws in Geography. The study makes contributions in both the research and application fields. By utilizing Average Nearest Neighbor, multi-distance spatial cluster analysis, and cluster and outlier analysis, it introduces the comprehensive methodology for ex post analysis of geodetic control points’ spatial patterns as well as the quantification of geodetic networks’ uniformity to regularly dense and regularly thinned. Moreover, it serves as a methodological resource and reference for the Head Office of Geodesy and Cartography, not only the maintenance, but also the further densification or modernization the geodetic network in Poland. Furthermore, the results give surveyors the ability to quickly assess the availability of geodetic points, as well as identify environmental obstacles that may hamper measurements. The results show that the base geodetic control points are evenly dispersed (one point over 50 sq. km), however they tend to cluster slightly in urbanized areas and forests (1.3 and 1.4 points per sq. km, respectively).

1. Introduction

Geodetic control networks are widely used for investigating and monitoring any geographical features and phenomena. They are of great priority to many human actions, like geodetic and geophysical measurements, surveying, civil and environmental engineering, GIS development, and gathering spatial data. These networks provide a series of accurate positional measures that cover the analyzed area, thus delivering a time-continuous infrastructure enabling studying and tracking of crustal deformation [1], as well as seismic and geomorphologic landform activities (e.g., [2,3,4]). The ground stability of geodetic control points provides quality to the results of geometric measures, enables an accurate and detailed description of land topography, and changes relative to Earth processes, natural hazards, and climate change [5]. Control surveys serve as the basis for initiating or reviewing subordinate surveys for property boundary delineation [6], route and construction planning [7,8], engineering object displacements [9], cadastral and topographic mapping [10], as well as many environmental applications and cultural heritage preservation and monitoring [11,12]. They are indispensable as a reference framework for land administration systems, especially for georeferencing spatial objects [13,14,15]. Many scholars point out that the production and quality control of any geographical data and products requires references to the base geodetic network that realizes the CRS (Coordinate Reference System) and provides the consistent frame for GNSS users [5,16,17]. Hence, the base geodetic control points (thereinafter refed as BGCPs) are of utmost importance for georeferencing of manuscripts of historical maps [11,18,19] or images from unmanned aerial vehicles [20].

Geodetic control points should cover an area relatively evenly to enable accurate and cost-effective measurements [21,22,23]. Although there is a variety of studies considering the design and densification of geodetic control networks [16,24,25] as well as investigating the influence of topographic objects that hinder the visibility of the horizon, interfere with satellite signals, and, finally, affect the quality of geodetic control stations’ positioning [26,27,28], a profound analysis of the spatial pattern of geodetic control points still requires investigation. The problem of geospatial distribution of geodetic control points was previously discussed in several publications [29,30,31]. However, these studies were concerned with the detail (third-order) of geodetic network point analyses in relatively small (less than 200 sq. km) rural areas. The results, related to surveying units and then grouped according to land use types, showed that geodetic control points are scattered with significantly visible groupings along roads, railways, and built-up areas. Moreover, the number and density of geodetic control points depend on the development of the area in question, and 35% to 50% depend on the land cover, mainly in locations of built-up areas, roads, and railways [30,31].

The density of geodetic control points is specified by the National Mapping Agencies. In Poland, all requirements for geodetic control networks are included in the Regulation of the Ministry of Administration and Digitization [32] implementing the Geodetic and Cartographic Law of the Polish Parliament. The document, like many other national regulations [33,34,35,36], assumes that in order to ensure an appropriate accuracy of survey measurements, the country’s territory should be covered evenly by high-order geodetic control networks. The anticipated network density is, in general, expressed by the inter-point distance or the density of the geodetic control points, and varies between counties, depending on their level of development [37].

The objective of this study was to investigate the spatial pattern of the geodetic control points of the Base National Geodetic Network (BNGN) in Poland, in particular, by performing an ex-post analysis of the uniformity of spatial distribution and determination of factors that disturb the assumed evenness. This is the first comprehensive study covering the entire country using GIScience methods, i.e., incremental spatial statistics, cluster and density analysis, as well as first and second order statistics. This study makes contributions in both research and applications. Based on a GIScience theoretical background, it introduces a comprehensive methodology for ex post analysis of geodetic control points’ spatial pattern as well as the quantification of geodetic network uniformity into regularly dense and regularly thinned. The conducted tests of the first and second order statistics provide, in the form of graphs and maps, exhaustive information on the spatial distribution type (regular, random, or clustered) of the geodetic points. The study is intended to be a methodological resource and reference for the National Mapping Agency and the Head Office of Geodesy and Cartography for the maintenance, as well as the further densification or modernization, of the base geodetic network in Poland. Furthermore, the results give surveyors the ability to quickly assess the availability of geodetic points, as well as identify environmental obstacles that may hamper measurements. Moreover, the method of investigating the spatial pattern of the BGCPs could be applied in any country or region, and after a small modification, on any man-made point feature of which the location is assumed to not be distributed randomly.

The remainder of the article is structured as follows: Section 2 provides an overview of the National Geodetic Control Network in Poland, Section 3 presents the area and data used, the fundamentals of the point pattern analysis, and steps in data analyzing, Section 4 and Section 5 describe the obtained outcomes in detail and provide discussion with the results, and finally, Section 6 concludes this study.

2. Overview of the National Geodetic Control Network in Poland

The control networks in Poland, maintained by the Head Office of Geodesy and Cartography, are designed and established according to the Geodetic and Cartographic Law [38] and the pertinent implementation rules established by the Ministry of Administration and Development [32]. Information about geodetic control points is stored in spatial databases, and includes the point number and unique identifier, coordinates, survey uncertainty, and information about their source, type of monumentation, and topographic details (e.g., location, access details, sketch plans); enabling the recovery or reconstruction of the point in terrain, and allowing the maintenance of the basic control network.

The geodetic network in Poland consists of points of the horizontal and the vertical networks. The horizontal network is hierarchical and comprises three orders of geodetic control networks. The main criterion for assigning a geodetic control point to the appropriate network order is the accuracy of determining its position, expressed as the mean error. The fundamental (first-order) network consists of 127 GNSS stations that belong to the Active Geodetic Network–European Position Determination System (ASG-EUPOS) out of which 24 are located in the neighboring countries [39]. These stations meet the criteria established by the Subcommittee EUREF and provide an average density of one point per 20,000 km². The base geodetic control is the second-order geodetic network providing the density of one point per 50 km². The first and second-order geodetic networks, according to the Ministry Regulation [32], form the so called base (main) geodetic network in Poland [40]. All points belonging to these networks are monumented and have precisely measured positions (0.01 m for horizontal position and 0.02 m for geodesic height [32]) as they form the basis for determining the positions of other points. All BGCPs should be situated on stable ground, far from sources of electromagnetic interference or surface vibrations (e.g., caused by heavy traffic), objects that cause reflection of satellite signals (e.g., water reservoirs, flat metal surfaces, walls, wire fences, high antenna towers), as well as terrain barriers above 10 degrees over the horizon. The detailed geodetic network is a third-order control with 1,378,383 points and the average density of control points of ca. one point per 22.5 ha. The mean error of the third-order point’s location is in the range from 0.05 to 0.10 m, however for newly established points it should be less than 0.07 m. The control points are accurately tied to at least three reference geodetic controls. The geodetic control networks are the backbone of the Polish Spatial Information Infrastructure [40,41].

3. Materials and Methods

3.1. Area of Investigation and Data Used

The study covers Poland, a central European country with the total area of 312,679 km². Poland has the eighth largest and one of the most dynamic economies in the European Union [42]. It is highly developed both socially and economically, with a Gini coefficient equal to 26.8 [43]. Simultaneously, the country holds a very high (33rd) worldwide position in the Human Development Index (0.865) [44]. The country is generally agricultural and flat, with more than 76% of the territory below 200 m above mean sea level. Agricultural land covers 60.2% of its surface, followed by forests and bushes that take as much as 30.7%; water—2.1%; wastelands—about 1.5%; and built-up and urbanized areas—5.5%. Among the built-up and urbanized land, communication areas have the largest share of about 56% [45]. Data used in the study are summarized in the

Table 1. Spatial data characteristics.

Data Set	Abbreviation	CRS 1	Data Format	Date	Data Provider
Base National Geodetic Network	BNGN	ETRS89	txt	2016	Central Geodetic and Cartographic Documentation Centre
CORINE Land Cover	CLC	ETRS89	shape	2018	Copernicus Land Monitoring Service
Vector Smart Map level 2	VmapL2	WGS84	VPF	2018	Polish Military Geography Directorate

¹ Coordinate Reference System.

and described below.

Coordinates of the 6723 BGCPs stations were obtained from the Central Geodetic and Cartographic Documentation Centre, in Warsaw. The points’ coordinates are given both in ETRS89 and Cartesian coordinate system “PL1992” (the EPSG code 2180). The geographic locations of points are shown in Figure 1a.

The CORINE Land Cover (CLC) delivers information on the biophysical characteristics of the Earth’s surface and its changes derived on visual interpretation of satellite images. The land cover nomenclature is hierarchical, including three levels of thematic details grouped in five major land cover types: Artificial, agricultural, forests and semi-natural areas, wetlands, water bodies, with 44 classes at the third, most detailed level. The CLC data is available from COPERNICUS land monitoring services [46] free of charge. The CORINE Land Cover (Figure 1b) data was used to investigate a relationship between BGCPs’ density and land cover.

VmapL2 delivers vector-based geospatial data at medium resolution (1:50,000). Data are structured in nine thematic layers, i.e., boundaries, elevation, hydrography, industry, physiography, population, transportation, utilities, vegetation [47]. The study utilized data on road and rail networks, as subsets of the transportation layers.

3.2. Research Hypotesis and Methods

Based on an in-depth literature and national regulation review, we hypothesized that the spatial pattern of the base geodetic control points in Poland depends on the scale of investigation. Country-wide it could be perceived as regularly dispersed, although, locally it is regularly densified or regularly thinned. Moreover, based on CLC and VmapL2 data, we discover the BGCPs that are located in the signal interference zone due to environmental factors, such as water reservoirs, heavy traffic paved roads, electricity lines, and other geographical features. This is of particular importance for both surveyors and mapping authorities, because some advanced measurement and calculation techniques for measuring BGCP position should be recommended. CLC data allows also to investigate dependence between the Thiessen polygons’ area and shape and the land cover structure.

The hypothesis was proven using GIS tools based on Tobler’s laws of geography, in particular, the first law that states ‘everything is related to everything else, but near things are more related than distant things’. This law serves as the basic concept for spatial statistical methods that are vastly relevant for spatial analyses and modeling, e.g., global and local indicators of spatial autocorrelation (Moran’s I, Anselin Local Moran’s I), and spatial interpolation (kernel density estimation). The second law (‘the phenomenon external to an area of interest affects what goes on inside’) emphasizes spatial heterogeneity, particularly scale-dependent concentration and decentralization [48].

The BGCPs spatial pattern was investigated by testing complete spatial randomness using nearest-neighbor analysis and Ripley’s K-function, and global Moran’s I statistics. The regularity of geodetic point locations analysis was based on morphological parameters (area, shape and nearest distance between two adjacent polygons) of Thiessen polygons generated from the BGCPs. As a foundation for the quantification of the geodetic network evenness into regularly dense and regularly thinned cluster, the outlier analysis (Anselin Local Morans I) was implemented. The adopted methodological framework is briefly shown in Figure 2.

3.2.1. Point Pattern Analysis

Point pattern analysis (PPA) is a set of methods for the investigation of spatial arrangements of points in space, usually two-dimensional, and could be expressed as a set X = {x ∈ D} where D, the study region, is a subset of R², a two-dimensional Euclidean space [49]. Diggle [50] defined a spatial point pattern as a set of locations, irregularly distributed within a study region, which have been generated by random mechanisms. Moreover, he formulated the three basic simplified spatial point patterns: Complete spatial randomness (CSR), regularity, and clustering that are very useful at an early stage of spatial analysis. Literature presents many methods of spatial point pattern analysis, classified by Cressie [51] in five main categories: Quadrat methods, distance methods, kernel estimators of intensity, nearest neighbor distribution functions, and K-function analyses. As noticed by [52], each method has its pros and cons, and works well for certain datasets or regions.

Three methods, namely: Kernel density, nearest neighbor analysis, and Ripley’s K-function, were used to examine the spatial geodetic control point pattern in this study. They are briefly described below. Kernel density, a nonparametric estimation, belongs to the local density techniques based on a first-order property of the pattern, which concerns itself with the variation of the observations’ density across a given area [53]. The kernel method computes a localized density for subsets of the study area; the sub-regions overlap one another, providing a moving sub-region window, defined as a kernel. The density, λ$\left(s\right),$ at s (a vector location anywhere in R and $s_1,\dots , s_n$ are the vector locations of the n observed events) is estimated as [54]:

$${\widehat{\lambda }}_{\tau }\left(s\right)={\sum }_{i=1}^n\frac{1}{2}k\left(\frac{s-s_i}{\tau }\right)$$

(1)

where, k() represents the kernel weighting function; s, the center of moving window; and τ, which is greater than 0, the bandwidth.

Deng et al. [55] demonstrated that the choice of the bandwidth τ has a critical impact on the results of density analysis. As noticed by Gatrell et al. [54], increasing the bandwidth results in smoothening of the spatial variation in intensity, while its reduction results in an increasingly ‘spiky’ estimate.

The nearest neighbor analysis and Ripley’s K-function belong to the distance methods. They explore how the points are distributed relative to one another (a second-order property of the point pattern) as opposed to how the points are distributed relative to the study extent [56,57,58]. The PPA is tested against the hypothesis of complete spatial randomness, under which the spatial pattern is a realization of a Poisson point process. An average nearest neighbor (ANN) analysis measures the average distance from each point in a given area to its nearest point which compares to the distances between nearest points and distances that would be expected on the basis of chance. If the average distance is less than the average for a hypothetical random distribution, the distribution of the analyzed features is considered clustered. If the average distance is greater than a hypothetical random distribution, the features are considered dispersed [59]. The nearest neighbor analysis, developed in 1950 [54], is one of the oldest distance statistics. The null hypothesis, tested by the two-tail test of CRS, is rejected if and only if:

$$[z_{ANN}]>z_{p/2}$$

(2)

where, $z_{ANN}$—observed distribution of analyzed sample of points, $z_p$—standard normal distribution, and p—the significance level. A $z_{ANN}$ value lower than $z_p$ shows that the analyzed point pattern is more clustered, while a value larger than $z_{p/2}$ indicates a pattern more dispersed than random.

The K-function, known as Ripley’s K-function, examines all inter-point distances instead of computing separate nearest neighbors. It illustrates how the spatial clustering or dispersion of point change when the neighborhood size changes. For data analysis, the variance stabilized Ripley K-function, called the L function, is generally used (Equation (3)):

$$L\left(d\right)=\sqrt{\frac{A{\sum }_{i=1}^N{\sum }_{i=1,j\ne 1}^N{k}_{i.j}}{\pi N\left(N-1\right)}}$$

(3)

where, d—the distance, N—total number of points (features), A—the total area, and k_i,j—a weight connected with edge correction function.

Particularly, in this study, Ripley’s K function counts the number of neighboring BGCPs found within a given distance of each individual point. The number of observed neighboring BGCPs is then compared to the number of points that are expected based on a completely spatially randomness. The advantage of the function is that it allows for the investigation of how BGCPs point pattern distribution is changing with scale. However, it should be remembered that patterns are suspect at larger distances due to edge effects.

3.2.2. Thiessen Polygons

Measuring the evenness of the BGCPs’ spatial location is based on the analysis of the area and shape of Thiessen polygons formed over the geodetic points’ site. Thiessen polygons, also called Voronoi diagrams, define individual areas of influence around each of a set of points, hence any location within a Thiessen polygon is closer to its linked point than to any other [60]. This property has caused that Thiessen polygons are perceived as a proper construct to test the CSR hypothesis as well as the regularity of spatial distribution [61,62]. In considering the equal division characteristics of the Thiessen polygon in spatial tessellation, the method can be used to investigate such issues as nearest points, proximity, adjacency, and accessibility analyses [63].

The shape indices (SI) of Thiessen polygons were calculated according to the equation, firstly introduced by [64]:

$$SI=\frac{P}{2\sqrt{\pi A}}$$

(4)

where, P—denotes the perimeter of the Thiessen polygon, and A—the area of the polygon. The SI value 1.128 refers to the square-shaped polygon, while the value 1 denotes a circle-shaped.

The diversities in the shape, area, and land cover structure of Thiessen polygons were measured by the coefficient of variation, CV (coefficient of disparity), a dimensionless number that quantifies the degree of variability relative to the mean [65]:

$$CV=\frac{\sigma }{\mu }100\%$$

(5)

where, CV is the coefficient of variation, σ—the standard deviation, and μ—the mean.

Clusters of small and big Thiessen polygons were found by analyzing the statistical distribution of the area and shape of the polygons. The diversity of Thiessen polygons land cover structure, expressed as percentage of artificial areas, agriculture, forest and bushes, water bodies, and marshes, as well as land cover diversity (i.e., number of classes—NC) and fragmentation (i.e., number of patches—NP), was assessed by the statistical measures of central tendency, position, and dispersion (mean, median) and the coefficient of disparity (CV).

3.2.3. Spatial Autocorelation, Cluster and Outlier Analysis

Spatial autocorrelation on the entire study area was investigated by the Global Moran’s I index according to Equation (6) [66].

$$I=\frac{n{\sum }_{i=1}^n{\sum }_j^n{w}_{ij}\left(y_i-\overline{y}\right)\left(y_j-\overline{y}\right)}{n{\sum }_{i=1}^n{\sum }_j^n{w}_{ij}{\left(y_i-\overline{y}\right)}^2}$$

(6)

where, n is the number of attributes; y_i is the attribute value in region i; y_j is the attribute value in region j; $\overline{y}$ is the mean value of the attribute in the study site; and w_ij are weights according to the neighborhood of the attributes. Moran’s I values closer to −1 indicate negative or inverse spatial autocorrelation between the Thiessen polygons’ morphological parameters, i.e., adjacent Thiessen polygons are characterized by dissimilar parameters values. The values closer to +1 indicate positive or direct spatial autocorrelation and highlight the regions with similarly high or low attribute values. Moran’s I value 0 shows a global spatial randomness pattern of distribution.

The local indicator of spatial autocorrelation of Thiessen polygons was explored by the means of Anselin Local Moran’s I (LISA) statistics in line with Equation (7) [67]:

$$I_i=\frac{\left(y_i-\overline{y}\right){\sum }_{j=1}^n\left(y_i-\overline{y}\right)}{\raisebox{1ex}{${\sum }_{j=1}^n{\left(y_i-\overline{y}\right)}^2$}\!\left/ \!\raisebox{-1ex}{$n$}\right.}$$

(7)

where, n is the number of attributes; y_i is the attribute value in region i; y_j is the attribute value in region j; $\overline{y}$ is the mean value of the attribute in the study site.

LISA is generally used for the decomposition of Moran’s statistics, the global measures of autocorrelation. The results, presented as a choropleth map, depict (by the COType attribute) statistically significant (95 percent confidence level) cluster of high values (HH) and low values (LL), as well as outlier where a high/low value are surrounded by low/high values (HL/LH). Based on the COTtype attributes’ values, the regularly dense and regularly thinned BGCPs’ patterns were defined.

The regularly dense cluster (RD) is the sum of the Thiessen polygons that belong to the subset formed by the COType LL of the Thiessens’ area (LL_AREA) and nearest neighbor distance between geodetic points (LL_NND) minus the Thiessen polygons of shape index (SI) value lower than $SI<\overline{SI}$-${\sigma }^2$, i.e., SI outliers (SI_out), as seen in the Equation (8). The elimination of very irregular polygons makes it possible to exclude the polygons that are adjacent to the country border, thus eliminating the edge effect [54]. The regularly thinned cluster (RT) is composed of the Thiessen polygons which COType attribute of an area and a nearest neighbor fields take value HH (Equation (9)).

$$RD = ({\mathrm{LL}}_{\mathrm{AREA}}\cup {\mathrm{LL}}_{\mathrm{NND}})\backslash {\mathrm{SI}}_{\mathrm{out}}$$

(8)

$$RT=H{H}_{AREA}\cup H{H}_{NND)}$$

(9)

The RD regions with a closer location between BGCP (LL) and RT—more distant (HH)—are shown on a map. This map, available as a web application, could help surveyors to assess the cost and duration of measurement work.

ArcGIS 10.5 and Statistica 13.1 were used for conducting the spatial and statistical analysis.

4. Results

4.1. General Characteristic of BGCPs’ Location

The Base Geodetic Control Network in Poland comprises 6723 points, including 103 ASG-EUPOS stations. The geographical points location highly (with the coefficient of determination equals to 0.92) corresponds to the main land cover type, namely: Artificial areas, agriculture, forest, and bushes (see

Table 2. Global density of BGCPs over CORINE land cover types.

CORINE Land Cover Types	Number of BGCPs (Frequency)	Area [km2]	Density BGCPs/1km2	Density BGCPs/50 km2
Artificial area	452	17,605	0.026	1.284
Agriculture, including:	3417	179,307	0.019	0.953
Arable land	2931	128,966	0.023	1.136
Pastures and meadows	486	50,341	0.010	0.483
Forest	2692	96,209	0.028	1.399
Bushes	141	5669	0.025	1.243
Others (e.g., bare rock, dispersed vegetation)	21	6701	0.003	0.157
Poland	6723	312,544	0.022	1.076

and Figure 3a). As many as 57.37% of BGCPs are situated on agricultural land, 30.78% in forest, and 5.65% in urban areas. There are relatively few points on pastures and meadows. This could be explained by the difficulty of maintaining a point’s monumentation stability in lands with high flood risk and variable groundwater levels, which, according to [68], is typical for Polish meadows and pastures, located mainly in river valleys.

The average density of the BGCPs in Poland amounts to 1.08 point per 50 km² and is in line with the national regulation [32].

A total of 6256 points are situated no further than 200 m from the transportation network (

Table 3. The number of BGCPs in the vicinity of land cover objects.

The Buffer Zone [m]	Number of BGCPs Near Roads	Number of BGCPs Near Railways	Number of BGCPs Near Water Bodies	Number of BGCPs Near Electricity Lines
10	157	0	2	19
20	343	2	5	33
50	966	22	15	103
100	1853	69	60	215
200	3094	168	219	474
total	6256	261	301	844

, Figure 3b), of which as many as 343 (5.1%) are within a 20 m buffer around paved roads, and two points are closer than 20 m from the railway.

The relatively close (up to 200 m) distance of as much as 93.05% of geodetic control points to paved roads makes them easily accessible to surveyors, which is undoubtedly a great advantage of the Base Geodetic Control Network in terms of its use and maintenance.

Nevertheless, there is a set of 383 BGCPs located no further than 20 m from environmental objects that could provide significant disruption of an electromagnetic signal, and finally, decrease the accuracy of positioning (see Table 3 and Figure 4).

Local density of BGCPs varies from 0.2 to 3.65 points over 50 km² and shows places of higher and lower point density, which are underlined in red and orange in the kernel density map (Figure 5a). A significant increase in the density of geodetic points is observed in major metropolitan agglomerations and heavily industrialized regions, such as mining areas (Figure 5b).

The highest density is observed in the Upper Silesian Polycentric Metropolitan Area (from 3 to 3.5 points per 50 km²) followed by Warsaw, Lodz, Wroclaw, Szczecin, Bydgoszcz-Torun-Grudziadz, tri-cities metropolitan areas, where the density of BGCPs varies from 1.2 to 2.85 points per 50 km². The lowest values (less than 0.5 over 50 km²) are mainly found in the Southern Pomerania Lakelands, Sandomierz and Nida Basins, as well as the North Polesie Plain, which are mostly agricultural and forested regions with a high share of natural vegetation and relatively sparsely populated areas, and the percentage of urbanized areas lower than 1.17 (Figure 5b).

4.2. Point Pattern Analysis

The null hypothesis presuming Complete Spatial Randomness (CSR) for BGCPs’ pattern analysis was rejected, using Average Nearest Neighbor tools. The $z_{ANN}$ score amounts to 1.39, with a 99% confidence level, which indicates that there is less than 1% likelihood that the dispersed pattern of the BGCPs could be the result of random chance and that the points are dispersed over the territory of Poland. The observed mean nearest inter-point distance equals 5.583 km, while the expected amounts to 4.013 km. Moreover, the coefficient of variation CV amounts to 21.66%, which means that the nearest distances between BGCPs do not vary significantly (see Figure 6).

Ripley’s K function (Figure 7) is estimated for distances up to 20 km, in 1 km increments, and indicates two different departures from randomness.

At the distances from 2 up to 7 km, typical for base geodetic control networks, L(d) lies below the expected distance value (higher than expected approximately 2000 m) and indicates statistically significant (with the level of confidence of 99%) evidence of regular dispersion of geodetic points. At longer distances (greater than 7 km) L(d) lies above the expected distance, indicating spatial clustering. In general, at near distances, the BGCPs are expected to be dispersed, as the number of points within a given distance of each individual point is small, as the distance increases, each BGCP would typically have more neighbors, hence they are clustered.

4.3. Thiessen Polygons Morphology and Land Cover Structure

The analysis of the Thiessen polygons area and shape also emphasizes that the distribution of BGCPs in Poland is generally dispersed. The shape of polygons, expressed as the shape index (Eq.4), is almost equal. The SI values vary from 0.095 to 0.935 (range 0.840), with the mean of 0.826, and median of 0.840 (see

Table 4. Morphology and land cover characteristics of Thiessen polygons (confidence level 95%).

Variables	Mean	Median	Minimum	Maximum	Range	Quartile (Range)	Variance	Std. Dev.	Coef. Var. [%]
Thiessen polygons area in km2	46.480	46.636	3.685	120.806	117.120	16.881	227.271	15.076	32.434
Shape index of Thiessen polygons	0.826	0.840	0.095	0.935	0.840	0.004	0.005	0.070	8.450
Heavily urbanized [%]	6.65	3.32	0.00	99.96	99.96	5.23	122.22	11.05	166.16
Sparsely urbanized [%]	2.79	1.70	0.00	37.11	37.11	3.22	11.97	3.46	123.97
Agriculture [%]	53.40	57.51	0.00	99.90	99.95	39.10	650.53	25.51	47.76
Forest and bushes [%]	32.45	26.52	0.00	99.98	99.98	35.64	613.72	24.77	76.35
Water bodies [%]	2.13	0.00	0.00	97.54	97.54	1.88	36.78	6.06	285.18
NC	8.94	9.00	1.00	16.00	15.00	2.00	3.60	1.90	21.23
NP	33.90	32.00	1.00	109.00	108.00	19.00	232.67	15.25	44.99

, and Figure 8a). Moreover, 50% of polygons are characterized by the SI index that fluctuates from 0.812 (lower quartile) to 0.863 (upper quartile). H-spread (IQR) amounts to 0.051, showing very low dispersion, which is also visible in the CV index value, equal to 8.45. The exclusion of polygons neighboring the country border entails even smaller diversification of the SI, which underlines the almost double reduction of the CV index from 8.45% to 4.48%.

The analysis of the Thiessen polygons’ size distribution indicates little deviation from the normal distribution. The mean equals 46.48 km², and the median is 46.64km². Skewness (0.051) shows the fairly symmetric character, with a little longer right tail, while kurtosis (0.896) indicates the fewer and less extreme outliers than in the normal distribution. The area of 50% polygons ranges from 3.685 to 46.636 km², and 90% does not exceed 64.4 km² (Figure 8b). The coefficient of disparity (CV) of the polygons’ area amounts to 32.4% and emphasizes their medium differentiation.

The analysis of land cover structure in the Thiessen polygons reveals high diversity, emphasized by the CV coefficient, which varied from 285.18% for water bodies to 47.76% for agriculture areas. The land cover dispersion in dense urban areas (CV = 166.16) is higher than in sparsely urbanized ones (CV=123.97) and was followed by forested lands (CV = 76.35). On average, the number of CLC classes (NC=9) is nearly four times lower than the number of patches (NP = 34), which suggests that the landscape heterogeneity is lower than its fragmentation. Moreover, there were some polygons covered entirely by one land cover patch (NP amounts to one). On average, 53.4% of the area of Thiessen polygons is agricultural, followed by forests (32.4%). Urban areas occupy about 6.6%, and dispersed settlement 2.8%. Water bodies (lakes and rivers) cover only 2.1%.

Artificial land (heavily and sparsely urbanized), on average, occupies 9.44% of the Thiessen polygons. However, the smallest polygons, with an area of less than 16.5 km² (mean − 2.0 std. dev.), are covered by artificial surfaces in 32.06%, while for very large polygons (bigger than 77 km², mean − 2.0 std. dev.) artificial land covers only 4.64% (see Figure 9). The Pearson’s correlation coefficient (PCC), equal to −0.28 (significant at p < 0.05000), indicates the weak downhill linear correlation between the size of the polygons and acreage of urbanized areas. The level of urbanization, however, is moderately correlated with road density (PCC equals 0.50). Nevertheless, in small polygons, the average share of urbanized and agricultural areas amounted to 30% and the average share of forests amounted to 20%, which shows a fairly balanced land cover structure.

A moderate positive linear correlation (with the significant at p < 0.05) was observed between a polygon’s size and the number of patches (PCC = 0.56), while, the relationship between a polygon’s size and the area of agricultural land was very weak (PCC = 0.14). Likewise, the linear relationship between the shape index the land cover structure was similarly very weak (PCC lower than 0.15).

4.4. Quantification of Geodetic Network Uniformity

The Thiessen polygons’ variability in size and shape analysis supported by cluster and outlier analysis (Anselin Local Moran’s I), and preceded by a global autocorrelation assessment, made it possible to delineate the places of more or less densified geodetic networks. Global Moran’s I statistics, based on Thiessen polygons’ locations and attributes values, such as area, shape, and nearest distance to neighbor polygons’ centroids, returned positive I and z-score values (I = 1.06/z = 204.43, I = 0.24/z − 45.18, and I = 1.08/z = 207.78, correspondingly) with a 99% level of significance. This empowers rejection (with a probability greater than 1%) of the null hypothesis, assuming that the Thiessen polygons characterised by the attribute being analyzed are randomly distributed in a study area that covered all of Poland.. Therefore, the spatial distributions of high and low area, shape, and NN distance (NNd) values are clustered. Such characteristics of the spatial distribution of polygons was the basis for further analysis, namely determining the type of the polygon clusters, due to the high/low values of the attribute of surrounding polygons (see Figure 10a–d).

Out of 6723 polygons, 3914 (58.2%) created statistically significant clusters of autocorrelated polygons, characterized by the HH/LL or HL/LH area values. As much as 1234 Thiessen polygons (18.3%) create HH clusters (light red), and 1130 create LL clusters (light blue) (Figure 10a). The number of autocorrelated Thiessen polygons is less if the nearest distance between neighboring polygons is analyzed; a total of 2054, which is 30.6% of all polygons. As much 46.0% (945 points) form LL clusters and 870 (42.4%) form HH clusters. These Thiessen polygons comprise 14.1% and 12.9% of all polygons, respectively. If the shape (i.e., the SI) of polygons is considered (Figure 10c), polygons of LL and HL are mainly located along the country border. They constitute 43.4% of autocorrelated polygons, and just 10% of all of them. A further 746 polygons form HH clusters. However, LL and HL clusters are rather small, with an average of 11 polygons, so they were eventually excluded from further analysis aimed at finding local BGCP Thiessen polygon clusters.

The largest clusters of the large and regular Thiessen polygons (HH COType values of the attributes area and shape), named the regularly thinned clusters (RT), are located in the North-West (Pomerania Lakeland, Northern Podlasie Plain), East (east part of Northern Masovia Lowland), south-east (Lublin and Roztocze Upland), South (Sandomierz Basin, Nida Basin, Wester Beskidy Foothills), and the central part (Masovia Lowland and Hills) of the country.

The RT regions (shown in green Figure 10e) are mainly covered by forest, agriculture, or agricultural-forest lands, and occupy 92,644 km² (29.7% of Poland). The regions of regularly densified geodetic control points (RD, presented in red, Figure 10e) are created by small, regularly shaped Thiessen polygons. They are composed of a smaller number of polygons than the RT regions (

Table 5. Descriptive characteristics of regularly thinned cluster (RT) and regularly dense cluster (RD) regions.

BGN Regularity Type	RT	RD
Number of Thiessen polygons	1504	1130
Area in km2 (in %)	92,643.9 (29.7)	30,920.1 (9.9)
Thiessen polygons min area in km2	9.29	3.68
Thiessen polygons max area in km2	120.80	46.44
Thiessen polygons mean area in km2	61.60	27.26
Thiessen polygons area frequency distribution
NN distance mean value in m	6789	3770

), and are more dispersed throughout the country territory, covering mainly urbanized areas, i.e., large agglomerations and accompanying industrial areas. Tiny RD clusters are located along the country’s border, mainly West (Odra river bands), North (along Baltic coastal zone), and South (the Sudetes, the Tatra Mountains and the Western and Eastern Carpathians). The two largest clusters are located in Silesia Upland and Oświęcim Basin as well as in the central part of the Masovia and Podlasie Lowlands. RG regions cover just 9.9% of the country and comprise 16.8% of all Thiessen polygons. When it comes to the polygons’ size, RD regions are more homogeneous than RT clusters, because Thiessen polygons’ area ranges from 3.68 to 46.44 km² with an average value of 27.35 km², while, in RT, is takes values from 9.29 to 120.80 km² (see Table 5).

Base Geodetic Network densification and thinning also emphasize the shortest distances between adjacent (near neighbor) points, which ranges from 5840 to 20,580 m in RT clusters, and from 145 to 5560 in RD clusters. Distances below 500 m occur in Upper Silesia, in areas threatened by land surface deformation related to hard coal mining.

5. Discussion

New technological advances in GIScience have made it easier to justify Tobler’s laws of geography. Now, we are able to measure almost everything, but the choice of technology and methods has a significant impact on the result of the analysis and final conclusions. Tobler’s laws give the background for GIScience theory, while GIS software provides tools for spatial analysis. In our study, Tobler’s laws of geography are read as a statement about the form of geodetic control points’ spatial arrangement, particularly in the sense of similarity. Regardless of the nature and type of geodetic control networks, their task is to create a basis for determining the position of various objects on the Earth’s surface. Preferably, the entire area concerned should be covered at a uniform density with a simultaneously adjusted network of control survey stations. Both ANN and Ripley’s K-function are best suited to analyzing the type of spatial distributions. They help avoid the modifiable areal unit problem [69] as they treat space as continuous, without constraining data within any regions or units, and, contrary to the quadrat method, they are free from statistical bias caused by data aggregation [70]. Nevertheless, Ripley’s K and NN functions describe different aspects of a point process. In particular, processes with the same K(t) function may have different nearest-neighbor distribution functions, and vice versa.

As stated in many national regulations and standards, the spacing of the higher-order stations can vary from 5 to 16 km [23,33,71]. The conducted research proved that, in Poland, such requirements are fulfilled, with the mean observed interstation distance equal to 5.6 km. Moreover, each geodetic control point should be monumented with substantial stable survey monuments, and situated far from sources of electromagnetic interference and ground vibration [5,25,32]. The analysis revealed that these requirements are fulfilled, although, 5.7% of the base geodetic control points, which are located in the near vicinity of geographical features that could provide significant disruption of electromagnetic signal, have to be monitored on a regular, short-term basis in order to ensure a high accuracy of their position.

The review of world-wide scientific literature concerning geodetic points’ location, as well as international and national regulation on network design and monumentation, revealed that intertwining between global, high-, and low-order networks is clearly visible. This is due to the growing demand for high-precision location measurements of objects, especially engineering, which is ensured by GNSS-based measurement techniques and methods [71]. Therefore, there is a clear tendency to place detailed geodetic control points (low-order) in a position that maximizes the use of various measurement techniques, especially those that are intended for use in GNSS observations [21,28,72]. High-order geodetic stations are dispersed regularly [25,39,73,74], and this evenness is only slightly disturbed in densely urbanized areas [23] or those with challenging environmental conditions [75]. At the local level, based on Poland as a case study, the clustered distribution and its relation to land cover/land use is evident and well documented [22,29,30], and the analysis revealed high clustering along roads, water courses, and in dense urban areas. High-order geodetic control networks, with regularly dispersed reference stations, define the geodetic reference system in Poland [16]. To the contrary, local, detailed networks, that are a part of land administration systems, are tailor-made to the covered area, and are suitable for performing local measurements, including surveying deformations and the necessity to monitor engineering structure [28,72]. They are inevitably designed to match local purposes and scale.

Scale is, de facto, a fundamental concept in any analysis concerned with spatial pattern and processes [48,66] as the patterns and processes of spatial objects have an implicit scale of variation. As noticed by Goodchild [76,77], scale is often perceived as an alternative dimension to interpret the patterns and processes. Scale also implicite relates to land policy, especially in land administration domain systems and models [15,78]. Scale dependent variation of BGCPs’ spatial pattern is shown by Ripley’s K-function and spatial autocorrelation. Furthermore, the spatial pattern of geodetic control networks is locally hampered by environmental features or processes (e.g., surface deformation as a result of mining, vicinity of water reservoirs, sky obstruction by a forest canopy) and changes from regular to densified or thinned.

6. Conclusions

The novelty of this study relied on quantifying the observed dispersed spatial pattern of the base geodetic control points (BGCPs) into regularly dense and regularly thinned, depending on a cluster and outlier analysis. The rules and principles of designing geodetic control networks are well recognized and described in the literature, but mainly as theoretical, numerical examples. Our research examines the spatial distribution of geodetic control points monumented within the frame of the base geodetic network in Poland. Hence, it has a broader context and provides some results concerning the spatial pattern of the BGCPs in Poland, which slightly changes within scales. The spatial dispersion of the BGCPs was documented by an ANN analysis, one of the oldest spatial statistics. The z-score of 1.39 for the ANN analysis, placed in the tail of normal distribution, provides clear evidence for dispersion. This dispersion was confirmed by the second-order spatial statistic, namely the Ripley’s K-function, which shows that at the near distances from 2 to 7 km, which are typical lengths in base geodetic networks. This leads to the conclusion that the assumed uniform distribution of BGCPs over the country’s territory has been achieved, however, in regions where the inter-point distance equals 2 km on average, a considerable densification of BGCPs is observed.

Only 6% of the BGCPs could be influenced by environmental obstacles, such as close vicinity to express roads, water reservoirs, or electricity lines. Moreover 452 points are located in artificial areas, where the satellite signals may be disturbed in multiple ways. All these points should be monitored regularly to ensure that their coordinates are determined accurately.

The presented research faces challenges in examining the configuration of Polish base geodetic control points by means of advanced spatial pattern methods, especially in analyzing the distribution of basic control network points. The achieved results could serve the scientific community and also land surveyors. They could be helpful at the survey planning stage because they show the accessibility of points in terms of inter-point distance, geometric distribution, and environmental conditions in the vicinity of stations.

The originality of the research lies in its complex approach to analyzing the density of geodetic control points and the evenness of their spatial distribution. Although the methods used are well known and widely described in the literature, their combination is an added value, especially for the National Mapping Agency that is responsible for geodetic control networks maintenance as well as for surveyors. Ultimately, the goal of the present research is not only to quantify the uniformity of spatial distribution of the base geodetic control points in Poland, but also to find environmental obstacles that hamper the assumed regular dispersed distribution and, finally, to provide methodological fundamentals enabling to conduct the analytical procedure for geodetic control points spatial pattern assessment.