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

Authors: Simita Biswas; Tae J. Kwon
Journal: Journal of Sensors
 (2020)

---

## Abstract

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.


---

## Body

## 1. 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 disseminates 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 weather-related crashes. It is also used by travellers via RWIS-connected 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–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 Kwon and Fu (2017) 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–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–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.Figure 1
(a) TPI and (b) WSI maps of eight states in the United States (adopted from [9]).
(a)(b)


## 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–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.Figure 1
(a) TPI and (b) WSI maps of eight states in the United States (adopted from [9]).
(a)(b)


## 2. 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 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 classesFigure2 presents an overview of the proposed methodology.Figure 2
An overview of the proposed methodology.


### 2.1. Spatiotemporal Semivariogram Modelling

Geostatistical spatiotemporal semivariogram modelling was used in this study for spatial and temporal continuity analysis. In general, 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–22]. A typical semivariogram with its parameters is presented in Figure 3(a).Figure 3
A typical (a) semivariogram and (b) spatiotemporal semivariogram.
(a)(b)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 asz=zs,t∣sϵS,tϵT. Here, S is the spatial domain and T is the temporal domain. Thus, the general equation of a random field Z is −zi=Zs,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 Zs,t can be modelled as Zs,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 (hs) and temporal lag (ht) as
(1)γhs,ht=12nhs,ht∑k=1nhs,htzsk,tk−zsk+hs,tk+ht2.Here,γhs,ht is the estimated semivariance value, nhs,ht is the total number of pairs in the analysis domain, and zsk,tk is the measurement at spatial location sk and temporal location tk [23, 24]. A three-dimensional spatiotemporal semivariogram is presented in Figure 3(b).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.


### 2.2. 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–27].


#### 2.2.1. 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 asCseph,u=CshCtu. Thus, the equation of the variogram is γseph,u=sill·γsh+γtu−γshγtu. 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.


#### 2.2.2. 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 Cpsh,u=k·CshCtu+Csh+Ctu. The equation for the variogram can be written asγpsh,u=k·sillt+1γsh+k·sills+1γtu−kγshγtu. The expression of a joint sill is sillst=k·sills·sillt+sills+sillt. 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.


#### 2.2.3. 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,kStAni. The equation for the covariance function is Cmh,u=Cjointh2+k·u2. The equation for the metric variogram can be written as γmh,u=γjointh2+k·u2. Temporal distances are internally rescaled to an equivalent spatial distance to determine the equivalent factor in terms of dependence of 1 m separation in a second or a minute. The expression of spatiotemporal anisotropy is kStAni=spatialunit/temporalunit=m/sec/min.


#### 2.2.4. Sum-Metric Covariance Model

This model is a combination of spatial, temporal, and metric models. The equation for a covariance function isCsmh,u=Csh+Ctu+Cjointh2+k·u2. The equation for a sum-metric variogram can be written as γsmh,u=γsh+γtu+γjointh2+k·u2. Spatial, temporal, and joint nuggets are estimated separately in this model.


#### 2.2.5. Simple Sum-Metric Covariance Model

This model is the simplified version of the sum-metric model to restrict the spatial, temporal, and joint variograms to nugget-free models. A single spatiotemporal nugget is introduced in this model. The equation for a variogram isγssmh,u=nug·1h>0,u>0+γsh+γtu+γjointh2+k·u2. Here, spatial, temporal, and joint nuggets are set to 0; only a joint nugget is fitted.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].


### 2.3. 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 WSI-based 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 ith particle in the search space can be represented as xi=xi1,xi2,⋯,xin. 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 ith particle can be presented as pi=pi1,pi2,⋯,pin. The index of the best particle in the swarm is represented by the subscript g. The velocity of particle movement is represented by vi=vi1,vi2,⋯,vin. The particle is attracted by pbest and gbest during the search process according to
(2)vid=ωvid+c1ζpid−xid+c2ηpgd−xid,(3)xid=xid+vid,where d is the dimension, representing the total number of candidate RWIS sites, where 1≤d≤n; c1 and c2 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 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
(4)xid=11+e−vid.In the modified BPSO, a threshold probability,r, is set to control whether the xid 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 (c1) 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 (c2) 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 self-learning and enhancing social-learning. The steps associated with the modified BPSO algorithm are listed as follows [29–31]:Step 1. “m” particles are initialized with dimensions of velocity.Step 2. Velocities are converted to positions (probabilities) using the sigmoid function (Equation (4)).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ω, c1, and c2 using
(5)ωnew=ωold−ωmax−ωminnumberofiterations,(6)c1new=c1old−c1max−c1minnumberofiterations,(7)c2new=c2old−c2max−c2minnumberofiterations,Step 6. Update each particle’s velocity using(8)vid=ωvid+c1∗ζid∗Δxpid+c2∗ηid∗Δxgid,If thevid>vmax, then vid=vmax.Step 7. Update particles’ positions using Equation (4).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.


## 2.1. Spatiotemporal Semivariogram Modelling

Geostatistical spatiotemporal semivariogram modelling was used in this study for spatial and temporal continuity analysis. In general, 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–22]. A typical semivariogram with its parameters is presented in Figure 3(a).Figure 3
A typical (a) semivariogram and (b) spatiotemporal semivariogram.
(a)(b)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 asz=zs,t∣sϵS,tϵT. Here, S is the spatial domain and T is the temporal domain. Thus, the general equation of a random field Z is −zi=Zs,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 Zs,t can be modelled as Zs,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 (hs) and temporal lag (ht) as
(1)γhs,ht=12nhs,ht∑k=1nhs,htzsk,tk−zsk+hs,tk+ht2.Here,γhs,ht is the estimated semivariance value, nhs,ht is the total number of pairs in the analysis domain, and zsk,tk is the measurement at spatial location sk and temporal location tk [23, 24]. A three-dimensional spatiotemporal semivariogram is presented in Figure 3(b).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.


## 2.2. 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–27].


### 2.2.1. 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 asCseph,u=CshCtu. Thus, the equation of the variogram is γseph,u=sill·γsh+γtu−γshγtu. 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.


### 2.2.2. 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 Cpsh,u=k·CshCtu+Csh+Ctu. The equation for the variogram can be written asγpsh,u=k·sillt+1γsh+k·sills+1γtu−kγshγtu. The expression of a joint sill is sillst=k·sills·sillt+sills+sillt. 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.


### 2.2.3. 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,kStAni. The equation for the covariance function is Cmh,u=Cjointh2+k·u2. The equation for the metric variogram can be written as γmh,u=γjointh2+k·u2. Temporal distances are internally rescaled to an equivalent spatial distance to determine the equivalent factor in terms of dependence of 1 m separation in a second or a minute. The expression of spatiotemporal anisotropy is kStAni=spatialunit/temporalunit=m/sec/min.


### 2.2.4. Sum-Metric Covariance Model

This model is a combination of spatial, temporal, and metric models. The equation for a covariance function isCsmh,u=Csh+Ctu+Cjointh2+k·u2. The equation for a sum-metric variogram can be written as γsmh,u=γsh+γtu+γjointh2+k·u2. Spatial, temporal, and joint nuggets are estimated separately in this model.


### 2.2.5. Simple Sum-Metric Covariance Model

This model is the simplified version of the sum-metric model to restrict the spatial, temporal, and joint variograms to nugget-free models. A single spatiotemporal nugget is introduced in this model. The equation for a variogram isγssmh,u=nug·1h>0,u>0+γsh+γtu+γjointh2+k·u2. Here, spatial, temporal, and joint nuggets are set to 0; only a joint nugget is fitted.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].


## 2.2.1. 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 asCseph,u=CshCtu. Thus, the equation of the variogram is γseph,u=sill·γsh+γtu−γshγtu. 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.


## 2.2.2. 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 Cpsh,u=k·CshCtu+Csh+Ctu. The equation for the variogram can be written asγpsh,u=k·sillt+1γsh+k·sills+1γtu−kγshγtu. The expression of a joint sill is sillst=k·sills·sillt+sills+sillt. 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.


## 2.2.3. 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,kStAni. The equation for the covariance function is Cmh,u=Cjointh2+k·u2. The equation for the metric variogram can be written as γmh,u=γjointh2+k·u2. Temporal distances are internally rescaled to an equivalent spatial distance to determine the equivalent factor in terms of dependence of 1 m separation in a second or a minute. The expression of spatiotemporal anisotropy is kStAni=spatialunit/temporalunit=m/sec/min.


## 2.2.4. Sum-Metric Covariance Model

This model is a combination of spatial, temporal, and metric models. The equation for a covariance function isCsmh,u=Csh+Ctu+Cjointh2+k·u2. The equation for a sum-metric variogram can be written as γsmh,u=γsh+γtu+γjointh2+k·u2. Spatial, temporal, and joint nuggets are estimated separately in this model.


## 2.2.5. Simple Sum-Metric Covariance Model

This model is the simplified version of the sum-metric model to restrict the spatial, temporal, and joint variograms to nugget-free models. A single spatiotemporal nugget is introduced in this model. The equation for a variogram isγssmh,u=nug·1h>0,u>0+γsh+γtu+γjointh2+k·u2. Here, spatial, temporal, and joint nuggets are set to 0; only a joint nugget is fitted.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].


## 2.3. 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 WSI-based 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 ith particle in the search space can be represented as xi=xi1,xi2,⋯,xin. 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 ith particle can be presented as pi=pi1,pi2,⋯,pin. The index of the best particle in the swarm is represented by the subscript g. The velocity of particle movement is represented by vi=vi1,vi2,⋯,vin. The particle is attracted by pbest and gbest during the search process according to
(2)vid=ωvid+c1ζpid−xid+c2ηpgd−xid,(3)xid=xid+vid,where d is the dimension, representing the total number of candidate RWIS sites, where 1≤d≤n; c1 and c2 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 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
(4)xid=11+e−vid.In the modified BPSO, a threshold probability,r, is set to control whether the xid 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 (c1) 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 (c2) 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 self-learning and enhancing social-learning. The steps associated with the modified BPSO algorithm are listed as follows [29–31]:Step 1. “m” particles are initialized with dimensions of velocity.Step 2. Velocities are converted to positions (probabilities) using the sigmoid function (Equation (4)).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ω, c1, and c2 using
(5)ωnew=ωold−ωmax−ωminnumberofiterations,(6)c1new=c1old−c1max−c1minnumberofiterations,(7)c2new=c2old−c2max−c2minnumberofiterations,Step 6. Update each particle’s velocity using(8)vid=ωvid+c1∗ζid∗Δxpid+c2∗ηid∗Δxgid,If thevid>vmax, then vid=vmax.Step 7. Update particles’ positions using Equation (4).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.


## 3. Case Studies

### 3.1. Study Area and RWIS Network

The study area was selected based on RWIS data availability and distribution of RWIS stations to cover varying topography and weather conditions that can be challenging for driving. Eight states had sufficient data, RWIS station coverage, and number of stations for our analysis, including Colorado (CO), Iowa (IA), Kansas (KS), Minnesota (MN), Ohio (OH), Kentucky (KY), Nebraska (NE), and Wyoming (WY). These eight states provide a broad enough region, as well as a range of topographic and weather conditions, in which to explore our approach. The total available RWIS station count for the states are CO (147), IA (86), KS (56), KY (38), MN (98), NE (70), OH (182), and WY (81). The study period selected for this project included a winter season (October 2016 to March 2017) to best capture challenging winter driving conditions. Distribution of RWIS stations for the study area is presented in Figure4.Figure 4
Distribution of RWIS station for eight states in the United States.


### 3.2. Data Description

The description of various data sources which are used in this study is described in this section.


#### 3.2.1. Topography Data

A Digital Elevation Model (DEM) with resolution of 30 m data (28 GB) was downloaded from the USGS website (https://earthexplorer.usgs.gov/) for topographic characteristic analysis in ArcGIS 10.4.1 for twenty US DOTs. The resulting TPI map file has a size of 49 GB.


#### 3.2.2. Weather Severity Data

An ESRI shapefile containing the weather severity index used in this study was generated by Meridian Environmental Technology.


#### 3.2.3. RWIS Data

RWIS data was downloaded from Iowa State University (http://mesonet.agron.iastate.edu/RWIS/) and collected from state DOTs. The statewide RWIS data was downloaded as an excel file. Measurements from a typical RWIS station include, but are not limited to, air and surface temperature, visibility, wind speed, and road surface conditions, collected at 15- to 20-minute intervals. In total, 1026 stations are included in the analysis and 4368 hours of data was used.


### 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 asm=β0+f1x1+f2x2+⋯+fixi. Here, m is the variable of interest, β0 is the intercept, and fixi is the smooth function of predictor xi. The smooth function can be expressed as fixi=∑n=1msxn [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). Table1 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.Table 1
Descriptive statistics of RST for the study area (winter months of 2016-2017).
StatesCOIAKSMNOHKYNEWYOct-16Min-6.7-1.10.8-2.6-1.5-7.2-6.2-9.4Average15.516.920.812.618.320.016.711.9Max37.841.546.241.047.039.544.745.1StDev9.06.98.27.06.96.97.78.1Nov-16Min-29.8-8.6-7.6-9.9-24.0—-12.5-15.9Average7.69.512.56.010.1—8.75.0Max46.938.339.530.235.5—34.430.6StDev9.27.17.86.26.3—8.08.0Dec-16Min-22.6-29.5-24.9-29.9-19.6-10.4-28.1-29.5Average-1.5-2.70.8-6.61.74.5-3.3-5.4Max26.618.822.911.751.521.117.817.6StDev7.46.47.36.85.74.86.26.3Jan-17Min-26.9-24.0-20.9-30.3-19.1-17.2-26.9-29.1Average-0.6-2.12.7-6.82.96.4-2.5-5.5Max37.819.628.514.225.521.819.921.0StDev7.76.07.47.96.35.16.36.9Feb-17Min-17.7-15.9-10.8-26.0-22.7-9.6-19.4-23.4Average6.05.69.8-1.57.19.34.40.9Max37.732.438.627.234.530.436.130.7StDev8.98.29.08.17.35.98.07.7Mar-17Min-17.2-16.4-7.3-22.0-10.8-8.4-11.5-17.5Average11.07.013.92.59.012.19.47.2Max38.336.148.736.336.031.743.739.5StDev10.27.69.58.77.16.58.39.3Figure 5
Seasonal road surface temperature details for eight states.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.


### 3.4. Spatiotemporal Semivariogram Modelling Results

RWIS data for the study area were processed considering both space and time domains on a monthly basis, from October 2016 to March 2017. Road surface temperature (RST) data were aggregated using twenty-minute intervals for the time domain analysis. A space-time matrix was then formulated for each TPI and WSI zones. Spatiotemporal autocorrelation 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.Figure 6
A sample of a spatiotemporal semivariogram of TPI and WSI classes for RWIS data of November 2016.Figure 7
Spatial semivariogram ranges for (a) TPI class and (b) WSI classes.
(a)(b)According to Figure7, 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.Figure 8
Temporal semivariogram ranges for (a) TPI class and (b) WSI classes.
(a)(b)The effect of weather severity in the TPI 1 zone, which is the flatland area, is presented in Figure9. 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).Figure 9
Spatial and temporal ranges for a flatland area (TPI 1) with different weather severity.


### 3.5. 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 of100km×100km was used for density optimization. A 5km×5km 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 Figure10, 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 relationship with the impact than WSI that the spatiotemporal range has on optimal density.Figure 10
Normalized prediction error as a function of RWIS density for (a) TPI classes and (b) WSI classes.
(a)(b)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 Figure10 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.Figure 11
Number of RWIS stations for a benefit increment of 0.1 unit for (a) TPI classes and (b) WSI classes.
(a)(b)For determining the optimal RWIS density, marginal benefits (number of additional RWIS stations needed for a 0.1-unit incremental benefit) are calculated from Figure10 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 of additional stations decreases for further increments of marginal benefits for both TPI and WSI classes.Figure 12
Added number of RWIS stations for a marginal benefit increment of 0.1 unit for (a) TPI classes and (b) WSI classes.
(a)(b)Looking at Figures11 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.3-unit 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.Figure 13
RWIS density comparison for (a) TPI classes and (b) WSI classes.
(a)(b)The density optimization results are used to generate an RWIS density chart for TPI-WSI zones which is presented in Table2. 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 km2, 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 findings 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.Table 2
RWIS density for TPI-WSI zones for unit area (1/10000 sq km).
RWIS density for unit areaTPI classesTPI 1TPI 2TPI 3LBAvgUBLBAvgUBLBAvgUBWSI classesWSI 11.933.257.002.184.008.752.685.0010.50WSI 22.304.259.502.555.0011.253.056.0013.00WSI 31.983.507.502.234.259.252.735.2511.00WSI 42.755.2511.503.006.0013.253.507.0015.00


## 3.1. Study Area and RWIS Network

The study area was selected based on RWIS data availability and distribution of RWIS stations to cover varying topography and weather conditions that can be challenging for driving. Eight states had sufficient data, RWIS station coverage, and number of stations for our analysis, including Colorado (CO), Iowa (IA), Kansas (KS), Minnesota (MN), Ohio (OH), Kentucky (KY), Nebraska (NE), and Wyoming (WY). These eight states provide a broad enough region, as well as a range of topographic and weather conditions, in which to explore our approach. The total available RWIS station count for the states are CO (147), IA (86), KS (56), KY (38), MN (98), NE (70), OH (182), and WY (81). The study period selected for this project included a winter season (October 2016 to March 2017) to best capture challenging winter driving conditions. Distribution of RWIS stations for the study area is presented in Figure4.Figure 4
Distribution of RWIS station for eight states in the United States.


## 3.2. Data Description

The description of various data sources which are used in this study is described in this section.


### 3.2.1. Topography Data

A Digital Elevation Model (DEM) with resolution of 30 m data (28 GB) was downloaded from the USGS website (https://earthexplorer.usgs.gov/) for topographic characteristic analysis in ArcGIS 10.4.1 for twenty US DOTs. The resulting TPI map file has a size of 49 GB.


### 3.2.2. Weather Severity Data

An ESRI shapefile containing the weather severity index used in this study was generated by Meridian Environmental Technology.


### 3.2.3. RWIS Data

RWIS data was downloaded from Iowa State University (http://mesonet.agron.iastate.edu/RWIS/) and collected from state DOTs. The statewide RWIS data was downloaded as an excel file. Measurements from a typical RWIS station include, but are not limited to, air and surface temperature, visibility, wind speed, and road surface conditions, collected at 15- to 20-minute intervals. In total, 1026 stations are included in the analysis and 4368 hours of data was used.


## 3.2.1. Topography Data

A Digital Elevation Model (DEM) with resolution of 30 m data (28 GB) was downloaded from the USGS website (https://earthexplorer.usgs.gov/) for topographic characteristic analysis in ArcGIS 10.4.1 for twenty US DOTs. The resulting TPI map file has a size of 49 GB.


## 3.2.2. Weather Severity Data

An ESRI shapefile containing the weather severity index used in this study was generated by Meridian Environmental Technology.


## 3.2.3. RWIS Data

RWIS data was downloaded from Iowa State University (http://mesonet.agron.iastate.edu/RWIS/) and collected from state DOTs. The statewide RWIS data was downloaded as an excel file. Measurements from a typical RWIS station include, but are not limited to, air and surface temperature, visibility, wind speed, and road surface conditions, collected at 15- to 20-minute intervals. In total, 1026 stations are included in the analysis and 4368 hours of data was used.


## 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 asm=β0+f1x1+f2x2+⋯+fixi. Here, m is the variable of interest, β0 is the intercept, and fixi is the smooth function of predictor xi. The smooth function can be expressed as fixi=∑n=1msxn [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). Table1 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.Table 1
Descriptive statistics of RST for the study area (winter months of 2016-2017).
StatesCOIAKSMNOHKYNEWYOct-16Min-6.7-1.10.8-2.6-1.5-7.2-6.2-9.4Average15.516.920.812.618.320.016.711.9Max37.841.546.241.047.039.544.745.1StDev9.06.98.27.06.96.97.78.1Nov-16Min-29.8-8.6-7.6-9.9-24.0—-12.5-15.9Average7.69.512.56.010.1—8.75.0Max46.938.339.530.235.5—34.430.6StDev9.27.17.86.26.3—8.08.0Dec-16Min-22.6-29.5-24.9-29.9-19.6-10.4-28.1-29.5Average-1.5-2.70.8-6.61.74.5-3.3-5.4Max26.618.822.911.751.521.117.817.6StDev7.46.47.36.85.74.86.26.3Jan-17Min-26.9-24.0-20.9-30.3-19.1-17.2-26.9-29.1Average-0.6-2.12.7-6.82.96.4-2.5-5.5Max37.819.628.514.225.521.819.921.0StDev7.76.07.47.96.35.16.36.9Feb-17Min-17.7-15.9-10.8-26.0-22.7-9.6-19.4-23.4Average6.05.69.8-1.57.19.34.40.9Max37.732.438.627.234.530.436.130.7StDev8.98.29.08.17.35.98.07.7Mar-17Min-17.2-16.4-7.3-22.0-10.8-8.4-11.5-17.5Average11.07.013.92.59.012.19.47.2Max38.336.148.736.336.031.743.739.5StDev10.27.69.58.77.16.58.39.3Figure 5
Seasonal road surface temperature details for eight states.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.


## 3.4. Spatiotemporal Semivariogram Modelling Results

RWIS data for the study area were processed considering both space and time domains on a monthly basis, from October 2016 to March 2017. Road surface temperature (RST) data were aggregated using twenty-minute intervals for the time domain analysis. A space-time matrix was then formulated for each TPI and WSI zones. Spatiotemporal autocorrelation 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.Figure 6
A sample of a spatiotemporal semivariogram of TPI and WSI classes for RWIS data of November 2016.Figure 7
Spatial semivariogram ranges for (a) TPI class and (b) WSI classes.
(a)(b)According to Figure7, 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.Figure 8
Temporal semivariogram ranges for (a) TPI class and (b) WSI classes.
(a)(b)The effect of weather severity in the TPI 1 zone, which is the flatland area, is presented in Figure9. 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).Figure 9
Spatial and temporal ranges for a flatland area (TPI 1) with different weather severity.


## 3.5. 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 of100km×100km was used for density optimization. A 5km×5km 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 Figure10, 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 relationship with the impact than WSI that the spatiotemporal range has on optimal density.Figure 10
Normalized prediction error as a function of RWIS density for (a) TPI classes and (b) WSI classes.
(a)(b)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 Figure10 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.Figure 11
Number of RWIS stations for a benefit increment of 0.1 unit for (a) TPI classes and (b) WSI classes.
(a)(b)For determining the optimal RWIS density, marginal benefits (number of additional RWIS stations needed for a 0.1-unit incremental benefit) are calculated from Figure10 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 of additional stations decreases for further increments of marginal benefits for both TPI and WSI classes.Figure 12
Added number of RWIS stations for a marginal benefit increment of 0.1 unit for (a) TPI classes and (b) WSI classes.
(a)(b)Looking at Figures11 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.3-unit 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.Figure 13
RWIS density comparison for (a) TPI classes and (b) WSI classes.
(a)(b)The density optimization results are used to generate an RWIS density chart for TPI-WSI zones which is presented in Table2. 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 km2, 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 findings 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.Table 2
RWIS density for TPI-WSI zones for unit area (1/10000 sq km).
RWIS density for unit areaTPI classesTPI 1TPI 2TPI 3LBAvgUBLBAvgUBLBAvgUBWSI classesWSI 11.933.257.002.184.008.752.685.0010.50WSI 22.304.259.502.555.0011.253.056.0013.00WSI 31.983.507.502.234.259.252.735.2511.00WSI 42.755.2511.503.006.0013.253.507.0015.00


## 4. 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:(i)
Strong dependency of the spatiotemporal autocorrelation range of RWIS measurements with TPI and WSI values was observed(ii)
The zone with the highest topographic variation (TPI class 3 mountainous region) had a shorter range of spatiotemporal structure, whereas zones with lower TPI values (TPI class 1 flatland region) had a higher range(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 timeRecommendations 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



---

*Source: 1208692-2020-07-26.xml*