Spatial-Temporal Analysis of Point Distribution Pattern of Schools Using Spatial Autocorrelation Indices in Bojnourd City

: In recent years, attention has been given to the construction and development of new educational centers, but their spatial distribution across the cities has received less attention. In this study, the Average Nearest Neighbor (ANN) and the optimized hot spot analysis methods have been used to determine the general spatial distribution of the schools. Also, in order to investigate the spatial distribution of the schools based on the substructure variables, which include the school building area, the results of the general and local Moran and Getis Ord analyses have been investigated. A di ﬀ erential Moran index was also used to study the spatial-temporal variations of the schools’ distribution patterns based on the net per capita variable, which is the amount of school building area per student. The results of the Average Nearest Neighbor (ANN) analysis indicated that the general spatial patterns of the primary schools, the ﬁrst high schools, and the secondary high schools in the years 2011, 2016, 2018, and 2021 are clustered. Applying the optimized hot spot analysis method also identiﬁed the southern areas and the suburbs as cold polygons with less-density. Also, the results of the di ﬀ erential Moran analysis showed the positive trend of the net per capita changes for the primary schools and ﬁrst high schools. However, the result is di ﬀ erent for the secondary high schools. The results indicate that the schools are generally more densely populated in the central areas of the city, but it was found that the location of the hot and the cold polygons di ﬀ ered. At the neighborhood level, the schools within the neighborhood of Ghiam, Nader, Amiriyeh, Berengy, Bargh, Enghelab, Vosogh, Daneshsara, Behzisti, Park Shahr, Madrese Ferdowsi, Sareban Mahale, Sfa, Mofkham, Padgan Artesh, North 17 Shahrivar, Jafar Abad, Mirzakuchek khan, Manba Ab, Sayyidi, Hosseini Masoum, and Koi Imam Reza were more dense, and the density of schools was lower in the neighborhoods of Malaksh, Halghe Sang, Shahrak Alghadir, Koi Janbazan, Ghale Aziz, and Ahmadabad. At the statistical districts level, the density of the schools in the districts located in Ghiam, Hosseini Masoum, Sayyidi, Manba Ab, Mirzakuchekkhan, Bargh, Berengy, Amiriyeh, Nader, Padegan Artesh, Enghelab, Vosough, Daneshsara, Madrese Ferdowsi, Sarban Mahale, Safa, Mokhofam, and Sharak Golestan neighborhoods was higher, and in the statistical districts located in Bagh Motahari, Shahrak Hekmat, Imam Khomeini township, Koi Sadeghiyeh, Bagher khan 1, Koi Azadegan, and Mohammadabad Korah, the school density was lower. However, at the level of statistical blocks, the density of the schools in the blocks located in the neighborhoods of Imam Reza, Koi Behdari, Ferdowsi, Mosalla, Nader, Amiriyeh, Berengy, Mirazkochekkhan, Manba Ab, Daneshsara, Park Shahr, Vosough, Enghelab, and Mokhofam was higher.


Introduction
The rapid growth of the urban population due to the migration of people from the rural regions to urban centers and the lack of a well-planned and systematic planning system has caused many problems in most cities in Iran. Due to the rapid growth of the population and the physical condition of the cities, issues, which include the lack of proper spatial distribution of the facilities, have arisen. One of the most important problems is the reduction of the per capita urban services, which include the educational services. The public services have different types and roles regarding the support of society functions [1]. A problem that planners often deal with is determining the optimal distribution of the service centers. The equal distribution of the service centers, which include the educational centers, is especially important because of the critical role they play in the knowledge, development, and culture of each society.
The schools are crucial components of public city services [2]. The educational services facilitate learning and teaching, which are both physical and spatial. They are the resources that educational

Literature Review
The equity of education is a crucial factor to make sustainable education growth. The geographical literature has conducted abundant research about investigating the spatial patterns of different kinds of services and facilities, but little research about the educational spatial patterns has been conducted.
Talen [10] found significant inequalities with the access to primary schools in the West Virginia regions. They investigated whether or not the distribution of travel costs among the citizen's places (blocks) and the educational facilities is equally based on the density and demographics of the population. Thomas et al. [11] investigated the inequity in educational facilities using the Gini coefficient and the Theil coefficient. Cao et al. [12] evaluated the spatial distribution and the equality of educational facilities in Gansu. Ajala et al. [13] examined the inequity of access to primary education in a state in Ethiopia using different data, which included the number of schools, enrollments, and teachers. They discovered that accessibility to primary education services was not equitable.
Meschi and Scervini [14] presented a broad set of educational attainment and inequity measures. As a metric to measure inequality, they used standard deviation, the coefficient of variation, the Gini index, the Theil index, the mean logarithmic deviation, and the Atkinson indices. Yu et al. [15] used the Gini index to calculate the inequalities. Pengfei et al. [16] identified the regional variations in the accessibility of elementary schools in Xiantao, which indicated that the city centers have better accessibility to elementary schools than rural areas. Lobban [5] analyzed the spatial, the aspatial, and the demographic data in order to identify potential locations for secondary schools in order to guarantee an equal educational facilities distribution. Yang et al. [17] assessed elementary school accessibility considering the impact of roadway slopes in mountainous regions.
Scott et al. [18] analyzed public school spatial inequality based on public transit. Moreno-Monroya et al. [19] found that inequity in the central business district of Sao Paulo in Brazil could be worsened by the unequal accessibility of the secondary schools. They found that the students with low incomes often have limited accessibility to schools, because they live in places with inadequate transportation substructures. Xiang et al. [4] discovered disadvantaged groups and their distributions. They also identified possible differences in education accessibility and spatial patterns. Ye et al. [20] examined the accessibility of educational facilities at the neighborhood level in China among the various social groups. Their results showed strong spatial inequalities with accessibility to secondary school. Zhou et al. [9] used the Gini coefficient to calculate the inequity in educational facilities.
Some research used a spatial pattern analysis. Gao et al. [21] investigated the spatial accessibility based on the shortest distance from schools to the inhabitants. They also explored the inequity in compulsory education with the distribution of the spatial accessibility. They used the Theil index and the spatial autocorrelation analyses. They used global and local Moran's I and Getis-Ord Gi*. The study of the spatial patterns showed significant global autocorrelation and apparent clusters. Setyono and Cahyono [22] evaluated the number of school catchment areas in each unit, which included the villages and the sub-districts. They used Moran' I spatial autocorrelation analysis to discover the pattern of the primary school facilities. Agbabiaka et al. [3] examined the spatial distribution of the secondary education services. They investigated the educational facilities' availability, accessibility, and functionality. They used the Average Nearest Neighbor (ANN) to examine the spatial pattern of the secondary school facilities distribution. Bulti et al. [6] investigated the spatial distribution and the accessibility of the elementary school facilities. To analyze the educational facilities' spatial concentration and the spatial distribution pattern, they used the location quotient and the Average Nearest Neighbor (ANN) analysis. They indicated the inequity of the service delivery between the neighboring areas and showed a clustered pattern in the overall spatial distribution of the educational facilities. Sumari et al. [23] investigated the spatial distribution of the town of Abbottabad's primary and secondary schools in Pakistan. They found an unequal pattern distribution with the schools.
Even though educational inequity has attracted wide academic attention, little attention has been paid to the educational facilities' spatial-temporal patterns. Adams and Hannum [24] investigated variations in the accessibility of the educational facilities and the inequality over time in China using a Sustainability 2020, 12, 7755 4 of 26 generalized estimating equation method. Williams and Wang [25] assessed high school accessibility for the years 1990, 2000, and 2010. They investigated the percentage change in the educational accessibility in various regions. Relative to the suburban and the rural regions, they discovered that the urban areas had significantly lower accessibility levels. Yang et al. [26] used the Gini index for several years to measure China's educational inequity.
Wan et al. [27] investigated the regional inequity in elementary education. They also investigated the temporary inequity by calculating the number of elementary educational facilities for each year. They used the Theil index to assess the elementary schools' educational facilities' balances in all counties. Xu et al. [28] introduced a method for social-spatial accessibility in order to evaluate the various degrees of accessibility across the geographical, opportunity, and economic aspects for the various elementary educational facilities. For the years 2008 and 2018, they used a spatial-temporal analysis. Comprehensive studies examining the spatial-temporal patterns have been conducted in other types of facilities, such as parks, green spaces, and public libraries. Xing et al. [29] evaluated the spatial-temporal differences between the demand and the supply of two modes of traffic, which included walking and driving, in the city center in Wuhan between 2000 and 2014. They used the coefficient of variation to measure the alterations, and they used the analysis of spatial correlation to clarify the impacts of spatial differences by alterations in the accessibility of each mode. They concluded that the supply and the demand for green space both increased. They found the spatial disparities of the green park spaces were significantly different in walking and driving modes in various periods. They analyzed the factors of the changes in inequalities using a bivariate analysis, which confirmed that policies and the strategies were the critical factors relevant to green space access. They used the spatial autocorrelation method, which included Moran's I global univariate indices for spatial-temporal autocorrelations of disparity and Moran's I bivariate indices for correlations of spatially correlated accessibility disparity variables.
Ye et al. [20] investigated the variations in the distribution of green space accessibility among regions and among communities between the years 2010 and 2015. They concluded that access to green spaces in Macau is unevenly distributed across regions and population groups. However, green space access has become more equally distributed over time despite the decrease in the overall access to green space.
Page et al. [30] investigated the spatio-temporal differences in the accessibility of the public libraries. They indicated that between 2011 and 2018, the degree of spatial accessibility to public libraries in Wales decreased on average. This matter was correlated with a significant decrease in public library financing during the same period.
This type of research showed spatial inequality with the distribution of critical facilities and services, and the findings have implications regarding the decisions.

Basics of the Used Methods
The necessary concepts of the methods used in the research will be described in this section.

Average Nearest Neighbor (ANN)
In order to determine the type of distribution of the educational centers in this study, the Average Nearest Neighbor analysis method was used. The results are important, because they facilitate the identification of the causes that influence the formation of the different deployment systems. The results also provide basic services and essential services needed to implement the general policies and complete the suggested models for the delivery of educational services at the regional level. It should be noted that as a result of applying the different steps, this method yields an index called the Rate neighborhood (Rn), which ranges from zero to 2.15. The closer the value of Rn is to zero indicates the cluster distribution pattern. The closer the value of Rn is to 2.15 indicates a regular distribution pattern, and the number 1 expresses the random spatial distribution. In the Average Nearest Neighbor method, if the index is less than 1, the distribution is clustered, and if the index is greater than 1, there is a Sustainability 2020, 12, 7755 5 of 26 tendency to disperse. The Average Nearest Neighbor method is very sensitive, and a small change in the distribution of land use can cause many changes in the computations [31].

Optimized Hot Spot Analysis
Optimized hot spot analysis implements the Getis-Ord Gi* method use variables extracted from the input data features. This method uses the Getis-Ord Gi* statistics to generate a map of statistically meaningful hot and cold spots. To generate the optimum outcomes, it assesses the attributes of the input feature class. This method can be done without selecting a field or by selecting a numerical field. the optimized hot spot analysis is conducted by applying the following two approaches in this research [32]: • The events count method inside each cell: this approach aggregates the event point data and computes the number of events inside each cell.

•
The event count method within each polygon: to aggregate events into count parameters, there must be aggregation polygons to overlay the event point data in the polygons. The events are enumerated across each polygon. The method of counting point events within each polygon is commonly used to examine the spatial concentration of the data at different statistical levels, such as statistical districts and statistical blocks.

Moran's I
This is a general index that shows whether the pattern is clustered, scattered, or random. Each aggregated point is allocated an intensity amount and needs some variation of the values for the measurement of this metric. The points that have similar values are expressed in the high Moran I values, which are positive or negative. The Moran index values vary between +1 and −1. The values that approximate to +1 demonstrate clustering, and the values that approximate to −1 demonstrate random distribution of the data patterns. The significance of the values can be evaluated against a normal distribution. The output of this analysis consists of five parameters of the observed values, which include I, the expected values, the standard deviation, the z-score, and the p-value. In this index, if the observed value of I exceeds the expected value, it indicates a cluster. If the observed values of I are lower than the expected values, it indicates the absence of a cluster. In addition, if the p-value probability is low and the z-score is high, the null hypothesis is rejected and indicates that the data distribution is clustered [33]. Equation (1) illustrates how the Moran index is calculated [34]: In Equation (1), z i is the standard deviation of feature i from the mean (x i − X), is the location of the two features i and j, n is the total number of features, and S 0 is the sum of the weights of the features (in Equation (2)).
The z-score can be calculated according to Equation (3). In Equation (3), V [I] is the variance.

Getis Ord General G
The statistic identifies the clusters that are more likely than the randomly generated clusters and uses the aggregated data. For example, if an area has high values for the number of schools and its neighbors are the same, this test concludes that they are part of a hot spot. This analysis is most Sustainability 2020, 12, 7755 6 of 26 often used when there is a nearly uniform distribution of the values and we needed to search for the unexpected spatial values. The general G statistic is considered according to Equation (4) [35]: In Equation (4), x i and x j are descriptive values for i and j, and w i, j is the weight of the location among i and j. The z-score for this statistic is also calculated according to Equation (5). In Equation (5), 3.5. Differential Moran's I The differential Moran's I is an alternative to the bivariate Moran index between a variable at one point in time and its neighbors at the preceding time, and it controls for the locational fixed effects through the difference of variables or, in other words, by calculating the spatial autocorrelation of the variables (i) at two times t and t-1 (y i,t − y i,t−1 ). Thus, we consider the correlation when it can be the outcome of a fixed impact. For example, a series of parameters that determine the value of y i that stay fixed through the time period.
Particularly, if µ i is the fixed impact related to the position i, the value of any position at time t can be summed as the intrinsic value and the fixed value as y i,t = y * i,t + µ i . However, µ i does not vary over time. The existence of µ i causes a temporal correlation between y i,t and y i,t−1 , which exceeds the correlation between y * i,t and y * t−1 . By considering the first difference, the constant effect is removed and guarantees that any residual correlation is based solely on y * .
In particular, the differential Moran method is also a method of identifying the type of spatial data distribution that compares a similar variable's spatial autocorrelation at two various times. This method is used to detect the clustering or the randomness of a variable's variations in two dimensions: location and time.
The scatter plot of differential Moran is easily broken down into 4 segments. The upper-right segment and the lower-left segment align to the positive spatial autocorrelation, which are the equivalent values at adjacent locations. These parts are considered to be the high-high spatial autocorrelation and the low-low spatial autocorrelation, respectively. In comparison, the lower-right and the upper-left segments align to the negative spatial autocorrelation, which are different values at the adjacent locations. They are considered the high-low and the low-high spatial autocorrelations, respectively. The slope of the differential Moran scatter plot represents the correlation coefficient of changes of a variable over time [36].

Anselin Local Moran's I
In addition to the general Moran index, there is also the local Moran index, which is used to determine the location of clusters. Moran's local statistics explain the pattern of the spatial relation of a spatial parameter in the neighborhood. In other words, this statistic determines the degree of autocorrelation between the adjacent cell values in a geographical area and tests its significance. Its calculation formula is similar to the regression formula where there are two variables x and y, and the autocorrelation between them is calculated, so the difference of each variable from the mean is obtained. However, here x and y are a point and its neighbor; in fact, we use neighbors as well. The local Moran index is considered according to Equation (6). In this respect, the matrix w is a proximity matrix that is considered to be a neighbor and adjacent if the number is 1 and not adjacent if the number is 0 [33]. The local Moran index is distinct from its general index and is calculated for each cell according to Equation (6). Similar to the general Moran index, the results can be analyzed using the z-score [8].
In Equation (6), x i is a value for feature i, X is the mean of the descriptive values, w i,j is the spatial weight between two features i and j where the sum of weights is 1, and n is the total number of features. The z-score is also calculated for this statistic in Equation (7): In this statistic, the standard z-score is calculated and tested at the confidence level. In Equation (8) In this analysis, if the value of I i is positive and significant, it indicates that the existing cells are encircled by similar cells. If I i is a large positive number, it indicates that there are strong clustering zones, so that the similar values have formed clusters. On the other hand, the significant and negative values of I i indicate that the feature is encircled by the value-less features, which are not at all similar in values. These types of features are called scattered. The presence of these types of effects indicate a negative spatial correlation [37]: The local Moran index examines the relationship between the points and the neighbors, which can occur in the following four situations illustrated below.

•
High-High: When both its own spatial autocorrelation and its neighbors are positive. In this group, the high-value points are surrounded by the high-value points. This analysis can identify clusters that have the highest values and are encircled by high values. • High-Low: When it is a spatial autocorrelation that is positive, and its neighbors' autocorrelation is negative. In this group, the high values are encircled by the low values. • Low-High: When it is a spatial autocorrelation that is negative, and its neighbors' autocorrelation is positive. In this group, the points with low values are encircled by the high value points. • Low-Low: When both their spatial autocorrelation and their neighbors are negative. In this group, the low points are encircled by other points, which also have low values.
Points where their autocorrelation and their neighbors are both high or both low are important for the analysis.

Getis Ord Gi*
This tool measures the concentration of the high or low points in a given area. It can distinguish between the hot spots and the cold spots. In other words, when Gi* is a high number, it indicates a hot spot, and when Gi* is a low number, it indicates a cold spot. In other words, if the observed Gi* values exceed the expected values, it will form a hot spot, and if the observed Gi* values are below the expected values, it will form a cold spot [33]. Getis Ord Gi*local statistics are computed according to Equation (9) [34]: In Equation (9), x j is the descriptive value for j, w i,j is the spatial weight among i and j, and n is the whole number of features. X represents the mean value (Equation (10)), and S is the standard deviation (Equation (11)). The G*statistic is itself a z-score, so no further calculations are needed.

Materials and Methods
In this section, the study area will be introduced and the methodology will be examined.

Study Area
Bojnourd City is in Bojnourd County, located in the center of North Khorasan province. According to the administrative divisions, this city includes 3 districts, 11 regions, and 74 neighborhoods. The area of the Bojnourd city is 21,634,693 m 2 , and its population was 228,931 in the year 2016.
The statistics of the educational spaces in Bojnourd are discussed in the following section. There are 19 governmental educational units and ten non-governmental educational units for girls' primary schools, and the service area of the governmental primary schools for girls is 15.7% of the entire city, and 3.9% of the city is covered by the non-governmental girls' primary schools. There are 19 governmental educational units and eight non-governmental educational units for the boys' primary schools, and the service area of the governmental and the non-governmental primary schools for boys is 13.9% and 7.9%, respectively. In the boys' and the girls' mixed primary schools, there are seven educational units, and their service area is 11.6% of the area of the city.
The educational system in Iran is a set of programs, methods, educational systems, and basic regulations. These are applied by the education organization in the form of two secondary and primary education courses. The school is established based on the criteria of the education organization, and it is responsible for the approved education programs at a certain educational level. In terms of the school fees in Iran, they are divided into two categories, which include governmental schools and non-governmental schools. The students are free to attend governmental schools, and students can only attend these schools if they meet the age requirements and are within the service area of these schools. The non-governmental schools are schools that are established and managed through the participation of the people, which adhere to certain goals, rules, programs, and the general instructions of the ministry of education under the supervision of that ministry. These schools operate under the supervision of the ministry of education, and the school fees are collected from families in accordance with the rules and regulations.
There are 16 governmental educational facilities, seven non-governmental educational facilities, and one exceptional facility for the girls' junior high schools. The service area of the governmental, non-governmental, and exceptional junior high schools for girls is 30.4%, 10.8%, and 1.2%, respectively. For the boys' first junior schools, there are 15 governmental educational facilities and four non-governmental educational facilities, and the service area of the governmental and the non-governmental junior high schools for boys is 41% and 14.3%, respectively. There are 12 governmental and three non-governmental educational facilities for the boys' secondary high schools, and their service area is 47.6% and 18.11%, respectively. For the girls' secondary high schools, there are 12 governmental and five non-governmental educational facilities, and their service area is 43.8% and 16.6%, respectively. The study area is depicted in Figure 1. First, the location of Iran on the world map is shown in Figure 1a. Second, the location of the North Khorasan Province in Iran is shown in Figure 1b. Third, the location of Boujnourd County in the North Khorasan Province is highlighted in Figure 1c. Finally, the overall distribution of the schools in Bojnourd city is shown in Figure 1d.

Methodology
The main goal of this paper is to compare and evaluate the different types of point pattern analyses to investigate the spatial distribution of the schools. In this regard, for the spatial-temporal analysis of the school distribution patterns, the spatial and the descriptive information of the Bojnourd city schools, such as gender, educational level, region, district, neighborhood, number of students, and the number of classes was collected by the department of education of the North Khorasan Province and completed by field visit. In this regard, focusing on the point patterns, two types of analysis were used to investigate the general pattern of the spatial-temporal distribution of the schools. The first type of analysis investigated the general pattern of the spatial-temporal distribution of the schools without specific parameters. However, the second type of analysis was based on certain parameters, which included the substructure and the net per capita. For the first type of analysis, the Average Nearest Neighbor method and the optimized hot spot analysis were used. In the second type of analysis, the results of the general spatial autocorrelation analysis, such as Moran's I index, Getis Ord, and the differential Moran, were evaluated to detect the clustering or the randomness of the school distribution based on the different variables. Also, in order to determine the location of the clusters, Anselin local Moran's I and Getis Ord Gi* were used. The differential Moran method was also used to visualize the changes of the hot spots over time. The above analyses were performed using ArcGIS 10.4.1 and Geoda 1.12 software. Figure 2 shows the steps of the methodological flowchart that will be described below. In the Flowchart illustrated in Figure 2, there are three phases. In phase 1, the types of data used in this study, which include the covering point and statistical data, are displayed. In phase 1, all the required data is collected and prepared for its intended purpose. In this phase variables, descriptive data and statistical data are selected for analysis in the next phase. All types of analyses that were performed, which include spatial and temporal, are described in phase 2. In phase 2, all the analyses are performed for 2011, 2016, and 2018, and their results are compared. For the differential Moran, the 2016-2017, 2017-2018, and 216-2018 time intervals are compared. Due to the fact that the last statistical data was officially collected by the Statistics Center of Iran in 2011 and 2016, these two years were selected. Also, the latest descriptive information about the schools was related to 2018, which is why this year was selected. Finally, the

Methodology
The main goal of this paper is to compare and evaluate the different types of point pattern analyses to investigate the spatial distribution of the schools. In this regard, for the spatial-temporal analysis of the school distribution patterns, the spatial and the descriptive information of the Bojnourd city schools, such as gender, educational level, region, district, neighborhood, number of students, and the number of classes was collected by the department of education of the North Khorasan Province and completed by field visit. In this regard, focusing on the point patterns, two types of analysis were used to investigate the general pattern of the spatial-temporal distribution of the schools. The first type of analysis investigated the general pattern of the spatial-temporal distribution of the schools without specific parameters. However, the second type of analysis was based on certain parameters, which included the substructure and the net per capita. For the first type of analysis, the Average Nearest Neighbor method and the optimized hot spot analysis were used. In the second type of analysis, the results of the general spatial autocorrelation analysis, such as Moran's I index, Getis Ord, and the differential Moran, were evaluated to detect the clustering or the randomness of the school distribution based on the different variables. Also, in order to determine the location of the clusters, Anselin local Moran's I and Getis Ord Gi* were used. The differential Moran method was also used to visualize the changes of the hot spots over time. The above analyses were performed using ArcGIS 10.4.1 and Geoda 1.12 software. Figure 2 shows the steps of the methodological flowchart that will be described below. In the Flowchart illustrated in Figure 2, there are three phases. In phase 1, the types of data used in this study, which include the covering point and statistical data, are displayed. In phase 1, all the required data is collected and prepared for its intended purpose. In this phase variables, descriptive data and statistical data are selected for analysis in the next phase. All types of analyses that were performed, which include spatial and temporal, are described in phase 2. In phase 2, all the analyses are performed for 2011, 2016, and 2018, and their results are compared. For the differential Moran, the 2016-2017, 2017-2018, and 216-2018 time intervals are compared. Due to the fact that the last statistical data was officially collected by the Statistics Center of Iran in 2011 and 2016, these two years were selected. Also, the latest descriptive information about the schools was related to 2018, which is why this year was selected. Finally, the results of the analyses described in phase 2 are presented in phase 3. In phase 3, the results of the spatial distribution and the cluster discovery are investigated. In other words, both the type of distribution and the location of the clusters in this phase are determined.
Sustainability 2020, 12, x FOR PEER REVIEW 10 of 29 results of the analyses described in phase 2 are presented in phase 3. In phase 3, the results of the spatial distribution and the cluster discovery are investigated. In other words, both the type of distribution and the location of the clusters in this phase are determined.

Results
In the following, the results of the different methods used in this study are presented.

Results
In the following, the results of the different methods used in this study are presented.

Implementation of the Average Nearest Neighbor and the Optimized Hot Spot Analysis
To investigate the overall distribution of the schools for years 2011, 2016, 2018, and 2021, the Average Nearest Neighbor was used. In this analysis, which does not include any variables, the Rate neighborhood (Rn) index is calculated throughout the city area based on the average distance of each point to its nearest neighbors. As previously mentioned, the closer the value of the Rate neighborhood (Rn) index is to zero, the more it represents a clustered pattern, and the closer it is to 2.15, the more it shows a scattered pattern. An index value of 1 also indicates a random pattern. Based on the results shown in Table 1, because the Rate neighborhood (Rn) index is close to zero, the overall distributions of the schools in the years 2011, 2016, 2018, and 2021 are clustered. The Average Nearest Neighbor was also analyzed separately for the boys' and the girls' primary schools, first high schools, and secondary high schools in the years 2011, 2016, and 2018. In this part of the analysis, the data for the year 2021 was excluded due to the uncertainty of the gender of the schools to be exploited by 2021. Based on the results of the analyses presented in Tables 2-4    In the case of the girls' secondary schools, the overall distribution of the schools is random in the years 2016 and 2018. However, for the year 2011, the distribution of the schools is clustered. In the case of boys' schools, the overall distributions of the primary schools, first high schools, and secondary high schools for years 2011, 2016, and 2018 are clustered.
A subsequent analysis that was to investigate and show the overall distribution of the schools is the optimized hot spot analysis. As previously mentioned, the calculations of this analysis are based on the Getis Ord Gi* index and can be done in two ways. According to Figure  In the case of the girls' secondary schools, the overall distribution of the schools is random in the years 2016 and 2018. However, for the year 2011, the distribution of the schools is clustered. In the case of boys' schools, the overall distributions of the primary schools, first high schools, and secondary high schools for years 2011, 2016, and 2018 are clustered.
A subsequent analysis that was to investigate and show the overall distribution of the schools is the optimized hot spot analysis. As previously mentioned, the calculations of this analysis are based on the Getis Ord Gi* index and can be done in two ways. According to Figure    By applying optimized hot spot analysis using the point count method within each polygon, which is the school counting method within each polygon, the autocorrelation of the schools at different levels, such as neighborhoods, statistical districts, and statistical blocks was investigated for By applying optimized hot spot analysis using the point count method within each polygon, which is the school counting method within each polygon, the autocorrelation of the schools at different levels, such as neighborhoods, statistical districts, and statistical blocks was investigated for the year 2018. The results indicate that the schools are generally more densely populated in the central areas of the city, but it was found that the location of the hot and the cold polygons differed. At the neighborhood level, the schools within the neighborhood of Ghiam, Nader, Amiriyeh, Berengy, Bargh, Enghelab, Vosogh, Daneshsara, Behzisti, Park Shahr, Madrese Ferdowsi, Sareban Mahale, Sfa, Mofkham, Padgan Artesh, North 17 Shahrivar, Jafar Abad, Mirzakuchek khan, Manba Ab, Sayyidi, Hosseini Masoum, and Koi Imam Reza were more dense, and the density of schools was lower in the neighborhoods of Malaksh, Halghe Sang, Shahrak Alghadir, Koi Janbazan, Ghale Aziz, and Ahmadabad. At the statistical districts level, the density of the schools in the districts located in Ghiam, Hosseini Masoum, Sayyidi, Manba Ab, Mirzakuchekkhan, Bargh, Berengy, Amiriyeh, Nader, Padegan Artesh, Enghelab, Vosough, Daneshsara, Madrese Ferdowsi, Sarban Mahale, Safa, Mokhofam, and Sharak Golestan neighborhoods was higher, and in the statistical districts located in Bagh Motahari, Shahrak Hekmat, Imam Khomeini township, Koi Sadeghiyeh, Bagher khan 1, Koi Azadegan, and Mohammadabad Korah, the school density was lower. However, at the level of statistical blocks, the density of the schools in the blocks located in the neighborhoods of Imam Reza, Koi Behdari, Ferdowsi, Mosalla, Nader, Amiriyeh, Berengy, Mirazkochekkhan, Manba Ab, Daneshsara, Park Shahr, Vosough, Enghelab, and Mokhofam was higher.
In the next step, the overall distribution of the schools at the spatial scales of neighborhoods, the statistical districts, and the statistical blocks was examined in different years. For the neighborhood spatial units, by comparing 2011, 2016, 2018, and 2021, it was found that the density of the schools was higher in the downtown neighborhoods, and these neighborhoods play the role of hot polygons. Neighborhoods in parts of the south of the city were also identified as cold polygons in terms of the school density.
At the level of statistical districts, the results indicate that the statistical districts in the downtown areas for the year 2011 are known as hot polygons in terms of the school density. The school density was also lower in the statistical districts of the western and the northeastern parts of the city. However, for the years 2016, 2018 and 2021, due to the addition of new schools in the Golestan Shahr township, the density of the schools in these areas has also increased, and the statistical district related to this township is also a hot polygon. Also, with the construction of the new schools on the outskirts of the city, the hot spots of the schools in the downtown area has been gradually reduced. However, the statistical districts within the neighborhoods of Bagh Motahari, Zir Bagh Motahhari, Hekmat townships, Imam Khomeini townships, Koi Sadeghiyeh, Koi Azadegan, Bagherkhan 1, and Mohammadabad Korah are cold polygons.
At the level of the statistical blocks, the results indicate that the statistical blocks located in parts of the city center in the year 2018 are known as hot polygons in terms of the school density. However, for years 2011, 2016, and 2021, some of the statistical blocks related to Golestan city are also among the high-density polygons. For the year 2021, there are cold polygons in parts of the south of the city, which include fewer schools than the other blocks. The results of these analyses are shown in Figure 4.

Implementation of Moran's I and Getis Ord General G
In this section, in order to investigate the spatial autocorrelation based on the substructure variable and to determine the random or clustered distribution of the schools, Moran's I and Getis Ord General G indices were used. At first, Moran's I statistic was used to check for the random or the clustered distributions of the girls' primary schools, first high schools, and secondary high schools. It can be concluded that the overall distributions of the primary schools, first high schools, and the secondary high schools for girls based on the substructure variable in the years 2011, 2016, and 2018 are random, which is due to the closeness of the value of Moran's coefficient to zero.
Then, in order to confirm the results of the Moran's I analysis, the analysis of the Getis Ord General G statistic (which is in regard to the substructure variable) on the point layers of the girls' primary schools, first high schools, and secondary high schools was performed for the years 2011, 2016 and 2018. According to Tables 5-7, the results show that the overall distributions of the girls' primary schools, first high schools, and secondary high schools based on the substructure variable for the years 2011, 2016, and 2018 are random, because the z-score is close to zero. The Moran's I and Getis Ord General G analyses were conducted on the point data of the boys' primary schools, first high schools, and secondary high schools. The results are shown in Tables 8-10. For the boys' primary schools, first high schools, and secondary high schools, because the coefficient of Moran is close to zero for the years 2011, 2016, and 2018, the overall distribution of the boys' primary schools based on the substructure variable is random.
The result of the Getis Ord General G analysis on the boys' primary schools confirms the Moran's I analysis. This means that the Z-score is close to zero for the years 2011, 2016, and 2018, and the overall distribution of the boys' primary schools based on the substructure variable is random. For the boys' first high schools in the years 2011 and 2018, the distribution of the schools is random, because the Z-Score is close to zero. However, for the year 2016, the distribution of schools is clustered with a low concentration or cold spots, because the Z-score is below −1.65. The distributions of the boys' secondary high schools for the years 2011, 2016, and 2018 are also random. To investigate the spatio-temporal variations of the distribution pattern of the schools based on the net per capita variable, which includes the substructure area and the number of students, the differential Moran method was used. In the differential Moran analysis, which was performed using Geoda 1.10 software, the weight matrix was defined based on the existence of the common border method (0 and 1) according to the Equation (12). Also, in order to achieve a higher accuracy, 999 iterations were used for the simulation.
In this regard, using the calculations of the differential Moran index, the trend of the net per capita changes for all primary schools, first high schools, and secondary high schools were evaluated in the 2016-2017, 2017-2018, and 2016-2018 time intervals. The spatial-temporal autocorrelation results of the schools with respect to the net per capita variable are shown in the differential Moran scatter plot in Figure 5, which illustrates the pattern of the net per capita changes of the schools in the desired time intervals.

Implementation of Anselin Local Moran's I and Getis Ord Gi*
After confirming that the data distribution is clustered, the local indices were used to determine the location of the clusters. By performing the general autocorrelation analyses, it was found that in most cases the distribution pattern of the schools based on the substructure variable is random. However, with respect to the values of the Moran, the Getis Ord indices, and the value of Z-Score, there is a weak spatial autocorrelation among the data. In this regard, to show this kind of autocorrelation, Anselin Local Moran's I and Getis Ord Gi* analyses were used in this study.
Anselin Local Moran's I analysis was first performed given the substructure variable for the girls' primary schools in 2011, 2016, and 2018. The results show that the substructure of the girls' primary schools in the southwest of the city in the year 2018 was above average (high-high areas), and it was below average in a region of the downtown area (low-low areas), which indicates the presence of clusters in these areas. However, in another part of downtown, the high outliers that have a value higher than the neighbors were also found. In the year 2016, the high-high areas are seen in the southwestern part of the city. Furthermore, there are high outliers in another part of the city center. However, there are no clusters in the city for the year 2011, but there are high outliers in part of downtown.
Next, considering the substructure variable, a hot spot analysis was performed by calculating the Getis Ord Gi* index for the girls' primary schools in 2011, 2016, and 2018. Then, in order to create a continuous raster surface, the Inverse Distance Weighting (IDW) interpolation method was applied to the results of the Getis Ord Gi* method. The results of applying this hybrid method are as follows. Next, the inverse distance weighting interpolation method was used on the results of the Getis Ord Gi* analysis for the boys' primary schools in 2011, 2016, and 2018. The results of applying this hybrid approach were that for the substructure variable, which was in a small area east of the city in the Taher Gholam neighborhood, the distribution of the schools is clustered with low values in 2011 and 2016. For the year 2018 in the southeast of the city in the Mirzakuchakkhan neighborhood, the distribution of the schools is clustered with high values, and in the Koi Sadeghiyeh neighborhood west of the city, the distribution of the schools is clustered with low values. The results of applying the inverse distance weighting interpolation on the results of the Getis Ord Gi* analysis for the girls' and the boys' primary schools are shown in Figure 6.
the inverse distance weighting interpolation on the results of the Getis Ord Gi* analysis for the girls' and the boys' primary schools are shown in Figure 6.

Discussion
Spatial distribution of the schools in the city plays a big role in educational equity. Therefore, for educational equity, the spatial distribution of the schools should be evaluated. In addition, one of the most important issues in the field of educational spatial justice is evaluating the results of different scenarios and the policies in this field at different times. In this regard, this study examined the spatial-temporal distribution of education facilities. To fulfill the purpose of this research, two types of analyses were used. The first type of analysis investigated the general pattern of the spatial-temporal distribution of the schools without specific parameters. For the first type of analysis, the Average Nearest Neighbor and optimized hot spot analysis methods were used. However, the second type of analysis was based on certain parameters, which included the substructure and the net per capita. Using the Moran's I and the Getis Ord General G analyses, the distribution pattern of the schools based on the substructure variable was studied. The differential Moran method was also used to investigate the spatiotemporal variations of the distribution pattern of the schools based on the net per capita variables in the 2016-2017, 2017-2018, and 2016-2018 time intervals.
The results of the Average Nearest Neighbor analysis indicate that the overall distribution of the primary, first, and the secondary high schools is clustered in the years 2011, 2016, 2018, and 2021. The Average Nearest Neighbor analysis was also performed separately for the girls' and the boys' schools. The result is that the type of overall distribution of the boys' primary, first, and secondary high schools are clustered in 2011, 2016, and 2018. Regarding the girls' schools, the distribution pattern of the primary and the secondary high schools is clustered for the years 2011, 2016, and 2018. A cluster pattern was also obtained for the girls' first high schools in the year 2011, but in the years 2016 and 2018, the distribution of the schools is random.
The optimized hot spot analysis was also used to show the overall distribution of the schools and the location of the clusters. Applying the method of counting point events within each cell indicated the accumulation of hot spots in the central areas of the city and cold spots on the outskirts of the city and the Golestan township. Using the method of counting point events within each polygon, the overall distributions of the schools in different levels, such as neighborhoods, statistical districts, and statistical blocks were evaluated and compared. This overview shows the high density of the schools in the central areas of the city, but the difference between the location of the hot and the cold polygons can be detected in different spatial units. The distribution patterns of the schools at the spatial scales of neighborhoods, statistical districts, and statistical blocks were studied for the 2011, 2016, 2018, and 2021 years. The density of the schools in the city center is higher for all three spatial units of neighborhoods, statistical districts, and statistical blocks. At the statistical district level, with the addition of new schools in the Golestan township, the density of the schools in this area has increased for the years 2016, 2018, and 2021. Also, with the construction of new schools on the outskirts of the city, the density of the schools in the downtown area has gradually decreased. At the level of statistical blocks for years 2011, 2016, and 2021, in addition to the city center, some of the statistical blocks related to the Golestan township are also high-density polygons. In general, it can be said that with the construction of new schools in the years 2018 and 2021, hot spots of the schools in the downtown area has been reduced, but the density of schools in other parts of the city has increased. However, the southern area and the suburbs of the city are known as cold polygons with lower density and require more organization and management.
In most of the Moran and Getis Ord analyses performed in this study, considering the substructure variable, the distribution pattern of the schools is random. This indicates that the schools with a higher substructure value did not focus on specific areas, and the substructure of the schools is affected by other factors, such as population density, number of classes, and a number of students that have not been addressed in this study. Based on the results of the local spatial autocorrelation indices for the girls' primary schools, the clusters are located in the downtown and southern areas of the city. For the boys' primary schools, the clusters are located in the east and southeast part of the city. For the secondary girls' secondary high schools, the clusters are located in the southern part of the city. For the secondary boys' secondary high schools, the clusters are located in the west and the southwestern parts of the city Also, the differential Moran method was used to evaluate the spatial-temporal autocorrelation of the net per capita of the schools. For all primary schools and first high schools, the trend of the net per capita changes was positive over The balanced spatial distribution of services is one of the most important signs of social justice. The distribution of urban services and the establishment of justice in the distribution of the urban services at the neighborhood level play a very important role in taking steps towards sustainable development. The proper distribution of services prevents intra-city migration, vacancy of some neighborhoods, uneven growth of the city, and the imposition of additional costs. The results of this research can also be used as an indicator to the policy makers for allocating funds to different cities and sectors. The results of this research can help the policymakers to identify the needs and the shortcomings of the city and also improve quality of life standards. Achieving social justice and creating equal opportunities is one of the most important needs of human societies today. Failure to pay attention to this important principle will create deep inequalities in society, and its realization will greatly contribute to political stability and national authority. In addition, the results of this study can be effective in reducing the cost of households due to transportation to schools and reducing urban traffic by equalizing the spatial schools distribution.
Determining the optimal distribution of the service centers is an issue that planners often dealt with, and the balanced distribution of the service centers, which includes the educational centers, is important because of its crucial role in culture and education. In addition to the correct method of education, the desired quality of education also requires suitable educational spaces, which have desirable and appropriate spaces that meet needs, but they require a careful and thorough analysis. This is because the desired quality of education in addition to the correct method of education also require suitable educational spaces that have desirable and appropriate spaces that meet needs and are carefully studied and analyzed. Among the most underexamined critical subjects in the field of educational spatial justice is the study of the spatial-temporal patterns of educational amenities, which is used to compare the effects of various scenarios and policies and develop a strong understanding of how the administrations' strategies are aimed at fostering social justice and equalizing the allocation of the educational facilities within society.

Conclusions
Based on the different results for the different point pattern analysis, it can be concluded that for comprehensive results, different algorithms should be implemented. In this regard, different algorithms have been implemented in this study. Based on the results of the Average Nearest Neighbor analysis, it can be concluded that spatial justice has been not established in the spatial distribution of the schools. Generally, based on the results of the optimized hot spot analysis within the cells and the polygons, it can be concluded that the hot spots are located in the central areas of the city, and the suburbs of Bojnourd city are considered cold polygons with less-density. However, when we analyze the results spatiotemporally, we can conclude that the density of the schools in the downtown areas has gradually decreased, and has increased in the suburbs. This conclusion indicates that the authorities had/have special policies for the suburbs, but that it still wasn't/isn't enough. In general, by examining the spatiotemporal results of the local spatial autocorrelation índices, which include Moran and Getis Ord, we conclude that over time, the amount of clusters has declined. This shows that the policymakers have taken steps towards realizing spatial justice in society. Based on the results of the differential Moran method, we conclude that the highest concentration of city planning in Bojnourd between 2011 and 2018 has been on the primary schools and the first high schools, but that the planners paid less attention to the secondary high schools during this period. This may also be the result of demographic structure changes in Bojnourd, because the families have tended to have fewer children in recent years. As a result, the number of students in the primary schools and the first high schools has diminished, and the net per capita in the primary schools and the first high schools has increased. Therefore, based on the results of different algorithms and depending on the circumstances, different policies could be adopted.
Based on the results of the Getis Ord Gi*, decision-makers can prioritize educational resources and educational funding for neighborhoods in the hotspot regions. By creating a continuous raster from the results of the Getis Ord Gi* method, we can have a comprehensive view of different areas of the city. We also recommend that the authorities remove the appropriate educational policies and initiatives in the coldspot areas and then implement them in other neighborhoods. Thus, educational facilities would be transferred from areas with lower Getis Ord Gi* index to areas with higher Getis Ord Gi* index. Also, new schools should be located in areas with a lower Getis Ord Gi* index.
Finally, we suggest that in order to achieve more accurate results, the areal pattern of the school distribution should also be examined, and the parameters, such as population density, accessibility, number of students, and number of classes should also be included in the analyses. Also, based on site selection criteria, new schools can be located in the cold spots.