# V CONCLUSION

In our explorative analysis we prove the potentialities and limits of a vector-based approach for analyzing large and complex terrains. Further investigation will be also conducted in order to improve the usability of the SketchUp extension developed in this research. Thus, every territory covered by a DEM could be analyzed, by dynamically changing the terrain site.

Many parameters can also affect the modelling of a site. Closeness has to be defined according to the site location. We proposed a 300m buffer but this parameter could be subject of a sensibility analysis. This analysis suggests that there should be an optimal obstacle region radius for sunlight access assessment, depending on each site. Generally speaking, adaptive modelling should be the aim of all site analysis. Also, the accuracy of the skyline silhouette also alters the solar exposure. It relies on the distance between contour lines. A sensibility analysis should also be led on that parameter to determine how the level of detail of the terrain model affects the calculation of solar radiation indicators.

We proposed that the solar exposure could be defined efficiently thanks to a site-centric modelling. The relevant obstacle region around the observing point depends on the urban fabric. This approach could be extended to other case studies, for a wide set of density range: from the open field to the high medieval urban fabric.

# References

* Floriani, L. De, Magillo, P. (2003). Algorithms for visibility computation on terrains: a survey. *Environment and Planning B: Planning and Design*, 30(5), 709-728. http://doi.org/10.1068/b12979

* GDAL (2016). GDAL - Geospatial Data Abstraction Library, Version 1.10.1. GDAL Development Team, Open Source Geospatial Foundation. Retrieved from http://www.gdal.org

* Morello, E., Ratti, C. (2009). Sunscapes: “Solar envelopes” and the analysis of urban DEMs. *Computers, Environment and Urban Systems*, 33(1), 26-34. http://doi.org/10.1016/j.compenvurbsys.2008.09.005

* Morello E., Carneiro C., Desthieux G., (2010). The use of digital 3-D information to assess urban environmental quality indicators. In Schmitt G. et al. (eds.), *Future Cities: Proceedings of the 28th International Conference on Education in Computer Aided Architectural Design in Europe*, Zurich, Switzerland, pp. 388-395.

* Prévot, A., Rodriguez, D., Molines, N., Beckers, B. (2011). La modélisation 3D: une nouvelle voie pour les documents d'urbanisme ? Application à l'optimisation énergétique des bâtiments. *Revue Internationale de Géomatique*, 21(4), 557-583. http://doi.org/10.3166/rig.15.557-583

* Tandy, C. R. V. (1967). The isovist method of landscape survey. *Methods of Landscape Analysis*.
# Combination of error characterization and spatial error model to improve quality of digital elevation models

Tomaž Podobnikar*¹

¹ University of Ljubljana, Faculty of Civil and Geodetic Engineering, Slovenia

*Corresponding author: tomaz.podobnikar@fgg.uni-lj.si

**Abstract**

A processing chain for spatial data cleaning has been developed, in order to produce more reliable digital elevation model (DEM). The first step of this heuristic process is error characterization in a combination of (geo)statistical, empirical and visual exploratory analyses, which were experimentally realized with clustering analysis. The second step is developing a spatial error model that uses homogeneous regions derived from terrain parameters, such as roughness and land cover, in relation to reference information. To support the methods of error characterization and spatial error model in spatial datasets, a novel ‘locally systematic’ error type has been proposed. This more complex error type characterizes (spatial) positional error field, which is supplementary to traditional random, systematic and gross error types. It is supposed to be the most typical with its geomorphological and topographic components. It comprehends geomorphological shapes, patterns, structures, i.e. semantics of the terrain surface that is usually oversimply considered as a random error. The results show that a comprehensive understanding of the error properties through error characterization supports effectively processing of the individual variables and from them a reliable error surface, whereas both are useful for effective removal different types of error, with an exception of the random one.

**Keywords**

DEM, digital elevation model, DTM, digital terrain model, spatial data quality, statistical/empirical/visual quality assessment, uncertainty, error characterization

# I INTRODUCTION

The quantity and availability of digital spatial datasets rapidly increase, especially widely applicable earth observation (EO) data provided by a series of satellites. The digital elevation model (DEM), or often called as digital terrain model (DTM) was one of the first forms of digital geographical information which became available (Fisher and Tate 2006) and it is now a fundamental dataset to mapping our world. The main reason of recently improved DEM usability are new methods for data acquisition that have been introduced, especially from passive and active sensors on satellites, airborne, and unnamed aerial systems (UASs). These data are today considerably higher (spatial) resolution than before and therefore more comparable to the other EO. There is also a growing number of public accessible application platforms that additionally increase the usability of the DEM. However, the spatial characteristic of the DEMs quality is not always in accordance with their resolution and user's expectations. Additionally, there is still a gap between what the quality assessment experts can produce and the information that users can understand and use (Devillers et al. 2010).
The aim was to develop a solution that combines a number of different aspects of the DEM quality through error/uncertainty characterization and a spatial error model as an error surface or field, in order to get a better quality DEM. The proposed solution basically contributes in the DEM production chains with a spatial data cleaning procedure (detection, validation, correction; Kimball and Caserta 2004) as part of quality assurance. For both, the already produced DEM, or for implementation in the future DEM processing chain, a better quality DEM can be produced. An important background of the proposed solution is a conceptualization of the DEM taking into the account fitness for use – not only for multi-purpose demands where the main aim is to process a near-universal DEM – but also for specific demands to produce different DEMs for various applications. E.g., for the near-universal DEM we can define that the model excludes vegetation, snow cover, different buildings and bridges, but include glaciers. It is also good to know what kinds of errors (Leon et al. 2014) and their magnitude are people willing to accept.

## II COMBINED APPROACH

The proposed solution is based on a combination of statistical, empirical and visual (with eyes inspection) quality control and assurance, i.e. on quantitative and qualitative approaches. The most significant elements of the spatial data cleaning process is a combined approach for error characterization and a spatial error model building (Figure 1) in order to generate a better quality DEM. Order of data cleaning methods is important; therefore, a decision-making process is necessary to adopt final procedure among several alternative possibilities resulting from different methods.

<figure><img src="image_147.png" /><figcaption>Figure 1: A processing chain of spatial data cleaning to produce an improved DEM from the already well processed DEM.</figcaption></figure>

In order to effectively test the results of the proposed solution, we located the study areas in different parts of Slovenia, where DEM from contour lines, National DEM, lidar DEM, geodetic spot elevations database and other reference datasets, such as land cover and geological map were available; as well as the most typical landscapes were included.

### Terrain parameters and regions

The expected DEM error/uncertainty/quality characteristics and magnitude are related to terrain parameters, i.e. terrain complexity/morphology, sampling density, which have been observed by numerous studies (e.g. Carlisle 2005, Oksanen and Sarjakoski 2006, Erdogan 2010). Podobnikar (2005) suggests correlating the error field with terrain skeleton (peaks, pits, saddles, ridges, talwegs), roughness, slopes, and rate of vegetation cover. Guth (2006) proposed geomorphometric variables that appear to be most important in terrain description and which
access SRTM quality: steepness of roughness, elevation-relief ratio or coefficient of dissection, organization, heterogeneity, profile curvature, plan curvature, average elevation, fabric direction, massiveness, modality, and limitation. The high-resolution lidar data and derived DEMs are expected to vary with land cover types, such as open terrain, high grass weeds and crops, brush lands and low trees, fully forested and urban areas (Zandbergen 2011). Many other papers classify the elevation error of lidar DEM similarly as open terrain, weeds and crops, scrub and bushes, forested, build-up areas (Hodgson and Bresnahan 2004). They classify the open terrain data as highest and fundamental accuracy, where the assessment of the other terrain types is supplementary. The ISPRS (by Sithole and Vosselman) provides fifteen benchmark lidar point clouds on the arbitrary defined representative environments, on which selected error types of the produced DEMs can be compared since 2003.

The DEM error is also determined by geology (rate of karst – karstness), hydrological network; anthropogenic impacts on terrain (e.g. stone quarries, buildings, transport network infrastructure); larger areas of plains, hills and mountains; precipitation and temperature (e.g. in karst areas with the same geology the terrain is considerably rougher in mountains than in lower areas). Nevertheless, the predicted error very much depends on the data sources used (e.g. photogrammetric, lidar), scale/resolution and related information entropy or fractal dimension (Wise 2012), as well as on the interpolation/filtering methods used and on particular producers’ style of data processing.

Practically, the terrain parameters were derived from DEM, land cover data and geological maps. The terrain parameters were either continuous surfaces or discrete datasets. All these parameters were then separately classified to significant regions.

## Error characterization (for the 1st correction)

Bearing in mind traditional error theory there are three types of errors: random, systematic and gross. Typical systematic (e.g. wrong vertical geodetic datum) and gross errors (e.g. wrong altitude of particular locations) can be detected and removed. In reality, most of the DEM errors are not random, not stationary (spatially variable) and not from and identical normal distribution (Zandbergen 2011). Statistical methods are appropriate for error assessment only if a substantial number of checkpoints is used, what is rarely feasible in real applications (Liu et al. 2012). Therefore, the validity of the most applicable root-mean-square error (RMSE) metrics is questionable, particularly because the main error in DEMs is so called systematic-like (Oksanen and Sarjakoski 2006) or locally systematic (Podobnikar 2008); which is not normally distributed (e.g. as tested with the Q-Q-plot). A robust spatial statistics (median, quantiles) or qualitative methods are needed in this case (Höhle and Potuckova 2011). Consequently, we propose four main groups that characterize a positional error in DEM:

* random error – (not removable – vertical and planimetric); differently autocorrelated and scale dependent (Fisher and Tate 2006)

* locally systematic error – (potentially partly removable and partly inherent – vertical only); it is: (1) very complex and most typical error in DEM, but usually oversimply considered as a random error or even uncertainty; (2) similar to the concept of low and medium frequency errors (Imhof 1982); (3) spatially autocorrelated (e.g. in relation to sampling density), in the context of the certain area (mostly constrained) geographic neighbourhood (e.g. in relation to underestimated altitude of peaks); (4) topographic and semantic with specific shapes and patterns (e.g. in relation to geomorphological forms and structures); in this context the DEM is considered as a surface and not as independent individual points
* *systematic error* – (mostly removable – vertical and planimetric)

* *gross error* – (mostly removable automatically or manually – vertical and planimetric); wrong attributes of altitudes are owing to not appropriate filtering of lidar point cloud data or interpolation; it is often a semantically recognizable error due to misinterpretation in manual digitalization of maps or in stereo-photogrammetry

There are numerous different measures and indices used to characterize DEM error into main proposed or other groups, such as to a group of (geo)statistical, empirical and visual; or to a group with an error assessment on a DEM only, or using reference datasets. Here are some basic examples of this kind of exploratory analysis: RMSE calculation, using shaded DEM or contour lines calculated from DEM, relating independent hydrological network and DEM, applying Monte Carlo simulations and sensitivity analyses, computing histogram or swath profile, etc. The proposed measures and indices are used for error identification and consequently for removal gross errors (e.g. with correction/transformation attributes, interpolation/filtering, filling gaps with other datasets) and systematic errors (e.g. with transformation or georeferencing).

To practically demonstrate the error characterization, we developed statistical cluster analysis with visual interpretation, on the empirically processed terrain parameters: DEM 20 elevations, slope, aspect, curvatures, and an error field as difference between the DEM 20 and a lidar DEM, which was used as a reference ground truth (Figure 2). This result helps better understand patterns of different types of errors that occur in the target DEM, and removal gross and systematic errors (in our case local horizontal shifts).

<figure><img src="image_148.png" /><figcaption>Figure 2: Error characterization with clustering analysis (area of Koroška Bela, NW Slovenia, 2.2 x 2 km).</figcaption></figure>

## Spatial error model (for the 2nd correction)

Our philosophy for spatial error model building is similar to reverse engineering, where we need to extract a knowledge/semantics (in our case through the error characterization process) and
then reproducing it (in our case with processing the spatial error model). An advantage of the spatial error model calculation is to better understand, control and remove the errors in DEM. This model has been also used for learning purposes, in order to understand the key issues of potentially unsuccessful DEM processing design and subsequently to improve the design.

Various kinds of DEM error model are known, such as spatial conditional stochastic models with Monte Carlo simulations (Podobnikar 2008), or interpolated from checkpoints with kriging (Leon et al. 2014). Comparable is an approach with (linear) regression models that define a relationship between DEM error and morphometric character (distribution and scale of elevation; Carlisle 2005). Due to the proved non-stationarity of parameters that vary across space, a local GWR technique is a valuable solution, too (Erdogan 2010). An empirical error model that assumes correlation of the DEM errors with terrain morphology, sampling density, and interpolation method, has been developed based on approximation theory (Liu et al. 2012). This model doesn't require a reference dataset, i.e. any checkpoints, but only the source data and the DEM.

Our proposed spatial error model is empirical and predictive. It basically determines locally systematic error as a combination of different variables, where the procedure is following:

* significant terrain parameters are processed through defining indices that influence to a DEM error; this process requires knowledge acquisition about error sources and types

* homogeneous regions within each terrain parameters are calculated with classification (and employing segmentation methods, where the aim to establish these regions is similar to applying stratification used for better spatial data sampling): in each region a stationarity of the stochastic process is assumed, i.e. the error within regions needs to be homogeneous with a certain threshold; statistical tests using reference points and DEM, or operator's knowledge help in the region definition

* variables $V_i$ are calculated for corresponding terrain parameters with an assignment of the unique error values to each region

* empirical spatial error model (SEM) is processed as: $SEM = \sum_{i=1}^{n}(V_i \cdot w_i)$, where $V_i$ is variable (e.g. errors in regions of terrain parameter 'roughness'), and $w_i$ is weight on interval [0, 1], where $\sum_{i=1}^{n}w_i = 1$

* spatial data cleaning is applied using spatial error model and variables

For the practical application, we developed a spatial error model (SEM). We propose following significant terrain parameters: RR for main terrain characteristics, RN for slope of terrain, an RG for land cover information and RO for terrain skeleton. Homogeneous regions for all four terrain parameters were calculated. The RR was classified into following four classes: flat surface (plains), low hills, hills, and mountains, according to the slope, profile curvature, and elevation. The standard region RN was classified into three classes: 0°-5°/5°-20°/> 20°. The RG was classified into categories of bushes, deciduous forest, mixed forest, coniferous forest and open areas. The last one, RO was classified into regions of peaks, pits, saddles, ridges and talwegs. RMSE values were assigned to terrain parameter's regions to get variables $V_i$. We applied exploratory statistical tests, with regards to our previous knowledge and experiences. Finally, the SEM was developed. Weights were assigned empirically: $w_1 = 0.5$ for RR, $w_2 = 0.27$ for RN, and $w_3 = 0.23$ for RG. A sum up of all three variables bearing in mind the proposed weights was calculated using map algebra operators. Finally, the RO was overlaid in order to calculate a spatial error model (Figure 3a).
For assessment of the spatial error model quality, we calculated an *absolute error* (Figure 3b). It was calculated as an absolute difference between a DEM from contour lines digitised from topographic maps in scale 1:25,000, and considerably better National DEM 12.5 used as a reference, which was produced mostly with stereo-photogrammetric methods (Podobnikar 2005).

The result (Figure 3) demonstrates very similar spatial error pattern for both, the predicted and the calculated error. It is obvious that the prediction model considers the right terrain parameters, variables, and weights. However, the most obvious differences between both are due to interpolation problems (encircled in Figure 3b) and other not (yet) eliminated gross errors.

<figure><img src="image_149.png" /><figcaption>Figure 3: A DEM error field (darker is area, larger is error – up to ~8 m): a) spatial error model (prediction); b) absolute error, calculated as difference between independent reference and target DEM 12.5 m (area of surroundings of Ljubljana, central Slovenia, 20 x 15 km).</figcaption></figure>
The spatial error model (prediction) and the variables were separately used for removal particular locally systematic and gross errors. Figure 4 shows improvement of the DEM from contour lines with the variable of terrain skeleton (RO).

<figure><img src="image_150.png" /><figcaption>Figure 4: DEM from contour lines improvement with variable RO. The largest errors of underestimated surface (up to -3.5 m – white areas) were found around the sharpest mountain peaks (area of Julian Alps, NW Slovenia, 4.5 x 3 km).</figcaption></figure>

### III DISCUSSION AND CONCLUSIONS

The paper proposes a novel combined approach for spatial errors characterization and prediction with a spatial error model, in order to improve a quality of the digital elevation model (DEM). The method for error characterization with clustering analysis was developed, as well as a method for empirical spatial error model processing based on terrain parameters used to derive homogeneous regions and variables. Both methods use combined solutions with statistical, empirical and visual quality control, in order to describe spatial pattern of error and its magnitude. They are followed by DEM error removal and correction methods for attributes correction or filling gaps with other datasets. Furthermore, four main groups of error that characterize positional errors in DEM were proposed. Besides of the traditional random, systematic, gross, a locally systematic error has been suggested, with its very complex nature that includes geomorphological components. The latest is actually the most typical error in DEM, which is usually too simply considered as a random error.

The planimetric (horizontal) error of the DEM is actually more complicated to assess than of the vertical one. It is also difficult to separate between vertical and planimetric errors, especially on a flat surface. Considerably more challenging is to assess its geomorphological properties. Moreover, in the practice, it is needed a very good quality reference DEM for error characterization and removal, as demonstrated in this paper, which is usually not available. In the case when there is not available any reference data, it can considerably – but not absolutely helps to use a set of independently acquired DEMs, even if they are not perfect quality. Another problem is a complexity of the error field that cannot be easily distinguished from the roughness/rugosity of the terrain, what needs a special attention for the future research.

The resulted error model based on an absolute error is surprisingly high quality, and very nicely reflects the real situation. Besides of the support of error characterization methods, our expert knowledge was also implemented in this model. The solution with spatial error characterization and model considerably better explain DEM error field than the global measures, such as RMSE. An advantage of the proposed empirical spatial error model is that the error pattern has been better understood, controlled and the errors in DEM were removed/corrected.
The future design of the more automated solution needs implementation a mixture of complex experiences and methods. Our next step of research will be a reconstruction of full error field without only absolute values as in this study, applying more comprehensive sensitivity analysis and knowledge-based systems. We will step-by-step apply other terrain parameters, which may influence the quality of the error model. Since the error characters of InSAR or lidar DEMs are considerably different, the terrain parameters and corresponding classifications will consider these options with different scales, and not only the DEM from photogrammetric data sources as in this case. Moreover, our future error assessment methodology will better consider specific demands of particular users.

## Acknowledgements

Part of the study was supported by the Water Science and Technology, and Geotechnics Program [P2-0180] of the Slovenian Research Agency (ARRS). The author is also grateful to A. Vrečko for her contribution in the clustering analysis.

## References

* Carlisle B.H. (2005). Modelling the Spatial Distribution of DEM Error. *Transactions in GIS* 9(4), 521–540.

* Devillers R., Stein A., Bedard Y., Chrisman N., Fisher P., Wenzhong, S. (2010). Thirty years of research on spatial data quality: Achievements, failures, and opportunities. *Transactions in GIS* 14(4), 387–400.

* Erdogan S. (2010). Modelling the spatial distribution of DEM error with geographically weighted regression: An experimental study. *Computers & Geosciences* 36(2), 34–43.

* Fisher P.F., Tate N.J. (2006). Causes and consequences of error in digital elevation models. *Progress in Physical Geography* 30(4), 467–489.

* Guth, P.L. (2006). Geomorphometry from SRTM: Comparison to NED. *Photogrammetric Engineering & Remote Sensing* 72(3), 269–277.

* Hodgson E.M., Bresnahan P. (2004). Accuracy of Airborne Lidar-Derived Elevation: Empirical Assessment and Error Budget. *Photogrammetric Engineering & Remote Sensing* 70(3), 331–339.

* Höhle J., Potuckova M. (2011). *Assessment of the Quality of Digital Terrain Models*. EuroSDR.

* Imhof E. (1982). *Cartographic Relief Presentation*. Berlin and New York: Walter de Gruyter.

* Kimball R, Caserta J. (2004). *The Data Warehouse ETL Toolkit: Practical Techniques for Extracting, Cleaning, Conforming, and Delivering Data*. New York: John Wiley & Sons.

* Leon J.X., Heuvelink G.B.M., Phinn S.R. (2014). Incorporating DEM Uncertainty in Coastal Inundation Mapping. *PLoS ONE* 9(9), 12 p.

* Liu XH., Hu P., Hu H., Sherba J. (2012). Approximation Theory Applied to DEM Vertical Accuracy Assessment. *Transaction in GIS* 16(3), 397–410.

* Oksanen J., Sarjakoski T. (2006). Uncovering the statistical and spatial characteristics of fine toposcale DEM error. *International Journal of Geographical Information Science* 20(4), 345–369.

* Podobnikar T. (2005). Production of integrated digital terrain model from multiple datasets of different quality. *International Journal of Geographical Information Science* 19(1), 69–89.

* Podobnikar T. (2008). Simulation and representation of the positional errors of boundary and interior regions in maps. In *Geospatial vision, new dimensions in cartography*, (Lecture notes in geoinformation and cartography), Springer, pp. 141–169.

* Wise S. (2012). Information entropy as a measure of DEM quality. *Computers & Geosciences* 48, 102–110.

* Zandbergen P.A. (2011). Characterizing the error distribution of lidar elevation data for North Carolina. *International Journal of Remote Sensing* 32(2), 409–430.
<figure><img src="image_151.png" /></figure>

SPATIAL
ACCURACY
2016

5-8 JULY 2016

MONTPELLIER, FRANCE

# Spatial uncertainties
simulation and propagation
<figure><img src="image_152.png" /></figure>
# Exploring the uncertainty of soil water holding capacity information

Linda Lilburne*, Stephen McNeill, Tom Cuthill, Pierre Roudier

Landcare Research, Lincoln, New Zealand

*Corresponding author: lilburnel@landcareresearch.co.nz

## Abstract

Soil water holding capacity is an important soil property for understanding how much irrigation water is required and the quantity of nutrients that are likely to leach into groundwater. This soil profile level property is derived from horizon level data including soil water content, stone content, thickness and horizon type. Soil water content data is expensive to measure so is often estimated in a model that is based on more readily collected soil information. A model of the soil water content (or hydraulic response) has been developed and tested. Its inputs include sand and clay content, profile and horizon classifications. Thus, uncertainties in the derived estimate of soil water holding capacity are due to variability in the inputs to the soil hydraulic model, error in the model itself, and variability of the other key soil properties (stone content, thickness). The combined uncertainty is estimated in a new Soil Profile Simulator tool. It is based on a simulation approach that draws upon the statistical error model of the soil hydraulic model and expert information held in a national soil survey database in the form of probability distributions. This expert information characterises the variability of clay, sand, horizon thickness, stone content, and uncertainty in the classification information associated with the soil polygons. This paper describes the tool, reporting on the progress that has been made in deriving and visualising quantified estimates of uncertainty of soil water holding capacity. Recent advances in technology, including new R packages (aqp, VGAM) and Rserve have been behind this progress.

## Introduction

Information on soil hydraulic properties is essential for the sustainable management of irrigated agricultural land. Providing too much water is wasteful and contributes to contamination of water through leaching and runoff of nutrients. Too little water will impact negatively on yield. One of the key soil properties for managers of irrigated land is the soil's water holding capacity or profile available water (PAW), which represents the amount of water held in the soil that is readily accessible to plants. In New Zealand, this is taken to be the water content held between soil moisture tensions of 10 and 1500 kPa. Although this can be measured directly, it is an expensive process, so PAW is commonly predicted using other soil properties that are easier to measure or observe. These include texture, carbon, bulk density and soil morphology. A relationship between soil properties is known as a pedotransfer function (PTF).

The soil hydraulic response has a number of features that make model predictions difficult. First, the response is bounded (0–100%), so either a bounded-value regression (e.g. logistic) is required or a transformation is needed to an unbounded domain. Second, the response is monotonic with respect to the tension: that is, the response strictly decreases as the tension increases. Finally, the response at a given tension is correlated with the response at other tensions. This complexity also adds to the challenge of estimating the error of the predicted response.
The soil properties used as inputs by the PTF are also uncertain. This paper presents progress with estimating and communicating uncertainty of PAW, whereby the error model of the soil hydraulic PTF is combined with probabilistic information on variability of the key soil inputs of the PTF held in a soil survey database. A new visualisation R package 'aqp' is used to visualise the results.

## Methods

### Soil hydraulic error model

In forming an empirical model for soil hydraulic response, we use soil sample data available from the New Zealand National Soils Database (NSD), which provides the response at tensions of 0 (total porosity), 5, 10, 20, 40, 100, and 1500 kPa. For each sample, texture (sand, silt, clay fractions) data are available, as well as the soil classification, and other factors describing the soil sample.

Our methodology for response prediction uses a vector generalised linear model (VGLM). This, according to Yee (2015), can be thought of as a generalisation of the generalised linear model (GLM) with a vector of responses, which is free from many of the restrictions that the GLM method imposes. We use a vector of responses formed from the logit-transformed 1500 kPa response (transformed to give an unbounded range), as well as the difference between the logit-transformed 100 and 1500 kPa responses, the difference between the logit-transformed 40 and 100 kPa responses, and so on, up to the difference between the logit-transformed 5 and 0 kPa response. The response at one of the specified tensions is formed from the VGLM prediction for 1500 kPa, plus a succession of differences for lower tensions. The uncertainty of each marginal response (at 1500 kPa or a difference in response between tensions) is formed from the estimated VGLM model.

An error model for the difference between the 10 and 1500 kPa responses (i.e. estimates of total available water within a soil horizon or layer), or indeed for any other convenient combination of responses, has been verified using independent validation data (McNeill et al. in prep.). Uncertainty limits, calculated in terms of containment intervals, are estimated by simulation of the aggregated response.

### Soil Profile Generator tool

Lilburne et al. (2012) described how information on soil variability and uncertainty was being incorporated in a national-scale soil database for New Zealand called S-map (Landcare Research, 2015). The very limited amount of soil sample data meant that an expert knowledge approach was used to record information on the confidence of classification and base property attributes, the variation of soils in a polygon and their proportional reliability, and the range of values of key soil properties (Lilburne et al. 2008). Variability of quantitative soil properties is stored in the form of probability distribution functions (pdf). The lack of point soil sample data precluded a geostatistical approach to modelling and simulating PAW.

A new Profile Simulator tool has been developed that creates realisations of profiles based on the expert-derived information on soil variability and uncertainty. The key information used in this study on PAW is the PTF error model; expected variability of the stone, sand and clay content; type and thickness of functional horizons within a soil survey polygon; confidence in the Soil Order classification, rock type of the fine material and drainage class, as well as the reliability of the proportions of soil types within a polygon. Values for horizon stone, sand, and clay content,
and thickness were drawn from their respective pdfs. Each profile realisation is checked against a set of rules to ensure that it is still consistent with the soil definition. For example, the sum of the simulated horizon thicknesses must fit the pdf of the soil's depth. Realisations that do not fit the rules are discarded and regenerated. The PTF error for each soil profile realisation was simulated and the PAW calculated.

The tool is based on an architecture that links the S-map database in SQL Server with RServe – a server that responds to requests from clients by running a R script. Rserve enables R calculations to be performed on request without the need to start up an R session each time. A number of servlets have been developed that retrieve the stored uncertainty information from SQL Server, call a set of Monte Carlo functions to generate realisations of the key parameters required by the soil hydraulics error model, run this model on each set of parameter realisations, simulate model error, and finally, return the set of profile realisations with estimates of available water, in a range of useful formats (database tables, csv file, Rdata file). RServe was installed on a Linux-based server to allow for multi-threaded processing.

## Visualisation

The 'aqp' package (Beaudette et al. 2013) is an R package designed for working with soil information. It contains a SoilProfileCollection class to simplify the process of working with the collection of data associated with soil profiles. It also includes tools for plotting soil properties by depth and their associated variability.

## Results

Using a training dataset of 1007 points, a VLGM model was developed and tested on an independent dataset of 432 points (Figure 1). Histograms of the residuals for seven tension values, and the associated model diagnostics using training and validation data, showed similar patterns, indicating that the model was not overfitted. Simulated error for the 1500 kPa tension value and total available water for an example horizon is shown in Figure 2. Response distributions are formed by simulation, and can be processed to give containment intervals (e.g. 95%). While most response distributions appear symmetric, the response can in some cases be highly skewed where there is little training data (rare soil classes). The simulated water retention curves for the 500 realisations of a Barr_6a sibling (a moderately deep sandy loam) are shown in Figure 3.

The Soil Profile Simulator is fast, generating 10 000 realisations of a specified soil sibling as a Rdata file in 28 seconds. The tool can output the realisations as tabular data or as a SoilProfileCollection (the class used in aqp). Use of the R aqp package facilitates drawing graphs of the soil profile. Figure 4 compares the 500 realisations of estimated PAW with the PAW that is listed on the S-map fact sheets (123.17 mm). The fact sheet value is estimated using the mean soil hydraulic response with mean values for clay and sand content, along with mean horizon thickness, stone content, and modal horizon classification. The mean PAW of the 500 realisations is 122.96 mm with a standard error of 1.2356 and a standard deviation of 27.62. Figure 5 shows the 500 soil profile realisations in terms of stone, sand and clay content down the profile. The variability of PAW down the profile is shown in Figure 6, where the darker grey area is the 25% to 75% quantile, and the dashed lines show the 95% confidence interval.
<figure><img src="image_153.png" /><figcaption>Figure 1: Measured-versus-fitted plots of the soil hydrological response for seven tension values. The vertical lines are plus and minus one standard error of the prediction uncertainty.</figcaption></figure>

<figure><img src="image_154.png" /><figcaption>Figure 2: Histograms showing the simulated error (1500 kPa and total available water (TAW)) for an Allophanic soil with a loamy functional horizon (sand 15%, silt 57.5%, clay 27.5%).</figcaption></figure>
<figure><img src="image_155.png" /><figcaption>Figure 3: Simulations of the water retention curve for each of the four horizons of the Barr_6a sibling (a moderately deep sandy loam) as simulated by the soil hydraulic PTF error model.</figcaption></figure>
<figure><img src="image_156.png" /></figure>

Figure 4: Histogram of the profile available water (PAW) of 500 realisations of the Barr_6a sibling. The red line shows estimated PAW where the mean values are used for the soil hydraulic response; clay, sand, stone content; and horizon thickness. The circles along the x axis indicate the distribution of the PAW estimates.

<figure><img src="image_157.png" /></figure>

Figure 5: Soil profile realisations (n = 500) from the Soil Profile Simulator showing the variability of stone, sand and clay content down the soil profile of the Barr_6a sibling.
<figure><img src="image_158.png" /><figcaption>Figure 6: The variability of profile available water (PAW) down the Barr_6a soil sibling profile. The blue line is the median PAW from the 500 realisations, the solid and dashed grey lines indicate the 5% and 25% quantiles, respectively.</figcaption></figure>

The speed of the Soil Profile Simulator tool also allows many siblings to be simulated according to the uncertainty of the proportions of each soil type within the soil survey polygons. Thus map realisations of PAW can be generated.

## Discussion

### Spatial uncertainty

The use of an expert approach to recording information about uncertainty within S-map is currently limited by the lack of any quantitative information about the spatial accuracy of the polygon line work. In the future, as we build in knowledge of soil-landform relationships, either explicitly or via digital soil mapping approaches, we will be able to represent this spatial uncertainty. For example, in a case where a soil type is linked to concave slopes and gullies, the weaker the membership of a DEM pixel in the concave-slope/gully class, the more uncertain the association with the soil class. This uncertainty can be added to the simulation tool. An alternative approach could be to use an expert approach to simplistically quantify the possible offset of each polygon line boundary.

### Correlation

Another current limitation is the lack of information on both correlation between soil properties and spatial autocorrelation. The former was accounted for by using an approach whereby simulated profiles that did not meet known constraints were discarded. For example, S-map
contains pdfs of horizon thickness. Some horizons will be inversely correlated – if one horizon is thicker than the mean, then at least one other horizon must be thinner in order for the overall depth (up to 1 m) to be maintained. If the sum of the simulated thicknesses was not consistent with the range of the depth pdf, then that realisation was discarded and another regenerated. The tool stops after one thousand attempts to generate a valid realisation. However, lack of information on the spatial autocorrelation of soil properties means that this aspect cannot be addressed at this point.

## Conclusion

Conventional soil survey databases do not lend themselves to exploring the impact of uncertainty in soil information. However, progress has been made in deriving and visualising quantified estimates of uncertainty in key soil properties related to the soil water holding capacity. This has been enabled by recent advances in technology, including new R packages (aqp, VGAM) and Rserve.

## Acknowledgements

This work was funded by Landcare Research Core Funding. We thank Stella Beliss and Varvara Vetrova who kindly reviewed the manuscript, and Leah Kearns for editing it.

## References

* Beaudette, D.E., Roudier, P., O'Geen, A.T. (2013). Algorithms for quantitative pedology: A toolkit for soil scientists. *Computers & Geosciences* 52, 258–268.

* S-map - New Zealand's national soil layer. DOI: http://dx.doi.org/10.7931/L1WC7. Accessed: 2016-05-16.

* Lilburne, L., Hewitt, A., Ferriss, S. (2006). Progress with the design of a soil uncertainty database, and associated tools for simulating spatial realisations of soil properties. In: M. Caetano, M. Painho (Eds.), 7th International Symposium on Spatial Accuracy Assessment in Natural Resources and Environmental Sciences. Instituto Geografico Portugues, Lisbon, Portugal, pp. 510-519.

* Lilburne, L.R., Hewitt, A., Webb, T. (2012). Soil and informatics science combine to develop S-map: a new generation soil information system for New Zealand. *Geoderma* 170, 232–238.

* McNeill, S., Lilburne, L., Webb, T., Cuthill, T. (in prep.) Pedotransfer functions for hydrological properties of New Zealand soils using S-map information.

* Yee, T.H. (2015). *Vector generalised linear and additive models*. Springer.
# Simulation of realistic digitizing errors on geographical objects using constraints modeling the operator input process

Jean-François Girres*¹

¹ Université Paul-Valéry Montpellier 3 – UMR GRED, France

*Corresponding author: jean-francois.girres@univ-montp3.fr

The geometry of vector objects in a geographical database is mainly obtained through manual input processes performed by an operator. It is well known that this capture process can cause errors - called digitizing errors - that affect the positioning of geographical objects or the computation of geometric measurements (i.e. length and area). To model these impacts, simulation methods, based on Monte-Carlo approaches, are generally proposed. However, several problems can occur by simulating digitizing errors with simple random processes (e.g. topological errors or non-respect of the object's shapes). Indeed, these problems would have been avoided by an operator performing a classical manual input process. Thus, in order to simulate digitizing error with a higher degree of realism, it seems important to understand and model some operator's behaviors. In this context, this paper proposes an approach to simulate digitizing errors using a set of constraints that try to reproduce the capture mechanisms of the operator. To contextualize these issues, principles and limitations of digitizing error simulation methods are presented in section 1, before detailing in section 2 the proposed constraints to model the operator input process. Finally, an experiment of the proposed methods is performed in section 3, before concluding.

## I SIMULATION OF DIGITIZING ERROR : PRINCIPLES AND LIMITATIONS

The geometry of geographical vector objects is mainly produced through manual input process performed by an operator. This human capture process can generate accidental and systematic errors, called digitizing errors, that affect the positioning of the vertices of the geometry. These errors impact the positioning of objects, the computation of geometric measurements (length or area) or any analysis performed with the geometries of vector objects (e.g. buffers, overlaps). Thus, to model the positional uncertainty involved by digitizing error and assess its impact, several contributions (Keefer et al, 1988; Hunter and Goodchild, 1996) propose to use simulation methods, based on Monte-Carlo approaches. Using these approaches, simulated geometries are created by generating random errors in the positioning of each vertex of the geometry of the original objects. As proposed by Bolstad et al (1990), a normal law can be used to model the imprecision in the positioning of the vertices.

The simulation of errors on a geometry can be seen as the translation of the coordinates x and y of each vertex. As exposed on figure 1, two methods can be used to translate the vertices in a simulated geometry:

- in the first method (on the left), the coordinates of a vertex P are translated according to an angle $a$ (generated randomly), at a distance $d$ from the original vertex position ($d$ is generated using a normal law)

- in the second method (on the right), the coordinates x and y of a vertex P are translated at distances $dx$ and $dy$ from the original vertex position
<figure><img src="image_159.png" /><figcaption>Figure 1: Two methods to translate the coordinates x and y of a vertex P</figcaption></figure>

The first method, which is considered as more correct to represent accidental errors (Vauglin, 1997) is used in this study to simulate digitizing errors on a geometry. This operation is performed on each vertex of the original geometry in order to generate a simulated geometry. By repeating this operation on a large set of realizations (figure 2), it is then possible to study the sensibility of geometric measurements on simulated geometries (by comparing lengths or areas with the original geometry) and then assess the impact of digitizing error on measurements, as proposed by Goodchild et al (1999).

<figure><img src="image_160.png" /><figcaption>Figure 2: An original geometry (left) and 10 simulated geometries (right)</figcaption></figure>

However, several problems, which would have been avoided by a “human” operator, can occur by simulating digitizing errors with simple random processes. The first type of problem is related to topological errors (e.g. self-intersection of edges) as exposed in figure 3. Several contributions have already been proposed to avoid such topological inconsistencies (Hunter and Goodchild, 1996; Bonin, 2002).

<figure><img src="image_161.png" /><figcaption>Figure 3: Topological inconsistencies on a simulated geometry (in blue)</figcaption></figure>
The second type of problem deals with the non-respect of the original shape of the object after a simulation. For instance, as exposed in figure 4, the simulation of random digitizing errors on an original polyline with a straight shape (in red) can generate a distorted simulated geometry (in blue). As a consequence, the simulated polyline loses its original straight shape, and its length automatically increases (because a straight line minimizes the length). Thus, this simulation can be considered as unrealistic, because a human operator would have preserve the straight shape of the polyline.

<figure><img src="image_162.png" /><figcaption>Figure 4: Non-respect of the original shape of a straight line using random simulations</figcaption></figure>

These two examples illustrate a simple fact : the use of simple random models is not adapted to simulate human input processes. Indeed, the operator input process involves mobilizing cognitive mechanisms that are complex to model. So, in order to simulate digitizing errors with a higher degree of realism (and assess their impact on measurements), some operator's behaviors need to be understood and modeled in the simulation process. The next section of this paper will present propositions in order to model the operator input process.

## II SIMULATION CONSTRAINTS TO MODEL THE OPERATOR INPUT PROCESS

In order to simulate digitizing errors with more realism, this section presents a set of constraints that try to reproduce the capture mechanisms of the operator. Based on basic assumptions about the data entry process, three types of constraints are integrated in the digitizing error simulation model in order to preserve the realism of the generated objects : the weighting of errors on vertices according to their angularity (1), the weighting of errors on nodes according to their degree (2) and the correlation of errors on vertices using progressive translation from the extremities (3).

### 2.1 Weighting of errors on vertices according to their angularity

The first type of constraint on the simulation process deals with the weighting of errors on vertices according to their angularity. This constraint is integrated in order to avoid the distortion of straight lines, as exposed in figure 4, that a human operator would have preserved.

The basic underlying assumption is that the angle between successive vertices of a geometry is a representation of the reality determined with the consciousness of the operator. So, if successive vertices of a geometry are aligned, a realistic error simulation should preserve this alignment. On the other hand, if brutal turns in the direction between successive vertices are observed, the risk of digitizing error might be more important. This basic assumption involves to take into account the angularity between successive vertices as a constraint on the simulation of digitizing error.

To model this constraint, a weighting $q$ (between 0 and 1) of the distance $d$ of the error on each vertex $P$ of the geometry is determined according to the angle between successive vertices $P-1$, $P$ and $P+1$. For instance, $A$, $B$ and $C$ are three successive vertices. To determine the weighting $q_B$ of the error on the vertex $B$, its angularity $\alpha_B$ (in radian) is computed using equation 1.

$$ \alpha_B = \frac{|\widehat{ABC} - \pi|}{\pi} \qquad (1) $$
Once the angularity of each vertex is determined, a function can be used in order to determine the weighting $q$ according to the angularity $\alpha$. Indeed, based on our assumption, the more the angle is close to $\pi$, the more the error should be small. But it does not mean that the error is null (which would be the case with $\alpha = 0$). Moreover, it is extremely rare to get angles between successive vertices where the angularity $\alpha_B$ is close to 0 or $2\pi$.

As a consequence a function need to be defined in order to determine the weighting $q$ of the error according to the angularity $\alpha$. This function can take basic forms (e.g. $q=0.75\alpha+0.25$; $q=\alpha^{0.5}$) or more complex formulations, but the goal of this paper is not to define the best function to determine the appropriate weighting according to a given angularity. In the experiments (section 3), a basic form will be used.

Once the function is defined, the weighting of digitizing errors can be applied on each vertex of the geometry using a normal law $N(0, \sigma)$. As exposed on figure 4, which represents the simulation of errors using a weighting according to the angularity, the more the angularity is important, the more the error on the vertex is important.

<figure><img src="image_163.png" /><figcaption>Figure 4: Simulation of errors weighted by the angularity between successive vertices</figcaption></figure>

The example of simulation exposed on figure 4 shows that no errors is applied on the extremities of the polyline, where it is not possible to determine any angularity value. To determine the error applied on the extremities of a polyline, a second method, based on the degrees of the extremity nodes, is exposed in the following subsection.

## 2.2 Weighting of errors on extremity nodes according to their degree

Because it is not possible to apply a weighting based on the angularity between successive vertices on the extremities of a polyline, a method based on the degree of extremity nodes is proposed.

The underlying assumption formulated to determine this constraint is that the more a node is connected with other edges, the easier its localization should be for the operator, and the more precise its positioning should be. On the other hand, if an extremity is not connected with any other edge, we can expect that it is more difficult for the operator to locate the node precisely. As a consequence, the positioning of unconnected extremity nodes should be more imprecise, and the digitizing error should be more important.

Following this assumption and in order to simulate digitizing errors with more realism, we consider that the error applied on the extremities of a geometry have to be weighted according to the degree of the edges. Thus, for a node $N$ with a degree $d$ located at the extremity of a polyline, we propose to apply a weighting $q_N$ as defined in equation 2, using a constant $a$.

$$q_N = \frac{a}{d} \qquad (2)$$
As mentioned previously, the goal of this paper is not to define the most appropriate value of the constant *a* used to compute the weighting of errors on extremities. In this paper, we will use a value of *a*=1. Nevertheless, some other parameterizations of this constant could be experimented in further researches.

The weighting $q_N$ is then applied on the simulated errors on the initial and final nodes of the geometry using a normal law $N(0, \sigma)$. Figure 5 illustrates the application of the weighting of errors on extremity nodes according to their degrees, associated with the weighting of errors on vertices according to their angularity. It shows that the error applied on unconnected nodes is more important than on connected nodes.

<figure><img src="image_164.png" /><figcaption>Figure 5: Determination of the degrees of nodes (left) for the weighting of errors on extremities (right)</figcaption></figure>

This second constraint shows how to apply a weighting for the simulation of digitizing error on the extremities of a geometry according to its configuration. Finally, in order to take into account the dynamic dimension of the data capture process between the initial and the final node, a last constraint is proposed : the correlation of errors between successive geometries.

## 2.3 Correlation of errors between successive vertices

Several contributions (Heuvelink et al., 2007; De Bruin, 2008) consider that the simulation of errors on a geometry should integrate a correlation between successive vertices. It is admitted that an operator can generate a systematic error (i.e. a bias) during the input process of vector data. Nevertheless, it is difficult to simulate this error, because it depends of the operator and the way he works. But some other reasons (e.g. the dynamic dimension of the data entry process) can justify to integrate a correlation of errors between successive vertices.

To explain the integration of the correlation of errors between successive vertices, we can use the example of an operator capturing a road. If an error occurs at the extremity of the road (i.e. at the crossroad), we can assume that the operator will try to cushion progressively this error on the entire polyline, instead of correcting it brutally on a single vertex, what would generate a non-respect of the shape of the road. Indeed, the input process of vector data is a dynamic (and active) process. Thus, if the operator generates an error, he will try to correct the error without reducing his productivity. Following this example, we can consider that the error generated on the extremities (proposed in the previous subsection) is propagated progressively on the entire geometry, which supposes a correlation of errors between successive vertices.

To model the correlation of errors between successive vertices, we propose a method that generates a translation of each vertex, based on the errors on the extremities.

Using the method proposed in the precedent subsection, errors a generated at the extremities of a geometry according to the degree of the nodes, and following a normal law $N(0, \sigma)$. As a consequence a translation using vectors $\vec{v}$ and $\vec{u}$ is performed on the extremities A and C, and a translation is applied on an intermediate vertex B with a vector $\vec{w}$, computed with a weighted sum of vectors $\vec{v}$ and $\vec{u}$.
The weighting coefficients *a* and *b* are defined as the ratio between the length from the intermediary vertex to each extremity node with the total length of the polyline, as exposed on equation 3.

$$ \vec{w} = a * \vec{v} + b * \vec{u} \qquad (3) $$

where *a* is the ratio between the length of [AB] and the length of [AC] and *b* is the ratio between the length of [BC] and the length of [AC]

<figure><img src="image_165.png" /><figcaption>Figure 6: Translation of a vertex *B* according to the errors on the extremities *A* and *C*</figcaption></figure>

It is important to notice that in figure 6, the error on intermediate vertices (using the angularity between successive vertices, as proposed on subsection 2.1) is not applied.

## 2.4 Combination of constraints for digitizing error simulation

The combination of the three constraints varies according to the type of geometries and their characteristics. For polylines, the three different constraints proposed to generate realistic digitizing error simulations are combined for each realization of the simulation as follows: (1) error simulation on the extremity nodes according to their degree, (2) translation of intermediary vertices according to extremity nodes errors (i.e. correlation of errors), (3) error simulation on translated vertices according to their angularity.

For polygons, it is necessary to take into account the characteristics of the represented objects. We propose to differentiate the objects representing a partition of the real world (e.g. administrative units) and objects that don't represent a partition of the real world (e.g. water surfaces). For polygon objects representing a partition of the real world, the same methodology as the one proposed for polylines is applied. Indeed, these objects suppose the existence of common borders. We can then suppose that a correlation between successive vertices exist, because the operator will start to capture a border from an initial node (the border between three units) to a final node. As a consequence, he will cushion the error from the extremities on the entire border. Thus, these polygon objects can be simplified as a graph.

For polygon objects which don't represent a partition of the real world, only the method of error simulation on vertices according to their angularity is applied. Indeed, because these geometries are closed and isolated, the notion of initial and final nodes doesn't exist. As a consequence, the method proposed to generate errors on extremities (according to the degree of the nodes) and translation between successive vertices (according to errors on extremities) don't need to be applied.

The different methods proposed in this section are experimented in the following section on real objects.
### III EXPERIMENTS

To illustrate the functioning of the methods proposed to simulate digitizing errors, an experiment is proposed on a dataset of administrative units, extracted from the French BDCARTO database in the area of Handaye and Saint-Jean-de-Luz (south-west of France). The total area of these administrative units is about 133.43 sq. km.

<figure><img src="image_166.png" /><figcaption>Figure 7: The administrative units used for the experiment (left) and their conversion in a topological graph (right, with edges in red and nodes in black)</figcaption></figure>

To simulate digitizing errors, a normal law $N(0, \sigma)$ is used, with a value $\sigma = 4.97$ m. (determined using the automatic capture scale estimation method, see Girres, 2015 and Girres, 2012) and a function $q=0.75\alpha+0.25$ for the weighting of errors according to the angularity between successive vertices. Finally, 1000 realizations of the simulation are realized. To perform the simulation, administrative units, which are considered as a partition of the real world, are converted in a topological graph (figure 7). This method allows to extract nodes and to process each border between units as a polyline during the simulation. Finally, after the simulation, polygon units are re-built using the graph faces.

<figure><img src="image_167.png" /><figcaption>Figure 8: Distribution of area differences between original and simulated administrative units</figcaption></figure>

Results of the experiments show that the distribution of area differences between the original dataset and the simulated datasets presents a value of $\sigma_A = 0.008$ sq. km. ($= 0.8$ ha). Assuming that this distribution follows a normal law, we can consider that the area imprecision is about +/- $3\sigma_A$ (with a confidence of 99%) which means an area imprecision of +/- 2.4 ha (0.018% of the original area).
As a comparison, the same experiment performed with no constraints on the simulation of digitizing errors would have generated an imprecision of about +/- 4.8 ha (or 0.036% of the original area), which means an area imprecision multiplied per two. These results show that the use of a simple random digitizing error simulation process can generate important overestimations of the measurement imprecision, and that the integration of constraints in the simulation process allows to limit such exaggerations.

## CONCLUSION

As a conclusion, this paper presents original methods in order to simulate digitizing errors with a higher degree of realism, by integrating constraints modeling the operator input process. Observations show that classical random processes used to simulate digitizing errors can generate topological inconsistencies and non-respect of the original shape of objects, and these problems would have been avoided by a human operator. Based on a set of three assumptions about the operator input process, three methods are then proposed to integrate constraints in the simulation process. Results of the experiments show that the use of these constraints reduce the measurement imprecision generated by a classical random process, in addition with the suppression of potential topological inconstancies. However, several improvements need to be realized on this method, especially the question of the parameterization of error weightings (according to the angularity of successive vertices or the nodes' degrees). Further researches should be investigated in order to provide the most appropriate ways to parametrize the different constraints. In spite of these implementation problems, this approach proposes new materials in order to model digitizing error with more realism, and assess its impact on the positioning of geographical objects or on geometric measurements.

## References

* Bolstad, P. V., Gessler, P. and Lillesand, T. M. (1990). Positional uncertainty in manually digitized map data. *International Journal of Geographical Information Systems* 4(4), 399-412.

* Bonin, O. (2002). Modèles d'erreurs dans une base de données géographiques et grandes déviations pour des sommes pondérées; application à l'estimation d'erreurs sur un temps de parcours. PhD Thesis, Paris 6 University, France.

* De Bruin, S. (2008). Modelling positional uncertainty of line features by accounting for stochastic deviations from straight line segments. *Transactions in GIS* 12(2), 165-177

* Girres, J.-F. (2012). Modèle d'estimation de l'imprécision des mesures géométriques de données géographiques. *Application aux mesures de longueur et de surface*, PhD Thesis, Paris-Est University, France.

* Girres, J.-F. (2015). Estimation of geographical databases capture scale based on inter-vertices distances exploration, *ISPRS Annals of Photogrammetry, Remote Sensing and Spatial Information Sciences* 2(3), 305-310.

* Goodchild, M., Shortridge, A. and Fohl P. (1999). Encapsulating simulation models with geospatial data sets. In Lowell, K. and Jaton, A. (eds) : *Spatial Accuracy Assessment : Land Information Uncertainty in Natural Resources*, Ann Arbor Press, 123-129.

* Heuvelink, G. B. M., Brown, J. D. et van Loon, E. E. (2007). A probabilistic framework for representing and simulating uncertain environmental variables. *International Journal of Geographical Systems* 21, 497-513

* Hunter, G. J. et Goodchild, M. F. (1996). A new model for handling vector data uncertainty in geographic information systems. *URISA Journal* 8(1), 51-57.

* Keefer, B., Smith, J. and Gregoire, T. (1988). Simulating manual digitizing error with statistical models. In *Proceedings of GIS/LIS'88 Conference*, 475-483.

* Vauglin, F. (1997). Modèles statistiques des imprécisions géométriques des objets géographiques linéaires. PhD Thesis, Marne-La-Vallée University, France.
# How far spatial accuracy governs land-use changes monitoring frequency: the urban sprawl monitoring example

Jean-Pierre Chéry*<fnref n="1" />, Jean-Stéphane Bailly<fnref n="2" />, Valérie Laurent³, Nathalie Saint-Geours<fnref n="4" />

*Corresponding author: jean-pierre.chery@agroparistech.fr

## Abstract

In this paper, we illustrate how far spatial accuracy of a land-use map governs land-use changes monitoring frequency on a urban sprawl monitoring case study. From a specific Monte Carlo approach propagating uncertainties, confidence curves for minimal monitoring frequency to detect significant changes in urban sprawl indicators were built. Results showed that frequency decreased when uspscaling indicators but it also showed very low monitoring frequency for indicators at the lower level.

## INTRODUCTION

When setting up land-use monitoring systems, spatial uncertainties and their impact on the detection of indicator changes are usually ignored. To capture a significant change in land-use indicators is thus strongly related to its spatial resolution, the velocities of the process it represents and the accuracy of the used indicators. As a consequence, the required monitoring land-use change frequency, corresponding to the mimimum time step to ensure a significant change in indicators, also depends on these three factors: indicator spatial resolution, change process velocity and spatial indicator accuracy.

In this paper, we illustrate how far spatial accuracy of a map governs land-use changes monitor-
ing frequency focusing on urban sprawl monitoring system from satellite imagery in southern
France. Urban sprawl is a major challenge for land use planners: it causes the sealing of lands
closest to the urban centers into impervious areas, thereby transforming highly productive cul-
tivated soils, increasing flood hazards, fragmenting natural habitats, and raising complex issues
related to transportation and social diversity. To define new urban policies, integrated moni-
toring of urban sprawl is crucial. It is therefore important to choose an appropriate monitoring
frequency (Allen and Lu, 2003). Such monitoring relies on urban sprawl indicators which suffer
from uncertainties related to impervious area mapping process or from the population census
methodology. In most urban sprawl monitoring studies, however, these uncertainties effects
are not considered, and the monitoring frequency is chosen solely based on data availability.
In practice, time lag between successive indicator calculations may range from 5 years to 25
years, or even be irregular. As for monitoring biological diversity (Yoccoz et al., 2001), having
a carefully designed urban sprawl monitoring system in terms of choice of indicator, scale and
temporal frequency is essential.

To support land-use planning, this paper proposes a framework to infer indicator monitoring
frequencies, taking into account indicator uncertainties and scaling. This framework was ap-

<footnote n="1">¹AgroParisTech, TETIS, France</footnote>

<footnote n="2">²AgroParisTech, LISAH, France</footnote>

<footnote n="2">²Irstea, TETIS, France</footnote>

<footnote n="4">⁴ITK, France</footnote>
plied to 3 urban sprawl indicators at municipality, inter-municipality, department, region scales on the former Languedoc-Roussillon region, France.

# MATERIAL AND METHODS

## Study area

Languedoc-Roussillon is a former region of mainland France that covers 27,376 km² along the Mediterranean coast, and hosts about than 2.7 million inhabitants. Since the 1960s, intense population pressure in the region combined with poor urban planning led to population concentration in the major cities and their peripheral extensions in the coastal zone. Languedoc-Roussillon region contains 1,545 municipalities aggregated in 130 inter-municipalities which facilitate cooperation on practical urban management issues, and 5 departments. Therefore, four increasing spatial supports of urban sprawl indicator were used in the study: municipalities, inter-municipalities, departments, and entire region.

## Data

Dataset of impervious polygons covers the entire study area were obtained for 1997 and 2009 from the Dupuy et al. (2012) methodology, using satelitte image (RapidEye imagery for 2009) classification with numerous additive manual post-classification at 1:8,000 scale to reduce confusions.

Population data for along years were taken from censuses of the French Institute of Statistics and Economic Studies (INSEE), which was performed at municipality scale. INSEE also provides uncertainty measures for their data (INSEE, 2012), as a coefficient of variation (CV) along 5 municipality categories (Tab. 1). For small municipalities having less than 10,000 inhabitants, there was exhaustive sampling and therefore, no uncertainty.

<table><caption>Table 1: Median of the coefficient of variation of the number of inhabitants per municipality category (INSEE, 2012), and number of municipalities of the study area in each category</caption><thead><tr><th></th><th>CV median</th><th>Number of municipalities less than 10,000</th></tr></thead><tbody><tr><td>0 to 10,000 inhabitants</td><td>0</td><td>1517</td></tr><tr><td>10,000 to 19,999 inhabitants</td><td>1.02</td><td>17</td></tr><tr><td>20,000 to 49,999 inhabitants</td><td>0.78</td><td>6</td></tr><tr><td>50,000 to 99,999 inhabitants</td><td>0.56</td><td>2</td></tr><tr><td>More than 100,000 inhabitants</td><td>0.39</td><td>3</td></tr></tbody></table>

## Urban sprawl spatial indicators

Urban areas usually consist in namely Morphological Urban Area (MUA) obtained after morphological closing of a set of delineated impervious polygons. Closing consists in successive dilation and erosion (with a given radius *r*), allowing to merge impervious polygons. Based on the MUA, Balestrat et al. (2010) proposed a collection of spatial indicators of urban sprawl. For this study, we selected three indicators from this collection: the MUA area *A* (area-based), the dispersion coefficient *D* (shape-based) corresponding to the ratio between the cumulated areas of MUA polygons whose areas are smaller, respectively bigger than 3 hectares, and the population density *P* (mixing spatial and population data, i.e. number of inhabitants). Maps of these three indicators may be produced for each of the four spatial supports (municipality, inter-municipality, department, and region) (Laurent et al., 2014b).
## Urban sprawl indicators uncertainty simulation

Three uncertainty sources in the computation of urban sprawl indicators *A*, *D* and *P* were considered: i) the impervious polygons, ii) the inhabitant number per municipality, and iii) the radius *r* for dilatation-erosion used to build the MUA. Geometric uncertainties for both 1997 and 2009 impervious polygon layers were explicitely simulated using the approach developped Laurent et al. (2014b). Prior to the simulation, metrics on impervious polygons geometry errors were first computed based on the comparison of image-processed dataset with a reference dataset on 75 sampled municipalities. These metrics quantify both the shape of impervious polygons using the position error of polygon vertices (specific spatial covariance models) and the rates of omission and commission errors of the impervious polygons within municipality areas. Based on these metrics, random realizations of impervious polygons datasets were simulated for each municipality. First, the position of polygon vertices was randomly perturbed to account for geometric uncertainties. Then, small impervious polygons were randomly added or removed to account for polygon omission/commission. For more details, the reader is referred to Laurent et al. (2014a). Figure 1 presents a few examples of simulated impervious polygon simulations and resulting MUA.

<figure><img src="image_168.png" /><figcaption>Figure 1: Original impervious polygons from the automatic dataset and three examples of impervious polygons simulations for the municipality of Jonquieres</figcaption></figure>

For each municipality, the population was assumed to follow a Gaussian distribution $N(p, p * CV)$, where *p* is the municipalitys nominal population, and *CV* is the population coefficient of variation taken from Table1. Random realizations of population values were then drawn from $N(p, p * CV)$.

All sources of uncertainties were propagated through the urban sprawl indicator computation chain in a Monte Carlo framework. For each value of the radius *r*, a Monte Carlo experiment with N = 1,000 simulations was carried out while keeping total computational cost tractable.

## Infering the right monitoring frequency under uncertainties

The "right" monitoring frequency is defined here as the minimal time lag $\Delta_t$ between two successive indicator calculations required to detect a significant change in the indicator. It depends both on indicator uncertainty and on the magnitude of the indicator change over time. Calculating $\Delta_t$ is thus considered as a problem of significant differences between two independent random variables. Let $I_0$ and $I_1$ be the indicator random variables at times $t_0$ and $t_1$ respectively, following independent Gaussian distributions $N(i_0, \sigma)$ and $N(i_1, \sigma)$. $i_0$ and $i_1$ are the indicator mean values (nominal values) at $t_0$ and $t_1$ and $\sigma$ is the indicator standard deviation, assumed to be constant over time and infered from previous section 2.4. Assuming both unbiased indicator
means at a given time and a linear indicator trend with slope s, the random variable $\Delta_T = \frac{I_1-I_0}{s}$, can be inferred from the lowest significant difference ($I_1 - I_0$) at confidence level $1 - \alpha$ as:

$$ \Delta_t = \frac{q2\sigma}{s} \qquad (1) $$

where q denotes the quantile of the standard normal distribution. We can then build confidence curves for $\Delta_t$, as a function of s, for the usual confidence levels 90, 95, and 99% for instance. To choose the range of s values for the Languedoc-Roussillon region ($s_{LR}$), we calculated the average indicator slope over the study area between 1997 and 2009, assuming similar geometric and thematic uncertainties within the 1997 and 2009 impervious polygons datasets.

# RESULTS

## Urban sprawl uncertainties simulation

The nominal data (no uncertainty) and the simulation results at the municipality scale with 50 m radius r are presented in maps on figure 2. The nominal indicator maps (line 1 - Fig. 2) show the spatial patterns of the indicators. A has higher values along the coast and other plains than in the rugged inland areas. Conversely, D (dispersion coefficient) shows an opposite pattern, with higher values in the rugged inland areas, and lower values in the plains. The P (population density) map shows mostly low values with few municipalities located mostly in rugged areas having higher values. Comparing the simulated indicator maps to the nominal maps, one notices that, for all three indicators, the simulated spatial patterns is similar to the nominal ones. This good match can also be seen on the scatterplots, where most points are located close to the y = x line. For A, however, the simulated values were much higher than nominal values for municipalities having A smaller than 50,000 m2. This threshold of 50,000 m2 corresponds to a rectangle MUA of 500 * 100 m, which would be an unrealistically small village. We therefore infer that the nominal map underestimates A for those municipalities. This underestimation of A by the nominal maps may be due either to the presence of clouds in the satellite images, or to the omission of (parts of) polygons in the image classification. Both D and P present a general trend to overestimation. The CV (coefficient of variation) maps contain the uncertainty information: the higher the CV value, the higher the uncertainty. The three CV maps show similar spatial trends, with lower uncertainty along the coast and in other plains, and higher uncertainties in the rugged areas. The median CV value increases from A (25%), to P (31%), to D (58%), meaning that uncertainty increases from A, to P, to D. The maximum CV values were very high (2.09 for A, 6.66 for D, and 14.57 for P), indicating locally very high uncertainties for some municipalities.

For convenience, results for other radius values are not detailly shown. Indeed, for large radius r values, more polygons are joined into the MUA, increasing A and, thereby decreasing D and P. For all three indicators, the median CV values slightly increase with increasing radius, indicating that increasing the radius slightly increases the indicator uncertainty.

The impact of indicators spatial upscaling is shown on the boxplots of figure 3. As expected, A increases with increasing entity size, because the MUA areas are being summed over larger entities. On the contrary, D decreases with increasing entity size. The entity size does not impact the P median. As expected (Heuvelink, 1998; Saint-Geours et al., 2014), for all three variables, uncertainty decreases (decreasing CV values) when upscaling. For A and D, the CV
<figure><img src="image_169.png" /><figcaption>Figure 2: Results of Monte Carlo simulations, compared to the nominal results at municipality scale. The three columns present the results for the *A*, *D*, and *P* indicators. Line 1 shows the nominal maps, lines 2 and 3 the maps of the mean and *CV* respectively over the 1,000 simulations, and line 4 the scatterplots of the simulated versus nominal indicator values</figcaption></figure>

values are divided by a factor 3 when upscaling from municipality to inter-municipality, by a factor 7 from inter-municipality to department, and by a factor 2 from department to region. For *P*, the factors are smaller: 3, 1.4, and 1.7, in the same upscaling order as for *A* and *D*.

## Right monitoring frequencies

The confidence curves for the minimal monitoring time lag $\Delta_t$, as a function of the slope *s* of the urban sprawl indicator trends, are presented in figure 4. Looking at the *A* plot at municipality scale, we can see that, for all three confidence levels, the curves present a hyperbolic decreasing
<figure><img src="image_170.png" /><figcaption>Figure 3: Influence of the aggregation scale on the indicator mean values and coefficients of variation (CV). The sample sizes are: 1,545 for the municipality scale (mun), 130 for the inter-municipality scale (int), 5 for the department scale (dep), and 1 for the region scale (reg).</figcaption></figure>

shape, with $\Delta_t$ values starting at about 1,000 years for $s$ close to zero, and decrease to about 10 years for $s = 15\%$ per year. This reflects the fact that smaller changes take longer than bigger changes to be detected. The curve for 99% confidence level shows higher $\Delta_t$ values than the 95%, followed by the 90% confidence level curves, reflecting the fact that to increase confidence that the detected change is a true change, and not due to uncertainties, more time is required. At municipality scale, the curves for the three indicators have similar magnitude, decreasing from 1,000 to 10 years. The magnitude of the curves decreases with upscaling, down to decreasing from 100 to 0.1 - 1. This is because the uncertainties decrease with upscaling, therefore translating into smaller $\Delta_t$ values. The bigger the aggregation scale, the more $A$, $P$ and $D$ drift apart in terms of $\Delta_t$ values at high absolute slope values: a significant change in $A$ can be detected faster than a change in $D$, than a change in $P$. For all plots, even with a logarithmic y-axis, the confidence curves present a hyperbolic shape, indicating the highly non-linear relationship between $s$ and $\Delta_t$. Looking at the particular case of Languedoc-Roussillon for 1997-2009, the slope values are smallest for $P$, followed by $A$, and largest for $D$. The corresponding $\Delta_t$ are therefore ranked in the opposite order, with $P$ having the highest $\Delta_t$, followed by $A$, and $D$ having the smallest $\Delta_t$ (value ranges indicated on Figure 4). One should note that, because of the logarithmic y-axis scale, the roughly constant vertical space between the confidence curves means that the difference in $\Delta_t$ in years is actually much higher for high $\Delta_t$ values than for lower $\Delta_t$ values. Because of averaging, the slope values decrease when increasing the aggregation scale. Because uncertainties decrease when upscaling, $\Delta_t$ also decreases when upscaling. The 12-year time lag between the 1997 and 2009 automatic impervious polygons datasets is sufficient to detect significant differences in $A$ and $D$ at region and department scales, and $D$ at inter-municipality scale, but is not sufficient for $P$.
<figure><img src="image_171.png" /><figcaption>Figure 4: Confidence curves at 90, 95, an 99% confidence levels for the minimal monitoring time lag $\Delta_t$ as a function of the absolute value of the relative slope of the indicator changes $s$. The slope ($s_{LR}$) and time lag ($\Delta t_{LR}$) values at 90 and 99% confidence levels for the Languedoc-Roussillon case study are indicated on the plots in the axes units</figcaption></figure>

DISCUSSION

As revealed by figure 4, minimum monitoring time lags vary strongly depending on the considered indicator and scale. These results are based on the average level of uncertainties for the entire Languedoc-Roussillon, whereas the uncertainties may strongly vary locally. Another caveat is that constant slope over time was assumed to build the confidence curves, and they should therefore be used with caution for inferring prospective time lag values. They are not adapted to non-linear cases, such as instantaneous abrupt growth of urban areas. Despite these limitations, we can attempt to compare the order of magnitude of the monitoring frequencies recommended by figure 4 to operational urban monitoring systems. In Europe, the monitoring time lag of the CORINE land cover datasets, decreased from 10 to 6 years. In the United States of America, the National Land Cover Database is updated every 5 years. These two systems are therefore appropriate for monitoring urban areas at scales as fine as the French department scale, but not inter-municipalities and municipalities. The theoretical monitoring time lag, as calculated from equation 1 is, of course, very important for establishing a new operational urban monitoring system, but other factors come into play. First, the theoretical time lag is different for each indicator and varies depending on the considered administrative scale. In this case study, the level of indicator uncertainty of the *D* indicator is simply too high to detect slight
changes. This is reinforced by the cyclic behavior of *D* (succession of dispersion and absorption of impervious polygons in the MUA), violating the linear trend assumption in the analysis. In this case, a maximum monitoring time lag has to be included to allow detecting intra-cycle changes. The time lag required for obtaining a new impervious polygon dataset is also important. Vector databases are typically updated less often than new imagery is acquired, although parts of the study area might be missing from the imagery. In addition, any governments available budget and number of qualified employees may also cause practical limitations. Finally, trade-offs are required.

## CONCLUSION

A specific spatially explicit Monte Carlo simulation of urban sprawl indicator was developed according to the impervious areas geometric and thematic uncertainties, socio-cultural data uncertainties and closing operation uncertainty required to obtain urban areas. This simulation process allowed to propagate and scaling uncertainties up to three urban sprawl indicators calculation: MUA area (*A*), MUA dispersion coefficient (*D*), and population density (*P*), for the entire Languedoc-Roussillon region, France. As expected, upscaling the indicators to larger entities decreased their uncertainties. Confidence curves for minimal monitoring frequency to detect significant changes in indicator values showed that frequency decreased, for all four administrative levels. The highest frequency values were observed for the coarsest administrative scale (region), with values between 0.5 and 1 year for *D* and *A*, and 37 to 58 years for *P*.

## ACKNOWLEDGEMENTS

This work was funded by the French Research Agency in the framework of the program Investissements d'Avenir through the support of the GEOSUD project (ANR-10-EQPX-20).

## References

Allen J., Lu K. (2003). Modeling and prediction of future urban growth in the charleston region of south Carolina: a GIS-based integrated approach. *Ecology and Society* 8(2), 2.

Balestrat M., Chery J., Tonneau J. (2010). Construction of spatial indicators for decision making: interest in a participatory approach: the case of peri-urban languedoc. In *Proceedings of a symposium on Innovation and Sustainable Development in Agriculture and Food, Montpellier, France, 28 June to 1st July 2010.*, pp. hal-00539776. Centre de Coopération Internationale en Recherche Agronomique pour le Développement (CIRAD).

Dupuy S., Barbe E., Balestrat M. (2012). An object-based image analysis method for monitoring land conversion by artificial sprawl use of rapideye and irs data. *Remote Sensing* 4(2), 404-423.

Heuvelink G. B. (1998). Uncertainty analysis in environmental modelling under a change of spatial scale. *Nutrient cycling in Agroecosystems* 50(1-3), 255-264.

INSEE (2012). Recensement de la population: La précision du chiffre de population dans les grandes communes de mtropole. Technical report, Institut National de la Statistique et des Etudes Economiques.

Laurent V., Saint-Geours N., Bailly J.-S., Chery J.-P. (2014a, July 8-11). Local urban sprawl accuracy from image segmentation uncertainties simulation. In A. Shortridge, J. Messina, S. Kravchenko, and A. Finley (Eds.), *Accuracy 2014, 11th international Symposium on Spatial Accuracy Assessment in Natural Resources and Environmental Sciences*, ISARA, East Lansing, Michigan, USA, pp. 146-150.

Laurent V. C. E., Saint-Geours N., Bailly J.-S., Chery J.-P. (2014b, 21-23 may 2014). Simulating geometric uncertainties of impervious areas based on image segmentation accuracy metrics. *South-Eastern European Journal of Earth Observation and Geomatics* 3(2), 37-40. Special Issue: 5th GEOBIA.

Saint-Geours N., Bailly J.-S., Grelot F. (2014). Multi-scale spatial sensitivity analysis: application to a flood damage assessment model. *Environmental Modelling & Software* 60, 153-166.

Yoccoz N. G., Nichols J. D., Boulinier T. (2001). Monitoring of biological diversity in space and time. *Trends in Ecology & Evolution* 16(8), 446-453.
# ‘spup’ – a R package for uncertainty propagation in spatial environmental modelling

K. Sawicka*, G.B.M. Heuvelink

Soil Geography and Landscape Group, Wageningen University, PO Box 47, 6700 AA, Wageningen, The Netherlands

*Corresponding author: kasia.sawicka@wur.nl

## Abstract

Computer models are crucial tools in engineering and environmental sciences for simulating the behaviour of complex systems. While many models are deterministic, the uncertainty in their predictions needs to be estimated before they are used for decision support. Advances in uncertainty analysis have been paralleled by a growing number of software tools, but none has gained recognition for universal applicability, including case studies with spatial models and spatial model inputs. We develop an R package that facilitates uncertainty propagation analysis in spatial environmental modelling. The ‘spup’ package includes functions for uncertainty model specification, propagation of uncertainty using Monte Carlo (MC) techniques, and uncertainty visualization functions. Uncertain variables are represented as objects which uncertainty is described by probability distributions. Spatial auto-correlation within a variable and cross-correlation between variables is also accommodated for. The package has implemented the MC approach with efficient sampling algorithms, i.e. stratified random sampling and Latin hypercube sampling. The MC realizations may be used as an input to the environmental models called from R, or externally. Selected static and interactive visualization methods that are understandable by non-statisticians can be used to visualize uncertainty about the measured input, model parameters and output of the uncertainty propagation.

## Key words

R language, uncertainty analysis, uncertainty propagation, spatial models, Monte Carlo

# I INTRODUCTION

Computer models have become a crucial tool in engineering and environmental sciences for simulating the behaviour of complex static and dynamic systems. However, while many models are deterministic, the uncertainty in their predictions needs to be estimated before they are used for decision support. Currently, advances in uncertainty propagation and assessment have been paralleled by a growing number of software tools for uncertainty analysis, but none has gained recognition for a universal applicability, including case studies with spatial models and spatial model inputs. Due to the growing popularity and applicability of the open source R programming language we undertook a project to develop an R package that facilitates uncertainty propagation analysis in spatial environmental modelling. The tool is intended for researchers and practitioners who understand the problems of uncertainty in data and models, and are looking for a simple, accessible implementation of the universal
methodology for uncertainty assessment. At the same time, it is designed to enable more experienced users to easily understand, customise, and possibly further develop the code.

A number of computational tools are readily available to tackle the uncertainty quantification problem to different degrees. These include both free software, like OpenTURNS (Andrianov et al., 2007), DACOTA (Adams et al., 2009) and DUE (Brown and Heuvelink, 2007), commercial, like COSSAN (Schuëller and Pradlwarter, 2006), or free, but written for a licenced software, e.g. SAFE (Pianosi et al., 2015) or UQLab (Marelli and Sudret, 2014) toolboxes for MATLAB. A broad review of existing software packages is available in Bastin et al. (2013). To the best of our knowledge, however, none of the existent software is specifically designed to be extended by the environmental science community. The use of powerful but complex languages like C++ (e.g. Dakota), Python (e.g. OpenTURNS) or Java (e.g. DUE) often discourages relevant portions of the non-highly-IT trained scientific community from the adoption of otherwise powerful tools.

The R programming language is an important tool for development in numerical and statistical analysis. R has advantages through its advanced statistical capabilities and high-quality graphical output (Ripley, 2001), and is gaining widespread use in science and education. Furthermore, through the use of R packages, the software can be used for a variety of geoscience analyses and visualisations. It has grown tremendously over the last 20 years, with over 8000 packages at the time of preparation of this paper. There is a number of R packages invoking uncertainty analysis through sensitivity analysis or use of a Bayesian framework for model calibration. We have found only one package named 'propagate' that deals with uncertainty propagation explicitly, using similar approaches as described in this paper. The package 'propagate', however, does not provide functionality for spatial models and variables.

## II EMPLOYED (SPATIAL) UNCERTAINTY PROPAGATION ANALYSIS APPROACH

Uncertainty propagation aims to analyse how uncertainties in data (e.g. from measurement error, sampling, interpolation), combined with model uncertainties (e.g. in the model parameters and structure) propagate through the model (Heuvelink et al., 2007). Many environmental phenomena of interest are spatial, temporal or spatio-temporal in nature and often have strong correlations imposed by the physics and dynamics of the natural systems, making uncertainty evaluation difficult. The most frequently used approach represents uncertainty with probability distribution functions (pdfs). The pdf describes the relative likelihood for the random variable to take on a given value and typically it is viewed as a shape of the distribution, for example normal, uniform, lognormal or exponential. It is common for the pdf to be parametrized, i.e. to be characterized by distribution parameters. For example, the normal distribution is parametrized in terms of the mean and the variance, or uniform distribution is parametrized by minimum and maximum values. For situations in which pdfs can be estimated reliably, they have a number of advantages over non-probabilistic techniques. They include methods for describing cross- and auto- correlation between uncertainties, methods for propagating uncertainties through simple algebras or more complex environmental
models, and methods for tracing the sources of uncertainty in environmental data and models (Heuvelink 1998).

A frequently used method for the analysis of uncertainty propagation is the Monte Carlo (MC) method (Hammersley and Handscomb, 1979, Lewis and Orav, 1989). It is very flexible and can reach an arbitrary level of accuracy, and therefore it is generally preferred over analytical methods such as the Taylor series method (Heuvelink, 1998). The idea of the MC method is to compute the output of the model repeatedly, with input values that are randomly sampled from their marginal or joint pdf. The set of model outputs forms a random sample from the output pdf, so that the parameters of the distribution, such as the mean, variance and quantiles, can be estimated from the sample. The method thus consists of the following steps:

1. Characterise uncertain model inputs with pdfs.

2. Repeatedly sample from (spatial) pdfs of uncertain inputs.

3. Run model with sampled inputs and store model outputs.

4. Compute summary statistics of model outputs.

Note that the above ignores uncertainty in model parameters and model structure, but these can easily be included if available as pdfs. A random sample from the model inputs can be obtained using an appropriate pseudo-random number generator (Lewis and Orav, 1989). Note that a conditioning step will have to be included when the model inputs are correlated. Application of the MC method to uncertainty propagation with operations that involve spatial interactions requires the simultaneous generation of realisations from the spatially distributed inputs implying that spatial correlation will have to be accounted for (Heuvelink et al., 1989). For uncertain spatially distributed continuous variables, such as elevation, rainfall and soil organic carbon content, we assume the following geostatistical model:

$$Z(x) = \mu(x) + \sigma(x) \cdot \varepsilon(x) \quad (1)$$

where $\mu$ is the (deterministic) mean of Z, $\sigma$ is a spatially variable standard deviation of the prediction of $\mu$ (spatial variability of $\sigma$ reflects that in some parts of study area the uncertainty is greater than in other parts), and $\varepsilon$ is a standardized, zero-mean, spatially auto-correlated residual modelled with a semivariogram or a correlogram (Diggle and Ribeiro, 2007, Webster and Oliver, 2007, Plant, 2012). The random sample is drawn from the pdf of $\varepsilon$ to further calculate the realizations of Z.

The drawback of the MC method is that the accuracy of the uncertainty assessment is inversely related to the square root of the number of runs N. This means that to double the accuracy, four times as many runs are needed. In complex, multi-variable systems high accuracies are obtained only when the number of runs is very large (i.e. N≥1,000), which may cause the method to become extremely time consuming. The improvement on MC efficiency can be made by employing efficient sampling techniques (e.g. Latin hypercube sampling) and parallel computing.
### III 'spup' (SPATIAL UNCERTAINTY PROPAGATION) PACKAGE DESIGN

The adopted approach for uncertainty propagation analysis dictates the general package design. The 'spup' package provides functions for examining the uncertainty propagation starting from input data and model parameters, via the environmental model onto model outputs (Fig. 1). The functions include uncertainty model specification, stochastic simulation and propagation of uncertainty using MC techniques, as well as several uncertainty visualization functions.

<figure><img src="image_172.png" /><figcaption>**Figure 1.** The 'spup' package design. 'spup' comprises of functions for defining uncertainty model (I), quantifying uncertainty propagation (II) and storing output in a format of data or images.</figcaption></figure>

Uncertain environmental variables are represented in the package as objects whose attribute values may be uncertain and described by probability distributions. Uncertainty assumption may also be ignored, in which case, during the model run the user works with $\mu$ (Eq. 1) as the model input that best represent the reality. Both numerical (e.g. air humidity) and categorical data (e.g. land cover) types are handled. Spatial auto-correlation within an attribute and cross-correlation between attributes is also accommodated for. The attributes may be independent in space, for which a marginal probability density function (mpdf) is defined at each point in space, or may co-vary in space, for which a joint probability density function (jpdf) is defined. Different shapes of marginal pdfs are supported, whereas joint pdfs may be defined for groups of attributes characterized with the normal distribution only. The specification of correlations between errors in space and cross-correlations between objects or attributes is made under the assumption that the correlations depend only on the distance between locations.

For spatially correlated variables the package relies on the unconditional Gaussian simulation implemented in the 'gstat' package (Pebesma, 2004). For drawing realizations of uncertain variables without assumed correlations the package has implemented the MC approach with efficient sampling algorithms, i.e. stratified random sampling and Latin hypercube sampling. The design includes facilitation of parallel computing to speed up MC computation. The MC realizations for uncertainty propagation quantification may be used as an input to the environmental models called from R, or externally.
Selected static (adjacent maps and glyphs) and interactive visualization methods that are understandable by non-experts with limited background in statistics can be used to summarize and visualize uncertainty about the measured input, model parameters and output of the uncertainty propagation.

## IV APPLICATION EXAMPLE – MAPPING SOIL MOISTURE CONTENT FOR THE ALLIER CATCHMENT

As part of a research study in quantitative land evaluation, the World Food Studies (WOFOST) crop simulation model (van Diepen et al., 1989) was used to calculate potential crop yields for floodplain soils of the Allier river in the Limagne rift valley, central France. The moisture content at wilting point ($\Theta_{wp}$) is an important input attribute for the WOFOST model. Because $\Theta_{wp}$ varies considerably over the area in a way that is not linked directly with soil type, it was necessary to map its variation separately to see how moisture limitations affect the calculated crop yield.

Unfortunately, because $\Theta_{wp}$ must be measured on samples in the laboratory, it is expensive and time-consuming to determine it for a sufficiently large number of data points for creating the prediction map by kriging. An alternative and cheaper way is to calculate $\Theta_{wp}$ from other indicators which are cheaper to measure. Because the moisture content at wilting point is often strongly correlated with the moisture content at field capacity ($\Theta_{fc}$) and the soil porosity ($\Phi$), both of which can be measured more easily, it was decided to investigate how errors in measuring and mapping these would work through to a map of calculated $\Theta_{wp}$. Calculation of $\Theta_{wp}$ can be done using a pedo-transfer function, which in this case takes the form of multiple linear regression:

$$ \Theta'_{wp} = \beta_0 + \beta_1 \cdot \Theta_{fc} + \beta_2 \cdot \Phi + \delta \quad (2) $$

where $\Theta'_{wp}$ denotes measured moisture content at wilting point, $\beta_0$, $\beta_1$, and $\beta_2$ are the regression coefficients and $\delta$ denote residuals attributed to lack of model fit and measurement error. The regression coefficients were estimated using standard ordinary least squares regression, ignoring spatial correlation between the observations at the locations. The maps of $\Theta_{fc}$ and $\Phi$ were derived using co-kriging and accounted for spatial cross-correlation between $\Theta_{fc}$ and $\Phi$. Each component on the right hand side of Eq. 2 is subject to uncertainty, which will propagate to uncertainty about $\Theta_{wp}$. Following the adopted MC approach, for each variable and parameter the uncertainty model is defined and 1000 MC samples are drawn. For the spatial variables a linear model of co-regionalization (Wackernagel, 2003) is fitted with use of the ‘gstat’ package and possible realities are simulated. The joint pdf of the model parameters and structural error $\delta$ was estimated using Bayesian calibration (Van Oijen et al., 2005) (note, this is not included in the ‘spup’ package) and a random sample was drawn from their joint posterior distribution. 1000 realizations of $\Theta_{wp}$ was then calculated using Eq. 2 and summary statistics such as mean of prediction and standard deviation were derived.

If an uncertainty analysis with WOFOST would show that the errors in $\Theta_{wp}$ cause errors in the output of WOFOST that are unacceptably large, then the accuracy of the map of $\Theta_{wp}$
would have to be improved. In order to decide how to proceed in such a situation, the
contribution of each individual error source to the overall uncertainty in Θ<sub>wp</sub> was determined
as well. Figure 2 presents the results and these show that both Θ<sub>fc</sub> and Φ, rather than model
parameters and model structural error, form the main source of uncertainty. Thus, the main
source of error in Θ<sub>wp</sub> is the one associated with the kriging errors of Θ<sub>fc</sub> and Φ.

<figure><img src="image_173.png" /><figcaption>Figure 2 Results of uncertainty propagation for soil moisture prediction in the Allier catchment.</figcaption></figure>

V CONCLUSIONS AND FURTHER WORK

We present a tool for uncertainty propagation assessment based on the uncertainty quantification framework described in e.g. Heuvelink et al. (2007). As the theoretical framework and implementation of the package progress, its application to real cases will be necessary, both to test the algorithms and usability of the tool, and to demonstrate the importance of assessing uncertainty in environmental data. The ‘spup’ package is being developed and used within the project “Quantifying Uncertainty in Integrated Catchment Studies (QUICS)”. QUICS aims to carry out research in order to take the implementation of the Water Framework Directive (WFD) to the next level and improve water quality management by assessing the uncertainty of integrated catchment model water quality predictions. Currently, the potential case studies for the ‘spup’ application include uncertainty propagation analysis
with the LandscapeDNDC model (Haas et al., 2012) and German Schwingbach catchment data, and Metaldehyde Prediction model developed currently for the Severn Tent Water, water provider in the Midlands, UK. Finally, 'spup' will be introduced to the wider scientific community through CRAN (The Comprehensive R Archive Network), where many more challenges will be faced, including the time and resources required to implement an uncertainty assessment and the need to make uncertainty analyses understandable to non-statisticians.

## References

ADAMS, B. M., BAUMAN, L. E., BOHNHOFF, W. J., DALBEY, K. R., EBEIDA, M. S., EDDY, J. P., ELDRED, M. S., HOUGH, P. D., HU, K. T., JAKEMAN, J. D., SWILER, L. P. & VIGIL, D. M. 2009. DAKOTA, A Multilevel Parallel Object-Oriented Framework for Design Optimization, Parameter Estimation, Uncertainty Quantification, and Sensitivity Analysis: Version 5.4 User's Manual.

ANDRIANOV, G., BURRIEL, S., CAMBIER, S., DUTFOY, A., DUTKA-MALEN, I., DE ROCQUIGNY, E., SUDRET, B., BENJAMIN, P., LEBRUN, R., MANGEANT, F. & PENDOLA, M. OpenTURNS, an open source initiative to Treat Uncertainties, Risks'N Statistics in 520 a structured industrial approach. ESREL'2007 Safety and Reliability Conference, 2007 Stavenger, Norway.

BASTIN, L., CORNFORD, D., JONES, R., HEUVELINK, G. B. M., PEBESMA, E., STASCH, C., NATIVI, S., MAZZETTI, P. & WILLIAMS, M. 2013. Managing uncertainty in integrated environmental modelling: The UncertWeb framework. *Environmental Modelling & Software*, 39, 116-134.

BROWN, J. D. & HEUVELINK, G. B. M. 2007. The Data Uncertainty Engine (DUE): A software tool for assessing and simulating uncertain environmental variables. *Computers & Geosciences*, 33, 172-190.

DIGGLE, P. & RIBEIRO, P. J. 2007. *Model-based geostatistics*, Springer.

HAAS, E., KLATT, S., FRÖHLICH, A., KRAFT, P., WERNER, C., KIESE, R., GROTE, R., BREUER, L. & BUTTERBACH-BAHL, K. 2012. LandscapeDNDC: a process model for simulation of biosphere-atmosphere-hydrosphere exchange processes at site and regional scale. *Landscape Ecology*, 28, 615-636.

HAMMERSLEY, J. M. & HANDSCOMB, D. C. 1979. *Monte Carlo methods*, London, Chapman and Hall.

HEUVELINK, G. B. 1998. *Error Propagation in Environmental Modelling with GIS*, London, Taylor and Francis.

HEUVELINK, G. B. M., BROWN, J. D. & VAN LOON, E. E. 2007. A probabilistic framework for representing and simulating uncertain environmental variables. *International Journal of Geographical Information Science*, 21, 497-513.

HEUVELINK, G. B. M., BURROUGH, P. A. & STEIN, A. 1989. Propagation of errors in spatial modelling with GIS. *International Journal of Geographical Information Systems*, 3, 303-322.

LEWIS, P. A. W. & ORAV, E. J. 1989. *Simulation methodology for statisticians, operations analysts, and engineers* Pacific Grove, Wadsworth & Brooks/Cole.

MARELLI, S. & SUDRET, B. 2014. UQLab: A Framework for Uncertainty Quantification in Matlab. *Vulnerability, Uncertainty, and Risk*. American Society of Civil Engineers.

PEBESMA, E. J. 2004. Multivariable geostatistics in S: the gstat package. *Computers & Geosciences*, 30, 683-691.

PIANOSI, F., SARRAZIN, F. & WAGENER, T. 2015. A Matlab toolbox for Global Sensitivity Analysis. *Environmental Modelling & Software*, 70, 80-85.

PLANT, R. E. 2012. *Spatial data analysis in ecology and agriculture using R*, CRC Press.

RIPLEY, B. D. 2001. The R project in statistical computing. *MSOR Connections*. The newsletter of the LTSN Maths, Stats \textbackslash{}\&\textbackslash{} OR Network, 23--25.

SCHUËLLER, G. I. & PRADLWARTER, H. J. 2006. Computational stochastic structural analysis (COSSAN) - a software tool. *Structural Safety*, 28, 68-82.

VAN DIEPEN, C. A., WOLF, J., VAN KEULEN, H. & RAPPOLDT, C. 1989. WOFOST: a simulation model of crop production. *Soil Use and Management*, 5, 16-24.
VAN OIJEN, M., ROUGIER, J. & SMITH, R. 2005. Bayesian calibration of process-based forest models: bridging the gap between models and data. *Tree Physiol*, 25, 915-27.

WACKERNAGEL, H. 2003. *Multivariate Geostatistics: An Introduction with Applications*, Springer.

WEBSTER, R. & OLIVER, M. A. 2007. *Geostatistics for environmental scientists*, John Wiley & Sons.
# Combining spatial and thematic uncertainty and sensitivity analysis for mountain natural hazard assessment

Jean-Marc Tacnet*<fnref n="1" />, Guillaume Dupouy<fnref n="1" />, Franck Bourrier<fnref n="2" />, Frédéric Berger<fnref n="2" />, Dominique Laigle<fnref n="1" />, Nicolas Crimier<fnref n="3" />, Kamel Mekhnacha<fnref n="3" />, Michel Mémier<fnref n="4" />

*Corresponding author: jean-marc.tacnet@irstea.fr

## Abstract

This article presents an application of the Hybrid uncertainty propagation method to natural hazard assessment. It proposes a comparison with the usual probabilistic Monte-Carlo method, and propagates uncertainties related to both spatial and thematic variables.

## Keywords

Uncertainty, sensitivity, Monte-Carlo, Hybrid.

## I CONTEXT

Rockfalls and debris-flows are dangerous phenomena in mountains that cause severe damage to exposed assets and population (Figure 1). Risk assessment is based on both thematic and spatial information: physical phenomenon features such as height, speed, impact loads but also the extent of the phenomena is essential to be assessed (Figure 2).

<figure><img src="image_174.png" /></figure>

<footnote n="1">¹ Université Grenoble Alpes, Irstea, UR ETGR, Snow avalanches Engineering and Torrent Control Research Unit – 2 rue de la papeterie, BP 76, F-38402 Saint-Martin-d'Hères Cedex, France.</footnote>

<footnote n="2">² Université Grenoble Alpes, Irstea, UR EMGR, Mountain Ecosystems Research Unit – 2 rue de la papeterie, BP 76, F-38402 Saint-Martin-d'Hères Cedex, France.</footnote>

<footnote n="3">³ Probayes – 82 allée Galilée, F-38330 Montbonnot, France.</footnote>

<footnote n="4">⁴ Sintegra – 11 chemin des Prés, 38241 Meylan, France.</footnote>
For both rockfalls and debris-flows hazard assessment, numerical models are used. In our case:

- For rockfalls, the code RockyFor3D described in Dorren et al. (2006) and Bourrier et al. (2009) simulates the 3D propagation of the rocks as a succession of free flights through the air and rebounds on the soil, modelled by a Digital Elevation Model (DEM). Using input thematic variables related to the falling rock characteristics, it provides spatialized outputs, namely kinetic energy, height, speed of the boulder.

- For debris-flows, Laigle et al. (2003) proposes the lave2D model, dedicated to the computation of the unconfined free-surface spreading of materials with complex rheology. It is based upon the steep-slope-shallow-water-equations which are solved using a finite volume technique which requires first to mesh the domain of interest. Equations are solved taking into account the material rheological behavior and the field topography represented by a DEM. Using input thematic variables related to the input hydrograph and rheological parameters, it provides spatialized outputs, namely flow height and speed.

Both thematic data and DEM can be affected by imperfection (imprecision, uncertainty) depending on the terrain morphology, data acquisition and processing methods. The main issues are therefore:

- To evaluate the quality of available DEM depending on its nature and acquisition process (Lidar, satellite, commercial maps etc.).

- To analyse and represent the effect of spatial information uncertainty on models' results.

This paper describes how innovative methods and tools are used to assess both thematic and spatial data uncertainties and to characterize the impact of the DEM uncertainties on the results of numerical modelling simulation.

## II SPATIAL AND THEMATIC UNCERTAINTY PROPAGATION: SCIENTIFIC APPROACH AND MAIN RESULTS

The approach compares two uncertainty propagation approaches, to propagate both DEM and thematic uncertainty:

- The usual, probabilistic, Monte-Carlo method.

- A possibilistic, so-called Hybrid method proposed by Baudrit et al. (2006).

That last Hybrid approach considers the different aspects of information imperfection, especially its imprecision (lack of information, inaccuracy of measure...). It relies on the usual probability theory, the possibility theory as defined by Zadeh (1978) and Dubois et al. (2000), and the belief function theory detailed by Shafer (1976) and Smets et al. (1994). This method generalizes, under some restrictive conditions, the usual Monte Carlo method.

The DEM variability is modelled by a stochastic field using the new ModTer software developed by ProBayes under the ModTer project consortium, and detailed by Crimier et al. (2016) (Figure 3). It takes into account heterogeneity of DEM quality as an input of simulation model: it produces terrain simulations and confidence maps related to altimetric information.
Examples of results (quantile of rock passing heights) of a numerical simulation of rockfalls are proposed: they show the influence of data imperfection including those resulting from expert assessments on the simulation results (Figure 4). To demonstrate the effects of spatial data quality, a sensitivity analysis describes the contribution of the DEM uncertainty on the global uncertainty.

<figure><img src="image_175.png" /><figcaption>Figure 3: A global framework to evaluate and propagate both thematic and spatial uncertainties using new ModTer software</figcaption></figure>

<figure><img src="image_176.png" /><figcaption>Figure 4: Example of rockfall simulation results</figcaption></figure>

*Acknowledgements:* Those developments have been partially funded by the ModTer project: a RAPID project supported by the French Defence Procurement Agency (DGA) and the General Directorate for Competitiveness, Industry and Services (DGCIS).

## References

Baudrit C., Guyonnet D., Dubois D. (2006). Joint propagation and exploitation of probabilistic and possibilistic information in risk assessment. *IEEE Transactions* 14, 593–608.

Bourrier F., Dorren L., Nicot F., Berger F., Darve F. (2009). Toward objective rockfall trajectory simulation using a stochastic impact model. *Geomorphology* 110(3-4), 68–79.

Crimier N., Mekhnacha K., Tacnet J.-M., Memier M. (2016). Terrain uncertainties modelling software : Modter. In *submitted to ESREL2016 conference*, Glasgow, Scotland.

Dorren L., Berger F., Putters U. (2006). Real-size experiments and 3-d simulation of rockfall on forested and non-forested slopes. *Natural Hazards and Earth System Sciences* 6(1), 145–153.

Dubois D., NGuyen H., Prade H. (2000). Possibility theory, probability and fuzzy sets : Misunderstandings, bridges and gaps. In *Fundamental of Fuzzy Sets*, pp. 343–438.

Laigle D. Hector A.-F., Hübl J., Rickenmann D. (2003). Comparison of numerical simulation of muddy debris flow spreading to records of real events. In *Proceedings of the 3rd International Conference 'Debris-Flow Hazards Mitigation'*, Davos, Switzerland.

Shafer G. (1976). *A mathematical theory of evidence*. Princeton University Press.

Smets P., Kennes P. Kennes R. (1994). The transferable belief model. *Artificial Intelligence* 66(2), 191–234.

Zadeh L.A. (1978). Fuzzy sets as a basis for a theory of possibility. *Fuzzy Sets and Systems* 1, 3–28.
# Performing Multi-Temporal Spatial Data Analysis for Coastal Areas and Assessing Thematic Accuracy

Daniel Cohenca¹and Carlos Antonio Oliveira Vieira¹*

¹Federal University of Santa Catarina, Geoscience Department, Florianópolis, Trindade, SC 88040-900 Brazil

*Corresponding author: carlos.vieira@ufsc.br

## Abstract

It was used TM Landsat images from 1985, 1994, 2004 and 2014 and Object Based Oriented Analysis (OBIA) image classification techniques in order to study changes in land use and land cover for two cities: Passo de Torres and Balneário Gaivota at the south region of Brazil. Based on these analyses, four informational classes were chosen: Natural Vegetation, Water, Urban and Agriculture Use. Agriculture use includes croplands, livestock activities and reforestation areas. It was used post-classification multi-temporal analysis. The aim of this paper is to perform multi-temporal spatial data analysis in order to identify trajectories of changes on land use and land cover, during the last three decades. This paper also discusses an alternative methodology to assess the thematic accuracy in multi-temporal analysis. The multi-temporal analysis results show that this study area has considerably changed in the last three decades. The most frequent transition was due to the transformation of natural areas into cropland, livestock and reforestation ones, from the year of 1985 to 2014. It was observed that most of transition for urban class occurred from natural vegetated areas, while few transitions occurred from agricultural areas, even considering these classes as being more representatives in the whole study area. This transition was especially common in coastal urban division into plots of land. Besides understanding those processes of change, this study points out the importance of better planning, licensing and monitoring of urban coast division into plots of land to prevent the continuous loss of natural ecosystems and impoverishment of environmental quality. The accuracy assessment discussion shows the importance in choose a better classification strategy in order to perform the multi-temporal analysis in remote sensed data. Furthermore, it is observed that a considerable amount of research needs to be undertaken before the spatial characterization of thematic accuracy associated with multi-temporal analysis can be adequately reported in standardized format and legends.

## I INTRODUCTION

Brazilian coastal areas have experienced fast changes in recent years, on form and model of occupation driven by geographical, economic, social forces and public policies. The multi-temporal analysis of changes on land use and land cover is a method focused on understanding, measuring and analysing spatial and temporal changes, both qualitatively and quantitatively.

On the other hand, it is widely acknowledged that classification of remotely sensed imagery has variable and often poor quality. The cause and nature of these errors has been the subject of extensive researches (e.g., Lu and Weng, 2007) in order to improve the accuracy of remotely sensed products.

The aim of this paper is to perform multi-temporal spatial data analysis in order to identify trajectories of changes on land use and land cover, during the last three decades. This study also discusses an alternative methodology to assess the thematic accuracy in multi-temporal analysis.
## II STUDY AREA AND MATERIAL

The study area covers the municipalities of Passo de Torres and Balneário Gaivota in coastal southern state of Santa Catarina and covers an area of 25,030 hectares (Figure 1). The landscape can be considered homogeneous with elevation ranging from sea level to 30 meters. This study area is entirely composed of Holocene coastal sedimentary plain (Dieh and Horn Filho, 1996), originally covered by herbaceous marshes, typical of wetlands or fixative dunes, shrub or tree (Falkenberg, 1999; CONAMA, 1999). In addition to the native vegetation, the area has diverse types of land cover including family farms, monoculture timber species (*Pinus* and *Eucalyptus* cultivation), rice cultivation, livestock and urban areas at different stages of occupation.

<figure><img src="image_177.png" /><figcaption>Figure 1: Study area location - municipalities of Passo de Torres and Balneário Gaivota.</figcaption></figure>

## III MATERIAL AND METHODS

Landsat TM scenes were selected at intervals of about 10 years, consisting in six spectral bands for each scene, for the images classification process (Table 1). It is important to note that the precise geometric correction process is crucial for the multi-temporal analysis (Morisette and Khorram, 2000; Mundia and Aniya, 2005; Singh, 1989). In order to perform this geometric correction procedure, it was used the image Landsat 8 / OLI 01/30/2014 as reference image (available with geometric correction by the United States Geological Survey - USGS). It was used a plane coordinates UTM system and the Datum WGS84. The resulting average quadratic error was less than 15 meters (half pixel) for all scenes. The resampling method used was the nearest neighbour, in order to maintaining the radiometric properties of the original images. After geometric correction, the four images were clipped using the feature IBGE municipal boundaries of the study area.

In addition, high-resolution images (e.g., SPOT and RapidEye Images) and photogrammetric data (Flights years 1978 and 1996) were used in the accuracy assessment procedures as reference data.
<table><caption>Table 1: Sensors, dates and bands used (orbit 220, point 080).</caption><thead><tr><td>Satellite - Sensor</td><td>Date</td><td>Bands</td></tr></thead><tbody><tr><td>LandSat 5 - TM</td><td>09/07/1985</td><td>1,2,3,4,5,7</td></tr><tr><td>LandSat 5 - TM</td><td>18/07/1994</td><td>1,2,3,4,5,7</td></tr><tr><td>LandSat 5 - TM</td><td>14/08/2004</td><td>1,2,3,4,5,7</td></tr><tr><td>LandSat 8 - OLI</td><td>30/01/2014</td><td>2,3,4,5,6,7</td></tr></tbody></table>

### (i) Multi-temporal classification procedures

In the classification process, it was used an object oriented supervised classification method: OBIA (Object Based Image Analysis). According to Cohenca and Carvalho (2015) the OBIA method is based on image segmentation which consists of subdividing the image into homogeneous regions, called objects, which are the basic elements of the classification process (Benzet al., 2004). Each object has spectral characteristics defined by its average value of pixel's radiometry and also geometric characteristics such as: contextual information and texture. The segments are classified following a tree process, whose rules for the differentiation of the classes are defined by the interpreter (Francisco and Almeida, 2012).

For the whole set of bands were tested weights parameters, such as: scale, shape and compactness of the objects in each of the images, which define separability of objects and, consequently, the size and number of objects to be generated (Table 2).

<table><caption>Table 2: Parameters used for segmentation of Landsat images.</caption><thead><tr><td>Images</td><td>Scale</td><td>Shape</td><td>Compactness</td></tr></thead><tbody><tr><td>TM-09/07/1985</td><td>10</td><td>0.1</td><td>0.8</td></tr><tr><td>TM-18/07/1994</td><td>10</td><td>0.1</td><td>0.9</td></tr><tr><td>TM-14/08/2004</td><td>10</td><td>0.1</td><td>0.9</td></tr><tr><td>OLI-30/01/2014</td><td>90</td><td>0.1</td><td>0.5</td></tr></tbody></table>

The distinctive characteristics of atmospheric conditions on the acquisition of images and different sensors used, motivates each scene be classified separately. The evaluation of the quality of the segmentation was carried out through joint evaluation of the classification results (Darwish et al., 2003).

From field visits and analysis of high resolution images, were initially defined eleven thematic informational classes, considered representative of the diversity of types of coverage in the study area. In order to reduce source of imprecision and considering the objective of identifying the changes in coverage between native vegetation, areas of agriculture use and urban areas, regardless of the specific type of use in each segment, the initially defined eleven classes were grouped into four final classes: *Natural Vegetation*, *Water*, *Urban and Agriculture Use*. Agriculture use includes croplands, livestock activities and reforestation areas.

In order to assess the accuracy of the classification procedures, it was used the method proposed by Congalton and Green (1999). It was selected 200 pixels at random (50 samples
per class with a minimum distance of 100 meters between samples), for each of these scenes independently (Korting, 2007). Comparing classification results and reference samples were generated confusion (or error) matrices, from which were derived indexes, such as: overall accuracy, Kappa and Conditional Kappa coefficients for each scene.

## (ii) Multi-temporal analysis

It was chosen an overlay vector method for comparing the thematic classifications performed for each analysed stage, considering the complexity and diversity of land uses and types of land cover in the study area for the multi-temporal analysis of changes.

Data analysis was quantitatively performed globally as well as regionalized, using the limits of the different settlement pattern identified, based on the main type of use, property size and dynamics of occupation that have been mapped in historical survey of the occupation process. This mapping was crossed with maps of land use change, allowing identify settlement patterns that have undergone major change and what were these changes in different time intervals.

Thus, it was analysed the trajectory of multi-temporal land cover transitions in order to understand the changes detected between the studied decades, in both: global and regionalized manner.

## VI RESULTS

Overall accuracies, derived from the confusion matrix, ranged between 84% and 88%, and the Kappa coefficients ranged between 0.793 and 0.836 (Table 3). These accuracies values can be considered very satisfactory. According to the conditional Kappa coefficients, the class with the highest accuracy rating was the class water, followed by urban class. The classes that had greater confusion were Natural and Agriculture Uses. This kind of confusion seems to be due to a lower separability of these classes. Figure 2 shows thematic images for the years of 1985, 1994, 2004 and 2014.

<table><caption>Table 3: Overall accuracies, Kappa and Conditional Kappa Coefficients.</caption><thead><tr><th colspan="2">Classes</th><th>Images</th><th>1985</th><th>1994</th><th>2004</th><th>2014</th></tr></thead><tbody><tr><td rowspan="4">Conditional Kappa</td><td>Water</td><td></td><td>0,973</td><td>0,947</td><td>1,000</td><td>0,895</td></tr><tr><td>Natural</td><td></td><td>0,847</td><td>0,720</td><td>0,632</td><td>0,692</td></tr><tr><td>Agriculture</td><td></td><td>0,783</td><td>0,636</td><td>0,803</td><td>0,787</td></tr><tr><td>Urban</td><td></td><td>0,748</td><td>0,869</td><td>0,771</td><td>0,844</td></tr><tr><td colspan="2">Kappa</td><td></td><td>0,836</td><td>0,793</td><td>0,800</td><td>0,807</td></tr><tr><td colspan="2">Overall Accuracy</td><td></td><td>88%</td><td>84%</td><td>85%</td><td>85%</td></tr></tbody></table>

The multi-temporal analysis results show that 19.8% of the study area has changed in the last three decades. Between 1985 and 2014, the most frequent transition (57.1% of changes) was due to the transformation of natural areas into agriculture use. The expansion of the urban class over natural areas responded for 24.3% of changes. The conversion of agriculture use into urban use responded for 5.1% of change, while conversion of agricultural areas to vegetation regenerated accounted only for 2.8% of change (Table 4).
<figure><img src="image_178.png" /><figcaption>Figure 2: Image Classification results for the years of 1985, 1994, 2004 and 2014.</figcaption></figure>

It was observed that 82% of the transition to urban class occurred from natural vegetated areas, while 18% occurred from agricultural areas, even considering these classes as being more representatives in the whole study area. This transition was especially common in coastal urban division into plot of land.

Remnants of natural class accounted for 17% of the study area in 2014, while in 1985 it was 32%, been particularly concentrated in these coastal urban divisions into plots of land.

Performing multi-temporal spatial data analysis operations on data of unknown accuracy will result in a product with low reliability and restricted use in the decision-making process, while errors deriving from one source can propagate through the database via derived products. Moreover, current accuracy assessment methods are based on non-spatial statistics derived from the confusion or error matrix, which compares the output of a classifier against known reference data. Although these measures are in widespread use, none of them considers the spatial distribution of erroneously classified pixels, either implicitly or explicitly.
<table><caption>Table 4: Frequent transitions in land coverage.</caption><thead><tr><th rowspan="2">Transitions</th><th colspan="2">1985-1994</th><th colspan="2">1994-2004</th><th colspan="2">2004-2014</th><th colspan="2">1985-2014</th></tr><tr><th>Area (ha)</th><th>%</th><th>Area (ha)</th><th>%</th><th>Area (ha)</th><th>%</th><th>Area (ha)</th><th>%</th></tr></thead><tbody><tr><td>Natural → Agriculture Use</td><td>1090</td><td>49,2</td><td>815</td><td>63,5</td><td>1370</td><td>72,0</td><td>2822</td><td>57,1</td></tr><tr><td>Natural → Urban</td><td>522</td><td>23,5</td><td>322</td><td>25,1</td><td>346</td><td>18,2</td><td>1200</td><td>24,3</td></tr><tr><td>Agriculture Use → Urban</td><td>60</td><td>2,7</td><td>74</td><td>5,7</td><td>127</td><td>6,7</td><td>250</td><td>5,1</td></tr><tr><td>Agriculture Use → Natural</td><td>411</td><td>18,6</td><td>58</td><td>4,5</td><td>26</td><td>1,4</td><td>138</td><td>2,8</td></tr><tr><td>Other Transitions</td><td>134</td><td>6,0</td><td>14</td><td>1,1</td><td>33</td><td>1,7</td><td>532</td><td>10,8</td></tr><tr><td>TOTAL</td><td>2217</td><td></td><td>1282</td><td></td><td>1902</td><td></td><td>4942</td><td>100</td></tr></tbody></table>

Vieira and Mather (2001) presents one possible way to characterize the spatial distribution of the errors in a thematic classification is by generating a *distance image* showing the distance from individual pixels to the multivariate means of the classes to which they have been assigned (Figure 3c). Either the Euclidean distance or the Mahalanobis distance measure can be used. The former, however, implies spherical clusters in feature space, while the latter takes into account the covariance between the features on which the classification is based. The individual distances are scaled onto a 0-255 range, and displayed as a grey scale image. Darker pixels are spectrally "nearer" to their class centroid (in the sense of statistical distance), and are thus more likely to be classified correctly. On the other hand, pixels with higher distance values are spectrally further from the centroid of the class to which they were assigned, and are thus more likely to be misclassified. This process (presented in Figure 3b) allowed to get known the thematic uncertainties associated for each pixel of the classified image.

<figure><img src="image_179.png" /><figcaption>Figure 3: Spatial characterisation of classification errors using thematic image generated by Maximum Likelihood Classifier (100 by 100 pixels). A) Reference Image; B) Reliability Image; C) Distance Image.</figcaption></figure>
Another possible way to characterise the spatial distribution of the errors in a thematic classification is by directly comparing thematic images with their respective ground truth maps. One of the products of this comparison should be an *error binary image* (Figure 3b), in which each point takes the value 0 (correctly labelled) or 1 (erroneously labelled). On-going researches has been developed in order to propagate these uncertainties through the multi-temporal analyse model and a thematic uncertainty map could be generated, as suggested by Vieira and Mather (2001).

## V CONCLUSIONS

The multi-temporal analysis results show that has considerably changed in the last three decades. Between 1985 and 2014, the most frequent transition was due to the transformation of natural areas into agricultural, livestock and reforestation ones. The expansion of the urban class over natural areas responded for 24% of changes. The conversion of agricultural use into urban one responded on exceeding of 5% of change, while conversion of agricultural areas to vegetation regenerated accounted only for less than 3% of change. It was observed that most of transition for urban class occurred from natural vegetated areas, while few transitions occurred from agricultural areas, even considering these classes as being more representatives in the whole study area. This transition was especially common in coastal urban division into plots of land. In 2014, remnants of natural class accounted for 17% of the study area, while in 1985 it was 32%, been particularly concentrated in these coastal urban divisions into plots of land.

Besides understanding those processes of change, this study points out the importance of better planning, licensing and monitoring of urban coast division into plots of land to prevent the continuous loss of natural ecosystems and impoverishment of environmental quality.

The accuracy assessment results point out the importance in choose a better classification strategy in order to perform the multi-temporal analysis in remote sensed data. Furthermore, it is also shown that a considerable amount of research needs to be undertaken before the spatial characterization of thematic accuracy associated with multi-temporal analysis can be adequately reported in standardized format and legends.

## References

Benz, U.C.; Hofmann, P.; Willhauck, G.; Lingenfelder, I.; Heynen, M. (2004). Multi-resolution, object-oriented fuzzy analysis of remote sensing data for GIS-ready information. *ISPRS Journal of Photogrammetry and Remote Sensing*, v. 58, n. 3, p. 239-258.

Cohenca, D.; Carvalho, R. (2015). Comparação de métodos de classificação OBIA, Máxima Verossimilhança e Distância Mínima em imagem OLI/Landsat-8 em área de alta diversidade de uso do solo. In: *SIMPÓSIO BRASILEIRO DE SENSORIAMENTO REMOTO*, XVII, 2015, João Pessoa. Anais... São José dos Campos: INPE, p. 1035-1042.

CONAMA (1999). Resolução Nº 261 de 30 de junho de 1999. Aprova parâmetro básico para análise dos estágios sucessivos de vegetação de restinga para o Estado de Santa Catarina. Conselho Nacional do Meio Ambiente-CONAMA, Brasília: DOU, n° 146, 02/08/1999. p. 29-31.

Congalton, R.; Green, K. (1999). *Assessing the Accuracy of Remotely Sensed Data - Principles and Practices*. Boca Raton: CRC Press, Taylor & Francis Group, 183p.

Darwish, A.; Leukert, K.; Reinhardt, W. (2003). Image segmentation for the purpose of object-based classification. *International Geoscience and Remote Sensing Symposium*, v. 3, p. III: 2039-2041.

Dieh, F.L.; Horn Filho, N.O. (1996). Compartimentação geológico-geomorfológica da zona litorânea e planície costeira do estado de Santa Catarina. *Notas Técnicas*, v. 9, p. 39-50.
Falkenberg, D.B. (1999). Aspectos da flora e da vegetação secundária da restinga de Santa Catarina, sul do Brasil. *Insula*, Florianópolis, n. 28, p. 1-30.

Francisco, C. N.; Almeida, C. M. (2012). Interpretação de imagens orbitais por meio de sistema especialista para mapeamento e cobertura da terra em região montanhosa. *Sociedade e Natureza*, v. 24, n. 2, p. 283-302.

Korting, T.S.; Fonseca, L.M.G.; Castejon, E.F.; Namikawa, L.M. (2014). Improvements in sample selection methods for image classification. *Remote Sensing*, v. 6, n. 8, p. 7580-7591.

Lu, D.; Weng, Q. (2007). A survey of image classification methods and techniques for improving classification performance. *International Journal of Remote Sensing*, v. 28, n. 5, p. 823-870.

Morisette, J.T.; Khorram, S. (2000). Accuracy assessment curves for satellite-based change detection. *Photogrammetric Engineering and Remote Sensing*, v. 66, n. 7, p. 875-880.

Mundia, C.N.; Aniya, M. (2005). Analysis of land use/cover changes and urban expansion of Nairobi city using remote sensing and GIS. *International Journal of Remote Sensing*, v. 26, n. 13, p. 2831-2849.

Singh, A. (1989). Review article digital change detection techniques using remotely-sensed data. *International journal of remote sensing*, v. 10, n. 6, p. 989-1003.

Vieira, C. A. O., Mather, P. M. (2001). On the Assessment of Spatial Reliability of Thematic Images. In: HALLS, P., Peter. (Org.). *Innovation in GIS: Spatial Information and the Environment*. UK: Editora London, Cap.01, p. 120-135.
<figure><img src="image_180.png" /></figure>
<figure><img src="image_181.png" /></figure>

Spatial model sensitivity analyses
<figure><img src="image_182.png" /></figure>
# Multiway sensitivity analysis of the fusion of earth observation, topography and social media data for rapid flood mapping

Didier G. Leibovici *, Julian F. Rosser

Nottingham Geospatial Institute
University of Nottingham, Nottingham, UK

*Corresponding author: didier.leibovici@nottingham.ac.uk

## Abstract

Social media can be a valuable source of information for both detecting and mapping hazard events. Rapidly estimating flood inundation extent can aid decision-making during crises and help in damage assessment. A data fusion method, based on a Bayesian statistical model, which uses weights of evidence to calculate and combine the variables according to their influence in order to map the flooded areas, was previously developed for an experimental timeframe covering the flooding in Oxford, UK during January 2014. The method used three data sources: geotagged photographs, optical remote sensing and high resolution terrain mapping. In this work we aim at evaluating the sensitivity of this method. The multivariate sensitivity analysis proposed, uses an empirical non-parametric approach to generate the variations of the posterior probability of water presence from different levels of prior uncertainties of the input data sources. The sensitivity assessment is spatially depicted as a multiway correspondence analysis of the generated multi-entry table: uncertainty levels for each of the three data sources, the geographical space and the descriptive measures of the output variation, i.e. a 5-way array.

## Keywords

crowdsourcing, social media, inundation,

# I INTRODUCTION

Crowdsourcing has been used in various disaster management situations such as earthquakes (Barrington et al. 2011) and flood damage (Tomnod 2015). Typically, these crowdsourcing scenarios are task orientated. For example, the volunteers are assessing the aftermath of a disaster from satellite image interpretation. Increasingly, analysis may be undertaken on non-authoritative social media reports. For example, analysis of Tweets has been used to source data used in quasi-real time for flood extent estimation (Smith et al. 2015), where a hydrological flooding extent model with initialisations from the volunteered observed locations of flood. Rosser et al. (2016) proposed to use a data fusion method combining multiple evidence of water presence including social media data sources to estimate the inundation extent. This paper aims to perform a sensitivity analysis of that approach. The data fusion method utilises viewshed analysis (Panteras et al. 2015) and weights of evidence (Tehrany et al. 2014). The workflow adopted is shown in Figure 1 with the elements involved in the sensitivity analysis reported here highlighted in green (see section II). The three domains of data sources involved in the method are: Landsat-8 water mapping based on the Modified Normalised Difference Water Index (Xu, 2006), geotagged photographs sourced from Flickr to estimate the citizen's extent estimation based on viewshed analysis, and modelled topographic variables sourced from a high resolution Digital Elevation Model.
A classical approach to analyse the sensitivity of a model is to vary factors one-at-a-time and record the effects on the output, such as previously undertaken in multi-criteria mapping problems (Chen et al., 2013). Here, we suggest a multivariate sampling for each combination of the uncertainty levels for the three dimensions. After initial estimation of potential input uncertainties, we modulate these prior uncertainties at different levels to estimate the variation in output uncertainty from the Monte Carlo simulations. The output uncertainty of the posterior probability of water presence is summarised from the simulated distributions in each combination of the input uncertainties. The generated multi-entry table capturing the uncertainty associations (illustrating how sensitive the rapid inundation extend method is), is then analysed from multiway correspondence analysis (Leibovici 2010).

<figure><img src="image_183.png" /><figcaption>Figure 1: BPMN workflow of the rapid flood inundation data fusion method used in Rosser *et al.* (2016) – in green are the data input where uncertainty sampling takes place.</figcaption></figure>

The method is quite computationally challenging as for each of the 6 174 000 pixels of the studied area, 64 000 evaluations of the workflow are performed: 4 levels of variations in each of the three dimensions with 1000 simulations. On a smaller area of 759 980 pixels (12% of the total area) 1 simulation took 108 seconds (1 core 2.7Ghz), so potentially taking 666 days for the whole study area without parallelisation. For the smaller area and parallelising the code using 20 cores and with only 100 simulations per cell, this is reduced to ~22 hours computation time and will be our first experiment.

## II UNCERTAINTIES IN THE DATA FUSION WORKFLOW

Despite the use of quantile classes in the weight of evidence method (see details in Rosser *et al.* 2016) the uncertainty of the measurements or data capture could make the output posterior probability of water presence uncertain. The purpose of the sensitivity analysis is to evaluate
the effect of the uncertainties and to highlight how one data source alone, or in combination with another, can impact on the efficacy of the model.

The whole study area around the city of Oxford, UK, and its characteristics are described in Rosser et al. 2016; the data collection coverage for the January 2014 flood event are represented in Figure 2. Ground truth delineation of the flood extent was based on expert interpretation of colour infra-red aerial photography, provided by the Environment Agency Geomatics group (see Figure 6).

<figure><img src="image_184.png" /><figcaption>Figure 2: The number of Flickr flood images, Landsat-8 and aerial image (the latter used for ground-truth) acquisition times over the experimental time frame (shaded area). The gauged river level height is also shown.</figcaption></figure>

Besides the representivity and adequation of the data sources for the study, which can be discussed as potential bias, the basis uncertainties taken into account per data source dimension are:

* Landsat-8 dimension

The satellite image matching the peak river level is used for the whole study period, so some variations could be expected in the MNDWI. In the first instance MNDWI uncertainty was estimated from the overall spatial variation of the MNDWI itself but a locally dependent (say 25 pixels' neighbourhood) variance map could be used on the index or directly on the wavelength bands. This gave a standard deviation of 0.22.

* topographic dimension

The Lidar DTM has a given horizontal and vertical uncertainty which changed due to resampling from 1m pixel size to 5m pixel size. To simplify the experiment, slope and elevation uncertainties were estimated from the observed variations globally, respectively 3 and 16 as 1 standard deviation.

* Flickr dimension

The geolocation of the social media data was resampled using as basis a 68% circular error, usually given by the mobile phone (corresponding to 1 standard deviation). As no such information was available in the Flickr data, we took as basis a value of 5m.

The weight of evidence method assumes conditional independence of the 'theme variables' (the predictors) given the outcome:

$$p(V_1V_2...V_m | O) = p(V_1 | O)p(V_2 | O)...p(V_m | O) \quad (1)$$
Assuming normal distribution of the multivariate input, we will sample the multivariate inputs as univariates with the additional assumption that the dependency to the outcome lies in the mean parameter of the distribution.

### III MULTIVARIATE SENSITIVITY ANALYSIS

The multivariate sensitivity analysis applied here is based on error propagation through the data fusion workflow. Modulation of the basis uncertainties (section II) was undertaken with four levels of standard deviations: ½ *sd*, 1 *sd*, 2 *sd*, 3 *sd* where *sd* is the basis uncertainty as described in section II. Monte Carlo simulations were derived with normal distributions centred on the observed data values with standard deviations as the modulated uncertainties. A modified workflow taken from Figure 1 was used as the tasks in red ('MNDWI' and 'Slope & Elevation') did not need to run during the sensitivity analysis. Also, the other tasks ('Quantile classification', 'Otsu binarisation' and 'Weight of Evidence') were run using the initial estimated parameters or weights. Therefore the classes adopted the same cut-points and W+ and W- values for the weight of evidence associated to these classes for each simulation. Running the whole workflow with new class and weight estimations each time would encompass a sensitivity analysis, but also a reliability analysis of the workflow itself. Studying the variation of the W+ and W- values would nonetheless be an important reliability and model adequacy feature of the approach to investigate in future work.

The resulting uncertainty of the multivariate input variations has been summarised (using min, Q1, median, mean, Q3, max, var of the simulated distributions at each pixel) in the 5-way array shown in Figure 3.

<figure><img src="image_185.png" /><figcaption>Figure 3: Multivariate sensitivity analysis output stored in a 5-way array, t, s and e indexing each the levels of input uncertainties, g the pixel location and m the distribution indicator of the output uncertainty.</figcaption></figure>

To analyse this sensitivity captured in the 5-way array we can perform multiway correspondence analysis (see Leibovici 2010), allowing an association of the levels of input uncertainty to the intensity of the output uncertainty. The analysis can be performed either on the 3-way table when aggregating across the space and considering only the variance as a summary measure, or on the 4-way table keeping the space dimension. The 5-way analysis can be used to describe further potential associations within the output uncertainty distribution. Before doing so, it is worth looking at the overall variation of the output when aggregating the values over the spatial domain. Figure 4 depicts that overall the workflow is not sensitive to the earth observation
domain as the change in variation in these 4 levels is rather small. Note, also, that the maxim-
um variance observed on the map over the 100 simulations was 4.5 E-4, and the variance over
the pixels of the mean of the simulations was 1.12E-04! The patterns of variation between so-
sial media uncertainties and topographic uncertainties are very similar in each of the 4 levels of
earth observation (MNDWI). Nonetheless, the earth observation uncertainty level 4 shows a
larger spread over the topographic levels of uncertainty for social uncertainty level 3, therefore
indicating some presence of interaction here. Globally as expected as the uncertainty increases
(in social media and topographic data) so does the output variation: linear increasing trend at
all earth observation levels. There is not much difference in topographic uncertainty levels ex-
cept t1 lower and with an apparent bigger increase in interaction with social dimension. The in-
teraction topographic and social media is only convincing when looking at the difference s1-
s4 for t1 versus the others.

<figure><img src="image_186.png" /><figcaption>Figure 4: Spatially aggregated output variations (rescaled variance) for each uncertainty com-
bination level (s, t, e).</figcaption></figure>

When performing the multiway correspondence analysis (FCAk), the first principal tensor
fitted by the algorithm is the ‘independence’ principal tensor (vectors of the margins which can
be also automatically set), then the following fitted tensors relate to lack of independence (Leibovici 2010). It is interesting to note that the FCAk of the spatially aggregated data (3-way table) gives 99.98% as representing the independence whilst the analysis on the 4-way table, which includes spatial differentiations, gives 91.95%. So, there is some spatial variation (in the output uncertainty), even though both are high values showing a very small lack of independence between the input uncertainty dimensions (both chi-square values are very small). For the 3-way analysis, the lack of independence has its two-way component breakdowns expressing 2% associated to the social variation margin, 2.7% associated to the marginal topological variation and 84.2% associated to earth observation variation margin, confirming no effect of the latter. Nonetheless, the pure interaction captured in the following (in the order of the decomposition) represents 7% of the lack of independence: *e3* with *t1* and *s4* opposed to *e4* with *t1* and *s1*. The 4% left (100 - (84.2+2+2.7+7)) are concentrated in 3 principal tensors representing nearly 1%, two of these as two-way decompositions associated to the 7% principal tensor. The 4-way analysis puts most of the lack of independence (51%) into a gradient in social media variation levels and spatial differentiation associated with marginal topographic and earth observation variations (Figure 5). The next principal tensor represents 7%, then a few at 1.5%, then a multitude accounting for about 30%. In Figure 5, the blue areas associated to higher social media uncertainty levels (*s3* and *s4*) coincide to the edges of the viewsheds, the brown areas correspond to areas with limited Flickr data.

<figure><img src="image_187.png" /><figcaption>Figure 5: FCAk principal tensor n°70 associated margins to *t* and *e* representing 51% of lack of independence: gradient of *s* is *s1*=-1.2, *s2*=-0.7, *s3*=0.4 and *s4*=1.4 (Values are mapped using an equal intervals).</figcaption></figure>
IV DISCUSSION AND CONCLUSIONS

This paper illustrates a multivariate sensitivity analysis for a data fusion algorithm to estimate a flooding extent from a series of predictors. The approach defines each input variable as part of a sampling dimension, with a choice of prior levels of uncertainty per dimension. The multi-way correspondence analyses provide a characterisation of the potential multivariate associations sensitive to the outcome of the data fusion workflow. The example showed a small sensitivity globally and very small influence of the earth observation dimension (EO) without interaction with the other dimensions. If, overall, the increase in input uncertainty implied greater output uncertainty, the topographic dimension showed mostly a difference between its lowest uncertainty level (t1) against the other ones and in its interaction with the social media data dimension: bigger increase under t1 as the social dimension uncertainty than the other ones. Spatially, most of the sensitivity relates to the edges of the viewsheds when sampling with high uncertainty levels in the Flickr data, which can be expected.

<figure><img src="image_188.png" /><figcaption>Figure 6: Sensitivity analysis study area, ground-truth and estimated flood extent (derived from Receiver Operating Characteristic).</figcaption></figure>

In this sensitivity analysis, the spatial autocorrelation in the uncertainty sampling has been neglected or reduced simply to the difference in means of the distributions (normal distributions
centred on the observed values). A simple local uncertainty estimation discussed in section II could compensate some of this drawback as the variances could be locally correlated (even though this will still mean sampling at pixel level). Computationally, this would add a potentially expensive step in the initial uncertainty estimation of the inputs but would be done only once. The weight of evidence method producing the posterior probability of flood at a pixel given the predictors comes from the general form:

$$ \operatorname{logit}(p(O = 1 | V_1 V_2 \dots V_m)) = \operatorname{logit}(p(O = 1)) + \sum_j W^{\pm}(V_j) \quad (2) $$

using the prior and the weights of evidence (see Rosser et al. 2016 for more details). Our initial prior was uniform over the area and taken from Rosser et al. (2016) in order to allow direct comparison. Note that a Bayesian updating could be performed through deriving a Beta distribution from the Monte Carlo simulation (estimating the parameters from the mean and variance of the simulated sample) and using this to generate a new prior in the process. In the last 10 years Bayesian non-parametric methods (e.g. Gelfand et al. 2005) have developed which should allow flexibility in modelling spatial and multivariate dependency. However, for large datasets this could be still computationally prohibitive (Berrocal 2016).

## References

* Berrocal, V. J. (2016). Identifying Trends in the Spatial Errors of a Regional Climate Model via Clustering. *Environmetrics* 27 (2): 90–102.

* Chen Y., Yu J. Khan S., (2013). The spatial framework for weight sensitivity analysis in AHP-based multi-criteria decision making, *Environmental Modeling Software*. 48 129–140.

* Gelfand, A.E., Kottas, A., MacEachern, S.N. (2005) Bayesian nonparametric spatial modeling with Dirichlet process mixing. J. Am. Stat. Assoc. 100, 1021–1035

* Leibovici D.G., (2010). Spatio-temporal Multiway Decomposition using Principal Tensor Analysis on k-modes: the R package PTAk. *Journal of Statistical Software* 34(10), 1–34.

* Panteras G., Wise S., Lu X., Croitoru A., Crooks A., and Stefanidis A., (2015). Triangulating Social Multimedia Content for Event Localization Using Flickr and Twitter. *Transactions in GIS* 19 (5): 694–715. doi:10.1111/tgis.12122.

* Rosser J., Leibovici D., Jackson M., (2016). Rapid flood inundation mapping using social media, remote sensing and topographic data. Submitted to: *Natural Hazards*.

* Smith, L., Liang, Q., James, P. and Lin, W., (2015). Assessing the utility of social media as a data source for flood risk management using a real-time modelling framework. *Journal of Flood Risk Management, (Article first published online: 27 MAR 2015 DOI: 10.1111/jfr3.12154)*

* Tehrany MS, Pradhan B, Jebur MN (2014) Flood susceptibility mapping using a novel ensemble weights-of-evidence and support vector machine models in GIS. *J Hydrol* 512:332–343. doi: 10.1016/j.jhydrol.2014.03.008

* Xu H., (2006). Modification of normalised difference water index (NDWI) to enhance open water features in remotely sensed imagery. *International Journal of Remote Sensing*. 27. 3025-3033.
# A GPU-based Solution for Accelerating Spatially-Explicit Uncertainty- and Sensitivity Analysis in Multi-Criteria Decision Making

Christoph Erlacher*<fnref n="1" />, Seda Şalap-Ayça <fnref n="2" />, Piotr Jankowski<fnref n="2" />, Karl-Heinrich Anders<fnref n="1" />, Gernot Paulus<fnref n="1" />

*Corresponding author: c.erlacher@cuas.at

## Abstract

This study presents a GPU-based approach to accelerate the time consuming parts of the variance-based spatially-explicit Uncertainty and Sensitivity Analysis (U-SA). Performance comparisons were conducted in respect to a CPU-based NumPy and a GPU-based CUDA implementation. Preliminary results reported herein suggest that the proposed approach will provide a quantitative decision quality measure in complex and comprehensive spatial multi-criteria decision making processes and will allow reasonable computational times making spatially explicit variance-based SA applicable and attractive for large-size problems. Furthermore, it will be beneficial for different application domains like natural hazard risk assessment, landscape assessment, infrastructure planning, environmental impact assessment or identification of land use strategies for sustainable regional development.

## Keywords

Spatial Decision Support, GPU, CUDA, Parallelization, Spatially-Explicit Uncertainty- and Sensitivity Analysis, Multi-Criteria Decision Making, MCDM.

## 1. Introduction

Uncertainty and Sensitivity Analysis (U-SA) is a critical step in multi-criteria decision making. The major goal is to verify the robustness and stability of an implemented model with respect to the existing uncertainties. Saltelli et al. (2008) define sensitivity analysis as “*the study of how uncertainty in the output of a model (numerical or otherwise) can be apportioned to different sources in the model input*”. Uncertainties can refer to the criterion weights, the decision rules, the standardization process of the input criteria and the accuracy as well as the resolution of the input data (Hwang and Yoon, 1981; Ligmann-Zielinska and Jankowski, 2012; Ganji et al., 2016). Only a handful of studies focused on spatially-explicit U-SA of model input factors (Ligmann-Zielinska and Jankowski, 2012; Ligmann-Zielinska and Jankowski, 2014; Şalap-Ayça and Jankowski 2016). As stated by Feizizadeh et al. (2014), Ligman-Zielinska and Jankowski (2014) and Şalap-Ayça and Jankowski (2016) spatially-explicit uncertainty and sensitivity analysis can be a time consuming process, which depends on the project area, the number of criteria and models runs associated with the sample size.

<footnote n="1">¹ Carinthia University of Applied Sciences – Department of Geoinformation and Environmental Technologies, Austria</footnote>

<footnote n="2">² San Diego State University – Department of Geography, United States of America</footnote>
Therefore, this paper reports on the study of high-performance computing approach for spatially-explicit U-SA of the assessment of agricultural land units with a multiple-criteria decision making model. Two implementations, a CPU-based NumPy and a GPU-based CUDA (Compute Unified Device Architecture) solution, are compared regarding the computational time of performing the simulations. The Environmental Benefit Index (EBI) is used by the United States Department of Agriculture (USDA) to prioritize agricultural land units based on environmental benefits and select high-scoring units for USDA crop reserve program. This simple scoring model has been used by the Farm Service Agency (FSA) since 1990 to rank farmers' requests to enroll land into the Conservation Reserve Program during each general sign-up period (competitive bidding). The model includes five environmental factors at the top level ("Wildlife", "Water-Quality", "Soil-Erosion", "Enduring-Benefits" and "Air-Quality"). Detailed information concerning the EBI-Framework can be found in *Şalap-Ayça and Jankowski (2016).*

## 2. Theoretical Background

### 2.1. Multi-Criteria Decision Making

Multi-Criteria Decision Making (MCDM) techniques (Malczewski, 1999; Malczewski, 2006; Malczewski and Rinner, 2015) belong to the field of Spatial Decision Support Systems (SDSS) and consider a set of alternatives that are evaluated on the basis of conflicting criteria. MCDM problems include several objectives and each objective is expressed by a set of criteria where every criterion refers to a certain influence value that is indicated by the criterion weight. The combination of MCDM and Geographic Information Systems (GIS) facilitates the integration of the spatial aspects with MCDM. Geographical data are expressed as criterion maps, which can be differentiated into factor and constraint maps.

In short, the workflow of the decision making process includes the standardization of the criteria (transformation) regarding comparability and the application of decision rules like Weighted Linear Combination (WLC), Ordered Weighted Averaging (Yager, 1988; Malczewski, 2006) or Analytic Hierarchy Process (Saaty, 1980), which results in the monotonic rank-order of decision alternatives leading to final recommendation.

### 2.2. Spatially-Explicit Uncertainty and Sensitivity Analysis

As previously mentioned, sensitivity analysis is a critical part of MCDM in order to verify the robustness and stability of the implemented model. In spatial multi-criteria evaluation uncertainty can be potentially associated with the selection of decision criteria, criteria measurement (inaccuracy and errors) and expert's preferences that represent the criteria weights. Ligmann-Zielinska and Jankowski (2008) state, that a decision situation typically comprises as spatial as well as spatial aspects. However, the vast majority of reported sensitivity analysis studies and evaluation methods concern only the as spatial nature of decision situations while spatial characteristics are only included implicitly. Due to the fact that the spatial distribution of options and their criteria values may potentially influence the rank order of alternatives, the authors emphasize the consideration of spatial criteria (e.g. *proximity, compactness, contiguity, etc.*) and spatial weighting (varying the criterion importance over space and assigning different weights to different spatial units) in GIS-based multi criteria evaluation, which calls for new, spatially-explicit methods of U-SA. As described by Ligmann-Zielinska and Jankowski (2014) uncertainty analysis helps to quantify the variability of model outcomes, whereas sensitivity analysis focuses on the identification of decision criteria or criteria weights that are responsible for the variability. Ligmann-Zielinska and Jankowski (2014) examined the robustness of land suitability evaluation with help of Monte Carlo Simulation (MCS) and variance-based sensitivity analysis. Variance-based sensitivity
analysis is often recommended because it represents a model independent procedure and is applicable for spatially-explicit data. Furthermore, in contrast to one at a time SA (Si, first order sensitivity indices) variance-based SA also incorporates the interaction of input factors (ST, total order sensitivity indices). The (S,ST) pair offers a succinct yet comprehensive measure of input influence (Ligmann-Zielinska and Jankowski, 2014).

## 2.3. GPU-based Parallelization

In the course of this research project, a GPU-based prototype was developed focusing on an acceleration strategy for performing the MCS. MCS represents the most time consuming component of the variance-based spatially-explicit uncertainty and sensitivity analysis method, which depends on the number of criteria and model runs associated with the sample size. Tang and Jia (2014) present an approach that focuses on parallel Graphics Processing Units (GPUs) in order to enable and accelerate the sensitivity analysis of large agent- based modelling of spatial opinion exchange. A solution for accelerating time-consuming calculations is domain-specific and depends on the implemented mathematical functions. General-Purpose Computing on Graphics Processing Units (GPGPUs) are used for computer graphics and entertainment, or GPU co-processors that are constructed to accelerate massively parallel floating-point operations (Krömer et al, 2014). The GPU accelerated computing prototype utilizes the CUDA architecture for NVIDIA GPUs. CPUs (Central Processing Unit) consist of few cores and are optimized for sequential serial processing. GPUs have a parallel architecture that incorporates thousands of smaller but more efficient cores. The calculations are organized into threads, blocks and grids, and are distributed with respect to the streaming multiprocessors and cores. Detailed explanations on GPU programming can be found in Wilt (2013).

## 3. Methodology

The original spatially-explicit Sensitivity and Uncertainty Framework (Ligmann-Zielinska and Jankowski, 2014) refers to a Python implementation that incorporates two dimensional NumPy arrays for the weight samples (N) and the criterion maps (k). The weight samples are generated using the quasi-random Sobol's experimental design (Software Simlab) and the Ideal Point (IP) decision rule is used to compute the suitability surfaces (see Figure 1). The total number of model runs (R) is represented by the following formula:

$$R = (k + 2) * N \qquad (1)$$

In order to compare the performance of the original Framework with the CUDA-GPU version two-dimensional NumPy arrays representing evaluation criteria are flattened in order to calculate the suitability surfaces (see Figure 2). Consequently, the dimension of the final suitability surface is represented by the number of model runs (rows) and the number locations (raster cells). Each row of the final suitability surface represents one suitability map that is associated with the corresponding weight sample and the used decision rule (see Figure 3).
<figure><img src="image_189.png" /><figcaption>Figure 1: Displays the workflow to conduct a variance-based spatially-explicit U-SA.</figcaption></figure>

<figure><img src="image_190.png" /><figcaption>Figure 2: Represents the data structure to perform the GPU-based MCS.</figcaption></figure>
<table><caption>Figure 3: Illustrates the structure of the suitability surfaces, where each row represents one suitability raster map based on the decision rule and the criterion weights.</caption><thead><tr><th colspan="6">Weight Sample</th><th colspan="10">Output of Simulations</th></tr></thead><tbody><tr><td>#1</td><td>0.2</td><td>0.15</td><td>0.1</td><td>0.3</td><td>0.25</td><td>#1</td><td>0.45</td><td>0.1</td><td>0</td><td>0.72</td><td>0.35</td><td>0.61</td><td>0.79</td><td>0.55</td><td>0.38</td></tr><tr><td>#2</td><td>0.25</td><td>0.1</td><td>0.15</td><td>0.2</td><td>0.3</td><td>#2</td><td>0.35</td><td>0.245</td><td>0</td><td>0.495</td><td>0.37</td><td>0.63</td><td>0.775</td><td>0.525</td><td>0.375</td></tr><tr><td>#3</td><td>0.15</td><td>0.2</td><td>0.2</td><td>0.25</td><td>0.2</td><td>#3</td><td>0.315</td><td>0.23</td><td>0</td><td>0.545</td><td>0.38</td><td>0.635</td><td>0.775</td><td>0.535</td><td>0.35</td></tr><tr><td>#4</td><td>0.3</td><td>0.15</td><td>0.15</td><td>0.2</td><td>0.2</td><td>#4</td><td>0.375</td><td>0.215</td><td>0</td><td>0.52</td><td>0.355</td><td>0.625</td><td>0.765</td><td>0.605</td><td>0.35</td></tr><tr><td>#...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>#..</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr></tbody></table>

The algorithm to perform the MCS for the CUDA-GPU implementation is divided into the host part (CPU) and the device part (GPU). Parallel portions of the application are executed on the device as kernels, whereas the host represents the program logic for memory allocation, data partitioning and recombination (see Figure 4). The data partitioning depends on the computing and memory capabilities of the device in order to avoid memory overflow. Therefore, the first task of the controller is to identify the available computing capacity of the device. For example, GPU memory could be allocated by another software or task. According to the determined GPU resources the number of partitions is calculated. For each partition the device memory for a subset of the suitability surfaces has to be allocated, which is necessary for calling kernels. A kernel is launched as a grid of thread blocks and each thread block consists of a specified number of threads (e.g. warp size). The number of blocks per grid and the number of threads per block depend on the number of rows and columns of the suitability surface.

<figure><img src="image_191.png" /><figcaption>Figure 4: Overview about the CUDA-GPU based MSC computation.</figcaption></figure>

Figure 5 illustrates a more detailed representation of the CUDA-GPU based Monte Carlo Simulations to generate the suitability surfaces.
<figure><img src="image_192.png" /><figcaption>Figure 5: Detailed workflow of the CUDA-GPU based MSC computation.</figcaption></figure>

## 4. First Speed Up Results

In order to facilitate the performance tests, different NVIDIA GPUs were investigated. The results confirm the efficacy of employing GPU-based solutions for computationally demanding, spatially-explicit U-SA. For example, for the Tesla k40c device, a computational acceleration up to a factor of 150 was achieved. The case study tested model inputs with 5 input factors, 1,475,464 locations, and 2,664 simulation runs.

Tables 1-3 illustrate CPU- and GPU specification overview and the performance comparison of different GPUs as well as a performance comparison regarding the CPU (NumPy – a package for scientific computing with Python) and CUDA (Anaconda NumbaPro) implementation.

<table><caption>Table 1: Tesla GPU Specification Overview.</caption><thead><tr><th>Tesla - GPU Specifications</th><th>Tesla K40 Workstation</th><th>Tesla K20m Server</th></tr></thead><tbody><tr><td>Peak double-precision floating point performance (Tflops)</td><td>1.43</td><td>1.17</td></tr><tr><td>Peak single-precision floating point performance (Tflops)</td><td>4.29</td><td>3.52</td></tr><tr><td># of CUDA cores</td><td>2880</td><td>2496</td></tr><tr><td>Total Global Memory - GDDR5 (GB)</td><td>12</td><td>5</td></tr><tr><td>Memory bandwidth (GB/sec)</td><td>288</td><td>208</td></tr></tbody></table>
<table><caption>Table 2: Intel Xeon CPU Specification Overview.</caption><thead><tr><th>CPU Specifications</th><th>Intel Xeon E5-1620 v3 Workstation</th><th>Intel Xeon E5-2697 v2 Server</th></tr></thead><tbody><tr><th scope="row">Intel Smart Cache (MB)</th><td>10</td><td>30</td></tr><tr><th scope="row"># of Cores</th><td>4</td><td>12</td></tr><tr><th scope="row"># of Threads</th><td>8</td><td>24</td></tr><tr><th scope="row">Processor Base Frequency (GHz)</th><td>3.5</td><td>2.7</td></tr><tr><th scope="row">Max memory bandwidth (GB/sec)</th><td>68</td><td>59.7</td></tr></tbody></table>

Table 3 reports the relative speed ups for both Tesla GPU cards with respect to the computational demand of the CPUs. The Tesla K40c device (Workstation) is 1.809 times faster than the Tesla K20m device (Server).

<table><caption>Table 3: CUDA-NumPy Performance Comparison.</caption><thead><tr><th>Performance Comparison</th><th>Workstation</th><th>Server</th></tr></thead><tbody><tr><th scope="row">Avg. Elapsed Time for CUDA (min)</th><td>0.1215</td><td>0.2198</td></tr><tr><th scope="row">Avg. Elapsed Time for NumPy (min)</th><td>18.6712</td><td>25.7937</td></tr><tr><th scope="row">CUDA Speed Up</th><td>153.6724</td><td>117.3508</td></tr></tbody></table>

## 5. Discussion and Future Prospects

The implemented GPU-based prototype offers reasonable computation times to perform spatially-explicit U-SA, which makes this type of sensitivity analysis applicable and potentially attractive for distributed-output (spatially-explicit) models without using expensive servers or super computers. Furthermore, this solution allows the integration of further decision rules (WLC, OWA etc.) by changing the kernel function and can be applied for different application areas.

The proposed workflow incorporating the data preparation regarding the supported interface and the partitioning and recombination of the GPU-based computation is limited to local operations. For supporting the calculation of local variability (focal- and zonal operations) in respect to criteria values and weights, a conceptual definition for the data structure as well as partitioning- and recombination approach have to be adapted (Malczewski, 2011; Şalap-Ayça and Jankowski; 2016).

Further accelerations can be achieved by implementing additional CUDA kernels for the uncertainty and sensitivity computations. Additionally, the distribution of the CUDA calculations among many GPUs will be investigated too as the next step of reported here research.
# 6. Acknowledgement

Financial support for this work was provided by the National Science Foundation Geography and Spatial Sciences Program Grant No. BCS-1263071. Any opinion, findings, conclusions, and recommendations expressed in this material are those of the authors(s) and do not necessarily reflect the views of the National Science Foundation.

# References

Feizizadeh, B., Jankowski, P., Blaschke, T. (2014). A GIS based Spatially-explicit Sensitivity and Uncertainty Analysis Approach for Multi-Criteria Decision Analysis. *Computers and Geosciences*, Volume 64, pp. 81-95.

Ganji, A., Maier, H., Dandy, G. (2016). A modified Sobol' sensitivity analysis method for decision-making in environmental problems. *Environmental Modelling & Software*, Volume 75, pp. 15-27.

Hwang, C.L., Yoon, K., 1981. Multiple Attribute Decision Making Methods and Applications: A State of the Art Survey. Springer-Verlag, Berlin.

Krömer, P., Platoš, J., Snášel, V. (2014). Nature-Inspired Meta-Heuristics on Modern GPUs: State of the Art and Brief Survey of Selected Algorithms. *International Journal of Parallel Programming*, Volume 42, Issue 5, pp. 681-709.

Ligmann-Zielinska, A., Jankowski, P. (2014). Spatially-Explicit Integrated Uncertainty and Sensitivity Analysis of Criteria Weights in Multicriteria Land Suitability Evaluation. *Environmental Modelling & Software*, Volume 57, pp. 235-247.

Ligmann-Zielinska, A., Jankowski, P. (2012). Impact of proximity-adjusted preferences on rank-order stability in geographical multicriteria decision analysis. *J. Geogr. Syst.* Volume 14, Issue 2, pp. 167-187.

Ligmann-Zielinska, A., Jankowski, P. (2008). *A Framework for Sensitivity Analysis in Spatial Multiple Criteria Evaluation*. In: T.J. Civa et al., eds. GIScience 2008, LNCS 5266, Springer-Verlag Berlin, pp. 217-333.

Malczewski, J. (1999). GIS and Multicriteria Decision Analysis. New York, Wiley.

Malczewski, J. (2006). Ordered weighted averaging with fuzzy quantifiers: GIS-based multicriteria evaluation for land-use suitability analysis. *International Journal of Applied Earth Observation and Geoinformation*. Volume 8, pp. 270-277.

Malczewski, J. (2011). Local Weighted Linear Combination. *Transactions in GIS*, Volume 15, Issue 4, pp. 439-455.

Malczewski, J., Rinner, C. (2015). Multicriteria Decision Analysis in Geographic Information Science, Advances in Geographic Information Science. New York, Springer.

Saaty, T.-L., 1980. The analytic hierarchy process. New York, McGraw-Hill.

Şalap-Ayça, S., Jankowski, P. (2016). Integrating Local Multi-Criteria Evaluation with Spatially Explicit Uncertainty-Sensitivity Analysis, Spatial Cognition & Computation, DOI: 10.1080/13875868.2015.1137578

Saltelli, A., Ratto, M., Andres, T., Campolongo, F., Cariboni, J., Gatelli, D., Saisana, M., Tarantola, S. (2008). *Global Sensitivity Analysis: The Primer*. John While & Sons, Ltd.

Tang, W., Jia, M. (2014). Global Sensitivity Analysis of a Large Agent-Based Model of Spatial Opinion Exchange: A Heterogeneous Multi-GPU Acceleration Approach. *Annals of the Association of American Geographers*, Volume 104, Issue 3, pp. 485-509.

Wilt, N. (2013). The CUDA Handbook: A Comprehensive Guide to GPU programming. New York, Addison-Wesley.

Yager, R.-R. (1988). On ordered weighted averaging aggregation operators in multi-criteria decision making. *IEEE Transactions on Systems, Man and Cybernetics, Volume 18*, Issue 1, pp. 183-190, IEEE Press Piscataway, NJ, USA.
# Uncertainty propagation in urban hydrology water quality modelling

J.A. Torres-Matallana*<fnref n="1" />,<fnref n="2" />, U. Leopold<fnref n="1" />, G.B.M. Heuvelink<fnref n="2" />

## Abstract

Uncertainty is often ignored in urban hydrology modelling. Engineering practice typically ignores uncertainties and uncertainty propagation. This can have large impacts, such as the wrong dimensioning of urban drainage systems and the inaccurate estimation of pollution in the environment caused by combined sewer overflows. This paper presents an uncertainty propagation analysis in urban hydrology modelling. The case study was the Haute-Sûre catchment in Luxembourg for one yearly time series measured in 2010, and 10 individual rainfall events measured in 2011. The selection of model input variables for uncertainty quantification was based on their level of uncertainty and model sensitivity. Probability distribution functions were defined to represent the uncertainty of the input variables. We applied a Monte Carlo technique using a simplified model, EmiStatR, which simulates the volume and substance flows in urban drainage systems. We focus in loads and concentrations of chemical oxygen demand and ammonium, which are important variables for wastewater and surface water quality management.

**Keywords:** uncertainty propagation, urban hydrology modelling, EmiStatR, Monte Carlo

## I INTRODUCTION

Uncertainty is often ignored in urban hydrology modelling (Mitchell et al., 2007; Bach et al., 2014). Engineering practice typically ignores uncertainties and uncertainty propagation, among others because of lack of user-friendly implementations (Schellart et al., 2010). This can have large impacts, such as the wrong dimensioning of urban drainage systems and the inaccurate estimation of pollution in the environment caused by combined sewer overflows (CSOs).

Six approaches to address uncertainty in the context of urban water systems are identified (Walker et al., 2003; Refsgaard et al., 2007; van Keur et al., 2008; Bach et al., 2014): 1) determinism; 2) statistical uncertainty; 3) scenario uncertainty; 4) qualitative uncertainty; 5) recognised ignorance; and 6) total (unrecognised) ignorance. Following van Keur et al. (2008), *determinism* applies when we have knowledge with absolute certainty about the system under analysis. The *statistical* approach is useful when it is possible to describe in statistical terms the uncertainty. This occurs when errors and uncertainties can be quantified by probability distribution functions (pdfs). The *scenario* approach, in contrast, applies when quantitative probabilities cannot be determined. It is used when possible outcomes can be listed without pretending that the list is exhaustive and without attaching probabilities to each possible outcome (Brown, 2004). The *qualitative* approach is used when *uncertainty cannot be characterised statistically, and not all outcomes are known* (Brown, 2004). Recognised ignorance occurs in a situation of awareness of lack of knowledge (van Keur et al., 2008). Finally, *total ignorance* is the state of *complete lack of awareness about possible outcomes* (van Keur et al., 2008). Amongst the

<footnote n="1">¹Luxembourg Institute of Science and Technology, Luxembourg</footnote>

<footnote n="2">²Wageningen University, Netherlands</footnote>

<footnote>*Corresponding author: arturo.torres@list.lu</footnote>
approaches described above, in this paper we will limit ourselves to the statical approach to
characterise and propagate uncertainties.

Sources of uncertainty in the context of the evaluation and design of urban water infrastructure
are identified by Neumann (2007): model structure, model parameters, errors in input data, and
in numerical and computational procedures. Our focus is on model input as a primary source of
uncertainty.

Uncertainty propagation analysis can be done analytically by means of the Taylor series method, or numerically by means of Monte Carlo simulation. This paper presents an uncertainty propagation analysis in urban hydrology modelling. For this we apply a Monte Carlo analysis using the EmiStatR model, which simulates the volume and substance flows in urban drainage systems. We use a case study from the Haute-Sûre catchment in Luxembourg for one yearly time series measured in 2010, and 10 individual rainfall events measured in 2011. We focus on substances as chemical oxygen demand (COD) and ammonium (NH₄), which are important variables for wastewater and surface water quality management.

## II MATERIALS AND METHODS

### 2.1 The EmiStatR model

The EmiStat model (Klepiszewski and Seiffert, 2013) is a Microsoft Excel based model that was developed for the Luxembourgish water authority as a tool for the evaluation of planning scenarios of sewer systems. EmiStatR is the R implementation of EmiStat. It provides a fast estimation of combined wastewater emissions by imposing a strongly simplified representation of the real-world system. It can aid the planning and design of hydraulic structures and pollutant handling, without the requirement of extensive simulation tools. The EmiStatR model includes six main components to simulate combined sewage discharges of a catchment (Figure 1).

The sewer system under investigation includes tank structures to store first flush pollutant peaks. After filling of the storage volume a combined sewage overflow structure discharges subsequent volume and pollutant inflows exceeding the structures outflow from the wastewater treatment plant (WWTP) to the receiving water.

In EmiStatR a simple volume balance taking into account inflow volume, present storage capacity and outflow to the WWTP is implemented to simulate the volume in the tank structure. In case of an overflow the pollutant concentrations in the CSO are equivalent to the combined sewage inflow concentrations of the structure.

The pollutants typically taken into account are total COD and NH₄. The variable COD is the standard used in the framework of the dimensioning of CSO structures in Luxembourg. NH₄ represents a diluted substance which can have a significant impact on surface water quality due to possible transformation to ammonia (NH₃).

At the CSO tank structure a simple volume balancing takes place: 1) substance and volume flows are stored and discharged to the WWTP if the storage volume is not completely filled up; 2) if the storage volume is completely filled up the proportion of the volume inflow which is not discharged to the WWTP goes to the CSO.

### 2.2 Sewer systems of the Haute-Sûre sub-catchments

The study area is a sub-catchment of the Haute-Sûre catchment in the north-west of Luxembourg. The combined sewer system of the sub-catchment drains eight villages: Büderscheid
<figure><img src="image_193.png" /><figcaption>Figure 1: Main components of the EmiStatR model: 1) Dry Weather Flow (DWF) including Infiltration Flow (IF); 2) Pollution of DWF; 3) Rain Weather Flow (RWF); 4) Pollution of RWF; 5) Combined Sewage Flow (CSF) and pollution; and 6) Combined Sewer Overflow (CSO) and pollution. CSOT = CSO tank (Background adapted from: Sanitary-District (2015))</figcaption></figure>

(BUD), Dahl (DAH), Eschdorf/Ost (ESO), Goesdorf (GOE), Heiderscheid (HEI), Kaundorf (KAU), Nocher (NOC) and Nocher-Route (NOR). The local sewer system downstream each village has a CSO tank to store pollutant peaks in the first flush of combined sewage flows. Table 2 shows the general characteristics of each CSO tank for each village. Figure 2 depicts their locations and the delineation of the catchment. The main land use types in the villages are residential, smaller industries and farms. Outside of the villages forest as well as agricultural arable and grassland are the dominating land uses. The receiving water bodies at CSO structures are tributaries of the river Sûre.

The general input data of the EmiStatR model for simulating the eight CSO structures is presented in Tables 1 and 2.

## 2.3 Model input uncertainty assessment

Following recommendations from Nol et al. (2010), not all model inputs were taken into account in the Monte Carlo uncertainty propagation analysis. Only those inputs that have a large uncertainty and to which the model is sensitive were included. The selection of model input for uncertainty quantification was based on a list of all model inputs, their level of uncertainty (low or high) and the level of model sensitivity (low or high). The level of uncertainty of the inputs was defined by expert judgement, literature research, measurements of different model inputs in the sewer system, and interviews with experts. The developers of the EmiStatR model, which have detailed knowledge about the processes in the system, were also consulted. The level of model sensitivity was derived by interpreting the model structure and components, interviews with experts, and model runs with EmiStatR. The main task to quantify input uncertainty is
<figure><img src="image_194.png" /><figcaption>Figure 2: The Haute-Sûre catchment and villages of the sewer system. Location of CSO structures and precipitation measurement stations.</figcaption></figure>

to define the probability distribution function (pdf) that represents the uncertainty of the variable chosen. The uncertainties of selected model inputs were characterized with pdfs following Heuvelink et al. (2007).

Measurement campaigns were done in Goesdorf from 28th April to 24th June 2011, in Kaundorf from 22nd June to 18th August in 2010 and from 20th July to 5th August in 2011, and in Nocher-Route from 18th November 2010 to 27th April 2011.

## 2.4 Uncertainty propagation

After the definition of model input uncertainties a Monte Carlo simulation was done to propagate input uncertainty to model output. The Monte Carlo method runs the EmiStatR model repeatedly, each time using different model input values, sampled from their probability distribution. The method thus consists of the following steps:

1. Repeat N times:

(a) Generate a set of realisations of the uncertain model inputs

(b) For this set of realisations, run the model and store the output

2. Compute and store sample statistics from the *N* model outputs.

Here, *N* is the number of Monte Carlo runs, i.e. the Monte Carlo sample size. Common sample statistics that measure the uncertainty are the standard deviation and the width of prediction intervals, which can be easily calculated from the *N* Monte Carlo outputs.
<table><caption>Table 1: General input data of the EmiStatR test model.</caption><thead><tr><th>Variable</th><th>Value</th></tr></thead><tbody><tr><td><strong>Wastewater</strong></td><td></td></tr><tr><td>Water consumption, <strong>qs</strong>, [l/(PE · d)]</td><td>150</td></tr><tr><td>Pollution COD, <strong>CODs</strong>, [g/(PE · d)]</td><td>104</td></tr><tr><td>Pollution NH<sub>4</sub>, <strong>NH4s</strong>, [g/(PE · d)]</td><td>4.7</td></tr><tr><td><strong>Infiltration water</strong></td><td></td></tr><tr><td>Inflow, <strong>qf</strong>, [l/(s · ha)]</td><td>0.05</td></tr><tr><td>Pollution COD, <strong>CODf</strong>, [g/(PE · d)]</td><td>0</td></tr><tr><td>Pollution NH<sub>4</sub>, <strong>NH4f</strong>, [g/(PE · d)]</td><td>0</td></tr><tr><td><strong>Rainwater</strong></td><td></td></tr><tr><td>Precipitation time series, <strong>P</strong> [mm/min]</td><td>P1</td></tr><tr><td>Pollution COD, <strong>CODr</strong>, [mg/l]</td><td>80</td></tr><tr><td>Pollution NH<sub>4</sub>, <strong>NH4r</strong> [mg/l]</td><td>0</td></tr><tr><td><strong>Storm water runoff</strong></td><td></td></tr><tr><td>Real flow time in the catchment, <strong>tf</strong> [min]</td><td>20</td></tr></tbody></table>

# III RESULTS

## 3.1 Simulations of the CSO structures

Deterministic simulations of the EmiStatR model for the eight CSO structures were performed to compare these with the outcomes of the uncertainty analysis. The deterministic input variables of the simulations are given in Table 1. The individual deterministic input variables for each CSO structure are presented in Table 2. Nine output variables were evaluated. Three out of nine are related with water quantity: volume of filling of the CSO tank, VTank [m³]; overflow volume, V<sub>Ov</sub> [m³]; and overflow flow, Q<sub>Ov</sub> [l/s]. Three other output variables are related with water quality in terms of COD: overflow COD load B<sub>COD,Ov</sub> [kg]; overflow COD average concentration, C<sub>COD,Ov,Av</sub> [mg/l]; overflow COD maximum concentration, C<sub>COD,Ov,Max</sub> [mg/l]. The remaining three output variables are also related with water quality but in terms of NH<sub>4</sub>: overflow NH<sub>4</sub> load, B<sub>NH4,Ov</sub> [kg]; overflow NH<sub>4</sub> average concentration, C<sub>NH4,Ov,Av</sub> [mg/l]; and overflow NH<sub>4</sub> maximum concentration C<sub>NH4,Ov,Max</sub> [mg/l].

## 3.2 Model input uncertainty assessment

To assess the sensitivity of the model output to small changes in the model input, simulations with EmiStatR for each input variable with an increment and decrement of 10% with relation to the deterministic values (Tables 1 and 2) were performed. Basically, the variables Aimp, Qd, and V are the most sensitive variables of water quantity variables (VTank, V<sub>Ov</sub> and Q<sub>Ov</sub>). Regarding water quality in terms of COD the input variables CODs, CODr, Aimp, Qd, V, and P have greatest impact on output COD variables, whereas input variables qs, NH4s, qf, NH4r, Aimp, pe, Qd, V and P have the greatest impact on output NH<sub>4</sub> variables.

After evaluation of the model output sensitivity to input variables, and taking into account the degree of uncertainties of each input, we selected three input variables to be included in the uncertainty analysis: CODs, NH4s, and CODr. We leave the uncertainty propagation
<table><thead><tr><th>Variable</th><th>BUD</th><th>DAH</th><th>ESO</th><th>GOE</th><th>HEI</th><th>KAU</th><th>NOC</th><th>NOR</th></tr></thead><tbody><tr><td colspan="9"><em>Sub-catchment data</em></td></tr><tr><td>Land use [-]<sup>1</sup></td><td>R/I</td><td>R/I</td><td>R/I</td><td>R/I</td><td>R/I</td><td>R/I</td><td>R/I</td><td>R/I</td></tr><tr><td>Total area [ha]</td><td>16.5</td><td>16.2</td><td>22.2</td><td>16.5</td><td>30.0</td><td>22.0</td><td>16.2</td><td>18.6</td></tr><tr><td>Impervious area [ha]</td><td>4.9</td><td>5.0</td><td>11.1</td><td>7.6</td><td>11.0</td><td>11.0</td><td>3.4</td><td>4.3</td></tr><tr><td>Population equivalents [PE]</td><td>289</td><td>460</td><td>705</td><td>611</td><td>676</td><td>358</td><td>260</td><td>326</td></tr><tr><td>Flow time structure [min]</td><td>10</td><td>10</td><td>7</td><td>10</td><td>6</td><td>10</td><td>10</td><td>10</td></tr><tr><td colspan="9"><em>Structure data</em></td></tr><tr><td>Throttled outflow [l/s]</td><td>6</td><td>4</td><td>12</td><td>9</td><td>11</td><td>9</td><td>2</td><td>4</td></tr><tr><td>Volume [m<sup>3</sup>]</td><td>90</td><td>270</td><td>330</td><td>190</td><td>220</td><td>180</td><td>166</td><td>157</td></tr></tbody></table>

¹R = residential, I = industrial

<caption>Table 2: General input variables of the CSO structures of the EmiStatR test model.</caption>

analysis of *Aimp* and *P*, which also score high on the sensitivity-uncertainty ranking, for future work.

The field measurements were the basis to characterise input uncertainty of *CODs* and *NH4s*. Samples of COD and NH₄ in mg/l (91 in total for each variable) were analysed in the dry weather flow produced in the villages of Goesdorf, Kaundorf and Nocher-Route. An average wastewater amount was calculated for Goesdorf (153 l/(PE·d)), Kaundorf (112 l/(PE·d)) and Nocher-Route (94.3 l/(PE·d)). Table 3 presents the summary statistics of the dry weather flow measurements of COD and NH4 and the correspondent value of *CODs* and *NH4s*. We use italic font for *CODs* and *NH4s* to emphasize that these are input variables of the model.

<table><caption>Table 3: Summary statistics of dry weather flow measurements for *CODs* and *NH4s* characterization.</caption><thead><tr><th></th><th>COD<br>[mg/l]</th><th>CODs<br>[g/(PE·d)]</th><th>NH4<br>[mg/l]</th><th>NH4s<br>[g/(PE·d)]</th></tr></thead><tbody><tr><th>Min</th><td>62</td><td>7</td><td>16.1</td><td>1.7</td></tr><tr><th>Mean</th><td>926</td><td>104</td><td>44.4</td><td>4.7</td></tr><tr><th>Max</th><td>3454</td><td>528</td><td>81.2</td><td>10.8</td></tr><tr><th>St. deviation</th><td>632</td><td>88</td><td>18.6</td><td>1.9</td></tr></tbody></table>

Regarding *CODr*, no field measurements were available. Thus, expert judgement and values from the literature were used to characterise input uncertainty in *CODr*. A mean value of 80 mg/l and a standard deviation of 90 mg/l were used. For all three input variables, *CODs*, *NH4s*, *CODr*, we proposed a normal distribution to characterise input uncertainty with mean and standard deviation for each variable as defined above. In order to avoid negative values, we first transformed the variables by taking their natural logarithm.

## 3.3 Uncertainty propagation

We made a deterministic run of the model with values of the input variables as presented in Tables 1 and 2. Additionally, we performed 1,500 Monte Carlo simulations allowing *CODs*, *NH4s*, and *CODr* as stochastic input variables with characteristics as defined in the previous section. In this way the total uncertainty of output variables due to input uncertainty was calculated. Figure 3 shows an example of model output and the total uncertainty band of 5 and 95 percentile for overflow concentration of COD (centre inset) and overflow concentration of NH₄
(bottom inset). In case of a rain event that produces a CSO, the uncertainty in the model output is quite large. Also, there is a systematic difference between the deterministic and the median run. The latter is always slightly above the deterministic run.

Next the contributions of each one of the input variables to the total uncertainty was calculated as well, by calculating the difference between the total uncertainty and the uncertainty obtained in the stochastic simulation of the other two variables. For instance, the uncertainty contribution of CODs was calculated as the total uncertainty minus the uncertainty of the simulations running only NH4s and CODr in stochastic mode. Therefore, 4,500 additional Monte Carlo simulations were done to calculate the uncertainty contribution of the three input variables.

The results indicate that there is no uncertainty contribution of CODs, NH4s and CODr to output variables VTank, VOv, and QOv. CODs and CODr contribute to uncertainty in BCOD, Ov and CCOD, Ov, Av. CODr has the most important uncertainty contribution to BCOD, Ov (99.6%) and CCOD, Ov, Av (90.4%). Finally, NH4s contributes 100% of the uncertainty of BNH4, Ov and CNH4, Ov, Av.

## IV CONCLUSIONS AND FUTURE WORK

We applied an uncertainty analysis of the simplified urban drainage model EmiStatR. A characterisation of the input uncertainty of the three main input variables that control output uncertainty in water quality variables was done. We found that the uncertainty in the load of COD per capita per day in the sewage (CODs) and the concentration of COD in runoff (CODr) contribute to the uncertainty of the output variables: overflow COD load and overflow COD average concentration. CODr has the most important uncertainty contribution in the load and concentration of COD. The load of NH4 per capita per day in the sewage (NH4s) contributes totally in the uncertainty of overflow NH4 load and overflow COD average concentration.

Rainfall is one of the most important drivers in the definition of uncertainty of output variables as load and concentration of COD and NH4. Therefore a detailed analysis should be done to quantify the contributions of rainfall input uncertainty in the total uncertainty. The rainfall uncertainty will be addressed in follow-up research as well as uncertainty due to impervious areas.

A further development of EmiStatR allows to analyse spatial model inputs. In this sense, a semi-distributed modelling framework can be implemented and several sub-catchments of the system modelled simultaneously, taking into account the inherent spatial variability of the inputs as rainfall and impervious areas and propagate the uncertainty of such inputs through the model. This also will be worked out in future research.

## ACKNOWLEDGEMENTS

The work presented is part of the QUICS (Quantifying Uncertainty in Integrated Catchment Studies) project. This project has received funding from the European Unions Seventh Framework Programme for research, technological development and demonstration under grant agreement No. 607000 and the Luxembourg Institute of Science and Technology.
<figure><img src="image_195.png" /><figcaption>Figure 3: Precipitation, and overflow COD and NH₄ concentrations in deterministic and uncertainty outcomes. Location: Goesdorf. From 14 to 18 August 2010. Inset centre and bottom: the blue line indicates the median, the black dotted line indicates the deterministic run, and the gray area indicates the 90 per cent uncertainty band.</figcaption></figure>

## References

Bach P. M., Rauch W., Mikkelsen P. S., McCarthy D. T., Deletic A. (2014). A critical review of integrated urban water modelling - Urban drainage and beyond. *Environmental Modelling & Software* 54, 88–107.

Brown J. D. (2004). Knowledge, uncertainty and physical geography: Towards the development of methodologies for questioning belief. *Transactions of the Institute of British Geographers* 29(3), 367–381.

Deletic A., Dotto C., McCarthy D., Kleidorfer M., Freni G., Mannina G., Uhl M., Henrichs M., Fletcher T., Rauch W., Bertrand-Krajewski J., Tait S. (2012). Assessing uncertainties in urban drainage models. *Physics and Chemistry of the Earth* 42-44, 3–10.

Grigoriu M. (2012). *Stochastic Systems: Uncertainty Quantization and Propagation*. London: Springer-Verlag.

Heuvelink G. B. M., Brown J. D., van Loon E. E. (2007). A probabilistic framework for representing and simulating uncertain environmental variables. *International Journal of Geographical Information Science* 21(5), 497–513.

Klepiszewski K., Seiffert S. (2013). Statistische Erfassung von Entlastungsbauwerken der Mischwasserbehandlung im Einzugsgebiet der Chiers. MIGR EMISTAT-MW. Technical report, TUDOR Centre de Ressources des Technologies pour l'Environnement, Luxembourg.

Mckay M. D., Beckman R. J., Conover W. J. (1979). A Comparison of Three Methods for Selecting Values of