Optimizing the location of weather monitoring stations using estimation uncertainty

In this article, we address the problem of planning a network of weather monitoring stations observing average air temperature (AAT). Assuming the network planning scenario as a location problem, an optimization model and an operative methodology are proposed. The model uses the geostatistical uncertainty of estimation and the indicator formalism to consider in the location process a variable demand surface, depending on the spatial arrangement of the stations. This surface is also used to express a spatial representativeness value for each element in the network. It is then possible to locate such a network using optimization techniques, such as the used methods of simulated annealing (SA) and construction heuristics. This new approach was applied in the optimization of the Portuguese network of weather stations monitoring the AAT variable. In this case study, scenarios of reduction in the number of stations were generated and analysed: the uncertainty of estimation was computed, interpreted and applied to model the varying demand surface that is used in the optimization process. Along with the determination of spatial representativeness value of individual stations, SA was used to detect redundancies on the existing network and establish the base for its expansion. Using a greedy algorithm, a new network for monitoring average temperature in the selected study area is proposed and its effectiveness is compared with the current distribution of stations. For this proposed network distribution maps of the uncertainty of estimation and the temperature distribution were created. Copyright © 2011 Royal Meteorological Society


Weather monitoring planning as a location problem
Planning a network of monitoring stations is a topic that has attracted the attention of planners and researchers, as the advantages of proper planning and management of these facilities translate into both economic and scientific benefits. Two concurrent but incompatible main goals are identifiable in any planning action of such networks: the economic aspect favours a network with few stations, whereas an efficient network able to provide a good/accurate estimate of an observed variable is likely to require more stations. A balance between these two requirements must therefore be found. Even when this number is established, another important question is to locate an ideal set of stations that guarantee the best possible estimation results. These two questions, concerning the number and location of monitoring stations, have no immediate answer and can be difficult to resolve, given the large number of possibilities inherent to the combinatorial problems of this type, named normative location problems.
These problems, for which locational analysis deals with the formulation and solution methodologies, occur in many and varied contexts, addressing a number of questions such as the quantity, shape, size and internal relations with the facilities being located (Daskin, 1995;Church, 2002;Church and Murray, 2009). All normative location problems share some similarities, and diverse taxonomies have been suggested. These classifications explore the fundamental components that may be identified (Brandeau and Chiu, 1989;Laurini and Thompson, 1992;Hamacher and Nickel, 1998), and some of them extend the scope to include more than point location problems. In general, location problems are formulated in an economic context where demand and offer relate each other in a specific space and the location decision concerns the placement of offer while optimizing the cost and service function, called the objective function (a function describing the optimization objective in the problem). The most common objective functions are used 942 A. AMORIM et al. to designate the various problems (covering, median, centre and many more) and classify them accordingly. The elements of location problems can emerge in diverse formats. Traditional location problems include previously known sets of offer and demand elements that relate themselves by some specified distance function, turning the problem into a combinatorial search for best solutions.

Background
There are several studies mainly focusing on methods of achieving an optimized network for monitoring natural or anthropogenic phenomena, many of which deal with networks for monitoring groundwater quality: in Prakash and Singh (2000), the selection of optimal locations to expand a network for monitoring groundwater levels in a region of India is tested minimizing the kriging variance of estimation and determining where the greatest estimation errors occur. Cameron and Hunter (2002) present an algorithm for optimizing networks monitoring groundwater quality based on kriging geostatistical estimators and their estimation variance. Their aim was to find spatial and temporal redundancy scenarios that can justify the exclusion of sampling locations and the consequent reduction of the network. Also aiming at reducing the spatial density of sampling sites of a network for monitoring groundwater quality, Li and Hilton (2005) -once again for budgetary reasons -apply an ant colony optimization-type algorithm (inspired by the behaviour of ant colonies) to spot cases of spatial redundancy. A more recent study by Yeh et al. (2006) shows how the combination of multivariate geostatistical analysis with genetic algorithms is effective in defining an optimal network for monitoring groundwater quality. This minimizes the estimation variance of the spatial factor and, in the authors' opinion, provides enough information to fully understand the spatial phenomenon. The broader issue of environmental monitoring networks has also generated a large number of studies. Caeiro et al. (2002) and Nunes et al. (2004aNunes et al. ( , 2004b investigated the optimization of sediment sampling sites in the Sado River estuary, Portugal. They apply the simulated annealing (SA) optimization algorithm, searching for solutions that minimize both the variance and the average estimation error resulting from kriging interpolators.
Regarding the optimization of weather stations, several studies have been carried out, but with very distinct methodological approaches: Periago et al. (1997) develop a methodology to optimize the udometer network for rainfall observation in Catalonia. This methodology consists in building a model based on multiple linear regression (MLR) using several rainfall independent and explanatory variables. Regression residuals are analysed and the locations where they are higher are evaluated. Residuals are interpolated by kriging and the locations with the largest errors concerning estimation by MLR are obtained. The authors claim that some parts of the country are not well covered by the udometer network and there are few criteria for new locations. They also state that some of the criteria that point to a network with homogeneous spatial distribution become ineffective due to the complex topography of the country. Pardo-Igúzquiza (1998) also addressed this topic and presented a method that defines the optimal network for estimating rainfall events in river basins. By applying the SA optimization, he used two different approaches: first, he selected a subgroup of stations within a group of existing ones and then a second analysis was carried out to test the extension of that network. Each new solution was evaluated in terms of minimizing both the estimation variance by kriging and the costs of acquiring the data. The main purpose of the work was to select the location and number of pluviometers that provide the greatest estimation accuracy at the lowest cost. The author emphasizes the increase in the spatial variability for events lasting one day and its consequences when designing an optimal network. He concludes that the optimal configuration of a network for monthly observations is not identical to an optimal network for daily observations.
Geostatistical estimation is an issue that has attracted much attention. The uncertainty associated to spatial estimation can be derived from the estimation variance of estimates. However, this approach is based on two strong assumptions: the estimation errors are Gaussian distributed and are independent of the data values (see Goovaerts, 1997, pp. 259-262). Another approach is based on the estimation of the local probability function distribution conditioned to data values (see Goovaerts, 1997, pp. 262-264). These local distributions of the variables can be derived from stochastic simulations, which produce several realizations of the variable (see Goovaerts, 1997, pp. 262-264) or, simpler, estimated directly using an indicator approach.
The introduction of uncertainty of estimation by the indicator formalism has shown a high range of applications, from evaluating the distribution probability of a particular phenomenon to more complex calculations. To select optimal locations for environmental variables sampling, Pardo-Igúzquiza and Dowd (2004) apply the indicator co-kriging to quantify probability of occurrence simultaneously calculating estimation variances. Optimality is achieved by minimizing both the uncertainty of belonging to a cumulative distribution function class and the uncertainty of estimation variance. The authors emphasize the advantages of using the indicator formalism to enable uncertainty of estimation modelling as opposed to the traditional kriging, which only expresses the estimated values and the respective estimation variance.
A simpler approach involving less post-processing proposed by Bastante et al., 2005 evaluates the uncertainty of estimation in the context of decision support and makes use of information obtained from geological surveys in areas classified as exploitable/nonexploitable. Through indicator kriging, they estimated the mathematical expectancy and determined the priority areas for intervention by probability of occurrence. A similar approach is developed by Diodato et al. (2004). The authors design 943 a methodology based on indicator kriging for determining areas susceptible to soil degradation. They calculated and distributed probabilities through the indicator formalism, trying to determine the occurrence of erosive phenomena. The assessment of uncertainty allowed the delimitation of areas potentially affected by land degradation. Addressing a location problem where the goal is to establish equipment for toxic waste deposit, Rosenbaum et al. (1996) present a study on lithology classification based on samples collected from time to time. Following the indicator kriging methodology, the authors convert each lithological type into a binary variable, subsequently estimating the cumulative distribution function. The results obtained by these authors show not only the areas of uncertainty of estimation but also the probability of a particular type of lithology belonging to a determined class. Several authors are mainly concerned about the modelling of the probability cumulative distribution function and the aid it can provide when evaluating and solving several problems. The study by Amini et al. (2004) defines areas contaminated by heavy metals. This task generally involves an inherent uncertainty that must be taken into account when making decisions. This way, the authors estimate the probability of exceeding the limit for cadmium and lead concentration by the indicator formalism. The work by Fabbri and Trevisani (2005) raises the question of which methodology is appropriate for assessing the feasibility of exploitation for geothermal waters. These authors also propose an approach using indicator kriging describing its formality. They state that indicator kriging is a good estimator of the cumulative distribution function and afterwards simulates a calculation of the uncertainty of estimation. Although several references converge on the idea of this being a suitable approach for assessing the modelling uncertainty of various phenomena, there are others that draw attention to some constraints. Juang et al. (2003) designate the indicator kriging as one of the main methods of modelling the uncertainty of estimation. However, the authors stress the smoothing effect (overestimation of low values and underestimation of high values) it has in map production, emphasizing the particular case of pollutant modelling. In this work, the cumulative distribution function is calculated in the sample sites and then estimated by kriging and conditional simulation. When analysing the results, authors realized that kriging performs a much smoother distribution and the estimated values vary in the range of observed values. Lloyd and Atkinson (2000) evaluated three different geostatistical approaches -ordinary kriging, kriging with external drift (KED) and indicator kriging -to estimate the elevation. In this case, the authors noticed that indicator kriging consumed significantly more processing time when compared with other methods. However, despite a slower processing and a more elaborate formalism, this technique allowed the assessment of associated uncertainty of estimation, while the other methods could only be evaluated by their estimation variance.

Methodology
In this work, a methodology for planning a network of weather stations to monitor average air temperature (AAT) at soil level is presented. This methodology focuses on the optimization (by minimization) of the geostatistical uncertainty of estimation, modelling the network planning scenario as a location problem, but considering a changeable demand surface.

Solution methods for location problems
Solution methodologies for location problems are usually included in either exact methods or heuristic approaches (Daskin, 1995;Church and Murray, 2009). The choice of one of these top-level methodologies depends on the size and complexity of the problem. Exact methods look for and output the best solution for the problem (the global optimum or optima), sometimes exploring exhaustively the full solution space or using branch and bound techniques. Global optima are only achievable when the problem size is not too large, i.e. its computational complexity is affordable in terms of time and resources required for processing. On the other hand, heuristic approaches combine the search of good or fair quality solutions with the limited computation time or space available for solving complex and large-scale problems, dropping the requirement of achieving a global optimum. We will only discuss two of the latter methodologies that were used in this work, namely a greedy construction heuristic and the SA approach.
Greedy construction heuristic methods operate on a simple iterative process, starting from an empty initial solution and selecting, at each step, the best available alternative, taking it into to the solution set, and memorizing the best solution found throughout the process. The evaluation of each solution is based on the objective function. Two main strategies are developed: adding or dropping elements to or from the initial solution set.
In the particular case of point location problems, greedy algorithms present solutions based on the successive addition or reduction of elements. This is repeated until some halting conditions are met, such as reaching a preset number of points or crossing a threshold value in the objective function.
Some recent applications prove the efficiency of this type of heuristics when compared with other approaches. Suh et al. (2006) concluded that, in most cases, the greedy heuristics produced solutions quite close to optimal in the particular case of the location of traffic volume monitoring points on a communication network. In this study, the goal was to minimize the cost by reducing the number of monitoring points and maximizing the relevance of the information collected (less redundant).
A frequent concern in optimization methods is the possibility of converging and being trapped in local optima. To prevent this, distinct strategies have been proposed in the literature, either specific heuristics or generic metaheuristics. The use of the SA technique as an optimization metaheuristic follows the principles presented by Metropolis et al. (1953) on statistical thermodynamics. SA was proposed by Kirkpatrick et al. (1983) as an algorithm to solve the well-known combinatorial optimization problems. It reduces the risk of falling into local minima (or metastable solutions) common to iterative improvement methods, as unlike the latter, it accepts solutions that worsen the objective function. The authors proposed the use of the Metropolis (Metropolis et al., 1953) procedure from statistical mechanics, which generalizes iterative improvement by incorporating controlled uphill steps (to worse solutions). The procedure considers one small random change in the system at a certain temperature (T ), the name for a control parameter that decreases during the execution; the change in the objective function is OF; if OF ≤ 0, then the change in the system is accepted and the new configuration is used as the starting point in the next step; if OF > 0, then the probability that the change is accepted is determined by P ( OF) = exp(− OF/k b T ); a random number uniformly distributed in the interval (0, 1) is taken and compared with the former probability; if this number is lower than P ( OF), then the change is accepted. Hence, SA surpasses local optima and the dependence on values in the initial solution by considering lower quality solutions in a solution pool, and this may lead to different and better solutions in following iterations. The introduction of the decreasing temperature parameter reduces the occurrence of 'drastic' changes and makes them less probable in time, while a pattern for local search around the current solution is noted (D'Amico et al., 2002;Miller, 1996).

Modelling the problem
Modelling the choice of the best places for weather monitoring stations as a normative location problem, the offer is clearly composed by the set of possible places to locate stations, i.e. the location space. The objective function is the minimization of the maximum value in the uncertainty of estimation surface.
However, the demand is not fixed, as in the iterative process of locating the stations, the uncertainty surface changes because it is based on the uncertainty of estimation, which is different for each set of stations. In general, demand elements are fixed, i.e. their parameterized behaviour and configuration do not change with the location of a new facility, but when a new monitoring station is located, a new uncertainty of estimation must be calculated. Demand is then modelled as a variable surface depending directly on the offer distribution. After evaluating a spatial configuration of the existing stations and its effects on the temperature uncertainty of estimation, the optimization process starts and outputs the number and location of stations.
Following the evaluation of how stations are distributed and its effects on the uncertainty of estimation, the optimization stage determines how many and which places could host a station, enabling future estimations to achieve the least possible uncertainty. Not having a fixed number of stations, the goal is then to optimize the minimum quantity that makes the goal feasible. The spatial distribution of the offer elements is decisive in the problem. This distribution, together with the estimation variability, will determine the uncertainty. Therefore, the problem shows some comparable aspects with the covering model: adding a new station helps to reduce uncertainty of estimation in nearby areas, just like in covering problems a similar decision implies that all nearby demand elements are covered. Coverage is established for demand elements that are below a given uncertainty threshold value. Elements that demand new facilities are all those above the threshold. The problem, however, also shares some similarities with the typical centre model, since the objective function is a minimization of the maximum value exhibited by the uncertainty of estimation surface.

Producing solutions
To compare distinct solutions (i.e. the configurations of stations), a value must be assigned to each one. As described, uncertainty of estimation is used as a quality measure in this process.
The SA algorithm is used in the first part of the optimization procedure to test the permanence of each station under a reduction scenario. For this, the maximization of the average value of the uncertainty of estimation using the indicator kriging formalism is considered.
Using a reduced set of stations as a starting point, a greedy construction algorithm can be applied, providing solutions for the network expansion stage. The reason for choosing a distinct algorithm for the second stage lies in the fact that the solution space for this additive procedure is much larger than that in the previous step, as all the cells in the study area are candidates to locate new stations. The greedy algorithm pseudo-code is as follows: 1. Calculate the uncertainty for the current solution; 2. Find the highest value of uncertainty of estimation; 3. Get the estimated temperature value in the location found in the previous step; 4. Add a new element, simulating the location of a new station in the highest uncertainty of estimation place, considering its temperature value as the estimated value; 5. Repeat steps 1-4 until decision criteria are met.
In this process, it is important to evaluate how the uncertainty of estimation behaves when facing a progressive increase in the network, namely considering the newly included locations. These iterative processes require a repeated calculation of the uncertainty of estimation as they produce and evaluate the various solutions.

Case study
The methods were used to optimize a set of classical climatological stations operated by the Portuguese Institute of Meteorology (IM). AAT was selected as the study variable. Analyses focused on the value of the annual AAT.
The following steps are described for the case study: (i) data collection and analysis: in this step, a descriptive data analysis is performed. Using explanatory auxiliary information, the variables that describe the AAT are presented; (ii) geostatistical modelling of the variable and evaluation: in this step, models of geostatistical interpolation were essayed and evaluated for each of the three periods; (iii) execution of the optimization process: in this step, reduction and expansion scenarios were produced using the optimization algorithms of SA and greedy heuristic, respectively. Further detail on all the steps is provided in the following sections.

Data collection and analysis
A total of 89 stations in the Portuguese territory was selected and then complemented by 34 stations in the neighbouring country, Spain, located in an area close to the border. For these stations, data acquired from the Spanish Institute of Meteorology was available. This included the climatological normals for the period 1961-1990, for annual monthly average values. Not all the 123 observation stations had 30-year sets, and a minimum 20-year observation period was considered. After the compilation of data, the location of stations was carefully verified with the collaboration and technical expertise of IM technicians.
Air temperature in the Iberian Peninsula is known to have a reasonable correlation with a few spatial variables that might be used as auxiliary information in spatial modelling (Ribeiro et al., 1997, pp. 394-398;Mesquita and Sousa, 2009). The most used auxiliary variables in modelling climate-related characteristics, such as elevation, slope, aspect and distance from the sea, are directly obtained from a digital terrain model (DTM). After the delimitation of the study area, a hybrid DTM with 1km spatial resolution was produced combining shuttle radar topography mission (SRTM) data for the Spanish territory and altimetric information used in the 1 : 25 000 cartographic series from the Portuguese Army Geographic Institute described in Néry and Matos (2005). Derived data sets were produced and defined as auxiliary matrices. After verifying the absence of expressive colinearity between these independent variables, the condition to include them in a regression model was guaranteed. The statistical correlations between AAT and these variables are presented in Table I. AAT is significantly affected by the orography and distance to the ocean. Slope presents a slightly less expressive correlation, while aspect is not interesting for individual application in AAT modelling.
The analysis of the correlations between dependent and independent variables suggests a multiple combination of the independent variables, and an essay of multiple linear regression (MLR) was conducted with Statistica  software. The AAT variable was described through three of the four independent variables (E, S and D), presenting a correlation coefficient of 0.869 with the original variable.
As the new MLR variable presents good statistical correlation values with the AAT, this was consequently incorporated in the estimation models.

Geostatistical modelling
Local distributions of the variable were estimated based on the algorithms of the indicator KED, using the abovementioned variables obtained by MLR. To minimize order relation violations in the indicator kriging process, and to outline the difficulty in the variogram adjustment of ill-defined extreme indicators, an indicator median variogram was applied.
After the estimation of the local distribution, it was possible to locally calculate its average and variance and use them as estimate values and associated estimation uncertainty. Unlike kriging variance, this local estimation variance depends on the local values of the variable.
Once the geostatistical estimation was performed, its results were evaluated: external validation and study of errors, correlation calculation between estimated and observed values, and visual assessment of the produced maps.
The main goal of including auxiliary information in uncertainty of estimation through the indicator formalism was to evaluate its effect in the uncertainty spatial distribution. A clear example of such effect occurs in the area of Serra da Estrela, a mountain range in the centralnorthern interior region, where the inclusion of auxiliary information contributed to a decrease in uncertainty (Figure 1). The inclusion of elevation enabled a better modelling in the areas with high variability. Places where the sampling is scarce and where AAT displays a high variability explained by the behaviour of an associated variable are revealed.
The study on the variability of both uncertainty of estimation for AAT and its explanatory factors demonstrates that some areas require a higher density of stations than others. This happens because some areas present higher variability, while in others, the uncertainty of estimation is more homogeneous.

Optimization
The 89 stations were evaluated in terms of their spatial representativeness, using the SA algorithm to test the permanence of each station under distinct reduction scenarios. A maximization of the average value of the In the first stage of the optimization process, the SA algorithm was applied to evaluate the contribution of each station to the network total variability and to assess spatial redundancy. At this stage, several network scenarios were produced, analysing the uncertainty of estimation behaviour when facing this reduction. Runs for ten reduction scenarios of the current 89-element network were made for 85, 80, 75, 70, 65, 60, 55, 50, 45 and 40 stations, replicating a gradual decrease in the network. Each solution contains the stations that the algorithm considered the most important concerning spatial variability, removing the five most redundant stations in each scenario. Only the stations that add significant information (variability and uncertainty) should be included in the solution. In the SA process, for each of the ten scenarios, optimization runs were executed for three random initial solutions (S1, S2 and S3), and for each one of these, three cooling coefficients were applied (0.90, 0.95 and 0.97), totalizing 90 different runs. Each of the 30 initial solutions was independently generated from the 89-station starting set.
Tests were also made with the SA algorithm for optimization of the monitoring networks considering 2, 3 or 4-opt neighbourhoods (i.e. altering the number of stations allowed to be changed in each iteration). However, these tests did not improve the optimization results or made the convergence faster, on the contrary. In fact, these changes alter the objective function too much to allow good convergence. SA is known for being a very efficient algorithm in what regards convergence, while requiring long processing times. Given the complexity of the problem treated here, the trade-off between quality of the maxima and convergence time has been considered on par or better than competition. Table II presents the values of optima of the objective function obtained by running the SA algorithm with the parameters identified before and with the three random initial solutions. Using the best values for each solution Table II. Maximum values of the uncertainty of estimation, in°C, for each run of the simulated annealing algorithm using distinct random initial solutions S1, S2 and S3.  found by the SA, independently from the initial solutions or the used cooling schemas, it was possible to compile a single value for the objective function, for each of the ten evaluation scenarios. The same solution is frequently found using the same starting set of stations despite the distinct values of the cooling coefficient. Some differences are noticeable among the best values of the objective function in runs generated from different initial solutions. This behaviour reflects the intrinsic variability of the algorithm which is generally lower than 10%, leading to the conclusion that the results are consistent.

Network reduction
The application of the SA optimization algorithm has confirmed that decreasing the number of stations usually leads to a gradual increase in uncertainty of estimation ( Figure 2). Notice that in each scenario, the displayed solution is the most spatially representative and the least redundant. Figure 3 displays the scenarios for the reduction produced by the SA algorithm over the base map of the estimation uncertainty obtained for the initial set of 89 stations. In each suggestion to reduce the network, a range of combinations of sampling sites was considered. A tenstation reduction (using a ten-station decrease analysis) does not imply that stations that were previously removed (or a combination of ten stations that were removed from a previous scenario) will be repeated again. The network of 80 and 70 stations showed a removal of sampling sites distributed more or less evenly throughout the national territory. In some of the stations, the algorithm eliminates are close to each other, but this is not a decisive factor: when reducing from 89 stations to 80, not all the stations removed by the algorithm are located near other stations. It is noticeable, in the 80-station scenario, that some stations were removed in areas of lower uncertainty of estimation. In the 70-station scenario, some sampling sites were also removed in areas of lower uncertainty and in the vicinity of other sites. In the 60-station scenario, a significant reduction in stations in the centre-north area is noticeable, while the stations in the mountainous regions of Serra da Estrela and Douro Valley (northeastern zone) are preserved. In the South, although in the interior only a few stations are located, the Algarve coast (southern) retains most of its stations. In the scenario considering a reduction to 50 stations, the removed stations seem to be homogenously distributed in the study area, except in the flatter area in the centre-south where they disappear almost completely, and in the southern coast which has its stations preserved. These results are in accordance with the correlation found between AAT and elevation. Table III describes the statistics of the initial 89station set, of the set of stations in each scenario, and of the corresponding removed stations. In the analysis of these indicators what seems to be more significant is the difference in variance values. The remaining stations have higher variance values than the stations that were removed. This falls in line with the defined objective function, in which only stations with the highest uncertainty and variability were included in the optimal solution. The variance is noted to increase as the network drops stations. Both the mean and median AAT were similar in the two sets, suggesting that no relevant information about the average AAT values was lost. On the other hand, stations with more extreme AAT values (minimum and maximum) were preserved in the solution sets, but this did not force a significant bias in the calculation of central tendency statistics. These results indicate that the algorithm performed reasonably well, having produced monitoring networks with only marginal differences in the optimal parameters.
The solutions with 50 and 55 elements were excluded due to a more pronounced increase in the objective function ( Figure 2). Maps of estimations with the solutions obtained in each scenario were produced, and the relative difference maps between consecutive settings were evaluated. The most significant changes occurred when reducing from 70 to 60 stations, so the 70-station solution was chosen as a base network. This solution maintains most of the statistical characteristics of the initial set and was, therefore, considered a good candidate for the network expansion stage.

Network expansion
In the previous section, the identification of station representativeness and the output of a reduced network that eliminates redundant stations were presented.
In this section, we describe an essay on the network expansion, aiming to reduce the maximum value of the uncertainty of estimation. The algorithm developed for this particular application iteratively searches for locations where the highest uncertainty of estimation in the Portuguese territory occurs. While SA was used for reducing the dimension of the network, a greedy algorithm was chosen for the following step, where the size of the network increases. The choice of a distinct approach was that a faster -though less efficient -algorithm was needed for the increase, as the number of combinations made it impossible to use a more efficient but slow algorithm, such as the SA. In our case, the solution space (i.e.   . Optimal network of 70 weather stations proposed by simulated annealing and 10 new stations proposed by the greedy algorithm on estimation surface: respective uncertainty for a new network configuration (left) and estimated average temperature using KED with multiple regression as drift (right). This figure is available in colour online at wileyonlinelibrary.com/journal/joc number of combinations) when reducing was in the order of 10 25 solutions, while during the increasing stage, the solution space was much larger (100 stations could be located anywhere in a mesh with slightly less than 90000 cells with 1 km side). This greedy algorithm iteratively places stations until the uncertainty value (or the number of stations predefined by the user) is achieved. For each new proposed station, an estimated value of the annual AAT is adopted using the KED estimator with the MLR surface as drift, as referred previously. Figures 4 through 6 show the solutions proposed by the greedy algorithm for an increase in stations reflecting three new expansion scenarios with 80, 90 and 100 stations.
Up to 30 new stations were found scattered all over the Portuguese mainland, with a higher incidence in areas previously considered of greater uncertainty. It is interesting to verify that the new stations are mainly located in mountainous regions, thus demonstrating the lack of observation points in these areas. It is also observable that some locations are reinforced by the location of nearby temperature monitoring stations.
In the 70-station solution, the maximum uncertainty of estimation value was 4.630°C. As Figure 7 shows, the maximum uncertainty of estimation of the annual AAT decreases as new weather stations are added. There is a slightly stronger decrease in maximum uncertainty up to 80 stations. With the location of 30 new stations, there was a decrease of 2.097°C in the maximum value of the uncertainty of estimation, which represents a 45% reduction.
The implemented greedy algorithm suggested the location of new weather stations in areas with the larger values of the uncertainty of estimation. This provides not only noteworthy suggestions for future plans of new monitoring stations but also information concerning the uncertainty of estimation behaviour. It was also noted that mountainous regions, especially those located in the northern and central parts of the country, are the areas that require further densification of the observation network. An interesting result occurred in the 89-station expanded network, where the maximum uncertainty of estimation was significantly lower than that in the initial solution set. The greedy algorithm found locations coincident with areas with lack of sampling, mostly in mountainous regions. This fact might be utilized as a validator for the SA step.

Conclusions
This research has proposed a model and a methodology for addressing the problem of optimizing a network of weather monitoring stations. The methods explored in the methodology are usually applied to a series of location problems with unvarying demand surfaces that use concepts as the distance-based coverage. However, in the case of climatological data collected by weather monitoring stations, which generally is extended to an entire area using interpolation techniques, the quality of a network might be evaluated using other descriptors, such as the uncertainty of estimation that a particular network layout induces. In this case, models for optimizing the network will require the presence of a demand surface that changes with each modification in the network. The proposed approach models the optimization process for networks of weather monitoring stations using the uncertainty of estimation as a demand surface. The maximization of the average value of this surface, generated for each particular combination of stations, was used as a quantifiable element capable of expressing the individual spatial representativeness of the network features. Using SA, it was also possible to identify a reduced monitoring network that maintained most of the statistical characteristics of the initial set and was used as a base network of the most significant stations.
The greedy construction algorithm step was adequate to the location of new stations, considering a changeable demand surface and a very large solution space. New stations were located in regions with larger values of uncertainty of estimation, providing important suggestions for locating new monitoring points. An interesting result is the information on the behaviour of the uncertainty of estimation surface, which provides complementary data for the assessment of weather station locations.
The greedy and SA algorithms do not guarantee optimality. The greedy approach was used due to its properties of good performance and speed, although SA is generally considered to have better performance than the greedy. However, its applicability to very large problems, such as the expansion example presented here, is very limited because of prohibitive running times.
The findings in this article are related to three aspects. As the number of facilities decreases in the first part of the optimization, it is expectable that the average value of the uncertainty of estimation, which is the value under analysis, has a gradual increase, but it is not always the case. As this does not usually occur in the most current models of optimization, it is necessary to explore in each case to what extent a minimized network provides reasonable values of the uncertainty of estimation. It was also noticed in the optimized network that some of the new stations proposed by the greedy algorithm match existing stations previously eliminated by SA. So, applying the greedy algorithm is interesting in terms of not only the expansion and densification of the network by locating new stations but also for analysing proposals arising from the application of SA. Finally, it is possible to apply the proposed model in other location problems that require changeable demand surfaces, such as networks used as value providers in interpolation processes, since no particular assumptions on the variable intrinsic characteristics were imposed.