Input Variables in the Analysis of Output from a Computer Code. *American Statistical Association and the American Society for Quality*.

Minasny B., McBratney A. B. (2006). A conditioned Latin hypercube method for sampling in the presence of ancillary information. *Computers and Geosciences* 32(9), 1378–1388.

Mitchell V., Duncan H., Inman M., Rahilly M., Stewart J., Vieritz A., Holt P., Grant A., Fletcher T., Coleman J., Maheepala S., Sharma A., Deletic A., Breen P. (2007). State of the art review of integrated urban water models. In *Novatech 2007*, pp. 507–514.

Nash J. E., Sutcliffe J. E. (1970). River flow forecasting through conceptual models. Part 1 - a discussion of principles. *Journal of Hydrology (Amsterdam)* 10, 282–290.

Neumann M. B. (2007). *Uncertainty Analysis for Performance Evaluation and Design of Urban Water Infrastructure*. Ph. D. thesis, Swiss Federal Institute of Technology, ETH Zurich.

Nol L., Heuvelink G. B. M., Veldkamp a., de Vries W., Kros J. (2010). Uncertainty propagation analysis of an N2O emission model at the plot and landscape scale. *Geoderma* 159(1-2), 9–23.

Refsgaard J. C., van der Sluijs J. P., Højberg A. L., Vanrolleghem P. A. (2007). Uncertainty in the environmental modelling process - A framework and guidance. *Environmental Modelling and Software* 22(11), 1543–1556.

Sanitary-District (2015, May). *Combined Sewer Overflow*. Internet.

Schellart A. N. A., Tait S. J., Ashley R. M. (2010). Towards quantification of uncertainty in predicting water quality failures in integrated catchment model studies. *Water Research* 44(13), 3893–3904.

van Keur P., Henriksen H. J., Refsgaard J. C., Brugnach M., Pahl-Wostl C., Dewulf A., Buiteveld H. (2008). Identification of major sources of uncertainty in current IWRM practice. Illustrated for the Rhine Basin. *Water Resources Management* 22(11), 1677–1708.

Walker W., Harremoës P., Rotmans J., van der Sluijs J., van Asselt M., Janssen P., Krayer von Krauss M. (2003, March). Defining Uncertainty: A Conceptual Basis for Uncertainty Management in Model-Based Decision Support. *Integrated Assessment* 4(1), 5–17.

Zoppou C. (2001). Review of urban storm water models. *Environmental Modelling and Software* 16(3), 195–231.
# How can big data be used to reduce uncertainty in stormwater modelling?

Nanée Chahinian*<fnref n="1" />, Anne-Laure Piat-Marchand<fnref n="2" />, Sandra Bringay<fnref n="2" />, Maguelonne Teisseire<fnref n="3" />,
Elodie Boulogne<fnref n="4" />, Laurent Deruelle<fnref n="5" />, Mustapha Derras<fnref n="5" />, Carole Delenne<fnref n="4" />

*Corresponding author: nanee.chahinian@ird.fr

Buried utility networks are often not fully documented and inaccurate positioning of network elements may impact hydraulic modelling applications. This study aims to put forward a methodology to extract network characteristics from data that is posted on the web or is available through sector specific databases, using heterogeneous text scrapping methods.

## I INTRODUCTION

Urban growth is an ongoing trend and one of its direct consequences is the development of buried utility networks. With growing needs among consumers, new networks are being installed and more underground space is being occupied. Locating these networks is becoming a challenging task. Mispositioning of utility networks is an important problem for both industrialised and developing countries and will worsen as cities expand and their networks increase in size and complexity (Jamil et al., 2012; Metje et al., 2007). Over the past century, it was common practice for public service providers to install, operate and repair their networks separately (Rogers et al., 2012). Now local authorities are confronted with the task of combining data produced by different parties, having distinct formats, variable precision and granularity (Chen and Cohn, 2011). Although in certain countries contractors are now obliged by law to position all buried networks within set precision ranges, finding data related to older network branches is a cumbersome task. Once located these data are often unavailable at the desired precision or are prone to errors or omissions. Since the mid 90's, an increasing volume of data is posted on the web or is available through sector specific databases. An alternative and complementary approach to field surveys would be to track down archived data by using new methods of heterogeneous texts scrapping. Hence the objective of this work is to assess whether big data can be used to reduce uncertainty in stormwater modelling. The big data used in this study is original information scrapped from the web such as calls for tenders, newspaper articles, consumer complaints, etc. Information extracted with text mining techniques such as used in Kergosien et al. (2015) are particularly interesting to confirm or infirm the position, the depth, the material of buried network elements. Call for tenders often include additional technical descriptions (slope, local constraints, junctions, etc.) which may be used to build an attribute table, that contains the characteristics of the underlying pipes.

This study is a part of a global project which aims to recreate a stormwater and a sewage network in settings where no accurate information regarding the position or characteristics of buried utility networks is available. A previous study (Pasquet et al., 2016) put forward a methodology to detect manhole covers and inlet grates from aerial images. The results were encouraging but dealt only with the correct detection of the objects i.e. no information regarding the underlying

<footnote n="1">¹IRD, HSM (CNRS, IRD, UM), France</footnote>

<footnote n="2">²Univ. Montpellier 3, LIRMM (CNRS, UM), France</footnote>

<footnote n="3">³TETIS, IRSTEA, France</footnote>

<footnote n="4">⁴Univ. Montpellier, HSM (CNRS, IRD, UM), France</footnote>

<footnote n="5">⁵Berger-Levrault, France</footnote>
network was supplied. Positional errors of remotely sensed manhole covers and grates could
lead to errors on rim elevation, which in turn impacts the pipe length and slope. Detection errors
such as i) objects masked in the aerial photograph by trees or shade, ii) spots wrongly assumed
to be manhole covers or grates, or even iii) visible but out of service elements, may strongly
impact the shape of the network. The pipe characteristics (material, diameter and roughness)
may also affect the output hydrographs. Consequently, it is important to assess the sensitivity of
the hydraulic model's results to the network descriptors in order to determine the most relevant
parameters, and the precision required when determining their values.

The paper is structured as follows. Section II describes the materials and methods including
study site, the tools and methodologies used for the sensitivity analysis and the textual analysis.
Section III and section IV present the main results of our study and the new perspectives that
this preliminary analysis opens.

## II MATERIALS AND METHODS

### 2.1 Study site and data

The study site is located in the "Hôpitaux Facultés" district of the city of Montpellier (Southern France), where heavy rains often cause over flooding of stormwater drains. It is an extension of the University of Montpellier's main campus which has already been used to test the manhole cover detection techniques proposed in a previous study (Pasquet et al., 2016).

The database consists of digital orthophotos at 16 cm resolution available through the city's
open data platform; digital data on the wastewater network extracted from the GIS of Montpel-
lier Métropole, the inter-municipality agency of the metropolitan area of Montpellier and data
on the stormwater network supplied in digital format by the Municipality of Montpellier. A
quick display of the three data layers showed that the manhole covers and grates located on the
university campus were not reported in the city's GIS systems. A special request was made to
the university's assets' management department which provided PDF files of the network maps
that were georeferenced in QGis. Manhole covers and grated were digitised as point data, the
water pipes as lines and all the descriptive data written in the PDF were entered manually in the
attribute table.

## 2.2 Methodology

### 2.2.1 Hydraulic modelling and sensitivity analysis

The network built in the previous step was imported into the PCSWMM (Computational Hy-
draulics International) software for stormwater modelling. PCSWMM is built on the US EPA's
SWMM5 engine and uses the Saint-Venant (1871) equations to model water flow.

The upstream catchment characteristics are determined based on the Municipality's dataset. The catchment is divided into three sub-catchments. Two synthetic storms are generated to run the model: a 2-hour duration Chicago design storm (Keifer and Chu, 1957) is used as a reference rainfall event and a Desbordes (Desbordes, 1974) double-triangle design storm with a 5 year return period is used to account for more intense events. No model calibration is undertaken. In order to determine the most important information to be searched for by data mining, sensitivity analyses are undertaken for the following elements:

1. **Slope.** In the framework of this project, rim elevation may be easily determined via remote sensing data (DEMs, airborne Lidar, etc). However, this does not mean that the inlet elevation will be known. If no further information is provided, pipe slopes should
deduced from rim elevation. Two tests are thus carried out to determine slope values: i) a constant slope is applied on the entire network; ii) the pipe slope is assumed to be equal to the terrain slope with a zero value in case of a counter slope. Figure 1 shows the profiles corresponding to the PCSWMM simulations (see section III).

2. **Network shape.** To test the influence of the missing elements, such as grates that are invisible on orthophotos or that are not reported in databases, the network shape is oversimplified (Figure 2).

3. **Roughness.** Pipes are assumed to be made of either PVC or concrete.

4. **Pipe diameter.** The diameter values are either increased or decreased by 30 and 50%.

<figure><img src="image_196.png" /><figcaption>Figure 1: Recreating slope values.</figcaption></figure>

<figure><img src="image_197.png" /><figcaption>Figure 2: Network reconstruction. a) Real network b) Simplified network</figcaption></figure>

### 2.2.2 Textual Analytics

Our process explore textual resources scrapped from the web to confirm or infirm information obtained from hydrologic modelling and sensitivity analysis. Our process is divided into 3 steps as shown in Figure 3. The first step is to gather a set of relevant documents. Once this corpus is built, textual documents are analyzed in two phases. The second phase aims to identify thematic, spatial and temporal features. The third step is dedicated to the linkage of these features.

For the first step, we manually gathered 30 documents identified as relevant. This set of documents is composed of one local planning map, one intervention report, three annual reports of public service concession holders, three technical notes, two prevention plans against natural floods, one scientific paper, one news article, two commercial brochures, four sewerage master plans, two municipal newsletters, one public inquiry file, one file of administrative acts, three study reports, a list of service delivery's contracts concluded, one convention document, three technical specifications. Six documents were identified as irrelevant.

For the second step, we used two different methods so as to identify three types of expressions: thematic, temporal and spatial, in our set of documents.

Our chain of treatments which aims to recognize spatial and thematic information is based on Unitex<fnref n="1" />, which is a free collection of programs dedicated to corpus processing based on linguist-

<footnote n="1">¹Unitex has been developed at LADL (Laboratoire d'Automatique Documentaire et Linguistique) http://www-igm.univ-mlv.fr/ unitex/</footnote>
tic resources. Resources used by the software consist of: i) Dictionaries: description of simple and compound words with their lemmas and a set of grammatical and semantic codes; ii) Grammars: representations, as graphs, of linguistic phenomena on the basis of finite state automata; iii) Lexicon-grammar tables: matrices describing syntactical properties of some words in which rows correspond to verbs and columns to syntactical properties. We applied the transducers' cascade "CasEN_Quaero" which has been created in order to recognize named entities (Maurel et al., 2011). A transducer is a set of boxes which contains symbols (e.g. characters' sequences, a syntactical code) that we aim to recognize. Each graph in a cascade is applied, one after the other, to the corpus which has been modified by the previous graph (Friburger, 2002).

The recognition of temporal features was carried out using Heideltime<fnref n="2" />. Heideltime is a rule-based tagger that retrieves temporal informations (absolute and relative dates, durations) from texts. Expressions identified are normalized according to the TIMEX3 annotation standard<fnref n="3" /> which describe temporal expressions, dates...).

Finally, the second phase of the process is then implemented as follows (Fig. 3):

* **Thematic identification:** we established a list of words related to the field of interest, composed of 77 words which describe objects such as manhole covers and stormwater or wastewater networks, and 62 terms that depict their characteristics (e.g. diameter, depth, flow) and measurement units (e.g. m 3/s, Eq-hab, m). We completed this list with terms extracted from two lexicons (UNESCO, 1992; Andréassian et al., 2002) and one reference work (Bourrier, 2008). This lexicon covers 1% of our corpus (percentage of annotated documents in the corpus). The automatic extraction of synonyms from 7 websites<fnref n="4" /> provided a list of 1721 terms which had been projected on texts retrieved in step 1. Finally, 683 appeared in the corpus, thus expanding the initial list to 344%. The lexicon extended this way covers 4% of our set of documents.

* **Spatial identification:** we used resources embedded in the CasEN_Quaero cascade of the Unitex software in order to recognize spatial information (e.g. "au nord de la route R12 allant de Montpellier à Lunel", "le boulevard du Languedoc").

* **Temporal identification:** we chose to use Heideltime so as to detect time expressions (e.g. "l'appel d'offre signé du 12 mai", "du 03/06/14 au 21/07/14").

The chain of treatments was evaluated by one annotator on a sample of 3 documents which contains 21 904 tokens. Thus, 285 spatial features, 213 temporal expressions and 346 thematic features were identified.

The third phase, corresponding to the linkage operation, will link the spatial and temporal features obtained in the previous phase with the thematic marking (Fig. 3). Depending on their degree of proximity (in the document, in the section, in the paragraph or in the sentence), spatial and temporal markings will be retained or not.

## III RESULTS

When no information is available on the inlet elevations or pipe slopes their values should be deduced from rim elevation data. The results obtained with the two approaches for this estimation are presented in Fig. 4, for the moderate and intense rainfall event. It can be seen that when assuming similar ground and underground slopes, the model underestimates both flood

<footnote n="2">²http://dbs.ifi.uni-heidelberg.de/index.php?id=129</footnote>

<footnote n="3">³http://www.timeml.org/tempeval2/tempeval2-trial/guidelines/timex3guidelines-072009.pdf</footnote>

<footnote n="4">⁴ http://www.cnrtl.fr, http://www.synonymo.fr, http://www.crisco.unicaen.fr, http://www.babla.fr, http://dico.isc.cnrs.fr, http://dictionary.reverso.net, http://dictionary.sensagent.com</footnote>
<figure><img src="image_198.png" /><figcaption>Figure 3: Global process of textual analysis.</figcaption></figure>

volume and peakflow. Indeed, imposing a null gradient for counter slopes yields overflows and thus diminishes the routed volume. For these events using a mean, constant slope does not induce great differences on the simulated hydrographs: changing the slope by up to 43% yields comparable results. The two values tested in this case correspond to either the mean slopes given by the real database or to the difference between upstream and downstream rim elevations.

The results are further confirmed for the high intensity rainfall event. However, slight differences can be seen based on the selected slope value; one yields a better estimation of the flood volume while the other one's results are closer to the flood dynamics. This raises questions about the slope value that should be used for network reconstruction. Although French technical guides recommend using a 0.01 m/m slope when designing stormflow networks, it seems that common practice, especially when extending existing networks, is to set the slope according to the input and output gradients. For a given pipe burial depth, this would imply dividing the network into distinct sections to avoid counter slopes or unrealistic elevation values.

<figure><img src="image_199.png" /><figcaption>Figure 4: Sensitivity to slope for an average event and an intense event. The slope 0.015 correspond to the mean of all the real slope values, and 0.0085 is directly the slope between upstream and outlet.</figcaption></figure>
Regarding the network's shape (Fig. 5) no real influence can be seen when modelling moderate intensity rainfall events. However, for heavier rainfall events, missing out inlet grates and consequent branching yields faster transfer times and upstream overflowing (Fig. 2).

<figure><img src="image_200.png" /><figcaption>Figure 5: Sensitivity to missing inlet grates</figcaption></figure>

For moderate rainfall, the roughness of the pipe material does not impact the output hydrograph. In case of high intensity events, high roughness may reduce flow velocity and yield greater over-flooding. However, roughness is often a calibrated parameter as its value cannot be easily determined by direct measurement. Consequently, lack of information on roughness can easily be overcome.

Finally, pipe diameter has nearly no influence on the output discharge as long as there is no transition from open-channel to closed-conduit flow. A significant reduction in pipe diameter will lead to over flooding.

Table 1 gives a ranking of the tested parameters in terms of impact on simulated flow hydrographs. Elevation and slope are the parameters which need to be estimated with the best accuracy. The aim of the textual analysis would be to identify any mention of these parameters and corresponding values in big data.

<table><caption>Table 1: Parameter ranking based on the sensitivity analysis</caption><thead><tr><td>Rank</td><td>Parameter</td><td>Justification</td></tr></thead><tbody><tr><td>1</td><td>Pipe slope</td><td>Faster transfer times if too steep</td></tr><tr><td>2</td><td>Pipe diameter</td><td>May cause overflowing if too small</td></tr><tr><td>3</td><td>Pipe roughness</td><td>Should be calibrated.</td></tr></tbody></table>

The first phase applied on the corpus allows the extraction of spatial, temporal and thematic features as illustrated in Figure 6. It can be seen that the program can accurately identify the name of the town ("ville de Vendôme"), the timeframe of the document ("depuis 1994"; "l'automne 2014"). By using the thematic pattern, the type of wastewater network ("réseau séparatif"-separate wastewater network) and the capacity of the wastewater treatment plant ("35000 équivalents/habitants") are also identified. However, the text also indicates that 98% of the wastewater network is routed through a separate network but the information is not picked up by the program. To overcome this, a new pattern could be added but the selection rules cannot be extended endlessly. A rate of "acceptable loss" should be set with the specialists. The quantitative results of the feature extraction are summarised in Table 2.

Regarding the results of the sensitivity analysis we carried out, with the exception of the rim elevation, the parameters correspond to terms which are too generic, with a wide variety of uses. "Pipe diameter" ("diamètre de conduite") for instance, was mentioned in reference to karst
<table><caption>Table 2: Quantitative Evaluation of the feature extraction.</caption><thead><tr><th></th><th>Spatial features</th><th>Temporal features</th><th>Thematic features</th></tr></thead><tbody><tr><th scope="row">False positive</th><td>35</td><td>21</td><td>14</td></tr><tr><th scope="row">Missing</th><td>124</td><td>22</td><td>34</td></tr><tr><th scope="row">Correct</th><td>285</td><td>213</td><td>346</td></tr><tr><th scope="row">Recall</th><td>0.70</td><td><strong>0.91</strong></td><td>0.91</td></tr><tr><th scope="row">Precision</th><td><strong>0.89</strong></td><td><strong>0.91</strong></td><td><strong>0.96</strong></td></tr><tr><th scope="row">F-Measure</th><td>0.78</td><td><strong>0.91</strong></td><td>0.94</td></tr></tbody></table>

aquifers, water supply and flooded caves. "Pipe diameter" pointed also to supply catalogs. The term "slope" was the one that needed most disambiguation. "The general slope of the terrain", "Average catchment slope" and "sloping" ("en pente" in French) are very frequently used. Pipe slope values which are not positioned closely to the words "pipe" and "slope" are often missed. As a general rule, quantified information related to the network's slope can not be easily detected because often there are many intermediate words between them and many expressions are used when reporting this values. This also applies to "Rim elevation". In order to obtain quantitative information about these parameters, it is necessary to overcome this problem. It would be necessary to investigate other annotation techniques such as the one described in (Berrahou et al., 2013).

<figure><img src="image_201.png" /><figcaption>Figure 6: Example of an annotated text</figcaption></figure>

## IV CONCLUSION AND PERSPECTIVES

A methodological framework was put forward to use big data in order to improve the accuracy of stormwater modelling networks. A sensitivity analysis was carried out to determine the most sensitive parameters and text mining techniques were used to extract information about them in big data scrapped from the web. The primary results are encouraging as spatial, temporal and thematic patterns are correctly identified. However, quantitative information about the network's elements such as elevation and slope can be missed sometimes. Future work will look into overcoming this problem by either setting new patterns or creating an index which will give the expert an indication on the presence of possibly interesting numerical data in a document. The final step will consist in building a model based solely on the information gathered through the big data. The uncertainty associated with each characteristic will be assessed and hydraulic simulations using a classical modelling software will be carried out to assess the uncertainty
transmission in each step of the process until the simulation of the output hydrograph. Other types of documents could also be considered such as photos or scanned plans with captions. Automatic scrapping from the web based on research engines such as Google or Yahoo and contextual analysis of snippets will also be explored (Opitz et al. (2014)).

## Acknowledgment

This study is part of the project "Cart'Eaux" funded by the European Regional Development Fund (ERDF).

## References

Andréassian V., Sarkissian V., Chelmiki W., Stanescu V., Moussa R. (2002). Cemagref.

Berrahou S. L., Buche P., Dibie-Barthélemy J., Roche M. (2013). How to extract unit of measure in scientific documents? In K. Liu, A. L. N. Fred, and J. Filipe (Eds.), *KDIR/KMIS 2013 - Proceedings of the International Conference on Knowledge Discovery and Information Retrieval and the International Conference on Knowledge Management and Information Sharing*, Vilamoura, Algarve, Portugal, 19 - 22 September, 2013, pp. 249–256. SciTePress.

Bourrier R. (2008). *Les réseaux d'assainissement Calculs - Applications - Perspectives*. Paris: Editions Tec et Doc/Lavoisier.

Chen H., Cohn A. (2011). Buried utility pipeline mapping based on multiple spatial data sources: A bayesian data fusion approach. In *IJCAI-11, Barcelona, Spain*, pp. 2411–2417.

Desbordes M. (1974). *Réflexions sur les méthodes de calcul des réseaux urbains d'assainissement pluvial*. Montpellier: PhD Thesis USTL.

Friburger N. (2002). *Reconnaissance automatique des noms propres : application à la classification automatique des textes journalistiques*. Tours: PhD Thesis University of Tours.

Jamil H., Z. N., Mohid Yussof M. (2012). Underground utility mapping and its challenges in Malaysia. In *FIG working week 2012. Knowing to manage the territory, protect the environment, evaluate the cultural heritage*, Rome, Italy, pp. 15.

Keifer J., Chu H. (1957). Synthetic storm patterns for drainage design. *Journal of Hydraulic's division* (83), 1–25.

Kergosien E., Alatrista-Salas H., Gaio M., Güttler F., Roche M., Teisseire M. (2015). When textual information becomes spatial information compatible with satellite images. In *KDIR*, pp. 301–306.

Maurel D., Friburger N., Antoine J., Eshkol-Taravella I., Nouvel D. (2011). Cascades autour de la reconnaissance des entités nommées. *TAL* (52), 69–96.

Metje N., Atkins P., Brennan M., Champan D., Lim H., Machell J., Muggleton J., Pennock S., Ratcliffe J., Redfern M., Rogers C., Saul A., Shan Q., Swingler S., Thomas A. (2007). Mapping the underworld: State of the art review. *Tunnelling and underground space technology* (22), 568–586.

Opitz T., Azé J., Bringay S., Joutard C., Lavergne C., Mollevi C. (2014). Breast cancer and quality of life: Medical information extraction from health forums. In C. Lovis, B. Séroussi, A. Hasman, L. Pape-Haugaard, O. Saka, and S. K. Andersen (Eds.), *e-Health - For Continuity of Care - Proceedings of MIE2014, the 25th European Medical Informatics Conference, Istanbul, Turkey, August 31 - September 3, 2014*, Volume 205 of *Studies in Health Technology and Informatics*, pp. 1070–1074. IOS Press.

Pasquet J., Desert T., Bartoli O., Chaumont M., Delenne C., Subsol G., Derras M., Chahinian N. (2016). Detection of manhole covers in high-resolution aerial images of urban areas by combining two methods. *IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing* 9(5), 1802–1807.

Rogers C., Hao T., Costello S., Burrow M., Metje N., Chapman D., ... Saul A. (2012). Condition assessment of the buried utility service infrastructure: a proposal for integration. *Tunnelling and Underground Space Technology* 28, 202–211.

UNESCO (1992). *Glossaire international d'Hydrologie*.
# Sensitivity analysis of spatio-temporal models describing nitrogen transfers, transformations and losses at the landscape scale.

Jordi Ferrer Savall¹*, Cyril Benhamou¹, Pierre Barbillon²,
Patrick Durand³, Marie-Luce Taupin⁴, Hervé Monod², Jean-Louis Drouet¹

¹ UMR ECOSYS, INRA, AgroParisTech, Université Paris-Saclay, 78850, Thiverval-Grignon, France

² UMR MIA-Paris, AgroParisTech, INRA, Université Paris-Saclay, 75005, Paris, France.

³ UMR SAS, INRA, Agrocampus Ouest. 84215, Rennes, France

⁴ UMR MaIAGE, INRA, Université Paris-Saclay, 78350 Jouy-en-Josas, France

*Corresponding author: jordi.ferrersavall@agroparistech.fr, jordi.ferrer.savall@gmail.com

## Abstract

Modelling complex systems such as agroecosystems often requires the quantification of a large num-
ber of input factors. Sensitivity analyses are useful to assess the impact of imput factors on model out-
comes and to determine the appropriate spatial and temporal resolution of models. Comprehensive
spatial and dynamic sensitivity analyses were applied to the model Nitroscape, a deterministic spa-
tially distributed model describing nitrogen transfers and transformations in rural landscapes, simulat-
ing five years of farm management in an intensive rural area of 3 km2. 11 global imput factors were
explored, characterizing spatial resolution of the simulations, the physical parameters of the landscape
and the agricultural practices. Their impact on 29 spatially-distributed model outcomes was assessed
through standard analysis of variance. Cluster analyses were applied to summarize the results of the
sensitivity analysis on the ensemble of model outcomes. The purpose of this paper is to present the
methods and preliminary results that we developed to carry out and to visualize comprehensive spatial
and dynamic sensitivity analyses of multiple input factors on multiple outcomes.

## Keywords

Please insert a few keywords here

Modeling nitrogen transfers, transformations and losses from processes integration at the landscape scale requires using complex models which are spatially and temporally explicit. Model predictions may vary considerably, depending on the initial conditions and input parameters. Therefore, it is important to characterize their robustness and variability through sensitivity analyses.

A variety of techniques have been developed to perform sensitivity analysis in spatially and dynamic models at different stages. Methods for exploring model inputs may range from the simplest one-factor-at-a-time screening techniques proposed by Morris (1991) to the exploration of representative combinations of factors in factorial designs, as described by Chen (2011). Recently, Moreau et al. (2013), proposed a method to take into account the spatial distribution and the granularity of input parameters in sensitivity analyses. Besides, Marrel et al. (2011) proposed a technique to extend variance-based methods well suited to analyse scalar outputs to dynamic and spatially explicit outcomes, rendering spatial maps of the sensitivity of each output to each model input. Finally, Ligmann-Zielinska (2013) showed that spatially explicit sensitivity analyses and analyses that aggregate outputs at different scales of description are complementary, but at the same time can lead to different decisions regarding input factor prioritization.

Here, we present an analysis of sensitivity of the NitroScape model developed by Duretz et al. (2011), with the specific aim to compare the impact of spatial granularity of the model to the impact of global parameters describing physical characteristics of the landscape and agricultural practices on spatially-distributed model outcomes. A central concern of the presented work is to provide tools for integrating the results of spatially
explicit sensitivity analyses applied on multiple model outcomes.

Nitroscape is a deterministic, spatially explicit and dynamic model describing the transfers and transformations of reactive forms of nitrogen (Nr) in rural landscapes of a few km² to a few tenth of km². It couples four modules characterizing processes of farm management, biotransformations and transfers of Nr by the atmospheric and hydrological pathways. It simulates Nr flows and losses within and between several landscape compartments: the atmosphere, the hydro-pedosphere (groundwater, water table and streams) and the terrestrial agroecosystems (livestock buildings, croplands, grasslands and semi-natural areas).

In order to evaluate the impact of global model inputs on the model outcomes, 11 parameters were selected. They characterized the spatial resolution of the model (size of the horizontal grid and depth of surface layers), the physical features of the landscape (soil characteristics) and the agricultural management (fertilization). A complete fractional factorial design with 11 categorical factors and 3 levels per factor was generated. The resulting plan comprised 243 units and the design resolution was 5, which allowed estimating main effects and two-factor interactions unconfounded.

The range of explored input values represented both the range of values observed in real settings and the degree of uncertainty in the measured or estimated values. Landscape structure (e.g. topography, land use), initial conditions (e.g. soil mineral content, agroecosystem characteristics) and forcing inputs (e.g. agricultural practices, meteorology) were common for all units in the factorial design. Their values were set from real measurements to represent an intensive rural area with mixed crops, pig farming and unmanaged ecosystems characterized by humid climatic conditions and little temperature contrasts.

Each run of the plan simulated nitrogen transfers, transformations and losses in a virtual landscape representing 300 ha of a catchment area with agricultural and farming activity, during 5 years, with daily outputs at the outlet and monthly outputs on the whole grid after an initialization period of two years. The impact of model inputs was evaluated on 29 model outcomes: 5 variables describing the catchment outflow, 9 spatially-distributed variables describing inter-compartment fluxes and 15 spatially-distributed variables describing the local state of the system.

NitroScape spatially-distributed output variables were aggregated to obtain either time-series describing spatially aggregated outcomes or maps of temporally aggregated outcomes. A principal component analysis (PCA) was applied to time series and spatial distributions to reduce their dimensionality. The influence of factors on each model outcome (expressed as aggregated data or in terms of its principal components) was explored through a standard analysis of variance (ANOVA) considering up to second-order interactions: the variability of each component of each outcome (day in time series, pixel in maps, or principal component), spread over the experimental design, was adscribed to every factor and two-factor interaction by partitioning the sum of squared deviations. The relative importance of factors and interactions was measured with sensitivity indices based on this partitioning. Finally, a cluster analysis and a PCA were applied to summarize the ensemble of results of the sensitivity analysis on the whole set of model outcomes, in order to identify outcomes with similar response to model inputs and to discriminate the effect of each factor and two-factor interaction on the ensemble of outcomes.

The methods here presented may be used to perform dynamic and spatial sensitivity analyses of models with multiple outcomes. The novel contributions of our proposal can be split into three aspects. Firstly, we simultaneously considered the spatial resolution of the model, the physical characteristics of the landscape and the agricultural practices in a single experimental design. This allows evaluating not only the main effect of the spatial resolution, but also its interactions with the other factors. Secondly, we aggregated either spatially or temporally each outcome and we carried out a global sensitivity analysis on the aggregated variables. This allows identifying patterns in the temporal sequences and in the spatial distributions of the sensitivity indexes that may arise from the model structure (seasonality and spatial structure of the landscape, respectively). And thirdly, we analysed the ensemble of results of the sensitivity analyses on multiple outcomes. This allows grouping together model outcomes that are mainly sensitive to the same inputs and identifying the factors that impact the ensemble of outcomes in a similar way.
## References

Chen H., Cheng C. (2011). Fractional Factorial Designs. *Design and Analysis of Experiments, Special Designs and Applications* 3, 299.

Duretz S., Drouet J.L., Durand P., Hutchings N.J., Theobald M.R., Salmon-Monviola J., Dragosits U., Maury O., Sutton M.A., Cellier P. (2011). NitroScape: a model to integrate nitrogen transfers and transformations in rural landscapes. *Environmental Pollution* 159 (11), 3162–70.

Ligmann-Zielinska A. (2013). Spatially-explicit sensitivity analysis of an agent-based model of land use change. *Int. J. Geographical Information Science* 27 (9), 1764-1781.

Marrel A., Iooss B., Jullien M., Laurent B., Volkova E. (2011). Global sensitivity analysis for models with spatially dependent outputs. *Environmetrics* 22 (3), 383-397.

Moreau P., Viaud V., Parnaudeau V., Salmon-Monviola J., Durand P. (2013). An approach for global sensitivity analysis of a complex environmental model to spatial inputs and parameters: a case study of an agro-hydrological model. *Environmental Modelling & Software* 47, 74-87.

Morris M.D. (1991). Factorial sampling plans for preliminary computational experiments. *Technometrics* 33 (2), 161-174.
# The transformed optimal transportation problem: sensitivity and segregation of the children-to-school constrained assignment in Lausanne

Théophile Emmanouilidis, Guillaume Guex, François Bavaud

University of Lausanne, Switzerland

*Corresponding author: Theophile.Emmanouilidis@unil.ch

**Abstract:** We present a fuzzy approach aiming at assigning origins (children home locations) to destinations (schools) while considering capacity and segregation constraints. Solving the classic instance of the optimal transportation problem generally leads to many equivalent solutions, and can be computationally demanding for large networks. Introducing an additional origin-destination mutual information term, measuring the assignment cohesiveness and penalizing hard assignments, allows faster computing of an unique solution. The algorithm is built on two freely adjustable parameters : the temperature T controlling the assignment cohesiveness, and the mixing parameter c, favoring social diversity at the destinations. The approach yields relevant results in school organisation and is currently used by the public administration of Lausanne to plan children-to-school assignment as well as to evaluate and simulate the development of the schooling infrastructures.

**Key-words:** optimal transport, regularization, fuzzy clustering, fuzzy logic, social mixing modelling, location-allocation uncertainty, spatial decision making, Lagrangian multipliers

## I INTRODUCTION

May it be from an urbanistic, a social or from a governance point of view, the evolution of cities is a major challenge of our contemporary societies. By giving the opportunity to analyze spatial and social configurations or attempting to simulate future ones, GIS cannot be overlooked in urban planning and management. In five years, the population of the city of Lausanne has grown from 134'700 to 140'570 inhabitants while the numbers in public schools have increased from 12'200 to 13'500 schoolchildren. Demographic rise and dynamics constitute spatially heterogeneous processes that have a direct impact on school organization and planning. They require the adaptation of the schooling infrastructures and may generate uncertainty in schoolchildren location-allocation. Nowadays most cities have divided their territory into multiple school districts to assign children to school. But sooner or later such districts become inaccurate as they define fixed boundaries in an evolving urban landscape. Beyond the demographic aspects, assigning children to public schools is a complex task as it has to satisfy political and public demands, as well as to meet legal requirements. To help to the administrative services in performing this task we have developed a fuzzy-logic approach for children to school assignment.

Models for schools and schoolchildren location-allocation have been investigated in operational research ever since the 60's (Clarke and Surkis (1968), Koenigsberg (1968), Heckman and Taylor (1969)): see the reviews of Caro et al. (2004), as well as Castillo-Lopez and Lopez-Ospina (2015) for recent years. All the models seek to allocate geography units to facilities, while considering spatial (journey to school distance threshold, closest assignment), logistic (capacity of schools) (Delmelle et al. (2014)), segregation and or economic constraints (Antunes and Peeters (2000)). Minimizing the average distance for given children locations and school capacities amounts to the celebrated optimal transportation problem (see e.g. Villani (2009)), possessing in general many equivalent solutions, computationally demanding for large networks.
Minimizing the social segregation at schools also possesses a multitude of solutions whenever the number of children locations exceeds the number of social categories, as in real situations.

The intrinsic degeneracy of the minimal distance and the minimal segregation problems prevents the existence of a unique and stable spatial partition, made of the union of school assignment basins. By contrast, introducing an additional origin-destination mutual information term, measuring the assignment *cohesiveness* and penalizing hard assignments, permits to make the solution unique. It furthermore allows the quick iterative computing of the minimal distance and the minimal segregation problems, as well as convex combinations of the latter, provided the penalization term is large enough.

The formalism and issues presently addressed bear obvious connections with iterative fitting, doubly constrained gravity flows, regularized optimal transportation, and model-based clustering with fixed mixture proportions. It involves in addition a presumably original term, the *social segregation at destination*, and appears to be *analytically tractable* (iterative determination of the unique solution, duality theory and sensitivity analysis). All the ingredients can be combined into an unique objective functional, the *free energy*, depending upon two freely adjustable parameters, the *temperature T* controlling the penalization of hard assignments, and the *mixing parameter c*, favoring social diversity at the destinations (section II). Toy and real applications (section III) illustrate the formalism together with its arguably numerous virtues, and the conclusion (section IV) emphasizes its connections to sensitivity analysis and robustness.

## II FORMALISM

Formalizing the above considerations necessitates the introduction of three categorical variables, namely $O = “child domicile”$ with $n$ modalities, $G = “school locations”$ with $m$ modalities, and $A = “social characteristics at origin”$ with $p$ modalities.

### 2.1 Notations and definitions

Consider the $n \times m$ matrix $N = (n_{ig})$ counting the number of children living at $i$ and assigned at school $g$. Row margins $n_{i\bullet}$ (where “$\bullet$” denotes the sum over the replaced index) count the number of children at $i$. Column margins $n_{\bullet g}$ yield the absolute capacity of school $g$.

Relative $n \times m$ couplings $P = (p_{ig})$ with $p_{ig} = n_{ig}/n_{\bullet\bullet} = p(i, g)$ constitute the joint probability distribution of $O$ and $G$. Its margins are denoted by $f_i = n_{i\bullet}/n_{\bullet\bullet} = p(i)$ (origin weights) and $\rho_g = n_{\bullet g}/n_{\bullet\bullet} = p(g)$ (destination weights). The assignment *cohesiveness* is measured by the origin-destination mutual entropy, that is

$$I[P] := I(O: G) = \sum_{ig} p_{ig} \ln \frac{p_{ig}}{f_i \rho_g} = H(O) + H(G) - H(O, G) = H(G) - H(G|O) \quad (1)$$

where $H()$ denotes the corresponding entropies. By construction, $I[P] \ge 0$, with equality iff $p_{ig} = f_i \rho_g$, that is iff $O$ and $G$ are independent.

Let $D = (d_{ig})$ denote the $n \times m$ matrix of domicile to school pedestrian distances. The transportation cost or *energy* is the average journey to school distance:

$$U[P] := \sum_{ig} p_{ig} d_{ig} \qquad (2)$$
Finally, consider a “social variable” A (such as nationality, ethnicity, social class or sex) with p modalities a = 1, ..., p, and let Y = (yᵢₐ) denote the given n × p matrix specifying the proportion of social type a at origin i, i.e. yᵢₐ = p(a|i) with yᵢₐ ≥ 0 and yᵢ· = 1. The joint distribution of social types and destinations is given by Q = (qᵣₐg) with

$$q_{ag} = \sum_{i} y_{ia} p_{ig} = p(a, g) \quad \text{with} \quad q_{a\bullet} = \sum_{i} f_{i} y_{ia} =: r_{a} \quad \text{and} \quad q_{\bullet g} = \rho_{g}.$$

$r_a = p(a)$ is the overall proportion of type a. Let us measure the social segregation at destinations by the mutual information between A and G, that is

$$S[P] := I(A:G) = \sum_{ag} q_{ag} \ln \frac{q_{ag}}{r_a \rho_g} \qquad (3)$$

By construction, $S[P] \ge 0$, with equality iff $q_{ag} = r_a \rho_g$, that is iff A and G are independent (absence of segregation at destinations).

## 2.2 Minimizing the free energy

The simultaneous control of the average distance, segregation and cohesiveness can be achieved by minimizing the *free energy*

$$F[P] := U[P] + c S[P] + T I[P] \qquad (4)$$

among all couplings P with given margins f and ρ. Here $T \ge 0$ is a freely adjustable *temperature* parameter, and $c \ge 0$ a freely adjustable *mixing* parameter. Optimal transportation problem corresponds to $T = c = 0$. Setting $T > 0$ favors soft assignment, and setting $c > 0$ diminishes social segregation at destinations. $c < 0$ can also be adopted, *favoring social segregation*, provided $|c|$ is not too large: further analysis demonstrates that $F[P]$ is convex for $c \ge 0$ (and hence possesses an unique minimizer), but ceases to be convex for $c < -T$.

The minimization of $F[P]$ is a standard exercise in calculus. Setting to zero its first derivative (including the Lagrange multipliers term insuring that marginal constraints are satisfied) yields the identities

$$p_{ig} = f_i \rho_g \exp(-\beta d_{ig}) \exp(-\beta c \ell_{ig}) \exp(\beta(\lambda_i + \mu_g)) \qquad (5a)$$

$$\exp(-\beta\lambda_i) = \sum_g \rho_g \exp(\beta\mu_g) \exp(-\beta d_{ig}) \exp(-\beta c \ell_{ig}) \qquad (5b)$$

$$\exp(-\beta\mu_g) = \sum_i f_i \exp(\beta\lambda_i) \exp(-\beta d_{ig}) \exp(-\beta c \ell_{ig}) \qquad (5c)$$

$$\text{where} \quad \ell_{ig}[P] := \sum_a y_{ia} \ln \frac{q_{ag}[P]}{r_a \rho_g} = \sum_a \frac{p(i a)}{p(i)} \ln \frac{p(ag)}{p(a)p(g)} \qquad (5d)$$

Here $\beta = 1/T$ is the inverse temperature or *coldness*, $\lambda_i$ and $\mu_g$ the Lagrange multipliers insuring the origin, respectively destination constraints. The quantity $\ell_{ig}$ (possibly negative) behaves as a “social distance” in the expression $d_{ig} + c \ell_{ig}$, penalizing the assignment $i \to g$ whenever the dominant social characteristics a at i tend to be overrepresented at g.

Direct substitution shows the minimum free energy to be $F = \sum_i f_i \lambda_i + \sum_g \rho_g \mu_g$, an expression generalizing Kantorovitch duality theory of optimal transportation (e.g. Villani (2009)), and justifying the interpretation of the multipliers $\lambda_i$ and $\mu_g$ as the unit embarkment cost at origin i, respectively disembarkment cost at destination g.
## 2.3 Iterative determination of the optimal coupling

Starting from some initial coupling such as $P^{(0)} = f\rho'$, and some destination multiplier value such as $\mu^{(0)} = 0$, a new coupling $P^{(1)}$ is obtained by successively applying (5d), (5b), (5c) and (5a). The process is then iterated. Convergence occurs provided the cohesiveness contribution is large enough compared to the contributions of both energy and social segregation. That is, the temperature $T$ must be large enough to allow efficient regularization. Numerical experiments with $\mathbb{R}$ (2.15.2) for the toy network of section 3.2 shows that convergence occurs (say after ca. 1'000 iterations) provided both the following conditions hold:

i) $T \ge 0.003 \cdot d_c$, where $d_c$ is the Hausdorff distance between the set (support) of sources $S = \{i | f_i > 0\}$ and the set of targets $T = \{g | \rho_g > 0\}$, that is $d_c = \max(d'_c, d''_c)$, where $d'_c = \max_i \min_g d_{ig}$ and $d''_c = \max_g \min_i d_{ig}$. Below that threshold, overflow occurs in (5b) or (5c).

ii) $T \ge |c|$, that is $c \in [-T, T]$. Outside that range, the iterative process may exhibit periodic or non-converging behavior (especially for small values of $T$), betraying instability of the fixed point (yet existing and unique for $T > 0$ and $c \ge 0$ since $F[P]$ is convex).

## 2.4 Boundaries and placement error

High temperature favors low cohesiveness, that is soft memberships (of children into schools) $z_{ig} = p_{ig}/f_i = p(g|i)$. The conditional entropy $H(G|i) = -\sum_g z_{ig} \ln z_{ig} \ge 0$ measures the uncertainty in the assignment of $i$. High values of $H(G|i)$ characterize origins $i$ lying at the boundary of two or more assignment basins, and help selecting domicile candidates for school reassignments, if needed (e.g. changes in $f$, $\rho$ or $Y$ from year to year).

This being said, even origins far from the assignment boundaries obey $H(G|i) > 0$ for $T > 0$ in view of the softness of $Z$. Attributing a given origin to different schools is not always possible nor desirable in real situations, so $Z$ might have to be hardened into $Z^* = (z_{ig}^*)$, where $z_{ig}^* = 1$ if $g = \arg \max_h z_{ih}$, and $z_{ig}^* = 0$ otherwise (ties are resolved at random). In general, $Z^*$ does not satisfy the capacity constraint, that is the quantity $\rho_g^* = \sum_i f_i z_{ig}^*$ generally differs from $\rho_g$: for $f$ and $\rho$ given, equation $\rho = fZ'$ possesses generally no hard solution $Z$. The relative placement error $\delta = \sum_g |\rho_g - \rho_g^*|$ is expected to increase with $T$ (figure 1). The absolute placement error $N\delta$ (where $N$ is the total number of children) counts the number of misattributed children (excess or deficit) by the hardened assignment.

# III ILLUSTRATIONS AND APPLICATIONS

## 3.1 Example A

$n = 600$ children of the same kind (no segregation considerations here) are uniformly placed on a $25 \times 24$ grid, and have to be assigned to $m = 6$ schools of identical capacity of 100 each, so that the average city-block distance ($L_1$ metric) is minimal (1). As expected, the fuzziness of the boundaries as well as the placement error diminishes with the coldness $\beta = 1/T$.

## 3.2 Example B

Consider the random placement, without overlapping, of $n = 40$ origins and $m = 5$ destinations, on a $13 \times 13$ grid endowed with the city-block distance, with uniform origin weights $f_i = 1/n$ and uniform destination weights $\rho_g = 1/m$. In addition, we focus on gender segregation ($p = 2$), and generate a gender gradient across the origins, such that the proportion of girls is maximum in the south-west corner, and minimum in the north-east corner. Figure (2) depicts the hardened memberships $Z^*$ for each origin $i$, and also the proportion of girls versus boys at domiciles and at schools.
<figure><img src="image_202.png" /><figcaption>Figure 1: Example A: entropic boundaries $H(G|i)$ between school assignment basins, home-to-school assignment based upon the hardened membership $z_{ig}^*$ (section 2.4), group sizes and resulting placement error $\delta$, for decreasing values of the temperatures $T = 1/\beta$.</figcaption></figure>

<figure><img src="image_203.png" /><figcaption>Figure 2: Example B. I: Hardened membership $z_{ig}^*$ (I), II: the proportion of girls versus boys at origins $p(a|i)$ (circles) and at destinations $p(a|g)$ (triangles). With $d_c = 7$, $T = 16$ and while penalizing gender segregation at the destinations ($c = 70 > 0$) or, on the opposite, favoring it ($c = -270 < 0$).</figcaption></figure>

<figure><img src="image_204.png" /><figcaption>Figure 3: Example B: at low temperature $T = 0.02$, the iterative solution converges for $c \in [-1, 4.2]$. Mixings $c \in [0, 3]$ increase the cohesiveness $I$ but leave $U$ and $S$ unchanged. Values $c > 3$ produce a sharp decrease in segregation $S$, and a sharp increase in energy $U$, depicting an effective anti-segregation regime together with an additional average distance cost.</figcaption></figure>
## 3.3 Children-to-school constrained assignment in Lausanne

In the optimal transportation limit, that is in absence of social segregation considerations (i.e. $c = 0$), and for low temperatures $T \to 0^+$, the procedure has been shown (Guex et al. (2016), Emmanouilidis (2016)) to permit mapping of school attendance areas together with their fuzzy boundaries. As shown in Figure 4 (left), large values of the (spatially smoothed<fnref n="1" />) assignment entropy $H(G|i)$ (section 2.4) highlight fuzzy boundaries between basins, where the origin-to-destination attribution is most uncertain, and where the potential for school organisation and planning is most promising (e.g. balancing class sizes, or splitting up neighborhoods between two or more schools in order to promote social mixing). Furthermore, groups of attendance areas can be aggregated to design a new school district pattern. Mapping Lagrangian multipliers provides an overview of embarkment costs $\lambda_i$ and disembarkment costs $\mu_g$ (Figure 4, right). Such a visualization can be especially useful for the comparison and evaluation of various potential scenarios involving changes of school capacities or new schooling infrastructures.

<figure><img src="image_205.png" /><figcaption>Figure 4: Resulting allocation plan for children entering secondary school in summer 2016 ($\beta = 0.2$). Left: hard assignment $Z^*$ generates fairly compact attendance areas around schools. Children located in dark areas (high entropy $H(G|i)$) are more likely to be assigned to two or more schools, than the ones living in white areas. Right: the lack of schools in the western and northern outskirts of the city generates high (red orange) values of $\lambda_i$ (spatially smoothed large embarkment costs). By contrast, children living in the center and in the south benefit from shorter journeys to school (green, yellow). The top-right school is too small to enroll all the surrounding children, which results in a large disembarkment cost $\mu_g$.</figcaption></figure>

### 3.3.1 Assignment with gender segregation constraint

Running the analysis on a similar dataset with adding the mixing parameter $c$ leads to the results exposed in Figure 5. We here consider the assignment of 973 children into 9 schools. Cases I and II relates respectively the results with high ($c = -700$) and low ($c = 300$) gender segregation, while case III is free of segregation constraint ($c = 0$). Despite the almost perfect balance between the number of boys (486) and girls (487), their repartition within the city is locally heterogeneous (background maps of line B). Allocating children with high segregation ($c = -700$) generates almost five unisex schools (maps I-B & C). Four of them are located in

<footnote n="1">¹in this paper, spatial interpolation is systematically performed through kriging (ArcGIS settings: type spherical, variable search radius, number of points:10).</footnote>
the city-center while the fifth, on the top, is surrounded by a majority of girls. Highest values of $c$ promotes, on the opposite, gender mix within schools (II-B & C). Results obtained with $c = 0$ are almost the same ((III-B & C), showing a good gender balance and suggesting that proportion of boys and girls around schools are almost equal. Accordingly, gender segregation here involves spatial mixing of overlapping hard assignments (I-A), while gender mixing creates a spatial partition of homogeneous connected groups (II & III-A).

Although the range of the assignment uncertainty ($H(G|i)$, line D) tends to narrow as $c$ decreases, highest values characterize the same city parts in the three maps. The Lagrangian multipliers (line E) exhibit similar trends in cases I and II, where children of the west, north and far east suburbs have the longest journey to school (highest $\lambda_i$). With $c = 300$ (II-E) embarkment costs in these areas are more smoothed. The disembarkment costs $\mu_g$ vary little between the three cases, suggesting that gender segregation has, for this dataset, little effect on the former.

<figure><img src="image_206.png" /><figcaption>Figure 5: Children to school assignment with gender segregation constraint. A large negative value of $c$ in (I) produces segregative allocation plan tending towards unisex schools. On the opposite, a large positive value of $c$ promotes gender mixing within the schools.</figcaption></figure>
# IV CONCLUSION

Introducing some amount of randomness and uncertainty in an otherwise deterministic problem (the optimal transportation) results in a unique solution, easily and quickly computable, and *robust* with respect to small variations of child locations and populations, or school locations and capacities: once more, the regularizing virtues of the origin-destination mutual entropy are confirmed. Regularization permits in addition the effective identification and visualization of borders, of costs, and of alternative near-optimal solutions: those benefits are essential (and much appreciated) for adapting assignments to changes (turn-over of children at the origins, modifications of the pedestrian network) and for future planning.

As demonstrated in the present paper, this approach can be extended to handle in addition the issue of *social segregation*, by the introduction of a new *mixing parameter* controlling for the segregation. The resulting thermodynamic-like formalism permits to perform tractable *sensitivity analyses* of various kinds, assessing the amount of change in the quantities of interest induced by a change in the network geometry, origin or destination distributions, or by a change in the temperature or mixing parameters. In particular, the distance versus segregation trade-off can be made fully explicit (figure 3) - a crucial prerequisite for informed political decisions.

Instead of using pedestrian distances, the algorithm can be also be fed by alternative transportation costs, such as the total displacement time on uni- or multimodal transportation unoriented networks. Also, beside gender mixing, the segregation issues related to child age or ethnic origin can be tackled by the same approach, to be addressed in a future work. Among further possible developments, the question of dealing with flexible school penalized capacities (instead of with fixed capacities) seems worth investigating, as is the issue of imposing a minimal social mixing at each school, rather than on average only.

# References

* Antunes A., Peeters D. (2000). A dynamic optimization model for school network planning. *Socio-Economic Planning Sciences* 34(2), 101–120.

* Caro F., Shirabe T., Guignard M., Weintraub A. (2004). School redistricting: Embedding gis tools with integer programming. *Journal of the Operational Research Society* 55(8), 836–849.

* Castillo-Lopez I., Lopez-Ospina H. A. (2015). School location and capacity modification considering the existence of externalities in students school choice. *Computers & Industrial Engineering* 80, 284–294.

* Clarke S., Surkis J. (1968). An operations research approach to racial desegregation of school systems. *Socio-Economic Planning Sciences* 1(3), 259–272.

* Delmelle E. M., Thill J.-C., Peeters D., Thomas I. (2014). A multi-period capacitated school location problem with modular equipment and closest assignment considerations. *Journal of Geographical Systems* 16(3), 263–286.

* Emmanouilidis T. (2016, 02). *Analyse de réseau piéton et gestion scolaire: sinuosité, centricité et transport optimal régularisé*. Ph. D. thesis, Université de Lausanne, faculté des Géosciences et de l'Environnement.

* Guex G., Emmanouilidis T., Bavaud F. (2016). Transportation clustering: a regularized version of the optimal transportation problem. *Submitted for publication*.

* Heckman L. B., Taylor H. M. (1969). School rezoning to achieve racial balance: a linear programming approach. *Socio-Economic Planning Sciences* 3(2), 127–133.

* Koenigsberg E. (1968). Mathematical analysis applied to school attendance areas. *Socio-Economic Planning Sciences* 1(4), 465–475.

* Villani C. (2009). *Optimal Transport, Old and New*. Springer.
<figure><img src="image_207.png" /></figure>

Poster Session
<figure><img src="image_208.png" /></figure>
# A method for testing the similarity of spatial samples

M. Virtudes alba-Fernández¹, Francisco J. Ariza-López², José Rodríguez-Avi¹

¹Dpt. Ststistics and O.R. University of Jaén, Sapin

²Dpt. Cartographic engineering, geodesics and Photogrammetry. University of Jaén, Spain

*Corresponding author: mvalba@ujaen.es

## Abstract

The purpose of this work is to propose a method to analyse the similarity between two spatial point patterns. Such similarity should be understood in the sense that both point patterns come from the same spatial point distribution. For this aim, first, we make use of a space-filling curve as a tool to linealize the space, with independence of its dimension, second, we model the count of points in the resultant grid by means of the multinomial law, and after that we test the homogeneity of both multinomial distributions by using the negative of the Matusita's affinity. Finally, we evaluate the performance of the procedure by means of a simulation study.

## I INTRODUCTION

The understanding of spatial point patterns is one of the major challenges of geographical analysis and has interest in many sciences (e.g. biogeography, crop sciences, ecology, geology, etc.). Spatial patterns and spatial statistical sampling is also a major issue in spatial data quality assessment because sampling is a very common procedure in order to derive estimates or perform tests (Ariza-Lpez, 2002). For these reasons the evaluation of the spatial similarity of two observed samples can be of interest in many cases.

It is usual that an estimate of an attribute or property (e.g. concentration of a mineral, positional accuracy, etc.), is derived from a sample of points under some spatial distributional hypothesis for such sample (e.g. following a theoretic or an observed spatial pattern). We think that, previous to any error or uncertainty consideration about the attribute or property estimation, we must confirm the underlying hypothesis about the position of samples, in our case: the spatial distribution of the taken sample in relation to other observed spatial pattern.

Two common methods in use to investigate discrete point events are the distance- and area-based tests. In the case of distance-based tests, we use whatever distance (e.g. Euclidean) between two events in order to determine a random, clustered, or uniform spatial pattern to the points (Bailey and Gatrell, 1995). On the other hand, in the case of area-based tests, we count the number of point events within a predefined spatial area (e.g. a quadrant, a census unit, etc.) (Andresen, 2009).

Our proposal is very different to other techniques such as the Ripley's K function (Ripley, 1976). The K function is a distance-based test that characterizes point processes at many distance scales and allows the detection of different behaviors (e.g. random, clustering, inhibition) (Freeman and Ford, 2002); and it can be used to test different specific patterns (e.g. homogeneous Poisson process -complete spatial randomness-, Matern hard-core process, Strauss process, etc.) (Dixon, 2012). For the correct application of these functions it is needed to accept several assumptions that not necessarily are fully met in reality, and if conditions are not met, the output may be incorrect (Bolibok, 2008).
<figure><img src="image_209.png" /><figcaption>Figure 1: Examples of some geological structures that can influence spatial patterns: a) batholiths, b) falt, c) Synclinal/anticlinal, d) shift, e) cataclasite/mylonite</figcaption></figure>

In this setting, we develop an area-based test centered on the counting of positional events (not attributes), which does not assume a theoretical model for the spatial distribution of the events (spatial distribution free). So, it can be applied for the comparison of two spatial samples (e.g. two control test samples, two field works, presences at two different times, etc.), with independence of sample sizes. For instance, the Figure 1 shows different geological structures that can determine the presence of point events in the space or the distribution of control samples in the space, and this method allows the comparison in such cases, and any other.

This method follows two steps; the first one is to make use of the space-filling curves (Sagan, 1994) as a tool to order the space, with independence of its dimension. The space-filling curves provide a partition of the space for a fixed level of neighborhood. The second one is to model the count of points in the resultant grid by means of the multinomial law. This way, the problem to test the similarity between sampled spatial distributions is equivalent to the problem of testing the homogeneity of two multinomial distributions. For doing this, we considered as a test statistic the negative of Matusita's affinity (Matusita, 1967).

On the other hand, the spatial distribution of samples may have an effect on the representative-ness of the sample. This is especially true in complex situations, although, at the end of the day, for many spatial data properties, the spatial pattern can be considered uniform. To evaluate the usefulness of the proposal, we carried out a simulation study in order to analyze the effect of the space-filling curves and the test statistics used, as well as, the effect of the sample sizes and the level of iterations in the curves. In particular, we took the Hilbert's curve and the Peano's curve, and we analyze the uniform pattern and for other spatial patterns, specially, those patterns associated with some well-known geological structures. Our findings suggest the methodology is able to detect the similarity of two spatial samples of points for all the tried cases.

This paper is organized as follows, after this introduction follows the description of the approach, where we present the statistical basis and the procedure for its application; next a simulation experiment is developed using some spatial patterns and levels of the space-filling curves. Finally conclusions are presented.

## II DESCRIPTION OF THE APPROACH

As said in the Introduction, the aim of this work is to propose a formal procedure to test whether two spatial distributions of points can be considered as similar or not. By Similar one should understand that coming from the same spatial pattern. To show how this approach works, let us consider the bivariate case, although for a general dimension, the procedure follows the same steps. So that, let $X_{11}, X_{12}, \dots, X_{1n_1}$ and $X_{21}, X_{22}, \dots, X_{2n_2}$ be two independent samples of points from two spatial patterns with size $n_i$, $i = 1, 2$; and without loss of generality, let us suppose that both samples take values in the unit-square $S = [0, 1] \times [0, 1]$. The application
of a particular space-filling curve induces a partition on S with $M = 2^{\nu} \times 2^{\nu}$ squares, where $\nu$ represents the number of iterations in the space-filling curve construction. This way, for a given space-filling curve and a fixed level $\nu$, the sampled points can be grouped into M classes, $C_1, C_2, \dots, C_M$, or equivalently, taking values in $\Upsilon_M = (1, 2, \dots, M)$. The order in M is induced by the space-filling curve considered. Note that the degree of neighborhood is different for each space-filling curve, specially, when the number of iterations increases.

Let $\pi_i^t = (\pi_{i1}, \pi_{2m}, \dots, \pi_{iM})^t$ be the cell probabilities associated with each multinomial distribution, that is, $\pi_{im} = P[X_i = m]$, for $i = 1, 2$, and $1 \le m \le M$, and where the superscript $t$ denotes transpose. The application of a particular space-filling curve to both samples of points, provide us the number of points falling into each class in both samples. That is to say, we obtain the observed frequencies on each cell and hence, the maximum likelihood estimator of $\hat{\pi}_i$, say, $\hat{\pi}_i = \frac{n_{im}}{n_i}$, $i = 1, 2$, $m = 1, \dots, M$. As a result, the problem of testing whether two spacial distributions are equal is equivalent to test whether two multinomial populations are equal. So that, our objective is to test the following null hypothesis

$$H_0: \pi_1 = \pi_2. \qquad (1)$$

For this purpose, several measures able to discriminate between multinomial populations can be used. Here, we consider a f-dissimilarity measure between multinomial populations because the smaller these measure is, the harder it is to discriminate between them. Specifically, it is considered the negative of Matusita's affinity (Matusita, 1967) defined by

$$T = - \sum_{m=1}^{M} \sqrt{\hat{\pi}_{1m} \hat{\pi}_{2m}},$$

which is a member of a class of test statistics based on f-dissimilarity (see Zoografos, 1998; Alba et al., 2009) and references therein for further theoretical results).

To decide when to reject $H_0$, we need to know the null distribution of $T$, or at least an approximation to it. Following Zoografos (1998), under $H_0$, the asymptotic null distribution of the test statistic is given by

$$8 \frac{n_1 n_2}{n_1 + n_2} (1 + T) \xrightarrow{\mathcal{L}} \chi_{M-1}^2. \qquad (2)$$

So, we reject $H_0$ if

$$8 \frac{n_1 n_2}{n_1 + n_2} (1 + T_{obs}) \ge \chi_{\alpha, M-1}^2,$$

where $T_{obs}$ represents the observed value of the test statistic $T$, and $\chi_{\alpha, M-1}^2$ denotes the $1 - \alpha$ percentil of the chi-square distribution with $M - 1$ degrees of freedom, $0 \le \alpha \le 1$; or equivalently, we reject $H_0$ if the corresponding $p$-value is less than or equal to $\alpha$, that is, if

$$p = P_{H_0} \left[ \chi_{M-1}^2 \ge 8 \frac{n_1 n_2}{n_1 + n_2} (1 + T_{obs}) \right] \le \alpha,$$

where $P_{H_0}$ stands for the probability law under the null hypothesis.

As observed among others in Kim, 2009, Alba et al., 2009 or Jiménez-Gamero et al., 2014, the $\chi^2$ approximation is rather poor for small and moderate sample sizes, and the approximation
of the null distribution by means of a parametric bootstrap estimator behaves better than the asymptotic one. From these results, we approximate the p-value by bootstrapping (see Alba-Fernández et al. 2009 for the theoretical properties of the bootstrap estimator).

The application of the procedure to study of similarity of two spatial samples of points follows the steps:

(1) Given the two spatial distributions of points {$X_{1j}$}$_{1\le j\le n_1}$ and {$X_{2j}$}$_{1\le j\le n_2}$, choose a space-filling curve and the number of the iterations in its construction, $\nu$. These decisions will determine the value of $M$.

(2) Apply the space-filling curve and obtain $\hat{\pi}_i$, $i = 1, 2$.

(3) Calculate $T_{obs}$, the observed values of $T$.

(4) Approximate the corresponding p-value by bootstrapping.

(5) Conclude if the sample spatial distributions can be considered equal or not.

In addition, the bootstrap algorithm to approximate the p-value for testing (1) can be assessed as follows:

(i) Calculate $T_{obs}$, the observed values of $T$.

(ii) For $b = 1, \dots, B$, generate $2B$ independent bootstrap samples, {$X_{1,j}^{*b}$}$_{1\le j\le n_1}$, {$X_{2,j}^{*b}$}$_{1\le j\le n_2}$ from the pooled multinomial distribution $M(n_1 + n_2; \hat{\pi}_{01}, \dots, \hat{\pi}_{0M})$ where

$$ \hat{\pi}_{0m} = \frac{n_{1m} + n_{2m}}{n_1 + n_2}, \quad m = 1, \dots, M. $$

(iii) Calculate the values of $T$ for each couple of bootstrap samples, say $T^{*b}$, $b = 1, \dots, B$.

(iv) Approximate the p-value by means of $\hat{p} = \text{card}\{b : T^{*b} \ge T_{obs}\}/B$, respectively.

### III SIMULATION EXPERIMENT

To evaluate the performance of this method, we have carried out a simulation study. The goal of this experiment is twofold, the first objective is to analyze the behaviour of the proposed methodology for small sample sizes with respect to the type I error for some space-filling curves and levels of sweep. In other words, if the procedure is able to conclude if two sample of points come from the same spatial pattern when they really do.

Together with this, the second task is to evaluate the power of the procedure, that is, if this method is able to detect two samples of point which have been generated by different spatial patterns. Next we briefly describe the simulation experiment and display the results obtained. All computations in this paper have been performed using scripts written in the R language (http://www.cran.r-project.org).

The methodology we propose is independent of the spatial pattern of sample of points, however, we are going to evaluate its behaviour by considering some spatial patterns related to the geological structures mentioned before (see figure 1). In particular, we will identify the uniform pattern with Pattern 1, the bivariate normal $N(\mu, \Sigma)$ with $\mu = (0.25, 0.25)$, $\Sigma = \begin{pmatrix} 0.25 & 0 \\ 0 & 0.25 \end{pmatrix}$ simulates a synclinal/anticlinal(Pattern 2) and the bivariate normal $N(\mu, \Sigma)$ with $\mu = (0.5, 0.5)$, $\Sigma = \begin{pmatrix} 0.25 & 0.8 \\ 0.8 & 0.25 \end{pmatrix}$ reproduces a Batholiths (Pattern 3) (see Figure 2) for a scatter plot of a sample of size 150).

So, we have generated two uniform samples on the unit-square of sizes $n_1 = n_2 = 25$, we have applied the method following the steps (1)-(5) described above with $B = 1000$ bootstrap
<figure><img src="image_210.png" /><figcaption>Figure 2: Example of Pattern 1 (left), Pattern 2 (center) and Pattern 3 (right).</figcaption></figure>

replications. We repeated this 1000 times and we calculated the fraction of *p*-values less than or equal to 0.05 and 0.10 (denoted as *f*05 and *f*10 in tables), which are the estimated type I error probabilities for $\alpha$ = 0.05, 0.10, respectively. We have repeated the whole experiment for the sample sizes $n_1$ = $n_2$ = 50, 150. Among the available space-filling curves, we have considered the Peano's and the Hilbert's curve (identified as *n* and *h* in tables). Note the value of *M* is related to the sample size and $\nu$, and we try to find an agreement between them. Here, we try the values $\nu$ = 1, 2. We have repeated this experiment for Patterns 2 and 3. The estimated type I error probabilities are shown in Table 1.

Looking at this table, we can conclude that the estimated type I error probabilities are quite close to the nominal ones in all the tried cases and the simulations results do not show differences between the space-filling curves. On the other hand, we have also studied if the methodology is able to distinguish between spatial patterns, and for this task we generated a sample of points following the Pattern 1 of size $n_1$ = 25 and other sample of points from Pattern 2 of size $n_2$ = 25. We applied the methodology in the same conditions as before. We repeated the whole experiment 1000 times and we calculated the fraction of *p*'s values less than or equal to 0.5, 0.10, which are the estimated power for $\alpha$ = 0.5, 0.10 (now, *f*05 and *f*10 in tables). We repeated the experiment for samples of points following the Patterns 3. The estimated powers are shown in Table 2. From these results, we can say the procedure is able to distinguish clearly between two different Patterns.

<table><caption>Table 1: Estimated type I error probabilities.</caption><thead><tr><th rowspan="2">n<sub>1</sub> = n<sub>2</sub></th><th rowspan="2">ν</th><th rowspan="2">curve</th><th colspan="2">Pattern 1</th><th colspan="2">Pattern 2</th><th colspan="2">Pattern 3</th></tr><tr><th>f05</th><th>f10</th><th>f05</th><th>f10</th><th>f05</th><th>f10</th></tr></thead><tbody><tr><td rowspan="2">25</td><td rowspan="2">1</td><td>n</td><td>0.042</td><td>0.099</td><td>0.052</td><td>0.100</td><td>0.048</td><td>0.096</td></tr><tr><td>h</td><td>0.047</td><td>0.102</td><td>0.052</td><td>0.096</td><td>0.050</td><td>0.096</td></tr><tr><td rowspan="2">50</td><td rowspan="2">1</td><td>n</td><td>0.055</td><td>0.103</td><td>0.053</td><td>0.105</td><td>0.054</td><td>0.106</td></tr><tr><td>h</td><td>0.053</td><td>0.101</td><td>0.053</td><td>0.101</td><td>0.053</td><td>0.107</td></tr><tr><td rowspan="2">150</td><td rowspan="2">1</td><td>n</td><td>0.054</td><td>0.104</td><td>0.058</td><td>0.107</td><td>0.046</td><td>0.093</td></tr><tr><td>h</td><td>0.057</td><td>0.099</td><td>0.059</td><td>0.104</td><td>0.045</td><td>0.098</td></tr><tr><td rowspan="2">150</td><td rowspan="2">2</td><td>n</td><td>0.056</td><td>0.104</td><td>0.048</td><td>0.094</td><td>0.060</td><td>0.108</td></tr><tr><td>h</td><td>0.057</td><td>0.103</td><td>0.048</td><td>0.093</td><td>0.057</td><td>0.107</td></tr></tbody></table>
<table><caption>Table 2: Estimated power.</caption><thead><tr><th>n<sub>1</sub> = n<sub>2</sub></th><th>ν</th><th>curve</th><th>Pattern 1 vs. Pattern 2</th><th>Pattern 1 vs. Pattern 3</th><th>Pattern 2 vs. Pattern 3</th></tr></thead><tbody><tr><td rowspan="2">25</td><td rowspan="2">1</td><td>n</td><td>f05</td><td>f10</td><td>f05</td><td>f10</td></tr><tr><td>h</td><td>0.429</td><td>0.548</td><td>0.409</td><td>0.534</td><td>0.570</td><td>0.665</td></tr><tr><td rowspan="2">50</td><td rowspan="2">1</td><td>n</td><td>f05</td><td>f10</td><td>f05</td><td>f10</td><td>f05</td><td>f10</td></tr><tr><td>h</td><td>0.427</td><td>0.557</td><td>0.366</td><td>0.501</td><td>0.537</td><td>0.635</td></tr><tr><td rowspan="2">150</td><td rowspan="2">1</td><td>n</td><td>f05</td><td>f10</td><td>f05</td><td>f10</td><td>f05</td><td>f10</td></tr><tr><td>h</td><td>0.747</td><td>0.837</td><td>0.719</td><td>0.802</td><td>0.886</td><td>0.938</td></tr><tr><td rowspan="2">150</td><td rowspan="2">2</td><td>n</td><td>f05</td><td>f10</td><td>f05</td><td>f10</td><td>f05</td><td>f10</td></tr><tr><td>h</td><td>0.994</td><td>0.995</td><td>0.998</td><td>0.999</td><td>0.999</td><td>1.000</td></tr><tr><td rowspan="2">150</td><td rowspan="2">2</td><td>n</td><td>f05</td><td>f10</td><td>f05</td><td>f10</td><td>f05</td><td>f10</td></tr><tr><td>h</td><td>0.987</td><td>0.994</td><td>1.000</td><td>1.000</td><td>1.000</td><td>1.000</td></tr></tbody></table>

## IV CONCLUSION

To sum up, the procedure introduced in Section 2 takes advantage of the use of space filling curves as a way to linearize spatial distributions, and following the order induced by the space filling curve for a fixed level of its construction, the count of the number of points falling into each grid may be modeled by a multinomial distribution. For testing the homogeneity of two multinomial laws, the negative of Matusita's affinity is considered. Finally, the results of the simulation experiment reveals that the proposed method provide us a statistical tool in order to decide about the similarity of spatial samples.

## ACKNOWLEDGMENTS

Research in this paper has been partially funded by grant CTM2015-68276-R of the Spanish Ministry on Science and Innovation.

## REFERENCES

Alba-Fernández, V., Jiménez-Gamero, M.D. (2009). Bootstrapping divergence statistics for testing homogeneity in multinomial populations. Mathematics and Computers in Simulations, 79, 3375–3384.

Andresen M.A. (2009). Testing for similarity in area-based spatial patterns: A nonparametric Monte Carlo approach. Applied Geography 29, 333345.

Ariza-Lpez F.J. (2002). Calidad en la produccin cartogrifica. Madrid, ES, Ra-Ma.

Bailey T.C., Gatrell, A.C. (1995). Interactive spatial data analysis. Harlow, UK: Prentice Hall.

Bolibok L. (2008). Limitations of Ripleys k(t) function use in the analysis of spatial patterns of tree stands with heterogeneous structure. Acta Sci. Pol. Silv. Colendar. Rat. Ind. Lignar 7(1), 5-18.

Dixon P.M. (2012). Ripley's K function. Encyclopedia of Environmetrics, 2nd ed. John Wiley and Sons, Inc.

Freeman E.A., Ford E.D. (2002). Effects of data quality on analysis of ecological pattern using the k(t) statistical function. Ecology 83(1), 3546.
Jiménez-Gamero, M.D., Alba-Fernández, V., Barranco-Chamorro, I., Muñoz-García, J. (2014).
Two classes of divergence statistics for testing uniform association. Statistics, 48 (2), 367–387.

Matusita, K. (1967). On the notion of affinity of several distributions and some of its applications. Annals of the Institute of Statistical Mathematics, 19, 181–192.

Sagan H. (1994). Space-Filling Curves. Springer-Verlag.

Ripley, B.D. (1976). The second-order analysis of stationary point processes. Journal of Applied Probability 13, 255-266.

Zografos, K. (1998). f-dissimilarity of several distributions in testing statistical hypotheses. Annals of the Institute of Statistical Mathematics, 50, 295-310.
# Combining punctual and ordinal contour data for accurate floodplain topography mapping

Carole Delenne*<fnref n="1" />,<fnref n="3" />, Jean-Stéphane Bailly<fnref n="2" />, Mathieu Dartelle<fnref n="3" />, Nelly Marcy<fnref n="1" />, Antoine Rousseau<fnref n="3" />

*Corresponding author: carole.delenne@umontpellier.fr

In the framework of hydrodynamic modelling, topography is classically obtained by an interpolation of punctual field elevation surveys. A methodology, based on block conditional simulations, is presented to enhance mapping accuracy, using contour lines extracted from flooded areas in remote sensing data.

## I INTRODUCTION

Hydrodynamic models in two dimensions require a precise knowledge of the domain topography. However, as far as the modelling of great rivers or lakes is concerned, few accurate topographic records are generally available to accurately model the floodplain topography. The usual way to acquire accurate topographic information for hydraulics over floodplain remains the ground surveys, that provide punctual values at a very high cost. Based on these points, usual interpolation schemes often yield a too coarse and inaccurate topographic map for a realistic hydrodynamic modelling. Remote sensing data such as Lidar or stereo-photogrammetry data on non vegetated areas, constitute a good alternative to obtain large scale information, but yield other issues of acquisition and data processing. In addition, topographic Lidar operating with near-infrared lasers are still not suitable for an exhaustive survey of floodplains where submerged areas remain.

Progress in remote sensing data repeatability and spatial resolution now allows the automatic monitoring of water surfaces delineation from areal or satellite images either in optical or radar domains (see e.g. an application to the Inner Niger Delta in Ogilvie et al. (2015)). As for the recently launched Sentinel sensors, these data are becoming widely available with increasing spatial and temporal resolutions, allowing in turn the spatio-temporal monitoring of flooded areas. Flood dynamics from remote sensing data are known to be informative on floodplain topography for long (see e.g. Schumann et al. (2007) or Hostache et al. (2010)). Indeed, the extracted flooded areas may be considered as iso-elevation contour lines, as the hypothesis can be made that all the points located on the water/soil limit have the same elevation. If contour line elevation remains unknown, rank between the detected contour lines as well as between every contour line and surrounding topographic data points are known. Assuming this information, the next challenge is thus to develop a spatial data fusion method of ordered contour lines data and usual elevation data points to better model the floodplain topography.

Mixing punctual ranked and continuous data in spatial estimation was already proposed in literature, mainly to deal with highly skewed field distribution (Yamamoto, 2000) or for very large spatial datasets (Cressie and Johannesson, 2008). These methods used the uniform score transform and back-transform based on the standardized rank estimation (Saito and Goovaerts, 2000) to prevent the right estimate distribution and limit the smoothing effect of kriging. Uncertainty

<footnote n="1">¹Univ. Montpellier, HSM, France</footnote>

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

<footnote n="3">³INRIA, Lemon, France</footnote>
assessment of resulting punctual kriging estimates was proposed thanks to the transformation
process (Yamamoto, 2007). Other methods also exist to quantify the uncertainty of contour
lines in random fields (Lindgren and Rychlik, 1995; Wameling and Saborowski, 2001). But
in these studies, the objective was to quantify the uncertainty in contour lines location and not
to quantify the uncertainty of contour lines level for a known location. Consequently, these
different studies only partly touched the problem we face.

The objective of this paper is to propose a method to estimate the elevation of the ranked contour
lines from data points. This method starts from block conditional simulations estimates filtered
on rank statistics. After contour line elevation estimate, a usual kriging is performed to recon-
struct a topographic map over the floodplain. Starting from a synoptic and theoretical example,
the principles, advantages, properties and limitations of the proposed method are exposed. The
spatial accuracies obtained in contour lines estimates and topographic map are thus compared
to the usual block kriging estimates.

II MATERIAL AND METHODS

2.1 Case study generation

A reference exhaustive elevation field Z(s) on a floodplain D was generated from a Gaussian
spatial covariance model with range 20, sill 1 and no nugget effect. The resulting field is
centered on μ = 10 and respects distribution N(μ, σ) (see Figure 1-left): the blue colours
represent lower floodplain elevations meanwhile the yellow and green stem for higher elevation
areas. D is a 100 m × 100 m area and the dataset was generated on gridded D with resolution
r = 1 m, r denoting both the remote sensing image spatial resolution and the desired end-user
topographic map resolution.

<figure><img src="image_211.png" /><figcaption>Figure 1: **Generated case study.** Left: Simulated elevation field over D; Middle: flooded area at time 1 (blue) and resulting rasterized contour (black crosses); Right: flooded area at time 2 (blue) and resulting rasterized contour (black crosses).</figcaption></figure>

To reproduce the remote sensing data support (Ogilvie et al., 2015), flooded areas were gen-
erated and delineated at resolution 0.5 m for two different times t₁ and t₂ during rising water:
when areas lower than 9.5 m (Fig. 1-middle) and 8.5 m (Fig. 1-right) respectively are flooded.
Contour polylines were generated from these rasterized flooded areas with 0.5 m spatial reso-
lution. The sets of two generated contour polylines C₁ and C₂ were thus reduced to a set of
polylines vertices sᵢⱼ (j ∈ (1, 2) denoting the polyline) regularly located along the lines (black
cross on figures 1). C₁ and C₂ contained n₁ = 201 and n₂ = 520 points respectively. Four
contour lines correspond to $C_1$ and five to $C_2$ (Fig. 1). When contour lines were connected to the $D$ boundaries, only vertices within $D$ were considered.

In addition, a sample of $n = 250$ topographic points $s_1...s_n$ was randomly selected on $D$, excluding points selection at location closer than $b = 2$ m from polylines (Fig. 2-left). This ensures that no data points are located on a polyline, which would make too obvious the estimation of its elevation. We further assume that contour polylines elevation is unknown but that ranks between the polylines (i.e. the polyline vertices $s_i^j$) and between polylines and surrounding topographic data points are known and may be computed (Fig. 2-right).

<figure><img src="image_212.png" /><figcaption>Figure 2: **Generated dataset.** Left: Simulated elevation field over $D$, contour polylines $C_1$ and $C_2$ (black lines) and data points (red points). Right: ranked data (points and polyline vertices) and according legend.</figcaption></figure>

## 2.2 Method

The method we developed may be seen as a block conditional simulation process filtered by rank statistics, where a polyline $C_j$ is seen as a block, i.e. a region corresponding to an unconnected set of vertices. Such simulation is known to be superior to kriging whenever interest lies in global statements for a region rather than inference on individual points. In the process described in the algorithm 1 hereafter, assuming a stationary random function, a spatial model (variogram) $\gamma(h)$ is first estimated and modelled from the 250 data points $s_1...s_n$.

Once the variogram $\gamma(h)$ is modelled, an ordinary kriging is performed on each polyline vertex. From the kriging, a first estimate of polyline elevation results from vertices averaging. For polyline $C_j$, we denote further $\hat{z}_{ko}^j$ this estimate (corresponding to the usual block kriging estimate).

Next, a set of $N = 100$ conditional Multi-Gaussian simulations using the fitted variogram model and pathing through the 250 $s_1...s_n$ is performed on each on the $s_i^j$ contour vertices (Eq. 1):

$$z_{cs}(s_i^j) = \hat{z}_{ko}(s_i^j) + (z_{us}(s_i^j) - \hat{z}_{ko}^*(s_i^j)) \quad (1)$$

At the end of the conditional simulation, an estimate at polyline scale is computed by averaging simulated values on vertices. This remains consistent with the usual block kriging estimate
**Algorithm 1 Contour polyline estimation process and field reconstruction**

1: **for** j = 1 to j= 2 **do**

2: Estimate and model the variogram $\gamma(h)$ from n data points

3: Estimate ordinary kriging $\hat{z}_{ko}(s_i^j)$ on polyline vertices

4: **for** n = 1 to N = 100 **do**

5: Draw unconditional multigaussian simulation $z_{us}(s_k)$ and $z_{us}(s_i^j)$ on data location (k) and polyline vertices (i).

6: Estimate ordinary kriging at vertices from simulated values at data locations $z_{ko}^*(s_i^j)$

7: Compute conditional simulation on vertices (Eq. 1)

8: Compute N polyline estimate by averaging vertices conditional simulations

9: Filter keeping only polyline estimate conform to the initial ranks

10: Compute polyline estimate $\hat{z}_j$ averaging the kept polyline estimates

11: Compute standard deviation estimate from the set of kept polyline estimates

12: Affect polyline estimate $\hat{z}_j$ to each polyline vertex $s_i^j$

13: Perform ordinary kriging on all gridded D points from data $z(s_1) \cdots z(s_n)$ and polyline vertices estimates $\hat{z}_j(s_i^j)$

since vertices are regularly located along the polyline.

$$ \hat{z}_k^j = \frac{1}{n_j} \sum_{i=1}^{n_j} z_{cs}(s_i^j) \qquad (2) $$

We thus obtain N elevation estimates $\hat{z}_k^j, k \in (1 \dots N)$ for a given contour polyline $C_j$. Only the estimates respecting the rank are kept in the following, yielding a set of M estimates ($M \le N$) for a given polyline. From this set, the final estimate of polyline $C_j$ becomes:

$$ \hat{z}_j = \frac{1}{M} \sum_{k=1}^{M} \hat{z}_k^j \qquad (3) $$

with uncertainty characterized by the standard deviation $\sigma_{C_j}$ of the $\hat{z}_k^j$ estimates.

In a final step, polyline estimate $\hat{z}_j$ is affected to each polyline vertex $s_i^j$ in order to perform, in addition to the 250 data points a gridded map using ordinary kriging.

### III RESULTS

Results obtained for this theoretical test case are given in Figures 3 and 4. Figure 3-left shows the estimated empirical variogram (dots) and fitted variogram model (lines). Figure 3-right shows for $C_1$ and $C_2$ the resulting distribution of the estimates using conditional simulation, the final estimation equal to the average of conditional simulation (black vertical line) compared to the true value (blue vertical line) or the ordinary kriging (interpolation) estimate (red vertical line). In this first application, rank was always respected in the conditional simulation and for the ordinary kriging but, as shown in the following, this will not always be the case. The estimated elevation obtained for $C_1$ is 8.53 m (8.45 using ordinary kriging) for a true value equals to 8.5 m, with a standard deviation of the estimate being only 5.37e-03 m. For $C_2$, the difference between the true elevation (9.5 m) and the estimated value is lower than one centimetre (3 cm using kriging) with standard deviation of the estimate equals to 3.99e-03 m.
<figure><img src="image_213.png" /><figcaption>Figure 3: Left: Re-estimated empirical variogram from the $n = 250$ data points; Right: estimated elevation for $C_1$ and $C_2$.</figcaption></figure>

Figure 4 shows the reconstructed fields from the method we proposed compared to i) the true field and ii) the field reconstructed only from the 250 data points without considering the contour polylines data. A visual comparison suggests better results for the proposed approach. However, artefacts can be seen around the pool located at the bottom-right corner of the domain. In this area, the elevation increases with the distance to the pool minimum but decreases again at three places. This behaviour can be explained by the fact that very few topographic points are available around this pool.

On the contrary, the field reconstructed with the ordinary kriging method shows greater difference to the true field especially on the top part of the domain, despite a high concentration of ground-truth surveys.

To better assess the precision given by the two approaches, the root mean square difference $e$ between results and true field is computed. The error value obtained without considering the contour polylines is equal to 0.232 m, and to 0.199 m for the proposed method.

<figure><img src="image_214.png" /><figcaption>Figure 4: **Reconstructed elevation fields:** Left: true field $D$; Middle: ordinary kriging reconstructed field; Right: method (polyline conditional simulation) reconstructed field</figcaption></figure>

This first test suggests that the sampling procedure may influence the results, at least through the location of acquisition points. To check the robustness of the proposed approach, the process described above was repeated 100 times, changing the 250 points data sampling.
Figure 5-left shows the contour polyline $C_1$ and $C_2$ error distribution obtained. Figure 5-right shows the mean distance to the true field distribution. Clearly, the proposed method using conditional simulation is more robust with few error variance. The standard deviation of the results obtained with this approach is 0.056 m (resp. 0.037 m) for the contour line 8.5 (resp. 9.5) compared to 0.183 m (resp. 0.1 m) using the ordinary kriging method.

Rank filtering occurred 10% of the times for the $C_1$ estimation and never for the $C_2$ estimation. This suggests that rank filtering is useful when contours are in extreme values.

<figure><img src="image_215.png" /><figcaption>Figure 5: Errors distribution on 100 random sample sets of 250 data points. Left: $C_j$ estimation error; Right: mean distance to the true field after reconstruction. The vertical lines represent the mean of the distribution.</figcaption></figure>

Field survey being very costly, the aim is to reduce to the minimum the required elevation data points. The test is thus performed again with different numbers of points. An example of result obtained with a random sample of $N = 50$ points is given in Figure 6. It can be seen that the use of contour line enables to retrieve the pools shapes with better accuracy than an ordinary kriging based only on data points.

The comparison between true field and computed elevation fields with the two approaches is done using the root mean square error. These results are summarized in Table 1 for random samples obtained with $N = 10$ to $N = 250$. As expected, the error decreases if the number of available data points is high and the proposed approach always gives better results than the ordinary kriging.
<figure><img src="image_216.png" /><figcaption>Figure 6: **Reconstructed elevation fields using only 50 data points:** Left: true field *D*; Middle: ordinary kriging reconstructed field; Right: method (polyline conditional simulation) reconstructed field.</figcaption></figure>

<table><caption>Table 1: Root mean square errors between true field and reconstructed elevation fields for different *N*.</caption><thead><tr><th></th><th>N = 10</th><th>N = 50</th><th>N = 100</th><th>N = 150</th><th>N = 200</th><th>N = 250</th></tr></thead><tbody><tr><th>e<sub>ko</sub></th><td>1.064</td><td>0.635</td><td>0.361</td><td>0.403</td><td>0.435</td><td>0.232</td></tr><tr><th>e<sub>cs</sub></th><td>0.639</td><td>0.469</td><td>0.329</td><td>0.227</td><td>0.189</td><td>0.199</td></tr></tbody></table>

## IV DISCUSSION AND CONCLUSION

The objective of this work was to put forward a methodology to enhance the accuracy of the topographic description required by numerical hydrodynamic modelling, by using information of level lines available from remote-sensing data.

The first results obtained on a totally theoretical example, show that the topographic estimation benefits from such additional data; and the repetition of the process indicates that this result is robust to sampling. If the gain in accuracy may seem limited at this stage (about 3 cm in mean for 250 points), the influence of the sampling rate should be assessed more precisely on the two approaches, since the benefit taken from additional sources of information is expected to increase as less ground-truth data are available.

Moreover, a classical random sampling is obviously not the best approach to conduct a field survey. In the first presented test case, random sampling excluding the direct neighbourhood of the contour lines yielded a lack of data near the pool located at bottom-right part of the domain. This produced irrelevant estimation of the topography in some areas. Further tests are thus planed to infer guidelines in the sampling definition in order to decrease the number of data points required while maintaining a high accuracy in the results, thus enhancing the cost/benefice ratio. In the framework of hydrodynamic modelling, this work will integrate results of uncertainty and sensitivity analyses such as given in Guinot and Cappelaere (2009) or Delenne et al. (2012).

The methodology developed on theoretical test cases, will be assessed in a real-case study of the hydrodynamic modelling of the Vaccares lagoon for which amount of data are available for validation.

## References

* Cressie N., Johannesson G. (2008). Fixed rank kriging for very large spatial data sets. *Journal of the Royal Statistical Society: Series B (Statistical Methodology)* 70(1), 209–226.
Delenne C., Cappelaere B., Guinot V. (2012). Uncertainty analysis of river flooding and dam failure risks using local sensitivity computations. *Reliability Engineering and System Safety* 107, 171–183.

Guinot V., Cappelaere B. (2009). Sensitivity equations for the one-dimensional shallow water equations: Practical application to model calibration. *Journal of Hydrologic Engineering* 14, 858–861.

Hostache R., Lai X., Monnier J., Puech C. (2010). Assimilation of spatially distributed water levels into a shallow-water flood model. part ii: Use of a remote sensing image of mosel river. *Journal of Hydrology* 390(3-4), 257 – 268.

Lindgren G., Rychlik I. (1995). How reliable are contour curves? confidence sets for level contours. *Bernoulli*, 301–319.

Ogilvie A., Belaud G., Delenne C., Bader J.-C., Oleksiak A., Bailly J. S., Ferry L., Martin D. (2015). Decadal monitoring of the Niger Inner Delta flood dynamics using MODIS optical data. *Journal of Hydrology* 523, 358–383.

Saito H., Goovaerts P. (2000). Geostatistical interpolation of positively skewed and censored data in dioxin-contaminated site. *Environmental Science Technology* 34, 4228–4235.

Schumann G., Matgen P., Hoffmann L., Hostache R., Pappenberger F., Pfister L. (2007). Deriving distributed roughness values from satellite radar data for flood inundation modelling. *Journal of Hydrology* 344(1-2), 96 – 111.

Wameling A., Saborowski J. (2001). Construction of local confidence intervals for contour lines. In *Proceedings of the IUFR4.11 Conference in forest biometry*.

Yamamoto J. K. (2000). An alternative measure of the reliability of ordinary kriging estimates. *Mathematical Geology* 32, 489–509.

Yamamoto J. K. (2007). On unbiased back-transform of lognormal kriging estimates. *Computacional Geosciences* 11, 219–234.
# Monitoring spatial accuracy of oil palm cultivation mapping in southern Cameroon from Landsat series images

Prune Christobelle Komba Mayossa¹*, Sébastien Gadal¹

¹Aix-Marseille Université, CNRS ESPACE UMR 7300, France

*Corresponding author: prune.komba-mayossa@etu.univ-amu.fr

## Abstract

Studying and mapping palm grove evolution allow understanding the impact related to its cultivation. Our study aims to map industrial palm grove using Landsat series images and measures the accuracy of the produced maps. It was carried out in SOCAPALM industrial plantation, located in southern of Cameroon. For the mapping and assessment of accuracy, per-pixel classification and confusion matrix method were used, respectively. We obtained high correlated maps (Kappa =0.92 in 2001 vs 0.86 in 2015). However, some confusions were observed between vegetation and oil palm classes for the two periods, affecting the maps accuracy. These confusions are caused by the presence of mixed pixels resulting from the spatial and spectral characteristics of palm groves, the method used to map and validate the map, and uncertainty related to dada. To increase the accuracy, we suggest (1) to use another mapping method such as super-resolution mapping, (2) develop a classification system of cartographic products.

## Keywords

Elaeis Guineensis.Jaq, oil palm, remote sensing, spatial accuracy, monitoring.

## I INTRODUCTION

Oil palm (Elaeis guineensis.jaq) is a perennial oleaginous plant from Central Africa. In Cameroon oil palm has a high economic importance with an industrialisation dating from colonial period (Elong, 2003). The incomes generated by oil palm cultivation have developed agro-industries such as SOCAPALM (Cameroonian Society of Palm groves). This activity in high-yield and low cost of returns (Riva, 2013), is behind of socio-environmental damages whom deforestation, loss of biodiversity, pollution, etc. Studying and mapping palm grove evolution allow understanding the impacts related to this culture. Our work is placed in the context of the long-term management of Oil palm resources in Congo basin. This paper is focused on the first part of this project: studying, mapping and monitoring SOCAPALM industrial palm groves using remote sensing. Several studies have shown the capability of remote sensing to map plant resources (Domaç and al, 2004; Li and al., 2015; Koh and al, 2011; etc.). The accuracy of the produced maps is main question for the study of evolution, mapping and characterization of palm grove, and for the reliability of our method.
This first part aimed to map palm grove and assess the spatial accuracy of produced map.

## II STUDY AREA AND DATA

### II.1 study area

Located in the Congo Basin, the Republic of Cameroon is a central African country with a surface area of 475,440 Km² and 20,000,000 inhabitants in 2015. Two main climatic zones are defined. In south the climate is equatorial and sub-equatorial, temperatures range between 25°C and 35°C, annual abundant precipitations up to 1000mm. In north the tropical climate is Sudanese (high temperatures and scarce rainfall) and Sahelian (irregular precipitations). SOCAPALM industrial plantations are located in south west of Cameroon, in the Ocean County, not far from the port city of Kribi (see figure 1).

<figure><img src="image_217.png" /><figcaption>Figure 2: Localisation of study area, from CNRS ESPACE UMR 7300.</figcaption></figure>

### II.2 Data

Landsat 7 ETM+ and Landsat 8 OLI images of same season, but with different acquisition date were used (see Table 1). These images were acquired on the 26th April 2001 (ETM+) and 25th April 2015 (OLI), respectively.

<table><caption>Table 1. Sensor specification of Landsat ETM+ and OLI.</caption><thead><tr><td>Parameters</td><td>ETM+</td><td>OLI</td></tr></thead><tbody><tr><td>Spectral range</td><td>0.45 -12.5µm</td><td>0.45-235µm</td></tr><tr><td>Spatial resolution</td><td>MS:30m/PAN:15m</td><td>MS:30m/PAN:15m</td></tr><tr><td>Swath width</td><td>183km</td><td>185km</td></tr><tr><td>Spatial coverage</td><td>Non-continuous</td><td>Non-continuous</td></tr><tr><td>Total number of bands</td><td>8</td><td>11</td></tr><tr><td>Mode</td><td>MS/PAN</td><td>MS/PAN</td></tr><tr><td>Date of acquisition of image</td><td>April 2001</td><td>April 2015</td></tr></tbody></table>
### III METHODOLOGY

Methodology adopted for this study is depicted in Figure 2. It involves mapping palm grove using spectral classification, evaluating the precision of produced maps for both Landsat images 2001 and 2015.

<figure><img src="image_218.png" /><figcaption>Figure 2: Flowchart of methodology adopted in this study.</figcaption></figure>

#### III.1 Pre-processing

The pre-processing phase is composed of a chain of 6 calculs:

(i) Pansharpening method to improve the spatial resolution of the multispectral image (30m), merging it with the panchromatic image (15m), by the Brovet transform (Lacombe, 2008). The function resamples automatically the seven (Landsat 7) or eleven (Landsat 8) channels in the maximal resolution by using several methods; in our case, the cubic convolution method was used, to obtain a multispectral image of 15mx15m of resolution.

(ii) Radiometric corrections, images were calibrated in radiance by applying the equation:

$$L=a*CN+b, \text{ where } CN \text{ is the digital count, } a \text{ is the gain and } b \text{ is the bias. The coefficients } a \text{ and } b \text{ of sensor calibration are given in the metadata files. The luminance}$$

$(W.m^{-2}.sr^{-1})$ was calculated for each band.

(iii) As the presence of cloud (30% of cloud cover) and water area, may disturb the analysis, they were masked.

(iv) The vegetation index was calculated according to the equation: $NDVI=NIR-RED/NIR+RED$. It allows characterizing the vegetation cover in terms of level of Chlorophyll (Pouchin and al., 2002).

(v) To improve the detectability of objects, (Gadal, 2003) convolution filters were used: morphological (maximum filter) and directional filter with a kernel size of 3x3 pixels.

(vi) Eight regions of interest (ROIs) were defined according to the spatial organisation of landscape: growing, young, and mature oil palm, low vegetation, forest, bare ground/built-up areas and waterway.
## III.2 PROCESSING

### III.2.1 Palm grove mapping

First, we created an image with B, G, R, NIR, SWIR and NDVI channels, masked water area, extracted spectral signatures from all channels and computed radiance statistics. Second, to estimate ROIs separability, Jeffries-Matusita distance (JM) was used. This average distance between two classes (Wacker and Landgrab, 1972) takes values in the range [0-2]. Value over 1.8 indicates a very good separability, a value under 1.8 indicates poor separability.

Third, supervised maximum likelihood classification was applied on filtered channel and on other channels of the image.

### III.2.2 Validation of classification

To validate produced maps, control areas were digitized with Arcgis GIS software in eleven as estimated for each image. The resulting manual classifications could then be crossed with the maximum likelihood classification result, to produce a confusion matrix and Kappa index.

### III.2.3 Analysis of maps accuracy

The conventional methods of accuracy assessment of thematic maps were employed: confusion matrix (Congalton, 1991). The confusion matrix gives an overall evaluation of map accuracy and for classification results of each thematic class.

Kappa index assesses the correlation between obtained results (maps produced) and the ground truth. Kappa takes values in the range [0-1] and it's divided into five categories: very low agreement from 0 to 0.20; weak correlation from 0.21 to 0.40; moderate correlation from 0.41 to 0.60; substantial correlation from 0.61 to 0.80; high correlation from 0.81 to 1. This index (equation 1) is expressed in terms of overall accuracy observed (equation 2) and expected (equation 3).

$$K = \frac{a - b}{1 - b} \qquad (1)$$

$$a = \frac{1}{N} \sum_{i=1}^{NC} x_{ii} \qquad (2)$$

$$b = \frac{1}{N^2} \sum_{i=1}^{NC} (x_{i+} \cdot x_{+i}) \qquad (3)$$

Where <NC> is the number of classes; <N> the total number of observations; <x<sub>ii</sub>> the number of observations in the row i; <x<sub>+i</sub>> and <x<sub>+i</sub>> total observations in the line i.

# IV RESULTS AND DISCUSSION

## IV.1 Classification by maximum likelihood and classes' separability

In general, the computed JM index to assess ROIs separability shows that, ROIs defined have good separability. Thus, JM distances are 1.99 for growing and young oil palms; but poor separability between forest class and mature oil palm. Confusions are expected between these last two classes (see Table 2).
## IV.2 Accuracy and thematic confusions analysis

For both Landsat 8 and ETM+ images, maps produced have very good accuracy (see Table 2) with 90% in 2001 and 80% in 2015, respectively. Different classes of palm grove are recognized: 99.67% for growing oil palm, 80% for mature oil palm and 93.42%, for mature oil palm in 2001. The same trend was obtained in 2015. Kappa values shows high correlation: $Kappa = 0.92$ (2001) and $Kappa= 0.86$ (2015).

However, despite of the high values of Kappa, some thematic confusions are observed between oil palm classes, or oil palm and vegetation classes. For example, we observed confusions between forest and mature oil palm (17%); or between young and growing oil palm (5.4%).

During the classification process, some pixels belonging to mature palm, where classified in to forest classes, for example.

The quality of classification is directly related to class separability, which itself depends on the variation among pixel from different classes as compared to within-class pixel variation. On one side, uncertainty to the data including several factors such as variable incidence, shading effects contribute to within class heterogeneity and alter the specificity of the signal.

On the other side, as some component may be common to different classes, (for example low vegetation pixels in the growing or young oil palm classes), distinct classes may share mixed pixels (Komba Mayossa, 2014).

On the validation map side, ground truth plays an important role in map accuracy. The knowledge of ground is an important key; the results were obtained from validation of digitized map and classification result. Landsat image, have 30mx30m of spatial resolution. The photo-interpretation is difficult especially in a heterogeneous landscape, as our study area where errors during the sampling step, cannot take to account within-plot heterogeneity. Thematic confusion caused by the presence of mixed pixels affects map accuracy (Chitroub, 2007). These mixed pixels result from spectral and spatial characteristics of the landscape studied objects, and from the method used to map.

## IV.3 Prospects

The maps obtained have a good accuracy, but some thematic confusions were observed.

One limitation of spatial accuracy of oil palm grove mapping is bound to class heterogeneity, attendant mixed pixels and the method used. Several ways for improvement are possible. At first, per-pixel classification assumes that each pixel represents a single class only. Maximum likely hood algorithm, used to map palm grove, ignores the mixed pixel problem. One of the solutions is to use a technic which allows mapping at the sub-pixel scale, such as super-resolution mapping (Priyaa and sanjeevi, 2013; Muad and foody, 2010). Second, as the map validation method improves or decreases the accuracy of produced map, the development of a classification system of cartographic products can be most interesting. Indeed, further studies should be focused on the evaluation of reliability of the produced map by this classification system (Chalifoux and al, 2006).
Map 2001

<table><thead><tr><th>GROUND TRUTH</th><th>Growing Oil palm</th><th>Young Oil palm</th><th>Mature Oil palm</th><th>Low Vegetation</th><th>Shady Forest</th><th>Forest</th><th>Ground Built1</th><th>Ground Built2</th><th>Ground Built3</th><th>Waterway</th><th>TOTAL</th></tr></thead><tbody><tr><td>Growing Oil palm</td><td>99.67</td><td>2.87</td><td></td><td>0.01</td><td></td><td></td><td></td><td></td><td></td><td></td><td>102.25</td></tr><tr><td>Young Oil palm</td><td>0.33</td><td>93.42</td><td></td><td>5.4</td><td></td><td></td><td></td><td></td><td></td><td></td><td>99.15</td></tr><tr><td>Mature Oil palm</td><td></td><td></td><td>80</td><td></td><td>10</td><td>2.22</td><td></td><td></td><td></td><td></td><td>92.32</td></tr><tr><td>Low Vegetation</td><td></td><td>3</td><td></td><td>94.59</td><td></td><td>1.70</td><td></td><td></td><td></td><td></td><td>100</td></tr><tr><td>Shady Forest</td><td></td><td></td><td>2.88</td><td></td><td>84.15</td><td>7.07</td><td></td><td></td><td></td><td>6</td><td>100.12</td></tr><tr><td>Forest</td><td></td><td></td><td>1.22</td><td></td><td></td><td>88.99</td><td></td><td></td><td></td><td></td><td>112.45</td></tr><tr><td>Ground Built1</td><td></td><td></td><td></td><td></td><td></td><td></td><td>100</td><td></td><td></td><td></td><td>106.45</td></tr><tr><td>Ground Built2</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>100</td><td></td><td></td><td>100</td></tr><tr><td>Ground Built3</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>93.55</td><td></td><td>93.55</td></tr><tr><td>Waterway</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>94</td><td>94</td></tr><tr><td>TOTAL</td><td>100</td><td>100</td><td>100</td><td>100</td><td>100</td><td>100</td><td>100</td><td>100</td><td>100</td><td>100</td><td>1000</td></tr></tbody></table>

Overall accuracy 90% / Kappa=0.92

Map 2015

<table><thead><tr><th>GROUND TRUTH</th><th>Growing Oil palm</th><th>Young Oil palm</th><th>Mature Oil palm</th><th>Low Vegetation</th><th>Shady Forest</th><th>Forest</th><th>Ground Built1</th><th>Ground Built2</th><th>Ground Built3</th><th>Waterway</th><th>TOTAL</th></tr></thead><tbody><tr><td>Growing Oil palm</td><td>100</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>100</td></tr><tr><td>Young Oil palm</td><td></td><td>100</td><td></td><td>0.04</td><td></td><td></td><td></td><td></td><td></td><td></td><td>100</td></tr><tr><td>Mature Oil palm</td><td></td><td></td><td>99.51</td><td></td><td>1.58</td><td></td><td></td><td></td><td></td><td></td><td>101.13</td></tr><tr><td>Low Vegetation</td><td></td><td></td><td></td><td>100</td><td></td><td></td><td></td><td></td><td></td><td></td><td>100</td></tr><tr><td>Shady Forest</td><td></td><td></td><td></td><td></td><td>97.55</td><td>12</td><td></td><td></td><td></td><td></td><td>99.15</td></tr><tr><td>Forest</td><td></td><td></td><td>0.49</td><td></td><td>2.01</td><td>97.22</td><td></td><td></td><td></td><td></td><td>99.63</td></tr><tr><td>Ground Built1</td><td></td><td></td><td></td><td></td><td></td><td></td><td>100</td><td></td><td></td><td></td><td>100</td></tr><tr><td>Ground Built2</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>100</td><td></td><td></td><td>100</td></tr><tr><td>Ground Built3</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>100</td><td></td><td>100</td></tr><tr><td>Waterway</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>100</td><td>94</td></tr><tr><td>TOTAL</td><td>100</td><td>100</td><td>100</td><td>100</td><td>100</td><td>100</td><td>100</td><td>100</td><td>100</td><td>1000</td><td></td></tr></tbody></table>

Overall accuracy 80% / Kappa=0.89

<caption>Table 2. Confusion matrix and Kappa index (Land cover map 2001 and 2015)</caption>
<figure><img src="image_219.png" /><figcaption>Figure 3: Landsat map in 2001(Landsat 7 ETM+), from CNRS ESPACE UMR 7300.</figcaption></figure>

<figure><img src="image_220.png" /><figcaption>Figure 4: Landsat map in 2015 (Landsat 8 OLI-TIRS), from CNRS ESPACE UMR 7300.</figcaption></figure>
## References

Elong J.G. (2003). Les plantations villageoises de palmier à huile de la Socapalm dans le bas-Moungo (Cameroun): un projet mal intégré aux préoccupations des paysans. *Les Cahiers d'Outre-Mer. Revue de géographie de Bordeaux* 56(224), 401-418.

Chitroub S. (2007). Annalyse des composantes indépendantes d'images multibandes: Faisabilité et perspectives. *Revue de télédétection* 7(1-4), 3-4.

Gadal S. (2003). Reconnaissances multi-niveaux d'unités paysagères par segmentation automatique d'images satellites. Télédétection des informations géographiques. Editions Anne-Elisabeth LAQUES.

Rival A., Levang P. (2013). *La palme de controverse : Palmier à huile et enjeux de développement*, édition Quae.

Priyaa A., Sanjeevi S. (2013). Super resolution mapping of multispectral and hyperspectral images of peechi reservoir, south india. *Image, 2010*.

Wacker A. G., Langreb D.A., 1972: Minimum distance classification in remote sensing. *LARS Technical Reports*, 25 pp.

Lacombe J.P. (2008). *Initiation au traitement d'images satellites: travaux dirigés, cahier2*. Ecole nationale Supérieure Agronomique de Toulouse.

Pouchin T., Debriej., Bourcier A. (2002). L'observation de la végétation de l'Afrique de l'Ouest par télédétection spatiale: l'apport de l'indice de végétation normalisé. *Science et changements planétaires/Sécheresse* 13(3), 187-94.

Koh L. P., Miettinen J., Liew S. C., Ghazoul, J. (2011). Remotely sensed evidence of tropical peatland conversion to oil palm. *Proceedings of the National Academy of Sciences*, 108(12), 5127-5132.

Li L., Dong J., Njeudeng Tenku, S., Xiao X. (2015). Mapping Oil Palm Plantations in Cameroon Using PALSAR 50-m Orthorectified Mosaïc Images. *Remote Sensing* 7(2), 1206-1224.

Domaç A., Zeydanli U., Yesilnacar E., Suzen M.L., 2004: Integration and usage of indices, feature components and topography in vegetation classification for regional, *20th ISPRS Congress, Istanbul, pp. 204-208*.