Developing Statewide Optimal RWIS Density Guidelines Using Space-Time Semivariogram Models

Preventing weather-related crashes is a significant part of maintaining the safety and mobility of the travelling public during winter months. To help mitigate detrimental effects of winter road conditions, transportation authorities rely on real-time and near-future road weather and surface condition information disseminated by road weather information systems (RWIS) to make more timely and accurate winter road maintenance-related decisions. However, the significant costs of these systems motivate governments to develop a framework for determining a region-specific optimal RWIS density. Building on our previous study to facilitate regional network optimization, this study is aimed at considering the nature of spatiotemporally varying RWIS measurements and integrating larger case studies comprising eight different US states. Space-time semivariogram models were developed to quantify the representativeness of RWIS measurements and examine their effects on regional topography and weather severity for improved generalization. The optimal RWIS density for different topographic and weather severity regions was then determined via one of the most successful combinatorial optimization techniques—particle swarm optimization. The findings of this study revealed a strong dependency of optimal RWIS density on varying environmental characteristics of the region under investigation. It is anticipated that the RWIS density guidelines developed in this study will provide decision makers with a tool they need to help design a long-term RWIS implementation plan.


Introduction
Intelligent Transportation Systems (ITS) are an important part of modern transportation engineering and have a significant impact on society and our everyday life. Their main contribution is the improvement of transportation mobility and safety, particularly in how it relates to the prevention of weather-related road crashes, a vital and challenging issue, especially for countries with cold regions. In the entire US, over 1.5 million road crashes, 800,000 injuries, and 7,000 fatalities occur annually due to adverse weather [1]. In Canada, about 3,000 fatalities result from weather-related accidents and over one in 135 people experience driving-related injuries annually [2]. A road weather information system (RWIS) consists of a set of weather and road condition sensors installed at a station along the roadside. This combination of advanced sensors gathers, processes, and dissemi-nates the road weather and condition information that is then used by road maintenance authorities to make operative decisions before and/or during inclement weather events in order to enhance overall road safety and decrease weatherrelated crashes. It is also used by travellers via RWISconnected dynamic message signs to help them make more informed decisions when traveling during inclement weather events. But investment cost in these technologies is significant. Indeed, many North American transportation agencies have invested millions of dollars to deploy RWIS stations to improve the monitoring coverage of winter road surface conditions [3,4]. As transportation authorities are constantly challenged by limited budgets for installing additional stations, it is of paramount importance to address the questions of how many stations are necessary (optimal RWIS density) and where to locate them (optimal RWIS location) to maximize the return on their investment.
There are numerous studies conducted by researchers around the world to quantify the spatial coverage of RWIS ata and determine the optimal RWIS density and location based on available RWIS data [4][5][6][7][8]. The Federal Highway Administration (FHWA) initiated extensive effort to provide RWIS sitting guidelines based on the knowledge and experience of field operators [8]. A more recent study by  evaluated the dependency of optimal RWIS density on topographic conditions of the regions under investigation by examining three US states (Iowa, Utah, and Minnesota) and one Canadian province (Ontario). Their findings indicated that more RWIS stations would be needed in mountainous areas than in flatland areas. The study also asserted that a region with a longer spatial autocorrelation range would require fewer stations than a region with a shorter range [4]. Our recent effort [9] further extends the former work by investigating the dependency of optimal RWIS density on two different measures, namely, topography and weather severity, in an effort to improve generalization potentials and design a long-term strategic RWIS deployment plan. Despite the uniqueness of the method developed and analyses undertaken therein, the study dealt solely with the spatial domain, which does not account for the inherent temporal correlation of road weather and surface conditions. Since, not only would road surface conditions vary over space but also over time, there is therefore a need to incorporate the temporal domain in the optimization of RWIS location and density.
Spatiotemporal data analysis is a well-known geostatistical analysis technique that is commonly used to analyse the variability of parameters that have a tendency to fluctuate over space and time. It combines spatial and time series analysis concurrently to preserve the interactive effect of time variation on spatial domain and vice-versa. Previous research reveals that spatiotemporal analysis is more accurate than spatial analysis as both time and space are included in spatiotemporal analysis [10]. The underlying techniques have been previously used to model air pollutants by quantifying the space-time variability of certain particles' concentrations [10][11][12][13]. Spatiotemporal regression kriging has been used to predict precipitation using Moderate-Resolution Imaging Spectroradiometer (MODIS) and Normalized Difference Vegetation Index (NDVI) data [14]. In a study more relevant to the topic of interest, Wang et al. (2019) investigated the spatiotemporal variation of road weather and surface conditions using RWIS data from Alberta, Canada. In this study, the authors developed space-time semivariogram models using RWIS measurements from Alberta in order to examine the applicability of the method. The study's findings provided the spatiotemporal feature of the RWIS database [15]. However, the dependency of the spatiotemporal feature of RWIS measurements on topography and weather severity has yet to be scrutinized.
As discussed, there exists a large gap in the current literature for determining spatiotemporally optimal density and location of an RWIS network. Equally important, there are no existing guidelines available for winter road maintenance agencies to use and decide its optimality, especially for regions with no RWIS data-the information that has been deemed prerequisite for undertaking any RWIS-related analyses including location and density. Since the study of finding an optimal location using space-time features is worthy of a separate paper, we will stay focused on answering another important question: how many RWIS are required in a region with varying environments to provide sufficient coverage over space and time? In particular, we are interested in (1) investigating spatiotemporal autocorrelation of RWIS measurements (i.e., road surface temperature (RST)) using large scale case studies, (2) examining the effect of topography and weather severity of regions on spatial and temporal continuity of RST data, and (3) developing an RWIS optimal density chart that can facilitate the decision-making process for planning a long-term RWIS deployment strategy. The outcome of this study can be used as a guideline for an RWIS network expansion plan that applies to different topographic and weather severity zones without having to gather any additional datasets.
The rest of the paper is organized as follows. The next section briefly discusses the characterizations of topography and weather based on our previous study, followed by the proposed methodology, including a brief discussion on spatiotemporal semivariogram modelling, covariance models, and optimization process. The case study section includes the data description and processing details. The next section includes spatiotemporal analysis results and density optimization results. The last section highlights the findings and suggests recommendations for future work.
1.1. Characterizations of Topography and Weather. The study area was characterized based on topography and weather severity taken from our recent study [9]. These two parameters were selected as zonal classification standards. Topography-based zonal classification was conducted using the topographic position index (TPI), which defines a relative topographical variation of an area of interest and its neighbourhood. On the other hand, weather severity-based zonal classification was conducted using the winter severity index (WSI) that provides a numerical estimate of weather severity using annual average snowfall accumulation and duration, annual duration of blowing snow, and annual duration of freezing rain [16][17][18]. The TPI and WSI maps for eight states in the United States (Colorado (CO), Iowa (IA), Kansas (KS), Minnesota (MN), Ohio (OH), Kentucky (KY), Nebraska (NE), and Wyoming (WY)) are presented in Figure 1. Flatland area is represented by TPI 1, followed by hilly and mountainous regions as TPI 2 and TPI 3, respectively. Similarly, WSI 1 represents less severe weather regions, followed by moderate, high, and extremely high severe weather zones as WSI 2, WSI 3, and WSI 4, respectively. These classification standards will be used in the following sections to further enhance generalization potentials pertaining to determining optimal RWIS densities for any given region under investigation.

Methodology
Spatiotemporal semivariogram modelling was used in this study to evaluate the spatiotemporal variability of RWIS measurements. RWIS data for a winter season was downloaded and processed, and a space-time matrix was formulated as an input of the spatiotemporal analysis. The dataset 2 Journal of Sensors was classified based on previously developed TPI and WSI classes, and a separate analysis was conducted for each month and zone, which were aggregated to generate a seasonal spatiotemporal autocorrelation range for the TPI and WSI classes. Optimal RWIS densities were then determined by relying on spatiotemporal semivariogram analysis results. Density per unit area was calculated and compared for different topographic and weather-based zones. The research methodology of this study can be summarized as follows: (1) Development of spatiotemporal semivariogram models to examine spatial and temporal autocorrelation of RWIS data (2) Evaluation of the effective spatial and temporal range of continuity under different topographic and weather settings (3) Examination of the hypothesis that the spatiotemporal variability of road weather data is dependent on the topographic variation and weather severity of the region (4) Determination of optimal density of RWIS stations using modified particle swarm optimization (PSO) for different topographic and weather classes (5) Comparison of the optimal RWIS density per unit area for different classes Figure 2 presents an overview of the proposed methodology.  Figure 1: (a) TPI and (b) WSI maps of eight states in the United States (adopted from [9]).

Spatiotemporal Semivariogram
3 Journal of Sensors a spatial semivariogram is a plot of semivariance versus lag separation distances between pairs of measurement points. Semivariance is a measure of dissimilarity between two measurements as a function of separation distance. A semivariogram can be defined using three basic parameters: range, nugget, and sill. The value at origin defines the dissimilarity with zero separation distance, which is theoretically zero. However, this value could be nonzero due to some measurement or sampling error; this error is known as the nugget effect. Sill can be defined as the semivariance value at which the semivariogram levels off. Partial sill is often encountered during semivariogram analysis, which is the difference between the actual sill and nugget effect. The spatial range of autocorrelation is the distance at which the semivariogram reaches the sill value. Spatial autocorrelation is considered zero beyond this spatial range [19][20][21][22]. A typical semivariogram with its parameters is presented in Figure 3(a).
The traditional spatial analysis is incorporated with temporal analysis to consider both the spatial and temporal effects. Spatiotemporal analysis is conducted for variables that vary over space and time. A set of variables in a spatiotemporal field can be defined as z = fzðs, tÞ | s ϵ S, t ϵ Tg. Here, S is the spatial domain and T is the temporal domain. Thus, the general equation of a random field Z is −z i = Zðs, tÞ, i = 1, 2, 3, ⋯n × T. Here, n is the number of stations and T is the number of time points. The random fields Zðs, tÞ can be modelled as Zðs, tÞ = μðs, tÞ + εðs, tÞ. Here, μðs, tÞ is the deterministic part and εðs, tÞ is the stochastic part. The deterministic part refers to the spatiotemporal trend while the stochastic part is the zero-mean second-order stationary spatiotemporal random field [13]. The stochastic part is used for spatiotemporal semivariogram modelling. Spatial and temporal variance is estimated as half of the mean squared difference between pairs of data separated by a user-defined spatial (h s ) and temporal lag (h t ) as Here, γðh s , h t Þ is the estimated semivariance value, nðhs, htÞ is the total number of pairs in the analysis domain, and z ðsk, tkÞ is the measurement at spatial location s k and temporal location t k [23,24]. A three-dimensional spatiotemporal semivariogram is presented in Figure 3 As the estimated model is irregular, a mathematical model was used for smoothening the empirical variogram. The covariance models generally used for spatiotemporal variogram modelling are discussed in the following section.

Covariance Models.
A number of covariance models are used for spatiotemporal semivariogram modelling. The most popular and widely used models are described as follows [10,[25][26][27].

Separable Covariance Model.
A separable covariance model assumes that the spatiotemporal covariance function can be represented as the product of a spatial and temporal term. The covariance function can be written as C sep ðh, uÞ = C s ðhÞC t ðuÞ. Thus, the equation of the variogram is γ sep ðh, uÞ = sill · ðγ s ðhÞ + γ t ðuÞ − γ s ðhÞγ t ðuÞÞ. Spatial and temporal sill is ignored in this model and kept constant at 1. A joint sill (=1) is used, which combines both spatial and temporal effects.

Product-Sum
Covariance Model. This model assumes a new parameter, k, as a weighting factor of the product ðk > 0Þ. The equation for the covariance function is C ps ðh, uÞ = k · C s ðhÞC t ðuÞ + C s ðhÞ + C t ðuÞ. The equation for the variogram can be written asγ ps ðh, uÞ = ðk · sill t + 1Þ γ s ðhÞ + ðk · sill s + 1Þγ t ðuÞ − kγ s ðhÞγ t ðuÞ. The expression of a joint sill is sill st = k · sill s · sill t + sill s + sill t . Here, the spatial and temporal nugget is ignored and kept constant at 0; the joint nugget is used to account for both spatial and temporal effects.

Metric Covariance Model.
Identical spatial and temporal covariance functions are assumed in this model, except for spatiotemporal anisotropy. Spatial, temporal, and spatiotemporal distances are treated equally in a joint covariance model by matching space and time by a spatiotemporal anisotropy parameter, kðStAniÞ. The equation for the covari-   As several previous studies attested to the superior performance of the sum-metric model in fitting the spatiotemporal variogram using environmental parameters (i.e., smallest mean squared errors), this model was selected and used in this analysis [12,14].

Density Optimization.
Density optimization in this study was conducted in order to compare the number of RWIS stations needed per unit area in different TPI-and WSIbased classes. Semivariogram parameters generated from different classes were used in the optimization process as representative of the spatial characteristics of the zones. A randomly selected square region (area of 10,000 square km) within the study area was used as the experimental boundary wherein the solution was limited for density optimization of all classes. Following the same procedure for generating the optimal location solutions, RWIS density curves for different TPI and WSI zones were generated based on a predefined number of RWIS stations. The marginal increment of benefit associated with each additional RWIS stations was calculated to determine the optimal RWIS density.
Density optimization was conducted using particle swarm optimization (PSO). PSO is an evolutionary computation technique and population-based global optimization method developed by Kennedy and Eberhart in 1995 [28]. PSO is widely used and popular in scientific computations because it is easily implemented and computationally inexpensive. In this optimization process, a number of n-dimensional candidate points (particle) are placed in the search space of a function, and each of the particles evaluates the objective function at its current location. Each particle can be assumed as a potential solution presented by velocity and position [29,30]. Movement of each particle is determined based on its best fit location with one or more swarms, and the algorithm searches for optima by updating the generations [31,32]. The i th particle in the search space can be represented as x i = ðx i1 , x i2 , ⋯, x in Þ. Each particle in the swarm flies to the previous best position and global best position; those are named as "pbest" and "gbest," respectively. The best previous position of the i th particle can be presented as p i = ðp i1 , p i2 , ⋯, p in Þ. The index of the best particle in the swarm is represented by the subscript g. The velocity of particle movement is represented by v i = ðv i1 , v i2 , ⋯, v in Þ. The particle is attracted by pbest and gbest during the search process according to where d is the dimension, representing the total number of candidate RWIS sites, where 1 ≤ d ≤ n; c 1 and c 2 are positive constants; ζ and η are random adjustment factors with a range of 0 to 1; and ω is the inertia weight. The performance of each particle is measured using a predefined fitness function. A Binary Particle Swarm Optimization (BPSO) was proposed by Kennedy and Eberhart in 1997 for solving  Journal of Sensors integer-programming problems, as the original PSO is not suitable in this case [33]. The basic difference between these two methods is in updating of the particles' position. The sigmoid function is utilized in BPSO where every dimension in the position becomes a number between 0 and 1. A modified BPSO was proposed by Gu et al. in 2019 to solve the RWIS location optimization problem [29]. The similar method was adapted for use in our study for region-wise RWIS density optimization. The position of the particle is updated using In the modified BPSO, a threshold probability, r, is set to control whether the x id becomes 1 or not, where 1 represents the selection of the element. During optimization, the total number of RWIS stations (m) is set as a constant and the algorithm is set to select best-fit "m" number of locations in the search space. The original BPSO shows premature convergence because of a quick loss of diversity. To treat this problem, more randomness is added into the internal mechanism of the modified BPSO in order to expand the search space and allow the particle to escape from any possible local minima. Another addition is that, if more than one location has the same probability for an RWIS station, a mechanism is set in the modified BPSO to randomly select one of them as the solution. Lastly, in addition to the maximum velocity set to control the speed of convergence, the inertia weight (ω) and self-learning factor (c 1 ) are set to decrease from 0.9 to 0.4 and 2 to 0, respectively, in the search process. On the other hand, the society-learning factor (c 2 ) is set to increase from 0 to 2. The parameters are chosen to ensure that the particles can fly slowly while eliminating their ability for selflearning and enhancing social-learning. The steps associated with the modified BPSO algorithm are listed as follows [29][30][31]: Step 1. "m" particles are initialized with dimensions of velocity.
Step 3. Two top probabilities are selected in each particle's position, and they are set to the selected candidate points for locating RWIS stations, then ordinary kriging variance is calculated for all the unknown points as the fitness value.
Step 4. Memorize the current individual best positions and the global best positions.
Step 5. Update ω, c 1 , and c 2 using Step 6. Update each particle's velocity using If the v id > v max , then v id = v max .
Step 8. Update the individual best position and global best position by comparing fitness values. If the updated fitness value is smaller than before, accept the new solution.
If not, then repeat the process from Step 2.

Weather Severity
Data. An ESRI shapefile containing the weather severity index used in this study was generated by Meridian Environmental Technology.  3.3. Data Processing. RWIS data were processed to remove the missing and erroneous data using five steps: (a) data completeness test to identify missing data, (b) reasonable range test to find erroneous data, (c) cross-checking RST data with air-temperature data, (d) RST data pattern analysis, and (e) detrending RST data with respect to time using a Generalized Additive Model (GAM). Data completeness was checked by identifying the total missing data for each sensor. If the total missing data is more than fifteen percent, the associated sensor ID was marked and the data from that sensor was not used for analysis. A reasonable range was tested based on historical data ranges for the associated region and month. Filtered RST data were then cross-checked with air-temperature data ranges for any possible outliers. An RST data pattern analysis was performed by plotting the day of the month versus the average daily temperature for all selected sensors, for each state and each month.
All selected sensors were expected to show a similar pattern throughout the month. If any unusual pattern was noticed, the RST data for the associated sensors were further investigated for the time period of the unusual pattern. In total, 48 sets of data (six months for eight states) were analysed using the above described process. Finally, RST data was detrended with respect to time using a Generalized Additive Model (GAM) to incorporate shorter scale variation in the temporal domain, where GAM worked as a generalized linear model with linear predictors. The GAM function was formulated as Here, m is the variable   [15,34]. Descriptive statistics (means and standard deviations) revealed relatively less variation in average monthly temperatures in the midwinter months than in the shoulder months. Overall minimum and maximum temperatures for the study area were -30.3°C to 51.5°C (-22.5°F to 124.7°F). Table 1 presents the maximum, minimum, average, and standard deviations of RST for the study area. Because of erroneous data in November 2016 in Kentucky, it was excluded from the analysis. Figure 5 shows the seasonal maximum, minimum, average, and standard deviations of RST for the study area.
Spatiotemporal semivariogram modelling was performed using the R statistical package version 3.2.5 [10,35,36]. RWIS density optimization was coded in Python. To improve the computational efficiency, all optimizations undertaken in this study were run on the supercomputer "beluga" from the University of Alberta, managed by Calcul Québec and Compute Canada (https://www.computecanada .ca/), with 32 CPUs, each of which runs on a 2.4 GHz CPU and 1 GB memory.    Journal of Sensors of RST for each study zone was then analysed using spatiotemporal variogram modelling methods using the gstat package in R [35,36]. A sample of spatiotemporal semivariograms of different TPI and WSI classes is presented in Figure 6. According to Figure 6, higher spatial and temporal variation is noticed with increase of topographic and weather variations. High cyclicity is noticed for WSI 1, possibly due to insufficient datasets used in spatiotemporal modelling. Seasonal spatiotemporal analysis results for TPI and WSI zones are presented in Figure 7. According to Figure 7, relatively higher spatial and temporal ranges are obtained for TPI class 1 (flatland area), followed by classes 2 and 3 representing hilly and mountainous areas, respectively. Similar results were obtained for weather-based classes, where regions with less topographic variation and less severe weather have a higher spatial and temporal range as depicted in Figure 8. The range of autocorrelation decreases with the increase of topographic variation and weather severity. In general, there is a trend of higher autocorrelation range during midwinter months compared to shoulder months, which is as expected.

Spatiotemporal
The effect of weather severity in the TPI 1 zone, which is the flatland area, is presented in Figure 9. In our study area, the flatland region consists of three weather severity regions and there is a trend for a higher autocorrelation range in areas with less weather severity than areas with more severe weather severity regions, especially for the temporal range.
The effect of weather for the spatial range is negligible as the total difference is only 1.6 km. From this, we can conclude that topography can serve as a more intuitive measure for RWIS network planning than weather severity. Similar comparisons for other TPI zones have not been made because the TPI 2 zone consists of WSI 2 and 3 zones, and the TPI 3 zone is identical to the WSI 4 zone (see Figure 1).

Development of Optimal RWIS Density Guidelines.
Density optimization is conducted using all three semivariogram parameters (range, nugget, and sill), while the dependency of the semivariogram range upon topographic-weather characteristics of the region is our main concern in this analysis. Spatial parameters of a spatiotemporal semivariogram were used as input for density optimization for topographic and weather severity zones (3 TPI classes and 4 WSI classes). A hypothetical network of 100 km × 100 km was used for density optimization. A 5 km × 5 km prediction grid was generated in ArcGIS [37] to create the candidate sites for RWIS station placement. The objective function was formulated to minimize the mean ordinary kriging estimation variance and solved via PSO as described earlier. The algorithm will optimize the location and density of RWIS stations in an iterative process where stations are added one by one into the study area and locations are selected based on heuristic attempts to minimize the objective function. The total number of RWIS stations allowed is arbitrarily limited to 100 in this study to ensure the variation trend of the estimation error as density changes can be fully displayed. The number of iterations is set to 5,000, after which the search process for new RWIS station locations is set to stop. Prediction errors were normalized to make a valid and fair comparison among zonal classes.
According to the density optimization results as shown in Figure 10, topographic and weather severity classes with larger spatial ranges required a lower number of RWIS stations, except weather zone 3, which shows a similar trend to weather zone 1. Such findings are in agreement with analyses using TPI classes: WSI 3 zone consists mainly of northern parts of MN, classified as a flatland area (TPI 1) and found to have a lower density of RWIS based on topographic analyses. From this, it can be stated that topographic measures (TPI) provide a more intuitive and direct   Journal of Sensors relationship with the impact than WSI that the spatiotemporal range has on optimal density. The number of RWIS stations needed for a 0.1-unit increment of benefit (number of stations needed to reach normalized fitness value of 0.1) is calculated from Figure 10 and presented in Figure 11. According to Figure 11, the initial 0.1 unit of incremental benefit requires the highest number of RWIS stations with each subsequent incremental benefit requiring fewer and fewer additional stations. This is as expected because the amount of marginal benefit from additional RWIS stations should decrease as station numbers increase. A similar trend is observed for both TPI and WSI classes.
For determining the optimal RWIS density, marginal benefits (number of additional RWIS stations needed for a 0.1-unit incremental benefit) are calculated from Figure 10 and presented in Figure 12. According to Figure 12, the number of added RWIS stations for an initial marginal incremental benefit of 0.1 unit is the highest, and then the number   Figures 11 and 12 , it can be noted that the marginal benefit decreases significantly after the 0.3-unit increment of benefit mark. This is determined as the median of the range and number of associated RWIS stations for a 0.3unit increment of benefit increase which is defined as the optimal RWIS density. The number of RWIS stations needed for an incremental benefit of 0.2 units and 0.4 units is selected as the upper and lower bounds, respectively. The number of RWIS stations needed for 0.2-, 0.3-, and 0.4-unit benefit increments is recorded from Figure 11 and plotted in Figure 13. According to this figure, RWIS density is the lowest for TPI class 1 and increases with an increase in topographic variation. Similarly, RWIS station numbers for less weather severe regions are the lowest and increase with an increase in weather severity, except WSI 3 for the above noted reasons.
The density optimization results are used to generate an RWIS density chart for TPI-WSI zones which is presented in Table 2. RWIS density for each TPI-WSI region is calculated by combining the topographic and weather effects with equal weightage. In our previous analysis, only the spatial domain was used for density optimization, and the output showed a minimum of one to a maximum of two RWIS stations per unit area [9]. Here, in this study, both spatial and temporal domains were used for density optimizations. As a result, three to seven RWIS stations are required for a unit area of 10,000 km 2 , whereby attesting that the number of RWIS required to provide enough coverage over space and time is about three times more than that over space only. Such find-ings can readily be used by winter road maintenance agencies for planning a regionwide RWIS network, especially for regions with limited or no available RWIS stations.

Conclusion and Recommendations
Winter road maintenance is one of the most critical activities for transportation maintenance agencies, especially for cold region countries. The significant and critical information needed for making winter road maintenance decisions is related to road condition and weather data, which is often collected, processed, and transmitted by a road weather information system (RWIS). Effective and efficient planning of RWIS networks is a must and is needed for maximizing the monitoring coverage and benefit of RWIS. The effectiveness of an RWIS network depends on its density and spatial distribution. In this study, we investigated the representativeness of RWIS measurements in two analysis domains-space and time. Spatial and temporal continuities of the variable of interest (RST) were investigated using geostatistical spatiotemporal semivariogram analysis and compared to different previously established topographic and weather regions. Lastly, optimal RWIS densities for three TPI and four WSI zones were estimated from the density optimization output.
The key findings of this study are listed as follows:   (iii) Weather severity was determined as another factor influencing the spatiotemporal continuity range of RWIS data. Areas with less severe weather tended to have a higher spatial range (e.g., WSI class 1, whereas areas with more severe weather had a lower range in spatial autocorrelation) (iv) The desired RWIS density showed a strong dependency on topography and weather characteristics of the region. Higher RWIS density is required for regions with high topographic variation and high incidence of severe weather, while lower RWIS density is needed for less varied topographic regions with less incidence of severe weather to achieve similar levels of monitoring coverage (v) Findings presented in this research provide an important basis for strategically locating regional RWIS stations, which are optimal in collecting measurements over space and time Recommendations for further research in this direction are listed as follows: (i) This study covers a wide region of flatland area, but the hilly and mountainous area covers a very small part of the total study area due to data availability issues. Hence, more case studies covering wider regions should be conducted for a better understanding of the relationship between the spatial range of autocorrelation in RST on the topographic-weather features to develop a more robust quantitative relation between these parameters (ii) The study period of this project is limited to one winter season, including six months from October 2016 to March 2017. Thus, large temporal ranges could be considered to improve the level of confidence in the outcomes (iii) Universal kriging or kriging with external drift could be applied considering meteorological parameters (wind speed and direction, precipitation, humidity, cloud cover, vegetation cover, etc.) to better capture the dependency of RST data (or other key parameters including road surface condition index) on local meteorological parameters. Additionally, universal kriging with height as an explanatory variable could be applied to determine the effect terrain has on the optimal RWIS density (iv) A sensitivity analysis could be conducted to investigate how the resulting optimal densities would change with respect to some of the factors considered in the analysis, especially those coefficients used to generate a winter severity index (or even a winter severity index model) and TPI classification schemes (v) Since many transportation agencies have adopted the AVL technologies and/or mobile RWIS sensors to improve the road surface condition monitoring capabilities, it will be worthwhile looking into synthesizing the best practices for integrating some of the latest mobile technologies to develop RWIS planning guidelines

Conflicts of Interest
The authors declare that there is no conflict of interest regarding the publication of this paper.