Proceedings of Spatial Accuracy 2016

Jean-Stéphane Bailly, Daniel A. Griffith, Didier Josselin

► To cite this version:

Jean-Stéphane Bailly, Daniel A. Griffith, Didier Josselin. Proceedings of Spatial Accuracy 2016. UMR 7300 ESPACE - Avignon. Spatial Accuracy, Jul 2016, Montpellier, France. , pp.366, 2016, 978-2-9105-4510-5. hal-01348458

HAL Id: hal-01348458

https://hal.science/hal-01348458

Submitted on 16 Sep 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
<figure><img src="image_2.png" /></figure>

# Proceedings of Spatial Accuracy 2016

Editors:
Jean-Stéphane Bailly, Daniel Griffith & Didier Josselin

[ISBN: 978-2-9105-4510-5]

organized by

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

INRA
SCIENCE & IMPACT

LISAH

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

AgroParisTech

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

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

ISARA
International Spatial Accuracy
Research Association
spatial-accuracy.org

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

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

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

<figure><img src="image_10.png" /></figure>
<figure><img src="image_11.png" /></figure>
<figure><img src="image_12.png" /></figure>

## Spatial uncertainty in knowledge-based systems

Is spatial information in ICT data reliable?

**Maxime Lenormand, Thomas Louail, Marc Barthelemy, José J. Ramasco** 9

Modelling, interpreting and visualizing uncertainties for the North Wyke Farm Platform baseline field surveys

**Paul Harris, Chris Brunsdon, Lex Comber, Hadewij Sint, Robert Orr, Michael Lee, Phil Murray** 18

Machine learning to deal with uncertainty in knowledge base for multivariate clustering applied to spatial analysis

**Stéphane Bourrelly, Antonino Marvuglia, Ian Vázquez-Rowe** 24

Generation of Plateau-Approximation Fuzzy Zones

**Hazaël Jones, Serge Guillaume, Patrice Loisel, Brigitte Charnomordic, Bruno Tisseyre** 31

Modeling process chain of SPOT images for resources uncertainty to monitor change in forest cover

**Aimé Richard Hajalalaina, Dominique Hérvé, Eric Delaitre, Thérèse Libourel** 38

## Spatial accuracy quantification in mapping

Uncertainty quantification of interpolated maps derived from observations with different accuracy levels

**Gerard B.M. Heuvelink, Dick Brus, Tom Hengl, Bas Kempen, Johan G.B. Leenaars, Maria Ruiperez-Gonzalez** 49

Survey designs which maximize efficiency gains in ALS-based forestry plot imputation

**Gavin Melville, Christine Stone, Jan Rombouts** 52

Validation of Copernicus High Resolution Layer on Imperviousness degree for 2006, 2009 and 2012

**Christophe Sannier, Javier Gallego, Jochen Dahmer, Geoff Smith, Hans Dufourmont, Alexandre Pennec** 60

The polygon overlay problem in electoral geography

**Romain Louvet, Jagannath Aryal, Didier Josselin, Christèle Marchand-Lagier, Cyrille Genre-Grandpierre** 67

Modelling Uncertainty using Geostatistics, a Case Study in Ecuador

**Elena Chicaiza, Buenaño Xavier** 74

A Hybrid approach for land cover mapping based on the combination of soft classifiers outputs and uncertainty information

**Luisa M.S. Gonçalves, Cidália C. Fonte** 80

Sampling to validate a global cropland map

**Javier Gallego, Anne Schucknecht, François Waldner** 88

Exploring the Usefulness of Land Parcel Data for Evaluating Multi-Temporal Built-up Land Layers

**Johanes Uhl, Stefan Leyk, Aneta J. Florczyk, Martino Pesaresi** 95

Urban objects recognition feasibilities by airborne hyperspectral and multispectral remote sensing

**Sébastien Gadal, Walid Ouerghemmi** 101

International Conference on Spatial Accuracy (2016): Space-Time Kriging of Temperature over Van-Turkey

**Pınar Aslantas, Okan Yeler** 109

## Spatial inference under uncertainty

Where (we think) the wild things are: comparing citizen generated and formal measures of wilderness

**Alexis Comber, Steve Carver** 117

Uncertain Clustering of Social Specialization in Metropolitan Areas

**Giovanni Fusco, Cristina Cao** 122
<figure><img src="image_13.png" /></figure>

<table><tr><td>Detection of outliers in crowdsourced GPS traces<br><b>Stefan Ivanovic, Ana-Maria Olteanu Raimond, Sébastien Mustière, Thomas Devogele</b></td><td>130</td></tr><tr><td>Locational error impacts on local spatial autocorrelation Indices: a Syracuse soil sample Pb-level data case study<br><b>Daniel Griffith, Yongwan Chun, Monghyeon Lee</b></td><td>136</td></tr><tr><td>Modeling spatial risk of the Foot-Mouth-Disease epidemic in South Korea<br><b>Eunhye Yoo, JiYoung Lee</b></td><td>144</td></tr><tr><td>Sensitivity of DBSCAN in identifying activity zones using online footprints<br><b>David W. S. Wong, Qunying Huang</b></td><td>151</td></tr><tr><td>On Reliability of Remote Sensing Data and Classification Methods for estimating transition rules of the land-use Cellular Automata<br><b>Yulia Grinblat, Michael Gilichinsky, Itzhak Benenson</b></td><td>157</td></tr><tr><td>Evaluating Performances of Spectral Indices for Burned Area Mapping Using object-based image analysis<br><b>Taskin Kavzoglu, Merve Yildiz Erdemir, Hasan Tonbul</b></td><td>162</td></tr><tr><td>Irregularly sampled data in space and time: using Poisson kriging to reduce the influence of uncertain observations in assessing the risk of flatoxin contamination of corn in Southern Georgia, USA<br><b>Ruth Kerry, Brenda Ortiz, Ben Ingram, Brian Scully, EunHye Yoo</b></td><td>169</td></tr><tr><td>A comparison of optimal map classification methods incorporating uncertainty information<br><b>Yongwan Chun, Hyeongmo Koo, Daniel A. Griffith</b></td><td>177</td></tr></table>

## Scale effects in spatial accuracy

<table><tr><td>Getting the right spatial mix: optimising the size, type and location of renewable energy facilities<br><b>Alexis Comber, Jen Dickiie</b></td><td>185</td></tr><tr><td>Impact of compositional and configurational data loss on downscaling accuracy<br><b>Amy E. Frazier, Peter Kedron</b></td><td>190</td></tr><tr><td>A Hierarchical Scale Setting Strategy for Improved Segmentation Performance Using Very High Resolution Images<br><b>Taskin Kavzoglu, Merve Yildiz Erdemir</b></td><td>195</td></tr><tr><td>Exploring Accurate Spatial Downscaling using Optimization<br><b>Michael Poss, Didier Josselin</b></td><td>202</td></tr><tr><td>Downscaling of soil moisture with area-to-point geographically weighted regression kriging and uncertainty analysis<br><b>Yan Jin, Yong Ge, Jianghao Wang, Yuehong Chen, Xudong Guanong Ge</b></td><td>211</td></tr></table>

## Uncertainty in 3D geodata

<table><tr><td>A proposal of a statistical test to control positional accuracy by means of 2 tolerances simultaneously<br><b>José Rodríguez-Avi, Francisco Javier Ariza-López</b></td><td>221</td></tr><tr><td>Practical global elevation data error simulation<br><b>Ashton Shortridge, Joe Messina, Xue Li, Nick Ronnei</b></td><td>229</td></tr><tr><td>How far is far enough? Towards an adaptive and "site-centric" modelling integrating co-visibility constraints for optimal land use<br><b>Valerio Signorelli, Thomas Leduc, Guillaume Chauvat</b></td><td>233</td></tr><tr><td>Combination of error characterization and spatial error model to improve quality of digital elevation models<br><b>Tomaž Podobnikar</b></td><td>241</td></tr></table>
<figure><img src="image_14.png" /></figure>

## Spatial uncertainties simulation and propagation

<table><tr><td>Exploring the uncertainty of soil water holding capacity information<br>Linda Lilburne, Stephen McNeill, Tom Cuthill, Pierre Roudier</td><td>251</td></tr><tr><td>Simulation of realistic digitizing errors on geographical objects using constraints modeling the operator input process<br>Jean-François Girres</td><td>259</td></tr><tr><td>How far spatial accuracy governs land-use changes monitoring frequency: the urban sprawl monitoring example<br>Jean-Pierre Chéry, Jean-Stéphane Bailly, Valérie Laurent, Nathalie Saint-Geours</td><td>267</td></tr><tr><td>'spup' – a R package for uncertainty propagation in spatial environmental modelling<br>Kasia Sawicka, G.B.M. Heuvelink</td><td>275</td></tr><tr><td>Combining spatial and thematic uncertainty and sensitivity analysis for mountain natural hazards assessment<br>Jean-Marc Tacnet, Guillaume Dupouy, Franck Bourrier, Frédéric Berger, Dominique Laigle, Nicolas Crimier, Kamel Mekhnacha, Michel Mémier</td><td>283</td></tr><tr><td>Performing Multi-Temporal Spatial Data Analysis for Coastal Areas and Assessing Thematic Accuracy<br>Daniel Cohenca, Carlos Antonio Oliveira Vieira</td><td>286</td></tr></table>

## Spatial model sensitivity analyses

<table><tr><td>Multiway sensitivity analysis of the fusion of earth observation, topography and social media data for rapid flood mapping<br>Didier G. Leibovici, Julian F. Rosser</td><td>297</td></tr><tr><td>A GPU-based Solution for Accelerating Spatially-Explicit Uncertainty- and Sensitivity Analysis in Multi-Criteria Decision Making<br>Christoph Erlacher, Seda Şalap-Ayça, Piotr Jankowski, Karl-Heinrich Anders, Gernot Paulus</td><td>305</td></tr><tr><td>Uncertainty propagation in urban hydrology water quality modelling<br>Jairo Arturo Torres Matallana, U. Leopold, G.B.M. Heuvelink</td><td>313</td></tr><tr><td>How can big data be used to reduce uncertainty in stormwater modelling?<br>Nanée Chahinian, Anne-Laure Piat-Marchand, Sandra Bringay, Maguelonne Teisseire, Elodie Boulogne, Laurent Deruelle, Mustapha Derras, Carole Delenne</td><td>322</td></tr><tr><td>Sensitivity analysis of spatio-temporal models describing nitrogen transfers, transformations and losses at the landscape scale<br>Jordi Ferrer Savall, Cyril Benhamou, Pierre Barbillon, Patrick Durand, Marie-Luce Taupin, Hervé Monod, Jean-Louis Drouet</td><td>330</td></tr><tr><td>The transformed optimal transportation problem: sensitivity and segregation of the children-to-school constrained assignment in Lausanne<br>Théophile Emmanouilidis, Guillaume Guex, François Bavaud</td><td>333</td></tr></table>

## Poster Session

<table><tr><td>A method for testing the similarity of spatial samples<br>Maria Virtudes Alba-Fernández, Francisco J. Ariza-López, José Rodríguez-Avi</td><td>343</td></tr><tr><td>Combining punctual and ordinal contour data for accurate floodplain topography mapping<br>Carole Delenne, Jean-Stéphane Bailly, Mathieu Dartelle, Nelly Marcy, Antoine Rousseau</td><td>350</td></tr><tr><td>Monitoring spatial accuracy of oil palm cultivation mapping in southern Cameroun from Landsat series images<br>Prune Christobelle Komba Mayossa, Sébastien Gadal</td><td>358</td></tr></table>
<figure><img src="image_15.png" /></figure>
<figure><img src="image_16.png" /></figure>

SPATIAL
ACCURACY
2016

5-8 JULY 2016

MONTPELLIER, FRANCE

# Spatial uncertainty
in knowledge-based systems
<figure><img src="image_17.png" /></figure>
# Is spatial information in ICT data reliable?

Maxime Lenormand¹,*, Thomas Louail²,³, Marc Barthelemy⁴,⁵, José J. Ramasco²

¹Iirstea, UMR TETIS, 500 rue François Breton, FR-34093 Montpellier, France

²IFISC (CSIC-UIB), Campus UIB, E-07122 Palma de Mallorca, Spain

³Géographie-Cités (CNRS - Paris 1 - Paris 7), 13 rue du four, FR-75006 Paris, France

⁴Institut de Physique Théorique, CEA-CNRS (URA 2306), F-91191, Gif-sur-Yvette, France

⁵CAMS, EHESS-CNRS (UMR 8557), 190-198 avenue de France, FR-75013 Paris, France

*Corresponding author: maxime.lenormand@irstea.fr

## Abstract

While an increasing number of human activities are studied using data produced by individuals' ICT devices, there have been relatively few contributions investigating the robustness of results against fluctuations of data characteristics. In particular, when ICT data contain spatial information, they represent an invaluable new source for analyzing urban phenomena. Here, we present a stability analysis of higher-level information extracted from mobile phone metadata passively produced during an entire year by 9 million individuals in Senegal. We focus on two specific information-retrieval tasks: (a) the identification of land use in the region of Dakar by analyzing the temporal rhythms of the communication activity; (b) the identification of home and work locations of anonymized individuals, allowing for the construction of the Origin-Destination (OD) matrices for commuting flows. Our analysis reveals that the spatial distributions of land use computed for different samples are remarkably robust, with on average 80% of shared surface area between the different spatial partitions. The OD matrix is less stable with a share of about 70% of commuters in common when considering all types of flows. Better results can be obtained at larger levels of aggregation. These different results confirm that ICT data are mostly a very useful source for the spatial analysis of urban systems, but that their reliability should be tested more thoroughly.

## I INTRODUCTION

Massive amounts of geolocated data are passively and continuously generated by individuals when they use their mobile and connected devices: smart phones, credit cards, GPSs, RFIDs or remote sensing devices. This deluge of digital footprints is growing at an extremely fast pace and represents an unprecedented opportunity for researchers, to address quantitatively challenging problems, in the hope of unveiling new insights on the dynamics of human societies. Many fields are concerned by the development of new techniques to handle these vast datasets, and range from applied mathematics, physics, to computer science, with plenty of applications to a variety of disciplines such as medicine, public health and social sciences for example.

Although big data have the advantage of large samples sizes (millions of observations), and high spatio-temporal resolution, they also raise new challenging issues. Some are technical and related to the storage, management and processing of these data (Kaisler et al., 2013), and others are methodological, such as the statistical validity of analysis performed on such data. For example, in the case of mobile phone data, researchers have often no control on the data collection that is usually made for other purposes. In many cases various hidden biases can affect the spatial behavior of anonymized individuals. Observing the world through the lenses
of data generated by information and communications technologies (ICT) may therefore lead to
possible distortions, and possibly to erroneous conclusions (Lewis, 2015). It is thus crucial to
perform statistical tests and to develop methods in order to assess the robustness of the results
obtained with ICT data. In the research community that studies human mobility (in particular
in cities, and related urban dynamics Ratti et al. (2006); Louail et al. (2014); Calabrese et al.
(2015); Louail et al. (2015), efforts in this sense have been made in recent years, notably by
comparing the results obtained with different data sources (Tizzoni et al., 2014; Lenormand
et al., 2014; Deville et al., 2014; Alexander et al., 2015; Toole et al., 2015). However, the
robustness of results to sample selection has, to our knowledge, never been studied yet.

In the following we present two examples of such uncertainty analysis on results obtained with
mobile phone data recorded in Senegal in 2013 (de Montjoye et al., 2014). We concentrate on
two information-retrieval tasks: first, we evaluate the uncertainty when inferring land use from
the rhythms of human activity (Soto and Frías-Martínez, 2011; Frías-Martínez et al., 2012;
Toole et al., 2012; Pei et al., 2014; Lenormand et al., 2015), and second, we study the un-
certainty when identifying individuals' most frequented locations (Ahas et al., 2010; Isaacman
et al., 2011; Lenormand et al., 2014; Toole et al., 2015). We conclude by mentioning possible
future steps to clearly assess the relevance of various ICT data sources for studying different
phenomena.

## II STUDY AREA AND DATA DESCRIPTION

We focus here on the region of Dakar, Senegal. The mobile phone data consists in call detail records (CDR) of phone calls and short messages exchanged by more than 9 million of anonymized Orange's customers. They were collected in Senegal in 2013, and were released to research teams in the framework of the 2014 Orange Data for Development challenge (de Montjoye et al., 2014). We will use for our study the second dataset (SET2) that was made available by Orange. It contains fine-grained location data on a rolling 2-week basis at the individual level. For each of the 25 two-week periods, a sample of about 300,000 mobile phone users is selected at the country scale. Whenever one of these individuals uses his/her mobile phone during the two-week period, the time and his/her position (at the level of serving cell tower) is recorded. This information can be used to study human activity and mobility patterns in the region of Dakar that is here divided into 457 spatial subunits. The partition is the Voronoi tessellation constructed from the location of antennas in the city chosen as nodes. Each Voronoi cell approximates thus the activity zone served by the antenna located at its center (see Figure 1a).

## III INFERRING LAND USE FROM MOBILE PHONE ACTIVITY

Geolocalized ICT data have been widely used to infer land use from human activity (Soto and Frías-Martínez, 2011; Frías-Martínez et al., 2012; Toole et al., 2012; Pei et al., 2014; Lenormand et al., 2015). The basic idea is to divide the region of interest into zones, then extract a temporal signal of activity for each of these zones, and finally cluster together zones that display similar signals. Each of these clusters corresponds to a certain type of activity (Residential, Commercial, ...). We use here the functional approach proposed in Lenormand et al. (2015). The method takes as input, for each cell, a signal composed of 168 points (24h×7days), each value corresponding to the number of users located in this cell, at this hour of the day and this day of the week. These signals are then normalized by the total hourly activity, in order to subtract trends introduced by circadian rhythms. A Pearson correlation matrix between cells is then computed. Two spatial units whose activity is strongly correlated in time will have a high
positive correlation value. This similarity matrix can be represented by a undirected weighted network, which is then clustered using the Infomap community detection algorithm (Rosvall and Bergstrom, 2008). This method has the advantage to be non-parametric (the number of clusters is not fixed a priori).

In order to extract temporal signals of activity in the region of Dakar, we first need to estimate the number of people in each zone, per hour, for each of the 350 days of our original sample. To do so, we rely on the following criteria: each person counts only once per hour. If a user is detected in *k* different zones within a certain 1-hour time period, each registered position will count as $(1/k)$ 'units of activity' for each of these *k* cells. The average number of users per hour and per day is 30,500, which represents about 1% of the total population of the region of Dakar.

To assess the robustness of land use identification from mobile phone activity to sample selection, we apply the functional approach described above to 100 weeks drawn at random from our original sample (after removing outliers). Note that the days of the weeks are drawn separately, and they are therefore not necessarily chronologically ordered.

<figure><img src="image_18.png" /><figcaption>Figure 1: (a) Map of the region of Dakar displaying the three clusters. Colors vary from white to the most recurrent cluster identified in the 100 random sample. The color saturation depends on the number of times the zone was classified as the most recurrent cluster. The color code is red for *Residential*, blue for *Business*, green for *Nightlife* and orange for other types of land use. (b) Temporal patterns associated with the three clusters. The solid lines represent the average temporal profile computed over 100 random extractions, while the dashed lines represent the standard deviations.</figcaption></figure>

First, we observe that three clusters emerge systematically, and represent on average 95% of the total surface. The remaining 5% correspond to other clusters with no clear patterns, probably associated with some local one-time events. We show on Figure 1b the average temporal profiles and the variability around this average for each of these three clusters. Each of the clusters can be roughly associated with a dominant land use:

* A Residential activity corresponds to a high probability of mobile phone use during early mornings, evenings and week end days.

* A Business cluster displays a significantly higher activity from 9am to 5-6pm during weekdays.

* A Nightlife activity profile is characterized by a high activity during night hours (1am-4am).

The Nightlife cluster (in green) covers the area of the international airport, and also the neighborhood of "la Pointe des Almadies", where live mainly wealthy people – most of the rich/sophisticated nightclubs are also located there. The Business cluster covers Dakar's CBD ("Le
Plateau”), where are located most of the companies headquarters, and also the port. Finally the residential cluster cover the rapidly growing parts of Dakar peninsula, who profit from the highway construction. It is worth noting that the results are consistent with the ones obtained with another mobile phone dataset in Spain (Lenormand et al., 2015), except that in the case of Dakar, the method is not able to distinguish between industrial (or logistic) and leisure nightlife activities (see Lenormand et al. (2015) for more details).

<figure><img src="image_19.png" /><figcaption>Figure 2: Uncertainty when inferring land use from mobile phone activity. (a) Area covered by the different land use types, expressed as a percentage of the total surface. The values have been averaged over 100 random samples, and the error bars represent the standard deviation. (b) Probability density function of the shared surface area between each pair of spatial distributions according to the type of land use.</figcaption></figure>

As can be observed in Figure 2a, the area covered by the different types of land use is quite stable over the 100 samples, with the residential land use type representing on average about 55% of the total surface, while we observe about 20% for the Business and Nightlife clusters. Nevertheless, the stability of the proportion does not imply that they follow the same spatial distribution from one sample to another. To check for the stability, we compute the surface area shared by two spatial distributions $S_l$ and $S'_l$ of a given type $l$ and obtained with two different samples. The expression for this quantity is

$$SSA(S_l, S'_l) = 2 \frac{A_{S_l \cap S'_l}}{A_{S_l} + A_{S'_l}}, \qquad (1)$$

where $A_S$ denotes the surface area of spatial distribution $S$. Note that in our case $A_{S_l} \approx A_{S'_l}$ (Figure 2a). Similarly, we can defined the total surface area shared by two spatial partitions $P$ and $P'$ (with the same number and types of land use) of the region of interest,

$$SSA(P, P') = \frac{\sum_l A_{S_l \cap S'_l}}{\sum_l A_{S_l}}. \qquad (2)$$

The results are displayed in Figure 2b. The similarities between the 100 different spatial distributions is globally high, with on average 80% of shared surface area between the spatial partitions. The agreement is larger for residential and business clusters with an average shared surface area around 90%. For the nightlife land use type we find a result about 75%.
A map of the region of Dakar displaying the uncertainty associated with the land use identification is shown in Figure 1a. The colors represent the different land use types, where a zone has been assigned to the most recurrent cluster type over the 100 land use identifications. The color saturation is related to the uncertainty quantified with the number of times the zone was classified as given recurrent cluster: the color is darker if the uncertainty is low, paler otherwise. Most of the zones have been assigned the same clusters more than 80% of the time. This leads us to the conclusion that the identification of land use from mobile phone activity shows a high level of robustness to sample selection.

## IV IDENTIFYING HOME-WORK LOCATIONS FROM MOBILE PHONE ACTIVITY

Geolocalized ICT data are also widely used to identify the most visited locations of an individual during his/her daily life trajectory, allowing to extract the origin-destination (OD) matrices of commuting flows, a fundamental object in mobility studies. A simple heuristic is that the most frequented place of a user in the late afternoon/evening and in the early morning can be identified as a proxy for his/her place of residence, while the most frequented area during working hours can be a proxy to his/her work (or activity) place. This simple assumption allows the accurate determination of mobility flows at intermediate geographical scales (see for example Tizzoni et al. (2014); Lenormand et al. (2014); Alexander et al. (2015); Toole et al. (2015)). However, the robustness of the results to sample selection has never been investigated.

For each of the 25 two-week periods and for each user, we apply the following home and work location extraction procedure:

* For each hour of the two weeks period (weekends excepted) during which an individual used his/her mobile phone, we identify the most frequently visited zone during this hour (based on his/her geolocalized mobile phone activity).

* Hours of activity are then divided into two groups, daytime hours (between 8am and 5pm) and nighttime hours (between 7pm and 7am).

* **Filter 1:** We keep only the users who have been 'active' at least ten 1-hour periods during daytime, and ten 1-hour periods during nighttime (spread over at least half of the days of the two-week period).

* For both groups of hours (daytime and nighttime), we identify the spatial unit in which the user has been localized the highest number of hours.

* **Filter 2:** We keep only users with a fraction of hours spent at "home" location and "work" location larger than one third of the total number of locations visited during nighttime and daytime, respectively.

The two filters allow us to discard users not showing enough regularity in order to estimate their main nighttime ('Home') and daytime activity ('Work') locations. The source code of this method is available online<fnref n="1" />.

At the end of the process, after filtering out users living and/or working outside the region of Dakar, the remaining number of users is on average ≈ 65,000. This number is quite stable over the 25 two-weeks periods, varying at most by 15% around this average. The resulting 25 two-weeks commuting networks can then be compared using several similarity metrics, such as the one described in Lenormand et al. (2016). We consider 2 commuting networks *T* and *T*', where *T*<sub>*ij*</sub> is the number of users living in zone *i* and working in zone *j*, and we will use three different metrics, that encode different network properties. First, the common fraction of

<footnote n="1">¹https://github.com/maximelenormand/Most-frequented-locations</footnote>
commuters (CPC), varying from 0, when no agreement is found, to 1, when the two networks are identical, is estimated as

$$CPC(T, T') = \frac{2 \sum_{i,j=1}^{n} \min(T_{ij}, T'_{ij})}{\sum_{i,j=1}^{n} T_{ij} + \sum_{i,j=1}^{n} T'_{ij}} \quad (3)$$

Second, we will consider the common proportion of links (CPL) that measures similarity in the networks' topological structure, and is calculated as

$$CPL(T, T') = \frac{2 \sum_{i,j=1}^{n} 1_{T_{ij}>0} \cdot 1_{T'_{ij}>0}}{\sum_{i,j=1}^{n} 1_{T_{ij}>0} + \sum_{i,j=1}^{n} 1_{T'_{ij}>0}} \quad (4)$$

Third, the common share of commuters according to the distance (CPCd), assessing the similarity between commuting distance distributions is given by

$$CPC_d(T, T') = \frac{\sum_{k=1}^{\infty} \min(N_k, N'_k)}{N} \quad (5)$$

where $N_k$ stands for the number of users with a commuting distance ranging between $2k-2$ and $2k$ kms.

<figure><img src="image_20.png" /><figcaption>Figure 3: Uncertainty in the identification of the most frequently visited locations from mobile phone activity. (a) Boxplots of the CPC (blue), CPL (red) and CPCd (green) between the 25 two-weeks commuting networks (1x2), 12 one-month commuting networks (2x2) and 6 two-months commuting networks (4x2). One-month and two-months commuting networks have been obtained by aggregation of consecutive two-weeks commuting networks. (b) Boxplots of the CPC (blue), CPL (red) and CPCd (green) between the 25 two-weeks commuting networks, according to the grid cell side length (1 km, 2 km and 3 km). The CPC values obtained when considering only inter-zonal flows are also displayed.</figcaption></figure>

The boxplots for the CPC, CPL and CPCd values obtained by comparing the 25 two-weeks commuting networks are displayed in Figure 3a (1x2). The results are not completely conclusive, with different commuting networks showing a good agreement around 70% of commuters in common (considering both inter- and intra-zonal flows), but with a value falling down to 35% when only inter-zonal flows are considered. The CPL values are also quite low, around 35% of links are in common between the different networks. An encouraging result is that the common part of commuters according to the distance is very high, showing around 90% of similarity
between the 25 commuting distance distributions. However, it is important to keep in mind
that these mixed results are obtained with a few thousand users for each two-weeks commuting
network, drawn at a high spatial resolution with an average surface area equal to 0.5 km². In-
deed, as it can be observed in Figure 3, both temporal and spatial aggregations of the networks
greatly improve the results. The aggregation of consecutive two-weeks commuting networks
allow us to perform pairwise comparisons of 12 one-month commuting networks (2 × 2) and
6 two-months commuting networks (4), and we clearly observe that the temporal aggregation
improves the results. The two-weeks commuting networks can also be aggregated spatially,
by projecting the data on a regular lattice (see Lenormand et al. (2014) for more details about
the aggregation method), and results of the comparison of the two-weeks commuting networks
according to the size of the grid cells (resp. 1 km, 2 km and 3 km) are represented in Figure 3b.
Here again, we observe a significant improvement, with CPC values almost always larger than
0.75.

In order to go further, for each of the 25 two-week periods and for each user, we identify the
home and work locations for each of the two weeks considered independently, by following
the procedure described above. This allows us to assess the influence of the sampling of points
along individual trajectories when identifying the home and work locations. We then compare
the locations identified for each of the two weeks. Considering the high spatial resolution and
the small time window, a good agreement is obtained, with an accuracy of 85.3 ± 1.3% for
home (average ± standard deviation over the 25 two-week periods) and 79.1 ± 2.7% for work.
Moreover, 60% of the inaccurate locations are less than 2 kms distant from each other. We
can therefore conclude that the identification of users' home-work locations from mobile phone
activity also shows a high level of robustness to sample selection.

V DISCUSSION

Data passively produced through information and communications technologies have been in-
creasingly used by researchers since the middle of the 2000's to analyze a variety of human
processes and activities. In particular, our understanding of human mobility has been deeply re-
newed thanks to these new sources. The longitudinal tracking of anonymized individuals opens
the door to an enhanced understanding of human and social phenomena that could not be stud-
ied empirically with such a level of detail. However, these data may obviously suffer from a
number of biases (Lewis, 2015), which include in particular sample selection. Systematic tests
are then required for ensuring statistical validity, along with cross-checks between various data
sources.

With this in mind, we performed two uncertainty analysis of results obtained with mobile phone
data produced by millions of anonymized individuals and collected during an entire year. In the
first part of the analysis, we assessed the uncertainty when inferring land use from human activ-
ity, estimated from the number of mobile phone users at different moments of the week. A good
agreement was obtained between the land uses identified from 100 randomly selected samples
of individuals, with on average 80% of shared surface area between land uses in the resulting
maps. In the second part of the analysis, we investigated the influence of sample selection on
the identification of users' home and work locations. We first examined the impact of the se-
lection of users on the journey-to-work commuting networks extracted at the city scale. In our
case-study of the city of Dakar, we showed that the level of uncertainty was highly dependent
of the spatio-temporal resolution, and that good results were reachable with a reasonable level
of aggregation. We then analyzed the effect of the sampling of locations along mobile phone
users' trajectories on the identification of their home and work locations. Most of the locations identified with different samples were the same, or very close to one another.

For these two spatial information retrieval tasks, our results suggest that the level of uncertainty associated with sample selection is low. Further work in this direction include the reproduction of such uncertainty analysis with other datasets coming from different countries and data sources. An important aspect of the rapidly growing 'new science of cities' (Batty, 2013), which heavily relies on new data sources, is to be able to reproduce results with different datasets, and to characterize and control to what extent the information provided by different sources are biased in a particular direction.

More studies in this spirit need to be done to strengthen the foundations of the field dedicated to the understanding of urban mobility and urban dynamics through ICT data. From a publication point of view, trying to reproduce previous results with different data sources, or to estimate the robustness of previously published results, might not be as appealing as proposing new measures and models, but is crucially important as well.

## ACKNOWLEDGMENTS

Partial financial support has been received from the Spanish Ministry of Economy (MINECO) and FEDER (EU) under project ESOTECOS (FIS2015-63628-C2-2-R), and from the EU Commission through project INSIGHT. The work of TL has been funded under the PD/004/2015, from the Conselleria de Educación, Cultura y Universidades of the Government of the Balearic Islands and from the European Social Fund through the Balearic Islands ESF operational program for 2013-2017. JJR acknowledges funding from the Ramón y Cajal program of MINECO.

## References

* Ahas R., Silm, S. and J. O., Saluveer E., Tiru M. (2010). Using Mobile Positioning Data to Model Locations Meaningful to Users of Mobile Phones. *Journal of Urban Technology* 17(1), 3–27.

* Alexander L., Jiang S., Murga M., González M. C. (2015). Origin–destination trips by purpose and time of day inferred from mobile phone data. *Transportation Research Part C: Emerging Technologies* 58, Part B, 240–250.

* Batty M. (2013). *The New Science of Cities*. MIT Press.

* Calabrese F., Ferrari L., Blondel V. D. (2015). Urban sensing using mobile phone network data: a survey of research. *ACM Computing Surveys (CSUR)* 47(2), 25.

* de Montjoye Y.-A., Smoreda Z., Trinquart R., Ziemlicki C., Blondel V. D. (2014). D4D-Senegal: The Second Mobile Phone Data for Development Challenge. *arXiv preprint arXiv:1407.4885*.

* Deville P., Linard C., Martin S., Gilbert M., Stevens F. R., Gaughan A. E., Blondel V. D., Tatem A. J. (2014). Dynamic population mapping using mobile phone data. *Proceedings of the National Academy of Sciences* 111(45), 15888–15893.

* Frías-Martínez V., Soto V., Hohwald H., Frías-Martínez E. (2012). Characterizing urban landscapes using geolocated tweets. In *SocialCom/PASSAT*, pp. 239–248. IEEE.

* Isaacman S., Becker R., Cáceres R., Kobourov S., Martonosi M., Rowland J., Varshavsky A. (2011). Identifying Important Places in People's Lives from Cellular Network Data. In K. Lyons, J. Hightower, and E. M. Huang (Eds.), *Pervasive Computing*. Springer Berlin Heidelberg.

* Kaisler S., Armour F., Espinosa J. A., Money W. (2013). Big Data: Issues and Challenges Moving Forward. In *2014 47th Hawaii International Conference on System Sciences*, Volume 0, pp. 995–1004.

* Lenormand M., Bassolas A., Ramasco J. J. (2016). Systematic comparison of trip distribution laws and models. *Journal of Transport Geography* 51, 158–169.

* Lenormand M., Picornell M., Cantú-Ros O. G., Tugores A., Louail T., Herranz R., Barthelemy M., Frías-Martínez E., Ramasco J. J. (2014). Cross-Checking Different Sources of Mobility Information. *PLoS ONE* 9(8), e105184.
Lenormand M., Picornell M., Garcia Cantú O., Tugores A., Louail T., Herranz R., Barthelemy M., Frías-Martínez E., Ramasco J. J. (2015). Comparing and modeling land use organization in cities. *Royal Society Open Science* 2, 150459.

Lewis K. (2015). Three fallacies of digital footprints. *Big Data & Society* 2(2).

Louail T., Lenormand M., Picornell M., Cantú O. G., Herranz R., Frias-Martinez E., Ramasco J. J., Barthelemy M. (2015). Uncovering the spatial structure of mobility networks. *Nature Communications* 6.

Louail T., Lenormand M., Ros O. G. C., Picornell M., Herranz R., Frias-Martinez E., Ramasco J. J., Barthelemy M. (2014). From mobile phone data to the spatial structure of cities. *Scientific reports* 4.

Pei T., Sobolevsky S., Ratti C., Shaw S. L., Zhou C. (2014). A new insight into land use classification based on aggregated mobile phone data. *International Journal of Geographical Information Science* 28, 1988–2007.

Ratti C., Frenchman D., Pulselli R. M., Williams S. (2006). Mobile landscapes: using location data from cell phones for urban analysis. *Environment and Planning B: Planning and Design* 33(5), 727–748.

Rosvall M., Bergstrom C. T. (2008). Maps of random walks on complex networks reveal community structure. *Proceedings of the National Academy of Sciences* 105(4), 1118–1123.

Soto V., Frías-Martínez E. (2011). Automated land use identification using cell-phone records. In *Proceedings of the 3rd ACM international workshop on MobiArch*, HotPlanet ’11, New York, NY, USA, pp. 17–22. ACM.

Tizzoni M., Bajardi P., Decuyper A., King G. K. K., Schneider C. M., Blondel V., Smoreda Z., González M. C., Colizza V. (2014). On the Use of Human Mobility Proxies for Modeling Epidemics. *PLOS Comput Biol* 10(7), e1003716.

Toole J. L., Colak S., Sturt B., Alexander L. P., Evsukoff A., González M. C. (2015). The path most traveled: Travel demand estimation using big data resources. *Transportation Research Part C: Emerging Technologies* 58, Part B, 162–177.

Toole J. L., Ulm M., González M. C., Bauer D. (2012). Inferring Land Use from Mobile Phone Activity. In *Proceedings of the ACM SIGKDD International Workshop on Urban Computing*, UrbComp ’12, New York, NY, USA, pp. 1–8. ACM.
# Modelling, interpreting and visualizing uncertainties for the North Wyke Farm Platform baseline field surveys

Paul Harris*<fnref n="1" />, Chris Brunsdon<fnref n="2" />, Lex Comber<fnref n="3" />, Hadewij Sint<fnref n="1" />, Robert Orr<fnref n="1" />,

Michael Lee<fnref n="1" />, Phil Murray<fnref n="1" />

*Corresponding author: paul.harris@rothamsted.ac.uk

**Abstract**

This study demonstrates a new approach for the visualization of spatial uncertainty using data from three agricultural field surveys.

**Keywords**

caricRture; kriging; geographically weighted models; grasslands research

## I INTRODUCTION

The North Wyke Farm Platform (NWFP) is a systems-based, farm-scale experiment with the aim of addressing grassland agricultural productivity and ecosystem responses to different management practices. The 63 ha site captures the data necessary to develop a better understanding of the dynamic processes and underlying mechanisms that can be used to model how agricultural grassland systems respond to different management inputs.

Via cattle beef and sheep production, the underlying principle is to manage each of three farmlets in three contrasting ways: (i) improvement of permanent pasture (i.e. ‘business as usual’ green farmlet); (ii) improvement through the use of legumes (blue farmlet); and (iii) improvement through innovation (red farmlet). The connectivity between the timing and intensity of the different management operations, together with the transport of nutrients and potential pollutants from the NWFP is evaluated using numerous inter-linked data collection exercises, operating at various spatial and/or temporal scales. Fig. 1a maps the NWFP experiment, where its 15 hydrologically-isolated sub-catchments are shown; some of which consist of multiple fields. In this study, we introduce some of the modelling and visualization opportunities that are possible with this rich data resource, with respect to baseline field survey data only.

## II METHODS

We study three surveys that were conducted in the summers of 2012 and 2013. Baseline data entails that the *blue* and *red* treatments of NWFP farmlets are not as yet in place (thus the map in Fig. 1a is effectively all *green* during these surveys). Surveys sampled over a mixture of 25m and 50m grids for: (a) plant nutrients in 2013 (n = 544); (b) plant species in 2013 (n = 294); and (c) soil nutrients in 2012 (n = 250). These data sets (and many more) are freely available at www.rothamsted.ac.uk/farmplatform.

For plant nutrients, bulk stable isotopes of Carbon (δ¹³C) and Nitrogen (δ¹⁵N) along with Total Carbon (C), Total Nitrogen (N) and sward height were measured. For plant species, 18 species

<footnote n="1">¹ Rothamsted Research, North Wyke, Devon, EX20 2SB, UK</footnote>

<footnote n="2">² National Centre for Geocomputation, Maynooth Unuveristy, Ireland</footnote>

<footnote n="3">³ University of Leeds, Leeds, West Yorkshire, LS2 9JT, UK</footnote>
were observed, where two species, Agrostis stolonifera and Loium perenne, clearly dominated. For soils nutrients, δ¹³C, δ¹⁵N, C and N, together with Bulk Density, Soil Organic Matter and pH were measured (Fig. 1b). Surveys can be analysed as a self-contained entity or related to each other.

Our study implements geostatistical and geographically weighted (GW) models to these surveys. Outputs are visualized to reflect different aspects of uncertainty, showcasing functions provided in the caricRture R package (Brunsdon 2016). Functions extend the ‘sketchy’ rendering techniques of Wood et al. (2012) to use with spatial data.

Using caricRture, we map a model output of interest (predicted/estimated to a finer 15m grid) via short lines that are coloured to demark the output range. To reflect uncertainty, the lines are drawn with varying degrees of ‘sketchiness’, where straight lines suggest relatively low levels of uncertainty, whilst highly sketchy (‘hand-drawn’) lines suggest relatively high levels of uncertainty. The technique provides a useful way of communicating uncertainty, where information that is commonly presented as two maps, are presented using only one. We present four case studies, each with a different visualization theme.

<figure><img src="image_21.png" /><figcaption>Figure 1:(a) The North Wyke Farm Platform and (b) Map depicting the C/N ratio for soils.</figcaption></figure>

## III CASE STUDIES

### Visualization of spatial prediction uncertainty

In our first visualization, we predict C from the plant nutrients survey using universal kriging (UK) with a linear drift. Instead of using the UK variance as a measure of prediction uncertainty, we use the UK interpolation variance of Yamamoto (2000) which, unlike the UK variance, accounts for local changes in sample variance. The method is sensitive to the kriging neighbourhood size, and we choose N = 20. Fig.2a provides the UK prediction map, Fig. 2b provides the corresponding UK interpolation variance map, and Fig. 6a. provides the resultant ‘sketchy’ map - combining both UK predictions and their uncertainty. Thus the prediction levels
remain the same for maps in Fig 2a and Fig. 6a, whereas the high interpolation variances coloured dark green in Fig. 2b correspond to areas of high ‘sketchiness’ in Fig. 6a. Thus plant C predictions should be viewed cautiously in sub-catchments 2 and 7 (see also Fig. 1a).

## Visualization of species metric uncertainty

Next, we suggest how to visualize species richness together with species diversity (for the plant species survey). Species diversity is a function the number of species present (i.e. species richness) and the evenness of how the individuals are spread (Hurlbert 1971). We choose to view species richness as our main output (Fig. 3a) and species diversity (Shannon’s index) as its uncertainty or variance (Fig. 3b). To aid these visualizations we have smoothed both indices to the 15m grid. Fig. 6b provides the resultant ‘sketchy’ map - combining both species metrics. Species richness levels remain the same for maps in Fig 3a and Fig. 6b, whereas high levels of species diversity coloured dark green in Fig. 3b correspond to areas of high ‘sketchiness’ in Fig. 6b. Thus plant species are both rich and diverse, in for example, sub-catchments 14 and 15.

## Visualization of local multivariate data structure uncertainty

Here, we apply a localised PCA (GWPCA, Harris et al. 2011) to all seven variables of the soils nutrients survey. We then map the local percentage of the total variance (PTV) accounted for by the first two components (PC1 and PC2) in Fig. 4a, where it appears that the seven soils variables tend to be relatively uncorrelated in the central sub-catchments of the NWFP, but relatively correlated elsewhere. From the usual PCA, the single global PTV value is 64.6% for the first two components, and in this instance, we map the absolute difference between this global value and each local PTV value in Fig. 4b. In Fig. 6c, we use that depicted in Fig. 4b, to reflect local structural uncertainty in our soils data with respect to how different the local PTV data is to that found globally. Thus data structure in sub-catchment 4 is most like that found, if we were to naively apply a non-spatial PCA to this data.

## Visualization of local regression coefficient uncertainty

For our final visualisation, we apply a bootstrap methodology (Harris et al. 2015) to test for non-stationarity in the coefficients of a GW regression (GWR) as an alternative to a multiple linear regression (MLR), the null hypothesis. Our response is plant species diversity and our predictor is the C/N ratio for soils. GWR reduces AIC by 22 units from the MLR fit, suggesting clear value in a localised relationship. Fig. 5a maps the local regression coefficients for the C/N ratio for soils, which compares to a global coefficient of 0.015. Fig 5b maps *p*-values that have been formatted such that low *p*-values indicate coefficients which are significantly different to that found globally. Fig. 6d displays the corresponding ‘sketchy’ map, where high levels of ‘sketchiness’ indicate areas of significant non-stationarity in the coefficient estimates.

# IV CONCLUDING COMMENTS

We have presented four case studies on the value of ‘sketchy’ maps to communicate uncertainty. Further ‘sketchy’ visualizations are possible when it comes to relaying uncertainty in the boundaries of spatial data, such as those of each NWFP sub-catchment, which we are currently working on.

## References

* Brunsdon (2016). Representing Uncertain Geographical Information with Algorithmic Map Caricatures. *Geophysical Research Abstracts Vol 18 EGU2016-2609*, EGU General Assembly 2016.

* Harris P, Brunsdon C, Charlton M (2011) Geographically weighted principal components analysis. *International Journal of Geographical Information Science* 25(10), 1717-1736.

* Harris P, Brunsdon C, Gollini I, Nakaya T, Charlton M (2015) Using Bootstrap Methods to Investigate Coefficient Non-stationarity in Regression Models: An Empirical Case Study. *Procedia Environmental Sciences* 27:112-115.
Hurlbert SH (1971). The Nonconcept of Species Diversity. *Ecology* 52(4), 577-586.

Wood J, Isenberg P, Isenberg T, Dykes J, Boukhelifa N, Slingsby A (2012). Sketchy rendering for information visualization. *IEEE Transactions on Visualization and Computer Graphics* 18(12), 2749-2758.

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

<figure><img src="image_22.png" /><figcaption>Figure 2: (a) UK predictions of C (plant nutrients) and (b) corresponding UK interpolation variances.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>Figure 3: (a) Predicted species richness and (b) predicted species diversity.</figcaption></figure>
<figure><img src="image_24.png" /><figcaption>Figure 4: (a) Local PTV data from a GWPCA and (b) absolute difference in PTV from local to global.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>Figure 5: (a) GWR coefficients for C/N ratio for soils and (b) corresponding 'p-values'.</figcaption></figure>
<figure><img src="image_26.png" /><figcaption>Figure 6: ‘Sketchy’ maps for (a) Fig. 2; (b) Fig. 3; (c) Fig. 4; and (d) Fig. 5.</figcaption></figure>
# Machine learning to deal with uncertainty in knowledge base for multivariate clustering applied to spatial analysis

Stephane Bourrelly<fnref n="1" />

Antonino Marvuglia<fnref n="2" />

Ian Vázquez-Rowe<fnref n="3" />

## Abstract

We present a Tailor Made Machine Learning (TMML) methodology combining different clustering algorithms, spatial statistical methods and cartographic tools. The methodology is currently being programmed in an R package, especially designed to handle multivariate spatial datasets. We highlight the strengths of unsupervised clustering for the management of environmental health phenomena, also pointing out several uncertainty sources affecting the results of our analysis. In particular, we acknowledge that the traditional hierarchical clustering is usually applied without performing dynamic reallocations, integrating spatial key-concepts or discussing the quality of outputs. Therefore we describe the foundations of the TMML methodology, which is applied to deal with these uncertainties, as well as with the variety of possible outputs. The R functions are applied to the spatial dataset included in the package so to illustrate the procedure to apply for identifying the most accurate clustering output, in the context of a sustainable agriculture example in Luxembourg.

## Keywords

Clustering, lattice, uncertainty, decision-making, agriculture.

# I INTRODUCTION

When dealing with multidimensional datasets, clustering algorithms are often used to produce a reduced knowledge-base to simplify the representation of phenomena. The aim of multivariate clustering applied to spatial contexts is to form groups as heterogeneous as possible, composed by spatial units as similar possible. When clustering performs on lattice; spatial typologies are created in order to classify the administrative areas in meaningful clusters and summarise several dimensions of phenomena. In this way knowledge can be extracted from the multidimensional complexity with the purpose to conceive local policies, e.g. to ensure the economic growth and prevent socio-ecological conflicts (Delgado and Romero, 2016). The Hierarchical Clustering Analysis (HCA) is the most applied unsupervised method by public stakeholders, e.g. for the sustainable management of the agricultural activities (Boyacioglu and Boyacioglu 2008).

Unsupervised classification methods provide objective representations, unlike supervised methods which constrain the results by a known response variable. Unfortunately, even though

<footnote n="1">¹University of Lyon 3 – UMR 7300 (ESPACE), France</footnote>

<footnote n="2">²Luxembourg Institute of Science and Technology (LIST), Luxembourg</footnote>

<footnote n="3">³Pontificia Universidad Católica del Perú, Peru</footnote>

<footnote>*Corresponding author: s.bourrelly@hotmail.fr</footnote>
HCA presents a few advantages, the results are not optimal. Therefore some dynamic reallocation methods have been conceived to overcome this drawback (Hennig and Liao, 2010). However Levine (2000) notes that they are scarcely used in social sciences as well as the quality criteria which are rarely presented to discuss the statistical significance and the meaning of typologies. Moreover Levine points out that unsupervised clustering is usually used without integrating the key concepts of the scope of application; e.g. when algorithms are applied without taking into account the interconnectivity of the spatial units. Obviously, these shortfalls contribute to increase results' uncertainty, which might give rise to controversial decisions.

In this paper we briefly describe the methodology of a so-called Tailor-Made Machine Learning (TMML), which combines several clustering and data-mining algorithms with spatial statistic operators as well as cartographic tools. The TMML provides a methodology for taking into account the interconnectivity of spatial units. We present some of the functions of the TMML R package currently under development, to compare and map the outputs of the implemented clustering. Functions are applied to the case study shapefile of the TMML package, representing the environmental concentrations of chemical substances generated by the application of agricultural fertilizers in the administrative areas in Luxembourg.

## II MATERIAL

The shapefile of the TMML package was derived from a cadastre spatial layer, representing the agricultural land use on the administrative areas (communes) in the Grand-duchy. This polygon layer describes the 108.328 georeferenced parcels contained in the cadastre records in 2009. These agricultural surfaces are differentiated from their main Agricultural Land Class (ALC) among the 23 official types, such as cereals, oilseed, legumes, etc. (Figure 1).

<figure><img src="image_27.png" /><figcaption>Figure 1: Luxembourg maps with 108.328 parcels named according to their main ALC.</figcaption></figure>

In order to know the agricultural harvests in 2009 we used the agricultural statistics of the national database<fnref n="1" />. However, national statistics differentiate the agricultural surfaces from 27

<footnote n="1">¹ http://www.statistiques public.lu/stat/ReportFolders/</footnote>
different categories of crops; which have been matched according to the ALC description of cadastral data.

Fertilizers dosages (and in some cases yields data) have been derived from experts' communications, as well as from KTBL (2006). Chemical emissions from fertilizers application have been estimated according to (Nemecek and Kägi 2007). Finally the ArcGIS software has been used for aggregating, at the scale of communes, all the environmental emissions and computing the spatial indicators. We note $x_k^j$, the concentration level of the chemical substance 'j' in the administrative area of the commune 'k' for an environmental compartment, i.e. air, ground water or surface water (Table 1).

<table><caption>Table 1: description of spatial indicators $x_k^j$.</caption><thead><tr><th rowspan="2">j<br>substance</th><th colspan="3">AIR</th><th colspan="2">GROUNDWATER</th><th colspan="2">SURFACE WATER</th></tr><tr><th>x_NH<sub>3</sub><br>Ammoniac</th><th>x_NO<sub>x</sub><br>Nitrogen oxide</th><th>x_N<sub>2</sub>O<br>Nitrous oxide</th><th>x_P<br>Phosphorus</th><th>x PO4<sup>3-</sup><br>Phosphate</th><th>x_NO<sub>3</sub><sup>-</sup><br>Nitrate</th><th>x PO4<sup>3-</sup><br>Phosphate</th></tr></thead><tbody><tr><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr></tbody></table>

The command

* shapeLux = data(ml.chemichal)

-loads polygon shapefile of the TMML package. It describes the 116 communes of Luxembourg in 2009 (names, surfaces, index) and provides the values of the $7 x_k^j$. The Figure 2 displays the spatial distributions of the concentrations of $NO_x$ in air, $NO_3^-$ in surface water and $PO4^{3-}$ in groundwater, expressed in kg per km² of administrative area.

<figure><img src="image_28.png" /><figcaption>Figure 2: maps of $x_k^{NO_x}$ (upper left), $x_k^{NO_3^-}$ (upper right) and $x_k^{PO_4^{3-}}$ (lower).</figcaption></figure>
### III METHOD

The HCA algorithm performs an unsupervised classification for creating objective typologies. At the first step it groups the two most similar individuals (spatial units) so to form the first cluster (group). Then it iteratively maximises the inter-clusters inertia to group the other spatial units into nested cluster. Nested clusters are organized as a hierarchical tree called dendrogram, which allows specifying a number of clusters C. However with hierarchical clustering methods the composition of each of the generic clusters $C_i$ is not optimal, due to the nested merging process. To overcome this drawback dynamic reallocation methods have been conceived in the literature for iteratively interchanging the cluster composition and minimising the intra-clusters inertia. The most popular algorithms are the Partitioning Around Medoids (PAM) and K-means (Kaufmann and Rousseeuw, 1990).

Since the model hypotheses are different they provide different typologies. In addition they not allow objectively choosing a number of clusters; that is why they are usually initialised with the HCA outputs. The interpretation of typologies is a controversial issue, as only their statistical significance can be discussed from several quality criteria such as the R², the C-INDEX, the rate of inertia or the average silhouette (Hennig and Liao, 2010).

Moreover these clustering methods assume that individuals are independent, but each spatial unit influences the other ones located in its neighbourhood. Here we propose a methodology to integrate the spatial connectivity in the clustering analysis.

Firstly, in the TMML guidelines a conceptual neighbourhood pattern among those defined by Cliff and Ord (1981) should be chosen (Figure 3). In this way a neighbourhood matrix *W* is created. *W* contains weights $w_{kl}$, where $w_{kl} = 1$ if the spatial unit ‘k’ is directly adjacent to the spatial unit ‘l’ and $w_{kl} = 0$ otherwise. In order to test the interdependency between the adjacent values $x_k^j$ and $x_l^j$ the Moran’s I statistic can be used (Gaetan and Guyon, 2010).

The tests can be replicated at several orders *h* of contiguity, i.e. for *h* > 1. In this case the neighbourhood matrix is termed *W*(*h*) and the strength of spatial dependencies can be read on a correlogram. The correlogram chart displays the values of *h* as a function of the Moran's statistics *I*(*h*) and their two-side confidence interval values, to denote the presence of a spatial autocorrelation (Gaetan and Guyon, 2010). In this way clustering algorithms can be applied to lagged values of $x_k^j$, noted $xs_k^j$ to take into account the spatial dependency concept.

<figure><img src="image_29.png" /><figcaption>Figure 3: conceptual neighbourhood patterns at the order h=1.</figcaption></figure>

### IV RESULTS

The function

* `ml.setC(x=X, update.C=4)`

-performs a HCA and provides the chart of inertia differences between-clusters in decreasing orders, so to help in setting the number of groups *C* (Figure.4).
The function

* `ml.clust( shape=shapeLux, x=X, C=4, maps=T, bar.charts=T, ordinal=T, labCi = c("WEAK","MIDDLE","HIGH","MAJOR"))`

-performs the HCA, K-means and PAM for merging the spatial units in *C* groups. It maps the distributions of class labels for the clustering outputs. With an ordinal typology the class labels 'labCi' are assigned according to the average centres of *C*ᵢ. When *bar.charts=T* bar charts are returned to summarise the characteristics of communes within the cluster *C*ᵢ. The horizontal bars give differences in number of standard deviations, between the means of *C*ᵢ and the overall average values for each chemical emission (Figure.4).

<figure><img src="image_30.png" /><figcaption>Figure 4: chart of inertia differences (upper left); maps of typologies for PAM (upper right) and HAC (lower left) with related bar chart (lower right).</figcaption></figure>

The function

* `ml pattern(shape=shapeLux, x=X, j.lab="x_NOx", type="QUEEN", map=T)
Moran.test=T, correlogram="Moran", h.max=7)`

-explores, for each spatial indicator *j.lab*, the interdependencies of values from a particular *type* of spatial pattern. The global autocorrelation of $x_k^j$ is forecasted through the *Moran test* and the lag order of spatial dependencies, such as $h \le hmax$, can be assessed from the *correlogram* (Figure.5).
<figure><img src="image_31.png" /><figcaption>Figure 4: map of “Queen” connectivity pattern at the first order and Morans’I correlogram.</figcaption></figure>

The function

• `ml.clust(..., ..., spatial patterns=rep("QUEEN", 7), h=c(1,1,1,2,2,3,5), mat.weight="W", quality=T, maps.xs=T)`

-performs the HCA, K-means and PAM on the spatially lagged values of indicators $xs_j^j$, according to the **spatial patterns**, j-order dependencies *h* and standardised neighbourhood matrix, when `mat.weight="W"` (Figure.5). When `quality=TRUE` the quality criteria of clustering, highlighted in the method section, are returned.

<figure><img src="image_32.png" /><figcaption>Figure 5: maps of $xs_k^{\text{Nox}}$ and the typologies for PAM and HAC, with related bar chart.</figcaption></figure>
# V DISCUSSION

The variety of spatial typologies returned from the different clustering algorithms applied to the original or lagged spatial indicators, highlights the strong uncertainty induced in decision-making when the quality of outputs is not questioned.

Even though the uncertainties might be explored through some quality criteria, from the analyst's perspective it is important knowing to what extent the dynamic reallocations and the integration of spatial patterns really improve the statistical significance of the results.

In this respect an important uncertainty source in unsupervised clustering methods lays into the fact that the best clustering is usually chosen through statistical criteria. Indeed, the model selection issue is currently a recognised challenge in the unsupervised classification field (Hennig and Liao, 2010).

Ideally, the best clustering representation should be identified through a trade-off between the statistical significance of typologies and their spatial correlations with a decision criterion, identified by stakeholders. In particular, we point out the importance of defining the most suitable spatial typology for reducing the impacts of chemical emissions on eco-systems and therefore the cumulative socio-ecological inequalities. In future work, these latter could finally be estimated in line with the index of Morello-Frosch et al. (2011).

# References

* Boyacioglu H., Boyacioglu H. (2008). Water pollution sources assessment by multivariate statistical methods in the Tahtali Basin, Turkey. Environ Geol, 54, 275-282.

* Cliff A., Ord J. (1981). The problem of spatial autocorrelation. In Scott A. London: Pion. 25-55.

* Delgado A., Romero I. (2016). Environmental conflict analysis using an integrated grey clustering and entropy-weight method: A case study of a mining project in Peru. Environmental Modelling & Software, 77, 108-121.

* Hennig C., Liao F. (2010). Comparing latent class and dissimilarity based. Department of Statistical Science, UCL, Department of Sociology. University of Illinois.

* Gaetan C., Guyon X. (2010). Spatial Statistics and Modeling. Berlin: Springer.

* Kaufman L., Rousseeuw P. (1990). Finding groups in data: an introduction to cluster analysis. New York: Wiley.

* KTBL - Kuratorium für Technik und Bauwesen in der Landwirtschaft (Ed.) (2006). Faustzahlen für die Landwirtschaft, Darmstadt, Germany (in German).

* Levine J. (2000). But What Have You Done For Us Lately. 29(1), 35-40

* Morello-Frosch R, Zuk M, Jerrett M, Shamasunder B, Kyle AD. (2011). Understanding the cumulative impacts of inequalities in environmental health: implications for policy. *Health Aff (Millwood)* 30, 879-87.

* Nemecek T., Kägi T. (2007). Life Cycle Inventories of Swiss and European Agricultural Production Systems. Final report ecoinvent V2.0 No. 15a. Agroscope Reckenholz-Taenikon Research Station ART, Swiss Centre for Life Cycle Inventories, Zurich and Dübendorf, CH, retrieved from: www.ecoinvent.ch.
# Generation of Plateau-Approximated Fuzzy Zones

Hazaël Jones<fnref n="1" />, Serge Guillaume<fnref n="2" />, Patrice Loisel<fnref n="3" />, Brigitte Charnomordic<fnref n="3" />, Bruno Tisseyre<fnref n="1" />

*Corresponding author: hazael.jones@supagro.fr

## Abstract

When modelling spatial real data in environmental context, the border between two zones is rarely a sharp edge. Often, the transition from one zone to another is a gradual process. Fuzzy zones intend to model this uncertainty in zones by allowing a spatial point to belong partially to different zones with different membership degrees. Generating fuzzy zones is not easy to perform and can be time-consuming and difficult to interpret. It is also difficult to ensure the convexity of the generated fuzzy zones.

Our idea is to generate plateau zones instead of fuzzy zones. Plateau zones are an approximation of fuzzy zones. The interest of using plateau zones lies in the following: easier spatial coherence and convexity of zones, preservation of uncertainties, level and number of plateaus can be chosen depending on the application, interpretation and easier operations between zones.

Our algorithm is based on quantiles of spatial data in order to produce some isocontours. According to the desired number of plateaus, it is possible to adjust the quantile values to find next plateaus for each resulting zones. The goal of this representation is to provide the user a simplification of the spatial representation and to preserve uncertainties in order to use them in the decision process.

## Keywords

Fuzzy Zones, Plateau Zones, Spatial Data, Zoning Algorithm, Spatial Coherence

## INTRODUCTION

When summarizing spatial data to consider zones (agronomy, geography, etc.), the border between two zones is rarely a sharp edge. In this context, there are no obstacles to detect but management zones to define. Often, the transition from one zone to another is a gradual process (e.g. between clay soil and sand soil, we can have a transition zone with more and more sand).

Fuzzy zones could be used to model this uncertainty/gradual change in zones by allowing a spatial point to belong to several zones with different membership degrees. Fuzzy zones keep the zone uncertainties until the final decision is taken, in order to be able to use this information to relax operational constraints.

Many methods exist to build zones based on segmentation (Pal and Pal (1993), Pedroso et al. (2010), Roudier et al. (2008)) but little work has been done on the generation of fuzzy zones (see Crane and Hall (1999), Philipp-Foliguet et al. (2009)). In Philipp-Foliguet et al., the mem-

<footnote n="1">¹Montpellier SupAgro, UMR ITAP, Montpellier, France</footnote>

<footnote n="2">²Irstea, UMR ITAP, Montpellier, France</footnote>

<footnote n="3">³INRA, UMR MISTEA, Montpellier, France</footnote>
bership degrees are assigned based on the distance to typical zones and attribute values, it is
then difficult to ensure convexity of the generated fuzzy zones.

Plateau zones (see Kanjilal et al. (2010)) are an approximation of fuzzy zones as shown in Figures 1 and 2. Each plateau corresponds to a level with a membership degree. Kanjilal et al. (2010) have studied the interest of plateau zones for computation, but there are no computation methods to generate these plateau zones. This article is a first step to provide a method for such generation.

<figure><img src="image_33.png" /><figcaption>Figure 1: An example of a fuzzy zone modelling an air-polluted area (a) and its representation as a plateau zone (b) source : Figures from article Kanjilal et al. (2010).</figcaption></figure>

<figure><img src="image_34.png" /><figcaption>Figure 2: A 2-dimensional illustration of a fuzzy zone and its approximation by plateau zones.</figcaption></figure>

Plateau zones present an interest for several reasons:

* Enhanced spatial coherence and convexity of membership degree for each zone,

* The level and the number of plateaus can be chosen depending on the application. Particularly in agriculture, it is not possible to apply a continuous treatment between one level and another.

* Interpretation: it is more efficient to use a map with plateau zones than continuous fuzzy zones in the decision making phase. Experts often reason using categories and not using a continuum of fuzzy values.

* Effective geometrical operations can be calculated between zones (difference, union, intersection, etc.) thanks to the simplification obtained by plateau (see Kanjilal et al. (2010)).

The goal of this representation is to provide the user a simplification of the spatial representation and to preserve uncertainties in order to use them within the decision process.
## THE ALGORITHM FOR BUILDING PLATEAU ZONES.

The algorithm is based on quantiles of spatial data. The use of quantiles on spatial data intends to produce some isocontours. An isocontour (level set) is a line with the same value all along its length. According to the desired number of plateaus, it is possible to slightly adjust the quantile to find the next plateaus for each resulting zone. Quantiles used on spatial data are an efficient way to ensure the imbrication of plateaus with a spatial coherence. If expert knowledge is available on the data (thresholds to divide the data, knowledge about data distribution, etc.), it is possible to use it to set the quantile values. When no *a priori* information is available, quantiles can be chosen in order to divide data into equal proportions. The choice of the best quantile is not the focus of this paper.

<figure><img src="image_35.png" /><figcaption>Figure 3: Initial zones based on a quantile.</figcaption></figure>

Figure 3 shows the initial zoning given by one quantile (0.4) on simulated mono-dimensional geo-referenced data (generated using a Gaussian random field). Four distinct initial zones are obtained, and our goal is to first obtain plateaus from them, and to generate approximated fuzzy zones from these plateaus. We will focus on Zone 2 to illustrate the algorithm.

Figures 4(a), 4(b), 4(c) and 4(d) shows how the shape of Zone 2 evolves when the quantile is slightly adjusted. When applied to Zone 2, the algorithm can provide a computation of the kernel, $\alpha$-cuts and support of this zone (see Figure 2 for an illustration of this concepts). Figure 4(a) shows the potential kernels of Zone 2.

Algorithm 1 is the proposition for building approximated fuzzy zones from plateau regions. For an initial zone and an initial quantile, our procedure *BUILDPLATEAUZONES* produces *nPlateau* zones associated to *nPlateau* quantile values. It is based on the following principle: as the initial zone (based on quantile 0.4) is selected to divide zones properly, it is granted a membership of 0.5. When we define zones thanks to quantiles inside this zone, the membership increases (see Algorithm 1, lines 6-7), and when we define them outside of this zone, it decreases (lines 8-9). The deeper (further) the zone is inside (outside) the initial zone, the more (less) it is characteristic of the zone. The algorithm allows a membership degree between 0 and 1 (line 7 between 0.5 and 1 - line 9 between 0 (0 excluded) and 0.5). A plateau with a membership degree of 0 is not relevant as it will not be part of the zone. Formula $0.5 - \frac{n}{nPlateau+1}$ (line 9) avoids this case.
<figure><img src="image_36.png" /><figcaption>Figure 4: Zoning for different quantile values</figcaption></figure>
**Algorithm 1** Algorithm for building plateau zones

1: **procedure** `BUILDPLATEAUZONES(initialZone, quantile, step, nPlateau, data)`
   ▷ `listPlateaus` is a list of zones with associated memberships degrees
2:   `listPlateaus` ← ∅
3:   `ADDZONEDEGREE(listPlateaus, initialZone, 0.5)`
   ▷ compute new zones
4:   `n` ← 1
5:   **while** `n` ≤ `nPlateau`/2 **do**
   ▷ compute new internal zones (memberships > 0.5)
6:     `newZoneInt` ← `COMPUTEZONEINT(initialZone, quantile - n * step, data)`
7:     `ADDZONEDEGREE(listPlateaus, newZoneInt, 0.5 + n/(nPlateau))`
   ▷ compute new external zones (memberships < 0.5)
8:     `newZoneExt` ← `COMPUTEZONEEXT(initialZone, quantile + n * step, data)`
9:     `ADDZONEDEGREE(listPlateaus, newZoneExt, 0.5 - n/(nPlateau+1))`
10:    `n` ← `n` + 1
11:
12:   **function** `COMPUTEZONEINT(initialZone, quantile, data)`
13:   `relevantNewZones` ← ∅
   ▷ return a list of zones based on `quantile`
14:   `newZones` = `COMPUTEZONES(quantile, data)`
15:   **for all** Z ∈ `newZones` **do**
16:     **if** Z ⊆ `initialZone` **then** `ADDZONE(relevantNewZones, Z)`
   ▷ return the largest zone that belongs to eligible zones
17:   **return** (`BIGGERZONE(relevantNewZones)`, `quantile`)
18:   **function** `COMPUTEZONEEXT(initialZone, quantile, data)`
19:   `relevantNewZones` ← ∅
   ▷ return a list of zones based on `quantile`
20:   `newZones` = `COMPUTEZONES(quantile, data)`
21:   **for all** Z ∈ `newZones` **do**
22:     **if** Z ⊇ `initialZone` **then** `ADDZONE(relevantNewZones, Z)`
   ▷ return the smallest zone that belongs to eligible zones
23:   **return** (`SMALLESTZONE(relevantNewZones)`, `quantile`)
24:   **function** `COMPUTEZONES(quantile, data)`
25:   `zonesFromQuantile` ← ∅
26:   `isoContours` ← `GENERATEISOCONTOURS(quantile, data)`
27:   `zonesFromQuantile` ← `CLOSEZONES(isoContours, data)`
28:   **return** (`zonesFromQuantile`)

When zones are generated thanks to a new quantile (functions `COMPUTEZONEINT` and `COMPUTEZONEEXT`), relevant new zones must be included in the initial zone for `COMPUTEZONEINT` (and must include the initial zone for `COMPUTEZONEEXT`). As a single quantile can generate multiple contour lines, if several zones are relevant, we choose the one with the largest area for `COMPUTEZONEINT` (smallest for `COMPUTEZONEEXT`). An illustration of this choice can be seen on Figure 4(a) where zone 3 will be chosen as it is bigger than zone 2. `COMPUTEZONES` returns the list of zones generated on a quantile isocontour basis, it builds isocontours based on data (line 26), and closes them respecting the field border (line 27). The initial zone is built following the same principle but initial zoning is edited by removing/expanding zones which
are too small.

Some functions are not specified in the algorithm to be more concise. `ADDZONEDEGREE` adds a zone and the associated membership degree to the final list of plateaus *listPlateaus*. `ADDZONE` adds a zone in the relevant zones. `BIGGERZONE` and `SMALLESTZONE` select the largest zone (respectively, the smallest) between a set of relevant zones. `GENERATEISOCONTOURS` produces isocontours from data values corresponding to a quantile and close them respecting the field border (`CLOSEZONES`).

The result of the computation on Zone 2 is given in Figure 5. The dark green line represents the border of the kernel (plateau of level 1) of Zone 2. 0.75 plateau is in green, 0.5 in black, 0.3 in blue, 0.17 in light blue. The same procedure is followed separately for each initial zone.

Thanks to the use of quantiles, the shape of Zone 2 plateaus is not just a dilatation-erosion of the initial zone because plateaus are directly based on data.

<figure><img src="image_37.png" /><figcaption>Figure 5: Final result with 5 plateaus for zone 2 (dark green (kernel), light blue (support)).</figcaption></figure>

The algorithm allows the computation of fuzzy regions thanks to a plateau approximation. The computation is efficient (2.15 seconds for 5 plateaus on Intel core i7, 8 cores 2.7GHz, data size = 968) and provides plateaus on spatial data. Plateau zones are relevant because they provide both a simplification of the spatial representation and also the preservation of uncertainties.

This preliminary work can be extended in the following directions: First, it will be interesting to analyse the sensitivity of parameters. The initial quantile choice and the quantile step will influence the results. The bigger the value of the step, the wider the plateaus, the larger the uncertainty of each zone. Furthermore, a study considering the number of plateaus should be done in order to decide how many plateaus are needed to give the most relevant approximation. Finally, we want to show the interest of this representation using real agronomical data.
## References

Crane S. E., Hall L. O. (1999). Learning to identify fuzzy regions in magnetic resonance images. In *Fuzzy Information Processing Society, 1999. NAFIPS. 18th International Conference of the North American*, pp. 352–356. IEEE.

Kanjilal V., Liu H., Schneider M. (2010). Plateau Regions : An Implementation Concept for Fuzzy Regions in Spatial Databases and GIS. In *Proceedings of the 13th International Conference on Information Processing and Management of Uncertainty in Knowledge-Based Systems*, pp. 624–633.

Pal N. R., Pal S. K. (1993). A review on image segmentation techniques. *Pattern recognition* 26(9), 1277–1294.

Pedroso M., Taylor J., Tisseyre B., Charnomordic B., Guillaume S. (2010). A segmentation algorithm for the delineation of agricultural management zones. *Computers and Electronics in Agriculture* 70(1), 199–208.

Philipp-Foliguet S., Gony J., Gosselin P. H. (2009). FReBIR: An image retrieval system based on fuzzy region matching. *Computer Vision and Image Understanding* 113, 693–707.

Roudier P., Tisseyre B., Poilvé H., Roger J.-M. (2008). Management zone delineation using a modified watershed algorithm. *Precision Agriculture* 9(5), 233–250.
# Modeling process chain of SPOT images for resources uncertainty to monitor change in forest cover

Aimé Richard Hajalalaina*<fnref n="1" />,<fnref n="2" />, Dominique Hérvé<fnref n="3" />, Eric Delaitre<fnref n="4" />, Thérèse Libourel<fnref n="4" />,<fnref n="5" />

*Corresponding author: arhajalalaina@yahoo.fr

**Abstract.** In this paper, process chain and knowledge-based models of SPOT satellite images are proposed to help scientists, working in the field of the environment and in particular of the forest, to solve uncertainty of spatial resources (data, process) to monitor the change in forest cover which usually results of deforestation. Indeed, deforestation mobilizes all research, various methods of satellite image processing on forest dynamics are proposed. The SPOT images available, at present, are voluminous and heterogeneous in terms of spatial, radiometric and temporal resolutions (sensors HRV, HRVIR, HRG, NAOMI). The use of these SPOT images may solve uncertainty of data. In addition, phenomena such as deforestation require the analysis of time series of satellite images and the development of automated and reusable processing chains for monitoring change of forest cover. We propose to formalize these processing chains from modeling an abstract and concrete models based on existing standards in terms of interoperability (International Standard Organisation ISO and OGC Open Geospatial Consortium). The use of these standards solves uncertainty of process. These processing chains modelled will be capitalized, and diffusible in operational environments. Our modeling approach uses work-context concepts. These concepts need organization of human resources, data, and process in order to establish a knowledge-based connecting the two latter. This knowledge-based will be used to solve uncertainty of SPOT images resources for monitoring change in forest cover.

**Keywords:** forest cover, resources uncertainty, modeling, knowledge base, processing chain, satellite image, SPOT, work-context.

## I INTRODUCTION

Our concerns are on the application of computer science in the field of environment remote sensing. Specialists in remote sensing produce hypotheses which they validate from experimental protocols. Now, remote sensing and computer experts work together to automatize these protocols in processing chains. This automatization poses many problems due, on the one hand, to the volumes of satellite images coming from different sensors and, in the other hand, to the proliferation of more or less complex image processing methods required by environment remote sensing.

In this context, it is necessary to put in place systems which help to store and manage important image-streams as well as their processing, by taking into account their various origins (different sensors, different spatial, spectral and temporal resolutions). The restitution and the exchange of these pieces of information are a real challenge in terms of

<footnote n="1">¹Centre Universitaire de Formation Professionnalisante, Université de Fianarantsoa, Madagascar</footnote>

<footnote n="2">²Laboratoire de Recherche Appliquée Multidisciplinaire, Université de Fianarantsoa, Madagascar</footnote>

<footnote n="3">³Institut de Recherche pour le Développement (IRD UMR 220), BP 64501, Montpellier, France</footnote>

<footnote n="4">⁴Espace DEV, 500 rue JF Breton, Montpellier, France</footnote>

<footnote n="5">⁵LIRMM, 161 rue Ada Université de Montpellier II, Montpellier, France</footnote>
interoperability. This later has an advantage in reducing the uncertainty of image-processing methods through the capitalization and the mutualization of experiments between remote sensing specialists.

A lot of research-work, which represent high accuracy of classification, has been published on the study of forest cover using SPOT images, namely Hajalalaina et al. (2013), Souza et al. (2005), Achard et al. (2002), Kimes et al. (1999), which has contributes to improving the knowledge on forest dynamics. But a data and processing formalization is necessary to ensure much larger spread of the research results to the scientific community. It helps to resolve the problem of uncertainty in the SPOT-image processing methods on the study of forest cover.

The formalization of data and processing is realized in a work-context acquired by the MDWeb platform put into practice by Desconnets et al. (2007). This platform proposes a view of human resource organization. It allows the plateform administrator to reference the future users from predefined roles and rights, and to reference the ressources (data/processings) within a metadata-base.

We propose to use this of this formalization to coordinate the collection of SPOT images and their processing, necessary for monitoring forest-cover. The objective is to capitalize, harmonize and spread the resources which allow a better understanding of forest dynamics at different spatial and temporal resolutions. It is part if several achievable actions within different period of time. The first one is the formalization of SPOT image processing to ensure their sharing, their re-use and their interoperability. The second one, which is a medium-term objective, is setting-up of a platform for the sharing a mutualization of experiments on SPOT image processing methods on the study of forest cover.

First, we have collected the current norms and formalization of spatial data and processing, then proposed work-context models for the formalization of SPOT image processing chains for monitoring of forest-cover.

## II RESOURCE FORMALIZATION

The formalization corresponds with the description of the resources used by respecting the norms and standards which are in force. They are in the size of metadata. Various formalizations respecting the norms in force deal with syntactic and semantic aspects for the metadata. This definition and the application of these norms help to resolve the heterogeneity of the resource used. Dublin Core in BNF (2008) and ISO 19115 (2003) are the two norms which are valid and most used for the formalization of georeferenced data, like SPOT image, while ISO 19119 (2005) allows us to describe their processing. The formalization of processing chains use notion of work-context. The work-context is closely linked to the field of expertise. It can then be envisaged that it is built, through consensus, by experts in the area. Every scientist (the final user) will have at their disposal this context which they can develop or enrich in their turn. The construction of any work-context corresponds with the progressive on gradual organization of useful data and processing references, following three steps: human resource organization, data organization, and processing organization according to Libourel et al. (2010). Moreover even if this approach has not been standardized yet, it uses a relatively simple graphic symbolic language, called Simple Workflow Model (SWM) suggested by Lin et al. (2008), which allows the user scientists to handle easy appropriation concepts in a simple language.

The work-context of SPOT images in field of monitoring forest-cover is made up of three models bellow: human resource organization (figure 1) which manages the description of the platform users as well as that of their different roles and access rights, data organization (figure
2) which manages the description of SPOT image category according to sensors, processing organization (figure 2), which manages the description of SPOT images processing categories, work-context knowledge-base (figure 3) and the abstract model of processing chain of SPOT images to monitor change in forest cover (figure 4). This abstract model can be implemented by designing concrete model by using Orfeo Toolbox (2015) image processing.

<figure><img src="image_38.png" /><figcaption>Figure 1: Human resource organization model</figcaption></figure>

<figure><img src="image_39.png" /><figcaption>Figure 2: SPOT images organization and their processing models</figcaption></figure>

<figure><img src="image_40.png" /><figcaption>Figure 3: Extract of Work-context knowledge-base of SPOT images model</figcaption></figure>
<figure><img src="image_41.png" /><figcaption>Figure 4: Abstract model of multi-resolution and multi-temporal SPOT images processing chain to monitoring change in forest cover</figcaption></figure>

### III ASSESSMENT OF PROCESSING CHAINS MODELS OF SPOT IMAGES

The manual processing chains of SPOT images will always remain a tough intellectual process. As a matter of fact, this type of chain requires the help of a remote sensing expert all along the realization process. In addition, new images entering the chain require the re-start of
the process from the beginning till the end. This makes their re-use difficult for users who are not specialists in the field and it does not make it possible to capitalize the experts' knowledge to assist the occasional users in the valorization of new images. We can say that the manual processing chains are not adapted to the SPOT images actualizing critical phenomena like the monitoring of forest cover.

Faced with this advantages of the manual processing chains, it interest of our proposals on formalizing processing chains lies in the automatization of different processing on SPOT satellite images.

**Processing chain re-use:** Once processing chain is formalized, its execution (implementation) is carried out automatically by taking into account new images. This is what makes the chain re-usable and also makes several repetitions of the processing possible without expert assistance for using the new images to carry out their analyses corresponding to their needs. Those properties are actually adapted to the monitoring of the forest cover, in view of changing and evolving character of SPOT images. It presents a major interest for capitalizing the experiments on SPOT images processing methods in order to overcome the problems of uncertainty in the processing of forest cover evolution.

**Interoperability of processing chains:** The formalized chain can therefore, a priori, be implemented in a distributed environment of Grid or Cloud type. Thus, it would then now be interesting to display these chains in a distributed manner in order to profit from this possible interoperability. As a matter of fact, in the contexts of forest cover monitoring, the distributed dimension will bring about an important advantage to this approach for it allow to divide up the experiments and spread on the whole set of the system to put in place, the contribution of separate elements to the knowledge-bases of work-context. In addition, normalized exchange protocols allow distributed requests which can be interesting for simultaneous research of SPOT image processing methods used in monitoring forest cover, which makes it possible to fills up the gaps on the uncertainty of SPOT image processing methods.

**Cataloging processing chains:** Making catalogs is the first step that will lead to an exchange or a spread of resources. Often the catalog is the origin of the initiatives for the rapprochement between the actors who wish to share their SPOT image processing chains for the monitoring of forest-cover. We would like to stress that it does not answer a move to share the processing chains but the one to share the metadata. We take profit from this possibility that the processing chains modeled and formalized in this way can be included into a catalog.

**Processing chain exchange:** We can take profit from the processing chain exchange by taking into account the different human resources which are proposed for monitoring change in forest cover. As a matter of fact, in contract to the spread which can be made through a downloading website, the exchange of SPOT image processing chains generally requires an interaction between the actors involved. This encounter is often initiated for legal reasons in order to draw up an exchange convention that respects the objectives of each party, or each of the parties. The notion of exchange is sometimes mixed with that of partnership, much broader. A partnership is a process of strategy cooperation between at least to actors and whose objective is to achieve a goal through the common use of material, intellectual, human, financial means. It often has a political dimension. The exchange of processing chains is interesting in setting up a process concerning the physical realization of one of the aspects of partnerships between the different actors. We stress that the notion of exchange implies a bilateral relation (the producer versus the receiver) which focuses on a flow while the notion of sharing also allows considering all the related aspects, in particular the appropriation and sharing of the knowledge which can result the exchange.
**Spread of processing chains:** More modest on the level of interactions, the objectives on the spread are not a challenge for the capacities of appropriation of processing chains by the actors with who the producer will have little or even no interaction. It can be pointed out here that there is a marked opposition with the co-production objective which, instead of relying on work groups, prefers to spread the produced processing chains without any discussion. It's the logic of the action of the researchers who wish to make their processing chains accessible to a great number of people.

**Mutualization of processing chains:** The interest in the formalization of chain-processing is to encourage the mutualization of experience experiments between the researchers and is justified by various reasons. The reason is first of all to avoid useless duplication of efforts in the constitution of chain processing, avoiding at the same time redundancy of the result fragmentation which comes out. Considering the compartmentalization related to the multiplication of autonomous and non-coordinated chain-processing, the formalization of these chains represents an adequate solution to these defects.

Finally, our proposals on the formalization of processing chains constitute necessary conditions for the mutualization of experiments for monitoring forest cover. In fact, the mutualization of data implies formalizing the production and spread supervision modalities, according to an interoperable mode of a piece of information ready for use and whose quality is checked, for reason of common usage. By reinforcing the cooperation between producer and users of processing chains, the mutualization also favours the exchange of experience and good practice through networking.

**Processing-chain sharing:** At last, it can be said that our proposals on the processing-chain formalization facilitate the sharing of experiments to monitoring change in forest cover. As a matter of fact, the data sharing aims at providing the community of experiment-users with coherent formats and structures, which will enable them to fulfill their missions better and produce their experimentations according to common and pre-defined requirements. Conversely, this interoperability and the improvement of processing-chain availability can contribute to reinforcing the relationships between the producer and user organizations and therefore can help to fill up some gaps on the use of SPOT satellite images for the monitoring of forest cover.

**Experiment-sharing platform:** Starting from the advantages offered by over proposals on the modelling and formalization of SPOT image processing chains for the monitoring of forest cover we are going to present other advantages which come out. We can state that the formalization of processing chains can be valued by the setting-up of a platform for sharing, mutualization re-use and the spread of experiments on the valuation of SPOT images for the monitoring of forest-cover.

The knowledge on the SPOT image processing methods for the monitoring of forest-cover are modelled then formalized in the form of processing-chains. The knowledge acquired on SPOT image-processing has to be made available to the actors. This capitalization is conceived so that everyone's experiment does not remain continued to individual level, but serves the collective in a move for knowledge sharing allowing to reduce the uncertainty in concerned field. The preservation and transmission of experiments and knowledge through formalized processing-chains facilities the implementation of new experiment protocols using SPOT images. More other, the capitalization and valuation of experiments using processing-chains is part of knowledge management. It means that the starting point is strong hypothesis that all experience or knowledge can by organized, referenced, enriched on supports which are adapted and exchanged as knowledge that other people can make their own. This puts the stress on the re-use, interoperability and sharing, the re-use, the spread of processing-chains to
valve the new SPOT images available as well as the new innovating methods for the
monitoring of forest-cover.

The interests in formalizing processing-chains can be valued by the proposal on setting up a platform while allows the sharing, the mutualisation the spread, the naming of SPOT image processing-chains through catalogs for the monitoring of forest-cover. In addition, the existing platforms at the world level propose functionality for knowledge management following this view; we have proposed platforms for sharing SPOT images to deal with the uncertainty due to the lack of image of this type.

Knowledge management consists in grousing, sharing, and updating knowledge. It requires not only setting up mechanisms and procedures for uniting, organizing, presenting and spreading SPOT image processing-chains which are modelled and partners, but also proceeding to the assessments of these operation-results.

We have demonstrated how to produce SPOT image processing-chains which are reusable, sharable, spreadable, interoperable, and mutualizable. The step which can follow is implementation of a platform for sharing experiments. In view the interest the valorisation of our proposals, we can draw up a deduction that the solutions we found represent a help in the formalization of SPOT image processing-chains to capitalize the knowledge of the monitoring of the evolution of forest-cover.

IV CONCLUSION AND PERPECTIVES

This article shows the possibility for sharing and mutualizing the experiments or existing SPOT image-processing to deal with the uncertainly of the monitoring of forest-cover using the technique of remote sensing. The proposal on resource modelling and formalization for their sharing is one of the efficient techniques for resolve the problems of heterogeneity of SPOT image-processing. In fact the experiments or the processing methods of these images used for the monitoring of forest-cover will be re-usable for non-remote-sensing-expert-users. This decreases the uncertainly of the use of SPOT image processing in the forest area.

Our wok is limited to modelling, formalization and design of a resource-knowledge base
related to the work-context. This latter has its advantage because if can be enriched by new
recently-produced resources. It helps integrate new satellite images and recent processing
methods for analysing forest cover.

Next, the perspective of this work is the integration of resources formalised in that way. This formalization makes the resources interoperable, which makes it possible to integrate the resources easily into the MDWeb tools for the sharing and mutualization and WPS (Web Processing Services), for the execution. For the execution of processing-chains via the web, we propose the web WPS Service server OGC standard (2012) since 2005.The specification of the WPS is in the form of a generic interface allowing describe and carry out or execute image processing chains of satellite images, according to Machet et al. (2008) and Eberle and Strobl (2012), which we produced in this work. The WPS is based on the http protocol and XML language. The processing-chains presented in this work will be converted into XML files, and will be integrated into WPS for execution.

References

Achard F., Eva H.D., Stibig H.-J., Mayaux P., Gallego J., Richards T., Malingreau, J.-P. (2002). Determination of deforestation rates of the world's humid tropical forests. Science 297:999-1002.

BnF, (2008). Guide d'utilisation du Dublin Core (DC) à la BnF : Dublin Core simple et Dublin Core qualifié, avec indications pour utiliser le profil d'application de TEL, Version 2.0. Bibliothèque nationale de France /Direction des Services et des Réseaux / Département de l'Information bibliographique et numérique, France.
Desconnets J-C., Libourel T., Clerc S., Granouillac B., (2007). Cataloguing for distribution of environnemental ressources. *10th AGILE, International Conference on Geographic Information Science*, Aalborg University, Danmark.

Eberle J. and Strobl C., (2012). WEB-Based Geoprocessing and Workflow Creation for Generating and Providing Remote Sensing Products. Geomatica, Vol.66(1), pp.13-26. Canadian Institute of Geomatics.

Hajalalaina A. R., Grizonnet M., Delaître E., Rakotondraompiana S., Hervé D., 2013. Discrimination des zones humides en forêt malgache, proposition d'une méthodologie multirésolution et multisource utilisant ORFEO ToolBox. *Revue Française de Photogrammétrie et de Télédétection*, n° 201, pp. 37-48.

ISO19115, (2003). Geographic Information Metadata, ISO 19115. International Organization for Standardization (ISO), Genève, Suisse.

ISO19119, (2005). Geographic Information Service, ISO 19119. International Organization for Standardization (ISO), Genève, Suisse.

Kimes D.S., Nelson R.F., Salas W.A. and Skole D.L., (1999). Mapping secondary tropical forest and forest age from SPOT HRV data. *in International Journal of Remote Sensing*, 20:3625-3640.

Libourel T., Lin Y., Mougenot I., Pierkot C., (2010). A platform dedicated to share and mutualize environmental applications. In J. Filipe, J. Cordeiro J. (eds.), *in Proceedings of the 12th international conference on enterprise systems. ICEIS, International Conference on Enterprise Systems, 12.*, Madere Funchal, 8-12 juin 2010. Setubal : SciTePress, p. 50-57.

Lin Y., Mougenot I., Libourel T. (2008). Un nouveau langage de workflow pour les sciences expérimentales. *In INFORSID'08 : Atelier ERTSI Evolution, Réutilisation et Traçabilité des Systèmes d'Information*, Fontainebleau, France.

Machet E., Kamhi M., Jacquin M., Le Page M., Dejoux J.-F. and Dedieu G., (2008). Web Processing Service pour le traitement des images satellites. *CESBIO, CNES, Toulouse, France*.

OGC. (2012). Web Processing Service 2.0 Standard Working Group. http://www.opengeospatial.org/projects/groups/wps2.0swg. Accessed: 2016-04-29

*Orfeo ToolBox*, (2015). The ORFEO Tool Box Software Guide Updated for OTB-5.2.1. https://www.orfeo-toolbox.org/packages/OTBSoftwareGuide.pdf. Accessed: 2016-04-29. Centre d'Etudes Spatiales (CNES), Toulouse, France

Souza C., Firestone L., Silva L.M.,., Roberts D., (2003). Mapping forest degradation in the Eastern Amazon from SPOT 4 through spectral mixture models. *in Remote Sensing of Environment*, 87:494-506.
<figure><img src="image_42.png" /></figure>
<figure><img src="image_43.png" /></figure>

SPATIAL
ACCURACY
2016

5-8 JULY 2016

MONTPELLIER, FRANCE

# Spatial accuracy quantification in mapping
<figure><img src="image_44.png" /></figure>
# Uncertainty quantification of interpolated maps derived from observations with different accuracy levels

Gerard B.M. Heuvelink*¹², Dick Brus<fnref n="2" />, Tom Hengl<fnref n="1" />, Bas Kempen<fnref n="1" />, Johan G.B. Leenaars<fnref n="1" /> and Maria Ruiperez-Gonzalez<fnref n="1" />

## Abstract

Most practical applications of spatial interpolation ignore that some measurements may be more accurate than others. As a result all measurements are treated equally important, while it is intuitively clear that more accurate measurements should carry more weight than less accurate measurements. Geostatistics provides the tools to perform spatial interpolation using measurements with different accuracy levels. In this short paper we use these tools to explore the sensitivity of interpolated maps to differences in measurement accuracy for a case study on mapping topsoil clay content in Namibia using kriging with external drift (KED). We also compare the kriging variance maps and show how incorporation of different measurement accuracy levels influences estimation of the KED model parameters.

## Keywords

Africa, geostatistics, interpolation, kriging, measurement error, regression, soil

Spatial interpolation errors are an important source of uncertainty in many spatial modelling applications and analyses. Geostatistics provides the tools to quantify interpolation error through the so-called kriging variance or by using spatial stochastic simulation, but in standard kriging the measurement error of individual observations is rarely addressed explicitly. It is usually represented as a component of the nugget variance of the semivariogram, but this implicitly assumes that all measurements are unbiased and have the same random measurement error variance. In reality, different measurement precisions and accuracies may occur because the data used may be a merge of field estimates and laboratory measurements, may be measured using different instruments and laboratory methods, or may be derived using the same methods but in different laboratories. Often data are also measured indirectly through proxies, such as when soil properties are estimated from soil spectroscopy signals that are converted to soil property values using statistical methods such as Partial Least Squares Regression (Brown et al. 2006, Leone et al. 2012). In recent years observations are also increasingly generated through crowd-sourcing and volunteered geographic information initiatives, which may suffer from large measurement errors (Goodchild and Li 2012). These initiatives can yield large volumes of data at cheap or zero cost, but their accuracy will usually be less than that of institutional data. In this work we extend kriging with external drift to the case in which each individual observation can have a different measurement error variance. As

<footnote n="1">¹ISRIC - World Soil Information, The Netherlands</footnote>

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

<footnote>*Corresponding author: gerard.heuvelink@wur.nl</footnote>
a result observations with small measurement error variance carry more weight than
observations with large measurement error variance, both in regression modelling and kriging.

The methodology builds on well-known approaches in geostatistics that go back as far as
Delhomme (1978) and is also presented in text books such as Chilès and Delfiner (1999,
Section 3.7.1). It boils down to modification of the kriging matrix by adding the measurement
error variances to the diagonal elements of the covariance matrix. In case of correlated errors,
the off-diagonal elements will also be affected. Solving the kriging system using the modified
kriging matrix automatically decreases the kriging weights of observations with larger
measurement error variances. Also, the influence of measurements on estimation of the trend
coefficients is reduced when the measurement error variance is larger. While this is all well
known, it is rarely applied in practice. This is unfortunate, because differences in measurement
errors may have a large impact on resulting maps and hence the prediction accuracy can be
markedly improved if these differences were taken into account. Differences in measurement
error variances are typically also not included in estimation of the semivariogram and trend
coefficients. In this presentation we show that measurement error variance can fairly easily be
included in parameter estimation (both for estimation of the variogram parameters and
regression coefficients) by taking a maximum likelihood estimation approach. Further to that,
we account for systematic measurement errors by representing the (unknown) systematic error
as a zero-mean random variable that is equal for all observations from the same source.
Uncertainty about the variogram parameters can be incorporated by taking a Markov Chain
Monte Carlo approach.

The statistical methodology is largely known and fairly straightforward but requires
adaptations of existing software implementations. We implemented the methodology as R
scripts, which in future will be extended to scalable R functions. We use a digital soil mapping
application using topsoil texture data from the Africa Soil Profiles database (Leenaars 2013)
and the LandPKS project (Herrick et al. 2013) to map topsoil clay content for Namibia. We
compare prediction maps and prediction error variance maps with those obtained when
measurement error is ignored. Results show marked differences and indicate that measurement
errors should not be ignored, particularly when there are large differences in accuracy levels
between observations within the conditioning dataset. We also explore the sensitivity of
mapping results for different degrees of spatial autocorrelation of measurement errors.

One important reason that the methodology has not often been applied in practice is that it
requires that the measurement error variances of all observations are known. In reality, this is
seldomly the case because data come from many sources and their accuracies are rarely
recorded. In our case study we used the texture triangle and expert judgement and expert
elicitation (O'Hagan et al., 2006) to quantify the measurement error variances, but these are
no substitute for real values and hence it is important that point data used for spatial
interpolation are routinely accompanied by measures of their accuracy.

Extension from linear multiple regression and regression kriging to non-linear machine-
learning regression methods, such as artificial neural networks, support vector machines and
random forests, is less obvious. One approach might be to duplicate more accurate
observations or assign weights to observations depending on their measurement error variance,
but this can only partly solve the problem and has a large ad hoc character. More satisfactory solutions should be found, but these should be sufficiently generic and work for the entire family of machine-learning methods, because in many practical machine-learning applications blends of multiple algorithms are used to optimise performance. This is important too, because in recent decades spatial interpolation makes use more and more of explanatory information contained in covariates and machine-learning algorithms are increasingly popular because they are more flexible than linear methods and usually produce more accurate predictions (e.g. Hengl et al, 2015).

## References

*   Brown, D.J., Shepherd, K.D., Walsh, M.G., Mays, M.D., Reinsch, T.G. (2006). Global soil characterization with VNIR diffuse reflectance spectroscopy. *Geoderma* 132(3-4), 273-290.

*   Chilès, J.-P., Delfiner, P. (1999). *Geostatistics. Modeling Spatial Uncertainty*. New York: Wiley.

*   Delhomme, J.P. (1978). Kriging in the hydrosciences. *Advances in Water Resources* 1(5), 251-266.

*   Goodchild, M.F., Li, L. (2012). Assuring the quality of volunteered geographic information. *Spatial Statistics* 1, 110-120.

*   Hengl, T., Heuvelink, G.B.M., Kempen, B., Leenaars, J.G.B., Walsh, M.G., Shepherd, K.D., Sila, A., MacMillan, R.A., Mendes de Jesus, J., Tamene, L., Tondoh, J.E. (2015). Mapping soil properties of Africa at 250 m resolution: random forests significantly improve current predictions. *PLoS ONE* 10(6): e0125814.

*   Herrick, J.E., Urama, K.C., Karl, J.W., Boos, J., Johnson, M.V.V., Shepherd, K.D., Hempel, J., Bestelmeyer, B.T., Davies, J., Guerra, J.L., Kosnik, C., Kimiti, D.W., Ekai, A.L., Muller, K., Norfleet, L., Ozor, N., Reinsch, T., Sarukhan, J., West, L.T. (2013). The global Land-Potential Knowledge System (LandPKS): Supporting evidence-based, site-specific land use and management through cloud computing, mobile applications, and crowdsourcing. *Journal of Soil and Water Conservation* 68(1), 5A-12A.

*   Leenaars, J.G.B. (2013). *Africa Soil Profiles Database, Version 1.1. A compilation of georeferenced and standardised legacy soil profile data for Sub-Saharan Africa (with dataset)*. ISRIC Report 2013/03. Africa Soil Information Service (AfSIS) project. ISRIC - World Soil Information, Wageningen, the Netherlands. 160 pp.

*   Leone, A.P., Viscarra-Rossel, R.A., Amenta, P., Buondonno, A. (2012). Prediction of soil properties with PLSR and vis-NIR spectroscopy: application to Mediterranean soils from Southern Italy. *Current Analytical Chemistry* 8(2), 283-299.

*   O'Hagan, A., Buck, C., Daneshkhah, A., Eiser, J., Garthwaite, P., Jenkinson, D., Oakley, J., Rakow, T. (2006). *Uncertain Judgements: Eliciting Experts' Probabilities*. Chichester: Wiley.
# Survey designs which maximize efficiency gains in ALS-based forestry plot imputation

Gavin Melville<fnref n="1" />, Christine Stone<fnref n="2" />, Jan Rombouts<fnref n="3" />

## ABSTRACT

The use of airborne laser scanner (ALS) data to estimate forest resource inventory variables is now becoming widespread in Australia (Rombouts et al., 2010; Stone et al., 2011). ALS data is combined with survey plot data to construct model-based estimates of timber volume. In particular, volumes of timber products, sourced from plantations of *Pinus radiata* have successfully been estimated using nearest neighbor methods. One of the challenges in this approach is to construct samples which capture the full efficiency gains which are achievable using the multiplicity of variables which can be derived from the ALS data. Some of the sampling design approaches that have been investigated include random sampling, grid sampling, stratification, systematic selection and balanced sampling. Estimates have been examined from both a design-based and model-based perspective (Melville et al., 2015). This talk will present results based on several approaches including a novel method specifically designed for imputation which optimizes the survey design by using the distance properties of the sample in the space defined by the auxiliary variables.

## I INTRODUCTION

The use of remote sensing to measure key inventory metrics over large areas of forest has become widespread in recent years (e.g. Breidenbach et al., 2010; Latifi et al., 2010; McRoberts, 2012; Hudak et al., 2014; Rombouts et al., 2014; Dash et al., 2015). One of the remote sensing techniques currently being employed in Australian softwood plantations is LiDAR, also referred to as airborne laser scanning (ALS). It has been established that the use of LiDAR data in conjunction with ground measurement leads to efficient prediction of commercially important attributes including timber volume, basal area and stems per hectare (Rombouts et al., 2010).

LiDAR and/or other auxiliary data have been used as covariates in a range of prediction methods including imputation, linear regression and machine learning methods such as random forest. With imputation, for example, forest plots which have been measured on the ground are linked to non-measured plots according to their similarity in the covariate space defined by the auxiliary variables. Specific information relating to key attributes, such as timber volume, is assigned to non-measured plots from the most closely related measured plots as defined by the similarity structure.

Measurement of field plots is a time consuming and costly process requiring specialized inventory crews who often need to access remote and challenging terrain. In order to be effective the field plots must cover the full range of variability in the key attributes. Therefore, appropriate

<footnote n="1">¹Trangie Agricultural Research Centre, Mitchell Hwy, Trangie 2823, Australia</footnote>

<footnote n="2">²NSW Department of Primary Industries Forest Research, Level 12, 10 Valentine Ave, Parramatta 2124, Australia</footnote>

<footnote n="3">³One FortyOne Plantations, 152 Jubilee Hwy E., Mt Gambier 5290 Australia</footnote>

<footnote>*Corresponding author: Gavin.Melville@dpi.nsw.gov.au</footnote>
methods of selecting the field plots are essential in constructing imputation estimates which are efficient, robust and economical.

There are a variety of sampling designs which can lead to a good spread in the attributes of interest. These include methods which are focused on structural attributes such as stratification, and spatially systematic designs such as grid sampling. Methods have also been developed which utilise the auxiliary data explicitly including stratification based on LiDAR variables (Hawbaker et al., 2009), multivariate methods such as sensor-directed response surfaces (SRDS - Lesch, 2005), and methods employing geographical coordinates such as generalized random tessellation stratified sampling (GRTS - Stevens and Olsen, 2004). More recently methods have become available which construct a sample which is simultaneously balanced with respect to multiple design variables and these have also been applied to forest inventory (Grafström et al., 2014). The present study investigates a new type of sample, termed nearest centroid (NC - Melville and Stone, 2016), which uses a multivariate clustering algorithm to select a sample of field plots.

In this paper the NC sampling method is presented, and data from a *Pinus radiata* plantation in South Australia are used to illustrate the technique. Comparisons are provided with other types of plot selection strategies.

## II MATERIALS AND METHODS

### 2.1 Inventory approach

We define the “area of interest” (AOI) as the finite region over which predictions will be made together with (if separate) the finite region over which ground-based measurements will be made. The AOI is tessellated into a set of non-overlapping contiguous pixels called “virtual” plots which are used to define the population. Virtual plots are constructed to have similar dimensions to the ground-based plots and each virtual plot is associated with a set of metrics which is calculated from the LiDAR data and used as auxiliary information. Virtual plots are further characterised as “reference” plots - the set of plots which are selected for ground-based measurement, “candidate” plots - the set of plots from which the reference plots are chosen, and “target” plots - the set of plots for which predictions are made. Therefore the approach which is used in an actual forest inventory can be separated into the following steps:-

1. Define the area of interest

2. Define the area containing the candidate plots

3. Construct a population of virtual plots in areas (1) and (2)

4. Select a sample of reference plots from area (2)

### 2.2 Imputation

The k-nearest neighbour approach involves calculating the distance or similarity in the auxiliary space between measured reference plots and target plots in order to determine which reference plots are most similar each target plot. The auxiliary variables are chosen because of their ability to predict the variable of interest which, in commercial forests, is typically timber volume. For every target plot, the k-nearest neighbour imputation estimate is the weighted mean of the variable of interest from the k most-similar reference plots, calculated as

$$ \tilde{y}_i = \frac{\sum_{j=1}^{k} w_{ij} y_j^i}{\sum_{j=1}^{k} w_{ij}} $$
where $\tilde{y}_i$ is the prediction for target plot $i$, $y_j^i$ is the $j$th nearest reference plot to plot $i$, and $w_{ij}$ is the weight assigned to plot $j$ (McRoberts et al., 2007). In this study $w_{ij} = 1$ was employed throughout.

## 2.3 NC sample

The NC sample operates by partitioning the target population into clusters of plots having minimum sums of squares of distances, from plots to the cluster centroids. The number of clusters is chosen to be equal to the required sample size, *n*, and the clusters themselves are not required to be spatially contiguous. The clustering procedure which was employed in this study is termed *k*-means clustering (Hartigan and Wong, 1979) and was performed using the standardized auxiliary variables and a Euclidean distance metric. *K*-means clustering is normally used for multivariate analysis of complex datasets. After calculating the cluster centroids one then finds the plots in the candidate set which are closest to the centroids in the auxiliary space and these become the reference plots.

An example is shown in Figure 1 which illustrates the selection of 8 reference plots by forming the target set into 8 clusters in the space defined by the auxiliary LiDAR variables p1m (proportion of heights greater than 1m) and mqh (mean quadratic height). The plots closest to the cluster centroids (centroids shown as black stars) are then selected as the reference plots. The term “nearest centroid” derives from the fact that the reference plots are the nearest plots to the target plot centroids. During imputation these plots become the nearest neighbours to target plots in the same cluster (provided they are sufficiently close to the plot centroids). Generally there would more than two auxiliary variables.

<figure><img src="image_45.png" /><figcaption>Figure 1: Illustration of NC sample using two auxiliary variables</figcaption></figure>
## 2.4 Data

The *P. radiata* plantation used in this study is located in South Australia and occupies around 3300 hectares. The plantation is surveyed periodically to provide key product and management information. This study was based on data from 304 ground plots measured in 2012 and a list of plot attributes is given in Table 5.1 of the Final Project Report presented by Rombouts et al. (2014). The essential characteristics of the ground plots are that they have an area of 1000 m² and contain trees with an age range of 14-32 years and a mean age of 26.1.

LiDAR data were acquired in January 2012 by flying an aircraft equipped with an ALTM Orion device at an altitude of 800 m across the estate. The LiDAR first return data had a mean point density of 5.9 pulses m⁻². A detailed description of the LiDAR data specifications is provided in Table 5.2 of Rombouts et al. (2014).

Around 120 separate variables were extracted from the LiDAR point cloud and are summarized in Table 3. These variables were available for each virtual plot. Ten of these variables were eventually selected for imputation according to how well they were able to predict the key attributes of interest using an approach described in Rombouts et al. (2014).

<table><caption>Table 1: Auxiliary variables available for imputation models</caption><thead><tr><th>LiDAR</th><th>first return</th><th>last return</th><th>description</th></tr></thead><tbody><tr><td>pground</td><td>Y</td><td>Y</td><td>proportion of ground returns</td></tr><tr><td>p>x</td><td>Y</td><td>Y</td><td>proportion of heights > x (x=1,2,5,10m)</td></tr><tr><td>sd</td><td>Y</td><td>Y</td><td>standard deviation of heights</td></tr><tr><td>skew</td><td>Y</td><td>Y</td><td>skewness of height distribution</td></tr><tr><td>kurtosis</td><td>Y</td><td>Y</td><td>kurtosis of height distribution</td></tr><tr><td>h</td><td>Y</td><td>Y</td><td>mean height</td></tr><tr><td>hmax</td><td>Y</td><td>Y</td><td>maximum height</td></tr><tr><td>hmax4</td><td>Y</td><td>Y</td><td>four highest in each plot quadrant</td></tr><tr><td>hx</td><td>Y</td><td></td><td>x% percentile height (x=10%~100%)</td></tr><tr><td>d0.x</td><td>Y</td><td>Y</td><td>proportion of heights between 0% and x% hmax (x=10~100%)</td></tr><tr><td>h>0</td><td>Y</td><td></td><td>mean height (heights > 0m i.e. vegetation)</td></tr><tr><td>mqh>0</td><td>Y</td><td></td><td>mean quadratic height (heights > 0m)</td></tr><tr><td>mqh>1</td><td>Y</td><td></td><td>mean quadratic height (heights > 1m)</td></tr><tr><td>scanangle</td><td>Y</td><td></td><td>mean scan angle in the plot</td></tr><tr><td>non LiDAR</td><td></td><td></td><td></td></tr><tr><td>lop</td><td></td><td></td><td>thinning status (last operation)</td></tr><tr><td>nsq</td><td></td><td></td><td>site quality index</td></tr><tr><td>age</td><td></td><td></td><td>plantation age</td></tr></tbody></table>

## 2.5 Simulation approach

Three sampling methods were examined as part of this study. An approach often used in forest inventories is to place a grid over the AOI and select plots at the grid intersection points. LiDAR data are not required to construct either a grid sample or a completely random sample. The use of LiDAR data as *a priori* information in the sample design is aimed at obtaining efficiency gains from the plot selection process. In this study the sampling methods which were evaluated, in addition to random sampling, were locally balanced sampling (Grafström et al., 2014) and the proposed NC method.
In the sampling simulations below, the 304 study plots were divided at random into a target set (200 plots) and a candidate set (104 plots). Reference plots were selected from the candidate plots using each of three sampling strategies i.e. random, locally balanced, and NC. The prediction method used throughout was Euclidean imputation with $k = 1$. The number of reference plots was fixed at either 10, 25, 50 or 75 and the variable of interest was the timber volume in each plot.

The various sampling strategies were evaluated in terms of how well the known variable of interest was predicted. The simulations were repeated 10,000 times with each realisation comprising a new set of target plots, candidate plots and reference plots. Comparisons were done in terms of the relative bias (RB %) and the relative root mean squared error (RMSE %). The relative bias was calculated at the AOI level and is defined as

$$RB = \frac{\sum_j (\hat{Y}_j - Y_j)}{\sum_j Y_j},$$

where $\hat{Y}_j$ is the estimate of total timber volume over the AOI for the j'th realisation and $Y_j$ is the actual total timber volume over the AOI for the j'th realisation. At the AOI level the relative root mean squared error is defined as

$$RRMSE = \frac{\sqrt{\frac{1}{B} \sum_j (\hat{Y}_j - Y_j)^2}}{\frac{1}{B} \sum_j Y_j},$$

where B is the number of realizations. Relative efficiency measures were also calculated for each sampling strategy, using the simple random sample as a benchmark. For any pair of sampling strategies the relative efficiency is calculated as the squared ratio of the AOI-level RMSE values.

### III RESULTS

Table 2 presents the simulation results for the three sampling strategies and the four sample sizes which were investigated.

<table><caption>Table 2: Summary comparisons of bias, accuracy and relative efficiency of three sampling strategies for predicting mean timber volume using field plot measures and LiDAR data</caption><thead><tr><th>Sampling method</th><th>Sample size</th><th>Relative Bias (%)</th><th>Relative RMSE (%)</th><th>Relative efficiency</th></tr></thead><tbody><tr><td>Random</td><td rowspan="3">10</td><td>-0.1</td><td>4.7</td><td>1.0</td></tr><tr><td>Balanced</td><td>0.2</td><td>3.8</td><td>1.5</td></tr><tr><td>NC</td><td>-0.1</td><td>3.1</td><td>2.2</td></tr><tr><td>Random</td><td rowspan="3">25</td><td>0.2</td><td>2.8</td><td>1.0</td></tr><tr><td>Balanced</td><td>0.1</td><td>2.3</td><td>1.4</td></tr><tr><td>NC</td><td>0.5</td><td>2.1</td><td>1.7</td></tr><tr><td>Random</td><td rowspan="3">50</td><td>0.3</td><td>1.8</td><td>1.0</td></tr><tr><td>Balanced</td><td>0.2</td><td>1.7</td><td>1.1</td></tr><tr><td>NC</td><td>0.7</td><td>1.6</td><td>1.2</td></tr><tr><td>Random</td><td rowspan="3">75</td><td>0.2</td><td>1.6</td><td>1.0</td></tr><tr><td>Balanced</td><td>0.3</td><td>1.6</td><td>1.1</td></tr><tr><td>NC</td><td>0.5</td><td>1.5</td><td>1.2</td></tr></tbody></table>
The relative bias is very small with these data, irrespective of the sampling method. With small sized samples the relative RMSE varied from 4.7% using a random sample, to 3.1% using the NC sample. Hence the relative efficiency of the NC sample is 2.2 times that of the random sample. The efficiency gains are most pronounced when the sample size is small. For larger samples, e.g. n= 75, the relative efficiency of the NC sample, compared to the random sample, is 1.2. The relative RMSE results are also presented graphically in Figure 2.

<figure><img src="image_46.png" /><figcaption>Figure 2: Relative RMSE % vs sample size for three sampling schemes</figcaption></figure>
# IV DISCUSSION

This paper describes an approach to inventory design which is specifically aimed at imputation and essentially different to other sampling methods. The results in Table 2 illustrate efficiency gains which are more than double that of random samples. The implications for inventory design are that surveys can be constructed with around half the number of plots that would be required using a conventional sampling approach. The method can be employed either within or across strata although the results in this paper were achieved without stratification.

Note that using remotely sensed data for survey design necessitates having the data available prior to plot selection. Therefore either the data need to be newly acquired or need to be available from a previous campaign. Where existing data are used, they need to be sufficiently recent. The key criteria in this respect is the extent of the correlation between the remotely sensed data and the variables of interest.

One of the primary advantages of the NC sampling method is its flexibility. Small area estimation (SAE) provides a good illustration of this. SAE uses reference plots which are mostly (or completely) outside the AOI to make predictions within the AOI. Most of the existing sampling strategies, including balanced sampling, are not suited to SAE because it is not possible to target the sample specifically to the AOI. There may or may not be reference plots in the sample which are similar enough to the target plots to provide acceptable imputation estimates. By partitioning the small area into defined clusters, the NC sample permits the selection of reference plots which are closely matched to the target plots.

It is proposed to make the clustering approach available as an R function (R Core Team, 2015) to enable sample selection in any application where auxiliary data are available for a population of discrete units such as forestry plots.

# V CONCLUSIONS

The sampling method presented in this paper is an intuitive approach to forest inventory where imputation methods are to be used. The method has been trialled on *P. Radiata* datasets from four separate locations in eastern Australia. In every case it has proved to be superior to other sampling methods. With small samples in particular, relative efficiencies are substantially higher than random sampling methods and moderately higher than balanced sampling strategies.

# References

* Breidenbach J., Nothdurft A., Kändler G. (2010). Comparison of nearest neighbour approaches for small area estimation of tree species-specific forest using airborne laser scanner data. *European Journal of Forest Research* 129, 833–846.

* Dash J. P., Marshall H. M., Rawley B. (2015). Methods for estimating multivariate stand yields and errors using k-nn and aerial laser scanning. *Forestry* 88, 237–247.

* Grafström A., Saarela S., Ene L. T. (2014). Efficient sampling strategies for forest inventories by spreading the sample in auxiliary space. *Canadian Journal of Forest Research* 44, 1156–1164.

* Hartigan J. A., Wong M. A. (1979). A k-means clustering algorithm. *Applied Statistics* 28, 104–108.
Hawbaker T. J., Keuler N. S., Lesak A. A., Gobakken T., Contrucci K., Radeloff V. C. (2009). Improved estimates of forest vegetation structure and biomass with a lidar-optimized sampling design. *Journal of Geophysical Research* 114, GE00E04.

Hudak A. T., Haren A. T., Crookston N. L., Liebermann R. J., Ohmann J. L. (2014). Imputing forest structure attributes from stand inventory and remotely sensed data in western Oregon, *Forest Science* 60, 253–269.

Latifi H., Nothdurft A., Koch B. (2010). Non-parametric prediction and mapping of standing timber volume and biomass in a temperate forest: apredictors. *Forestry* 83, 395–407.

Lesch S. M. (2005). Sensor-directed response surface sampling designs for characterizing spatial variation in soil p<sup>C</sup>omputers and Electronics in Agriculture 46, 153–179.

McRoberts R. E. (2012). Estimating forest attributes for small areas using nearest neighbors techniques. *Forest Ecology and Management* 272, 3–12.

McRoberts R. E., Tomppo E. O., Finley A. O., Heikkinen J. (2007). Estimating aerial means and variances of forest attributes using the k-nearest neighbors technique and satellite imagery. *Remote Sensing of Environment* 111, 466–480.

Melville G., Stone C. (2016). Optimizing nearest neighbour information - a simple, efficient sampling strategy for forestry plot imputation using remotely sensed data. *Australian Forestry*, to appear.

Melville G., Stone C., Turner R. (2015). Application of lidar data to maximize the efficiency of inventory plots in softwood plantations. *New Zealand Journal of Forestry Science* 45:9, 1–16.

R Core Team (2015). *R: A Language and Environment for Statistical Computing*. Vienna, Austria: R Foundation for Statistical Computing.

Rombouts J., Ferguson I., Leech J., Culvenor D. (2010, September 14–17). An evaluation of the field sampling design of the first operational lidar based site quality survey of radiata pine plantations in South Australia. In *Proceedings of the 2010 Silvilaser Conference*, Freiburg, Germany.

Rombouts J., Melville G., Kathuria A., Stone C. (2014). Operational deployment of lidar derived information into softwood resource systems. *Final Report for Project PNC305-1213*, http://www.fwpa.com.au/rd-and-e/resources/611-operational-deployment-of-lidar-derived-information-into-softwood-resource-systems.html.

Stevens D. L. J., Olsen A. R. (2004). Spatially balanced sampling of natural resources. *Journal of the American Statistical Association* 99, 262–278.

Stone C., Penman T., Turner R. (2011). Determining an optimal model for processing lidar data at the plot level: results for a *Pinus radiata* plantation in New South Wales, Australia. *New Zealand Journal of Forestry Science* 41, 191–205.
# Validation of Copernicus High Resolution Layer on Imperviousness degree for 2006, 2009 and 2012

Christophe Sannier<fnref n="1" />*, Javier Gallego<fnref n="2" />, Jochen Dahmer<fnref n="3" />, Geoff Smith<fnref n="4" />, Hans Dufourmont<fnref n="5" />, Alexandre Pennec<fnref n="1" />

*Corresponding author: christophe.sannier@sirs-fr.com

## Abstract

The validation of a dataset such as the Copernicus Pan-European imperviousness degree high resolution layer requires considerable effort. A stratified systematic sampling approach was developed based on the LUCAS sampling frame focusing on a 2 stage stratification approach. A two-stage stratified sample of 20,164 1ha square primary sampling units (PSU) was selected over EEA39 based on countries or groups of countries which area was greater than 90,000km<fnref n="2" /> and a series of omission and commission strata. In each PSU, a grid of 5 x 5 Secondary Sample units (SSUs) with a 20 m step was applied. These points were photo-interpreted on orthophotos with a resolution better than 2.5m.

Initial results based on the binary conversion of the map by applying the 30% threshold indicate a level of omission and commission errors substantially greater than the required maximum level of 15% set in the product specifications. However, this assumes that complete information is available for each PSU which is not the case. An alternative procedure was applied to the quantitative continuous data considering the sampling error due to the SSUs selection which is expected to exhibit a more realistic assessment of the amount of omission and commission.

## Keywords

Sratified systematic sampling approach, LUCAS, Binomial confidence interval, Omission error, Commission error

# I INTRODUCTION

Pan-European High Resolution Layers (HRL) provide information on specific land cover characteristics, and are complementary to land cover / land use mapping such as in the CORINE Land Cover (CLC) datasets (Büttner et al. 2012) as part of the Land Monitoring Service (land.copernicus.eu) of the Copernicus programme, managed by the EC. The HRLs are produced from 20 m spatial resolution satellite imagery through a combination of automatic processing and interactive rule-based classification.

Five themes have been identified so far, corresponding with the main themes from CLC, i.e. imperviousness (the level of sealed soil), tree cover density and forest type, permanent

<footnote n="1">¹Systèmes d'Information à Référence Spatiale (SIRS) SAS, France</footnote>

<footnote n="2">²Joint Research Centre (JRC), Italy</footnote>

<footnote n="3">³GAF AG, Germany</footnote>

<footnote n="4">⁴SpectoNatura, UK</footnote>

<footnote n="5">⁵European Environment Agency (EEA), Denmark</footnote>
grasslands, wetlands and water bodies. Products with an initial pixel size of 20 by 20 m are aggregated into 100 by 100 m grid cells for final pan-European mosaic products. The imperviousness layer was the first to be produced during 2006-2008 from multi-sensor, bi-temporal and ortho-rectified satellite imagery, the same as used for the CORINE Land Cover 2006 update. The production of IMD2006 covered 38 European countries (32 EEA Member States and 6 West-Balkan countries). Since the 2006 production, a time series of imperviousness has been produced for reference years 2009 and recently 2012 over the whole area covered by the 39 member and cooperating countries of the European Environment Agency (EEA) representing a total of 6 million km². For each year it is available as a raster layer with 20 m resolution. At the time of undertaking this study, the area delivered for the year 2012 was around 90% of the total area and the study is based on this area (Spain, Greece, Cyprus and the French overseas regions are missing).

Built-up areas are characterized by the substitution of the original (semi-) natural land cover or water surface with an artificial, often impervious cover. These artificial surfaces are usually maintained over long periods of time. The imperviousness HRL captures the spatial distribution of artificially sealed areas, including the level of sealing of the soil per area unit. The level of sealed soil (imperviousness degree 1-100%) is produced using an automatic algorithm based on a calibrated normalised difference vegetation index (NDVI). A description of the Copernicus Imperviousness layer methodology was described by Gangkofner et al. (2010) for the 2009 update and by Lefebvre et al (2013) for the 2012 update. Similar methods were also applied in the USA for the development of the National Land Cover database (Xian et al. 2011).

A density threshold of 30% was used to derive the built-up layer from the imperviousness layer. This was not intended to be a separate product, but instead was calculated for the verification process only, because density products cannot be verified.

The objectives of this study were to develop and implement an accuracy assessment exercise (i) capable to confirm the results obtained by participating countries during the verification phase based on the built-up mask (Büttner, 2012), (ii) suitable to assess the accuracy of the imperviousness layer at EE39 (iii) whilst ensuring that the results can be analysed at biogeographical region or large country level to ensure that there are no major regional differences.

## II METHODS

The validation of this dataset over such a large area requires considerable and mostly unprecedented effort. Much of the literature on the assessment of the accuracy of impervious surface delineation from remotely sensed data tend to focus on relatively small areas (Ji & Jensen 1999, Yang et al. 2003, Chabaeva et al. 2007). The only comparable example is that of the National Land Cover Database of the Conterminous United States (Wickham et al. 2013) and a recent study from Hansen et al. (2014). However, this study focuses primarily on imperviousness changes and the accuracy assessment is therefore combined with that of the other land cover classes. In fact, relatively low accuracies are reported for the imperviousness change class, which is partially linked to the fact that imperviousness areas although constantly increasing, still occupy a very small portion of the total land area. In Europe, it is currently estimated that artificial surfaces represent less than 5% of the total EEA39 (Büttner et al. 2012) and it should
be expected that impervious surfaces would only represent a subset of this area, which makes it a very rare class and a particular challenge to assess its accuracy with a high degree of precision.

Wickham et al. (2013) applied a stratified random sampling approach using the land cover classes as strata. This ensured that sampling intensity is adapted to the occurrence of each class to reach a sufficient level of user's accuracy even for rare classes (Wickham et al. 2013). For the European Imperviousness layer, a similar approach was applied but based on a stratified systematic sampling approach using the EUROSTAT Land Use / Cover Area fram statistical Survey (LUCAS) sampling frame (Gallego and Delincé 2010) and focusing on a 2 stage stratification approach. The main advantage of using a LUCAS based approach is that a systematic approach ensures full traceability and it is also possible that sampling units will be shared for assessing several products, thus providing potential economies of scale. The first stage was to identify countries or groups of countries with an area greater than 90,000km² (see figure 1 below). This step allowed the analyse of the results for different countries and biogeographical regions to assess whether any heterogeneity in the quality of the data across different regions did emerge.

<figure><img src="image_47.png" /><figcaption>Figure 1. Level of reporting by country or aggregated countries as applied in the external validation.</figcaption></figure>

Then as a second stage, for each country or groups of countries, omission, commission and commission change strata were determined as follows:

* Commission 2006-2009-2012-2015: Imperviousness Degree 30-100% in 2006-2009-2012
* Omission High Probability 2006-2009-2012: Imperviousness Degree 0-29% & CLC impervious classes 2006-2009-2012

* Omission Low Probability 2006-2009-2012-2015: Rest of the area 2006-2009-2012

* Commission Change 2006-2009: all changes [increased and decreased]

* Commission Change 2009-2012 : all changes [increased and decreased]

As indicated the commission strata were included to increase the precision of user's accuracy and the inclusion of a high probability omission error was to increase the precision of the producer's accuracy.

A minimum of 50-100 1ha square primary sampling units (PSU) were deemed sufficient to reach the required precision and were selected per stratum resulting in a total of 20,164 PSUs over the EEA39 area. Each 1ha PSU corresponds to a single 100m pixels in the aggregated imperviousness layer. A detailed procedure to map the imperviousness, involving field data or photo-interpretation, of all PSUs would be too time consuming and expensive. Therefore, in each PSU a grid of 5 x 5 secondary sample units (SSUs) with a 20 m step was defined (see figure 2). These SSUs were photo-interpreted against orthophotos with a spatial resolution better than 2.5m to determine if they were sealed. If a point falls on the boundary of an impervious element, a shifting rule is applied so that roughly half of the points in this situation are classified as impervious.

<figure><img src="image_48.png" /><figcaption>Figure 2: Example of SSUs organised in a 5x5 20m grid</figcaption></figure>

### III THEMATIC ACCURACY BASED ON 30% THRESHOLD

A 30% threshold was applied to the imperviousness degree product to convert the continuous density product to a binary mask (Büttner 2012).

Thematic accuracy is presented in the form of an error matrix made out of the results of the interpretation of the samples and their actual values in the impervious layer. As explained in (Selkowitz & Stehman, 2011), unequal sampling intensity resulting from the stratified systematic sampling approach should be accounted for by applying a weight factor (p) to each

$$ \hat{p}_{ij} = \left(\frac{1}{N}\right) \sum_{x \in (i,j)} \frac{1}{\pi_{uh}} $$
sample unit based on the ratio between the number of samples and the size of the stratum considered:

(1)

where *i* and *j* are the columns and rows in the matrix, *N* is the total number of possible units (population) and *π* is the sampling intensity for a given stratum.

This is because the samples from smaller strata show a higher sampling intensity than those from the larger strata. Therefore, a correction for the sampling intensity will be applied to the error matrices produced, following the procedure described by (Selkowitz & Stehman, 2011) and applied by (Olofson et al., 2013). This leads to a weighting factor inversely proportional to the inclusion probability of samples from a given stratum. Not applying this correction could result in underestimating or overestimating map accuracies.

Initial results based on the binary conversion of the map by applying the 30% threshold indicate a level of omission and commission errors substantially greater than the required maximum level of 15% set in the product specifications. There is also considerable variability across different bio-geographical regions and group of countries and it should be noted that the 2012 layer exhibits higher accuracy to that of the 2009 and 2006 layer, which can be explained by the fact that there are always improvements made during the production of an updated dataset. However, omission errors appear to be heavily influenced by the presence of omissions in the low probability stratum which carry substantial weight. This may suggest that the stratification approach selected for omission probably needs to be revisited.

## IV ANALYSIS OF IMPERVIOUSNESS DEGREE VALUES

In addition, the assessment based on the binary mask assumes that complete information is available for each PSU, which is not the case since the PSU level estimate of imperviousness is based on the SSUs. Therefore, an alternative procedure was applied to the quantitative continuous data considering the sampling error due to the SSUs selection, which is expected to exhibit a more realistic assessment of the amount of omission and commission errors.

If we had a complete information on the cell for our reference data, a reasonable measure of the commission φ and omission ψ errors would be:

$$ \varphi = \frac{\sum_i pos(m_i-r_i)}{\sum_i m_i} \qquad \psi = \frac{\sum_i pos(r_i-m_i)}{\sum_i r_i} \qquad (2) $$

where *pos(x)* is the positive part, i.e. *pos(x) = x* if *x*>0 and *pos(x) = 0* if *x*≤0.

If the map reports a proportion *m*<sub>i</sub> and the reference data give a proportion *r*<sub>i</sub>,

For each sampling unit of 100 m we have a quantitative value in the map (estimated % in the satellite image classification) and a reference value that is an estimation obtained from a sample of 25 points. The number of impervious points that we are using as reference value has a probability distribution due to the within-cell sampling. If the within-sampling is random, the number of points follows a binomial B(25, p). In our case the sampling scheme is systematic, but we use anyhow the binomial as an approximation.
<figure><img src="image_49.png" /><figcaption>Figure 3: Representation of the behaviour of the 95% confidence interval for a 5x5 SSU grid over the whole range of imperviousness degree values</figcaption></figure>

Therefore we cannot say that there is any significant disagreement if $m_i$ lies within ($r_{0i}$, $r_{1i}$), a confidence interval corresponding to B(25, $r_i$). Figure 3 represents the behaviour of the 95% confidence interval for B(25, $r_i$). Notice that only for proportions close to 0.5 we can apply the usual Gaussian approximation that leads to an interval approximately ($r_i \pm 2s_i$), while for proportions close to 0 or to 1 the intervals are strongly asymmetric.

A possible adaptation of the formulas above for the commission $\varphi$ and omission errors $\psi$ would be:

$$ \varphi = \frac{\sum_i \text{pos}(m_i - r_{1i})}{\sum_i m_i} \qquad \psi = \frac{\sum_i \text{pos}(r_{0i} - m_i)}{\sum_i r_i} \qquad (3) $$

The results obtained show that the amount of error is less than 10% for commission and just under 20% for omission, which exceeds the desired level of accuracy for this layer with regards to commission errors. But again this is heavily influenced by omissions included in the large area low omission strata. When removing just 39 points from the sample corresponding to these large strata (out of more than 20,000), the level of omission error is then reduced to 13%, which is within the acceptable range of accuracy expected for the imperviousness layer. This again suggests that the stratification approach for assessing omission errors needs to be revised to better target areas of potential omission.

# V CONCLUSIONS

The results from the accuracy assessment exercise seem to suggest that the Copernicus Imperviousness High Resolution Layer appear to reach or even exceed the required level of
accuracy (less than 10-15% error for both omission and commission errors) particularly with
respect to commission errors, but the stratification approach would need to be further improved
to better target omission errors.

References

Büttner G., Kosztra B., Maucha G. and Pataki R. (2012) *Implementation and achievements of CLC2006*, ETC-LUSI, EEA, 65p.

Büttner G (2012) *Guidelines for verification and enhancement of high resolution layers produced under GMES initial operations (GIO) Land monitoring 2011 – 2013*, EEA, 47p.

Chabaeva, A., Hurd, J., & Civco, D. (2007). Quantitative assessment of the accuracy of spatial estimation of impervious cover. In *ASPRS Annual Conference Proceedigns 2007*. Tampa, Florida.

Gallego, J., & Delincé, J. (2010). The European Land Use and Cover Area-Frame Statistical Survey. *Agricultural survey methods*, 149-168.

Gangkofner, U., Weichselbaum, J., Kuntz, S., Brodsky, L., Larsson, K., & De Pasquale, V. (2010). Update of the European High-resolution Layer of Built-up Areas and Soil Sealing 2006 with Image2009 Data. 30th *EARSeL Symposium 2010: Remote Sensing for Science, Education and Culture*

Hansen, M. C., Egorov, A., Potapov, P. V., Stehman, S. V., Tyukavina, A., Turubanova, S. A., ... & Kommareddy, A. (2014). Monitoring conterminous United States (CONUS) land cover change with web-enabled landsat data (WELD). *Remote sensing of Environment*, 140, 466-484.

Ji, M., & Jensen, J. R. (1999). Effectiveness of subpixel analysis in detecting and quantifying urban imperviousness from Landsat Thematic Mapper imagery. *Geocarto International*, *14*(4), 33-41.

Lefebvre, A., Beaugendre, N., Pennec, A., Sannier, C., Corpetti, T. (2013). Using data fusion to update built-up areas of the 2012 European High-Resolution Layer Imperviousness. In *Proceedings of the 33rd EARSeL Symposium Conference*, Matera, Italy, 3–6 June 2013; pp. 321–328.

Olofsson, P., Foody, G. M., Stehman, S. V., & Woodcock, C. E. (2013). Making better use of accuracy data in land change studies: Estimating accuracy and area and quantifying uncertainty using stratified estimation. *Remote Sensing of Environment*, 129, 122–131.

Selkowitz, D. J., & Stehman, S. V. (2011). Thematic accuracy of the National Land Cover Database (NLCD) 2001 land cover for Alaska. *Remote Sensing of Environment*, June 2011, 115(6), 1401–1407.

Wickham, J. D., Stehman, S. V., Gass, L., Dewitz, J., Fry, J. A., & Wade, T. G. (2013). Accuracy assessment of NLCD 2006 land cover and impervious surface. *Remote Sensing of Environment*, 130, 294-304.

Yang, L., Xian, G., Klaver, J. M., & Deal, B. (2003). Urban land-cover change detection through sub-pixel imperviousness mapping using remotely sensed data. *Photogrammetric Engineering & Remote Sensing*, 69(9), 1003-1010.

Xian, G., Homer, C., Dewitz, J., Fry, J., Hossain, N., & Wickham, J. (2011). Change of impervious surface area between 2001 and 2006 in the conterminous United States. *Photogrammetric Engineering and Remote Sensing*, 77(8), 758-762.
# The polygon overlay problem in electoral geography

Romain Louvet*<fnref n="1" />,<fnref n="2" />, Jagannath Aryal<fnref n="2" />, Didier Josselin<fnref n="1" />,<fnref n="3" />, Christèle Marchand-Lagier<fnref n="4" />,
Cyrille Genre-Grandpierre<fnref n="1" />

**Abstract**

We developed an algorithm for reducing geometric differences between a source and a target dataset. The algorithm tackles the polygon overlay problem in electoral geography before using areal interpolation methods. Our results show that improvement in matching between statistical areas and polling areas can reduce up to 40% of areal interpolation errors. This is applied to two case studies: the city of Avignon, France, and the city of Hobart, Australia.

**Keywords**

polygon overlay problem, areal interpolation, spatial disaggregation, spatial aggregation, COSP

# I INTRODUCTION

When we explain elections based on socio-spatial context, one main methodological issue is the need to compare variables from different sources, i.e. electoral results with sociological and economic variables. The areal units used for mapping these variables, being designed for different purposes, usually come with different spatial resolutions and different boundaries. Therefore, it is an issue when studying the relationship between variables or trying to compare data over time (Fotheringham and Rogerson, 2013). This issue has been defined as one example of the Change Of Support Problem (COSP) called the polygon overlay problem when dealing with incompatible area to area spatial data (Gotway and Young, 2002). Being confronted with this problem, it is necessary to reallocate data from a source dataset to a target dataset, or in other words from the areal units with available data to the areal units of interest, by using areal interpolation methods.

Areal interpolation is a widely known and studied problem in spatial science and encountered with many type of data (Carson, 2013; Fotheringham and Rogerson, 2013). Many methods already exist (Tobler, 1979; Mugglin et al., 1999; Eicher and Brewer, 2001; Mennis and Hultgren, 2006; Reibel and Agrawal, 2007; Krivoruchko et al., 2011; Zhang and Qiu, 2011; Qiu et al., 2012; Lin et al., 2011), have been compared (Goodchild et al., 1993; Lam, 1983; Carson, 2013; Fotheringham and Rogerson, 2013; Do, 2015), and have been implemented. This paper does not propose a new method of areal interpolation but a process to reduce areal interpolation error by improving matching between source and target data.

Areal interpolation methods model spatial distribution based on strong assumptions such as homogeneity, isotropy, and stationarity. These assumptions, more fitted to natural phenomenon,

<footnote n="1">¹UMR ESPACE 7300 CNRS, Université d'Avignon, France</footnote>

<footnote n="2">²University of Tasmania, Discipline of Geography and Spatial Sciences, School of Land and Food</footnote>

<footnote n="3">³LIA, Université d'Avignon, France</footnote>

<footnote n="4">⁴CHERPA, Aix-en-Provence, LBNC, Université d'Avignon, France</footnote>

<footnote>*Corresponding author: romain.louvet@alumni.univ-avignon.fr</footnote>
can fail to model human spatial distribution. Although intelligent and sophisticated methods
can improve greatly their fitness to the actual spatial distribution, they are still subjected of
generating interpolation errors. Therefore the polygon overlay problem is raising questions
about the actual accuracy of analyses of electoral behaviour in their spatial context.

Another less common approach to this problem would be trying to reduce the mismatch between
the source and target dataset. Such a method is based on aggregating algorithms which create
new areal units optimizing an objective. This idea of an automated zoning procedure (AZP)
was developed by Openshaw (1977) to solve the Modifiable Areal Unit Problem (MAUP) and
has been applied as a solution to the polygon overlay problem by Martin (2003). Creating
more fitted areal units could be indeed useful in the case of electoral geography because both
geometric and attributes accuracy of polling areas are actually questionable (Bernard et al.,
2015).

In this study, we combine areal interpolation with our algorithm, similar to AZP and sliver
polygons eliminating tools, in order to reduce first the differences between the source and target
dataset before using areal interpolation. Then it is expected that the results of areal interpolation
will have fewer errors. This idea was applied to two case studies, one in France, and the other
in Australia.

II MATERIAL AND METHODS

<figure><img src="image_50.png" /><figcaption>Figure 1: 2012 national election polling areas (left) and 2010 IRIS statistical areas (right), respectively target and source data of the city of Avignon, France</figcaption></figure>
<figure><img src="image_51.png" /><figcaption>Figure 2: 2013 federal election polling areas (left) and 2011 SA1 statistical areas (right), respectively target and source data of the city of Hobart, Australia</figcaption></figure>

We are comparing two study areas with similar populations but different densities and polling systems: Avignon (figure 1) and Hobart (figure 2).

We choose to use four areal interpolation methods, based on: area weighting, binary dasymetric, Kriging, and Geographically Weighted Regression. As ancillary data, we used building areas from land use for Avignon, and buildings polygons and points for Hobart. For both study areas, we used the number of dwellings as explanatory variable.

0. given
0.1. source as statistical area polygons
0.2. target as polling area polygons
1. intersect = geometric intersection of source and target
2. selected layer = selection from intersect
3. loop for each feature in intersect
3.1. if feature in selected layer do
3.1.1. feature src = feature identifier in source
3.1.2. feature trgt = feature identifier in target
3.1.3. neighbours = get feature neighbours
4. loop for each feature neighbour in neighbours
4.1. neighbour src = neighbour identifier in source
4.2. neighbours trgt = neighbour identifyer in target
4.3. if neighbours trgt == feature trgt
4.3.1. pass
4.4. else if neighbour src == feature src
4.4.1. lengths = list of shared line lengths between feature and neighbour
5. set feature new id in target as the identifier of its neighbour with maximum value in lengths

<caption>Figure 3: Our algorithm in pseudo code</caption>
<figure><img src="image_52.png" /><figcaption>Figure 4: Steps of the implemented algorithm improving matching between source and target data by eliminating sliver polygons and aggregating polygons from the target data according to the source dataset boundaries.</figcaption></figure>

In order to measure the accuracy of the areal interpolation methods we used population counts known at a disaggregated level. We used four error measure metrics. Three of them are central error values: the Root Mean Square Error (RMSE), the Mean Absolute Error (MAE), and the Median Absolute Error (MedAE). The fourth one we used is the Relative Absolute Value (RAV).

Our algorithm (see figure 3 & figure 4) works similarly to sliver polygons eliminating toolsby merging a selected intersection to the neighbour polygon sharing the longest arc. But this algorithm aggregates selected polygons only within defined boundaries, in our case the source data limits. In order to select the intersections we computed source and target area ratios and the thinness index of each intersecting feature.

### III RESULTS & DISCUSSION

Our results comparing Avignon (figure 5, table 1), to Hobart (figure 6, table 2) show that our method can indeed reduce areal interpolation errors, regardless of the areal interpolation method, and that it is linked to how well source and target data originally match. Considering the RAE, our results even show that the total of areal interpolation error can be reduced by 10 percentage points, from 21.5% to 12.6% for Avignon in the case of Areal Weighting, which is approximately 40% half fewer errors. The comparison between Avignon and Hobart show that our algorithm's performance decreased slightly when source and target data already match well and the target dataset has a broader scale than the source dataset, (table 1, table 2).
<figure><img src="image_53.png" /><figcaption>Figure 5: Avignon, Areal interpolation errors of population count, original polling areas and polling areas modified by our algorithm</figcaption></figure>

<figure><img src="image_54.png" /><figcaption>Figure 6: Hobart, Areal interpolation errors of population count, original polling areas and polling areas modified by our algorithm</figcaption></figure>
<table><caption>Table 1: Avignon, geometric precision and interpolation errors, original target data and the new target</caption><thead><tr><th colspan="4">Polling Area Geometric Precision</th></tr><tr><th></th><th>Original Target Data</th><th>New Target Data</th><th>Precision Variation (%)</th></tr></thead><tbody><tr><th scope="row">Smallest Unit (ha)</th><td>7.9</td><td>9.3</td><td>-15</td></tr><tr><th scope="row">Number of Units</th><td>57</td><td>43</td><td>-25</td></tr></tbody><thead><tr><th colspan="4">Intersects with Statistical Area Geometric Precision</th></tr><tr><th></th><th>Original Target Data</th><th>New Target Data</th><th>Precision Variation (%)</th></tr></thead><tbody><tr><th scope="row">Smallest Unit (ha)</th><td>0.001</td><td>0.9</td><td>-99.9</td></tr><tr><th scope="row">Number of Units</th><td>281</td><td>94</td><td>-66</td></tr></tbody><thead><tr><th colspan="8">Interpolation Errors Variation</th></tr><tr><th></th><th colspan="2">RMSE</th><th colspan="2">MAE</th><th colspan="2">MedAE</th><th colspan="2">RAE</th></tr><tr><th></th><th>%</th><th>people</th><th>%</th><th>people</th><th>%</th><th>people</th><th>%</th><th>people</th></tr></thead><tbody><tr><th scope="row">Area Weighting</th><td>-23</td><td>-101</td><td>-24</td><td>-85</td><td>-21</td><td>-66</td><td>-43</td><td>-8489</td></tr><tr><th scope="row">Binary Dasymetric</th><td>-20</td><td>-86</td><td>-23</td><td>-75</td><td>-34</td><td>-90</td><td>-42</td><td>-7954</td></tr><tr><th scope="row">GWR</th><td>+14</td><td>+50</td><td>+18</td><td>+48</td><td>+25</td><td>+56</td><td>-11</td><td>-1697</td></tr><tr><th scope="row">Kriging</th><td>-17</td><td>-66</td><td>-27</td><td>-81</td><td>-47</td><td>-126</td><td>-45</td><td>-7714</td></tr></tbody></table>

<table><caption>Table 2: Hobart, geometric precision and interpolation errors, original target data and new target data</caption><thead><tr><th colspan="4">Polling Area Geometric Precision</th></tr><tr><th></th><th>Original Target Data</th><th>New Target Data</th><th>Precision Variation (%)</th></tr></thead><tbody><tr><th scope="row">Smallest Unit (ha)</th><td>2.3</td><td>36.3</td><td>-94</td></tr><tr><th scope="row">Number of Units</th><td>35</td><td>34</td><td>-3</td></tr></tbody><thead><tr><th colspan="4">Intersects with Statistical Area Geometric Precision</th></tr><tr><th></th><th>Original Target Data</th><th>New Target Data</th><th>Precision Variation (%)</th></tr></thead><tbody><tr><th scope="row">Smallest Unit (ha)</th><td>0.001</td><td>5</td><td>-99.998</td></tr><tr><th scope="row">Number of Units</th><td>415</td><td>238</td><td>-43</td></tr></tbody><thead><tr><th colspan="8">Interpolation Errors Variation</th></tr><tr><th></th><th colspan="2">RMSE</th><th colspan="2">MAE</th><th colspan="2">MedAE</th><th colspan="2">RAE</th></tr><tr><th></th><th>%</th><th>people</th><th>%</th><th>people</th><th>%</th><th>people</th><th>%</th><th>people</th></tr></thead><tbody><tr><th scope="row">Area Weighting</th><td>-4</td><td>-9</td><td>-8</td><td>-13</td><td>-10</td><td>-10</td><td>-26</td><td>-1531</td></tr><tr><th scope="row">Binary Dasymetric</th><td>-10</td><td>-13</td><td>-16</td><td>-17</td><td>-45</td><td>-43</td><td>-33</td><td>-1204</td></tr><tr><th scope="row">GWR</th><td>-10</td><td>-17</td><td>-4</td><td>-5</td><td>+15</td><td>+12</td><td>-23</td><td>-1058</td></tr><tr><th scope="row">Kriging</th><td>-2</td><td>-5</td><td>+1</td><td>+1</td><td>+21</td><td>+21</td><td>-20</td><td>-1022</td></tr></tbody></table>
## References

Bernard L., Marchand-Lagier C., Josselin D., Louvet R. (2015, June). Some ways to estimate the effects of the bad voters registration on electoral participation. In *Congrès AFSP 2015*, Aix-en-Provence, France.

Carson B. D. (2013). *Testing Kriging-Based Areal Interpolation for Census-Based Socioeconomic Data*. Master's thesis, University of Redlands.

Do V. H. (2015). *Les méthodes d'interpolation pour données sur zones*. Ph. D. thesis, Toulouse 1.

Eicher C. L., Brewer C. A. (2001). Dasymetric mapping and areal interpolation: Implementation and evaluation. *Cartography and Geographic Information Science* 28(2), 125–138.

Fotheringham S., Rogerson P. (2013, April). *Spatial Analysis And GIS*. CRC Press.

Goodchild M. F., Anselin L., Deichmann U. (1993). A Framework for the Areal Interpolation of Socioeconomic Data. *Environment and Planning A* 25, 383–397.

Gotway C. A., Young L. J. (2002). Combining incompatible spatial data. *Journal of the American Statistical Association* 97(458), 632–648.

Krivoruchko K., Gribov A., Krause E. (2011). Multivariate Areal Interpolation for Continuous and Count Data. *Procedia Environmental Sciences* 3, 14–19.

Lam N. S.-N. (1983, January). Spatial Interpolation Methods: A Review. *The American Cartographer* 10(2), 129–150.

Lin J., Cromley R., Zhang C. (2011, March). Using geographically weighted regression to solve the areal interpolation problem. *Annals of GIS* 17(1), 1–14.

Martin D. (2003, March). Extending the automated zoning procedure to reconcile incompatible zoning systems. *International Journal of Geographical Information Science* 17(2), 181–196.

Mennis J., Hultgren T. (2006). Intelligent dasymetric mapping and its application to areal interpolation. *Cartography and Geographic Information Science* 33(3), 179–194.

Mugglin A. S., Carlin B. P., Zhu L., Conlon E. (1999). Bayesian areal interpolation, estimation, and smoothing: an inferential approach for geographic information systems. *Environment and Planning A* 31(8), 1337–1352.

Openshaw S. (1977). A Geographical Solution to Scale and Aggregation Problems in Region-Building, Partitioning and Spatial Modelling. *Transactions of the Institute of British Geographers* 2(4), 459–472.

Qiu F., Zhang C., Zhou Y. (2012, September). The Development of an Areal Interpolation ArcGIS Extension and a Comparative Study. *GIScience & Remote Sensing* 49(5), 644–663.

Reibel M., Agrawal A. A. (2007). Areal interpolation of population counts using pre-classified land cover data. *Population Research and Policy Review*, 10–1007.

Tobler W. R. (1979). Smooth pycnophylactic interpolation for geographical regions. *Journal of the American Statistical Association* 74(367), 519–530.

Zhang C., Qiu F. (2011, March). A Point-Based Intelligent Approach to Areal Interpolation. *The Professional Geographer* 63(2), 262–276.
# Modelling uncertainty using geostatistics, a case study in Ecuador

Chicaiza Elena*¹ and Buenaño Xavier¹

¹Technical University of Madrid (UPM), Spain

*Corresponding author: gachisbar@gmail.com

**Abstract**

Spatial data includes positional and thematic components. In this study, the uncertainty analysis has been focused in thematic aspect. A total of 200 soil samples with metal determinations were analyzed. Graphical and numerical tools were employed to summarize and detect critical metals. Copper was identified as the cation of major concern in the study area. Kriging and Co-located Co-Kriging (with Zinc as secondary variable) were employed to carry out predictions. The leave-one-out cross validation method was employed to assess the variogram model used in estimation process, specifically mean square deviation ratio (MSDR) statistic was considered. Finally, unconditional simulations (500 realizations) were computed in order to get mean value and uncertainty of the random field of variable analyzed. The mean of variances obtained from simulations is 1.6 times higher than the mean of kriging variances.

**Keywords**

Uncertainty, simulation, validation.

# I INTRODUCTION

According Griffith et al. (2015), spatial data comprise two components: attribute and location. In this work, we focused the analysis in attribute component.

Currently, the statement that descriptions of spatial phenomena are subject to uncertainty is now generally accepted as explained by Chiles and Delfiner (2009). In this context, Chiles and Delfiner (2009) argues Geostatistics can be defined as "the application of probabilistic methods to regionalized variables".

This study uses R Development Core Team (2013) open-source statistical software. In order to carry out the geostatistical analysis, two specific packages or libraries developed by Ribeiro Jr. and Diggle (2001) (geoR) and Pebesma (2004) (gstat) were employed.

The study area is located in the southeastern of Ecuador, in Zamora province. The area is inside a mining concession where informal mining activities are carried out. The accumulation of heavy metals in soils can bring human health problems, specially when these are exposed directly, for example in this kind of mining activities.

In this study, metals in 200 soil samples were analyzed. A multivariate data exploration analysis was carried out with some metals (Cu, Zn, Cd, and Pb). Kabacoff (2011) recommends the use of relatively new graphical statistical analysis tool named corrgram (see Figure 1) that summarize the relationships between variables very well.

The Figure 1 shows that these variables have a positive correlation (hatching in 45 degrees) with exception of Cd-Pb relation (hatching in -45 degrees); the color saturation in the figure repre-
<figure><img src="image_55.png" /><figcaption>Figure 1: A corrgram of heavy metals soil samples.</figcaption></figure>

sents magnitude of correlation between variables. Notice the significant correlation between Cu and Zn.

## II CALCULATIONS

Copper has been identified as an element of major concern. For this reason, the Cu variable was chosen for the spatial estimation by two methods: a) Ordinary Kriging and b) Co-located Co-Kriging. It's interesting to estimate this variable with a covariate, in this case the Zn metal, in order to define the uncertainty threshold between the two methods employed. The Figure 2 shows variograms of Cu and Zn and the cross-variogram.

<figure><img src="image_56.png" /><figcaption>Figure 2: Direct variograms of Cu, Zn, cross-variogram and fitted functions.</figcaption></figure>

A initial exploration data analysis (EDA) of Cu variable was carried out in order to identify outliers and a possible presence of trend that could limit the application of geostatistical theory.

Through directional clouds plotting, spatial trend was discarded. Directional variograms (0, 45, 90, 135 degrees) confirm this assumption.
An important aspect to notice is that exists a nugget effect in the variogram model of both elements analyzed (Cu and Zn). In this case, the nugget effect is originated in the presence of micro-structure.

The same approach was carried out to analyze Zn variable, and the results were similar.

Some variogram fitting techniques can be employed, but according Oliver and Webster (2014) residual maximum likelihood (REML) method is the current best practice. This approach was carried out in the present study.

According Goovaerts (2001), the choice of an approach for uncertainty modeling should be guided by the answer to these questions: What type of uncertainty model is being sought?, What is the support for the uncertainty assessment?, Why do we model uncertainty?.

In this context, it is important to carry out a validation procedure. There are some ways to validate the model decisions, one of them is the cross validation technique.

Bivand et al. (2008) denotes that cross validation splits the data set into two sets: a modeling set and a validation set. When the number of splits is equal to the number of observations, the procedure is called leave-one-out cross validation.

For the cross-validation, leave-one-out technique was used. There are some statistics to analyze the cross-validation. The most interesting is the mean squared deviation ratio (MSDR), which is the mean of the squared errors divided by the corresponding kriging variances according Oliver and Webster (2014). The formula is shown in Equation 1.

$$MSDR = \frac{1}{N} \sum_{i=1}^{N} \frac{\{z(x_i) - \hat{Z}(x_i)\}^2}{\hat{\sigma}_K^2(x_i)} \quad (1)$$

In these equation $z(x_i)$ is the ith datum at $x_i$, $\hat{Z}(x_i)$ is the kriged prediction there, and $\hat{\sigma}_K^2(x_i)$ is the kriging variance.

Also mean error (ME) and root mean squared error (RMSE) were computed. The results are shown in Table 1.

<table><caption>Table 1: Uncertainty results with two geostatistical methods.</caption><thead><tr><th></th><th>ME</th><th>RMSE</th><th>MSDR</th></tr></thead><tbody><tr><td><strong>Ordinary Kriging</strong></td><td>0.104</td><td>45.318</td><td>1.017</td></tr><tr><td><strong>Co-Kriging</strong></td><td>-0.085</td><td>40.688</td><td>1.099</td></tr></tbody></table>

Another important tool to define uncertainty is simulation. Bivand et al. (2008) defines geostatistical simulation as the simulation of possible realisations of a random field, given the specifications for that random field (e.g. mean structure, residual variogram, intrinsic stationarity) and possibly observation data.

In this context, the drawing of envelopes is an interesting technique. Envelopes were computed assuming a (transformed) Gaussian random field model. According Ribeiro Jr. and Diggle (2001), simulated values are generated at the data locations, given the model parameters obtained by REML. Empirical variogram is computed for each simulation using the same lags as for the original variogram of the data. The envelopes are computed by taking, at each lag, the
maximum and minimum values of the variograms for simulated data. In this case, 99 simulations were carried out. Figure 3 shows the envelopes based in this criteria. The figure shows the range in that we would expect to model the variogram of Copper variable.

<figure><img src="image_57.png" /><figcaption>Figure 3: Envelopes of simulations.</figcaption></figure>

Also, unconditional simulations that ignore observations and only reproduce means and pre-scribed variability, have been carried out. A mean value of 71, corresponding to the mean of prediction map was considering to carried out the simulations.

The results of four (4) first realizations are shown in Figure 4.

<figure><img src="image_58.png" /><figcaption>Figure 4: Unconditional simulation of 4 realizations with REML variogram parameters</figcaption></figure>

A total of five hundred (500) realizations were carried out. The mean of variances obtained in these was computed; this mean is 1.6 times higher than the mean of kriging variances. This is a reasonable value, considering that observations were not included in simulations. Nevertheless
the high values of kriging variances and simulations's mean variances could be explained for nugget effect inclusion in the models. The mean of random field generated by simulations was 70.38 and the mean of kriging predictions was 70.96. The proximity of these values reflects that the random field was modeled correctly.

### III RESULTS

ME of two methods is almost equal. This situation demonstrates that kriging (in general) is an unbiased estimator. RMSE of ordinary kriging is higher than co-kriging due to the auxiliary variable provides more information to the model prediction.

On the other hand, MSDR of Ordinary Kriging is slightly lower than MSDR of co-kriging, probably due to the characteristics of this statistic. Lloyd (2010) in his study concerned with mapping precipitation amount, found that OK outperformed CK (with elevation as the secondary variable) and this was due, at least in part, to the weak global relationship between the two variables.

The map of kriging prediction is shown in Figure 5. The gradient color represents the Copper variable estimation and the contour lines represent variance prediction. The highest variances are located in the corner of southwestern sector. This zone exhibits a lower sampling density in relation with the rest of study area.

<figure><img src="image_59.png" /><figcaption>Figure 5: Map of kriging prediction and variance.</figcaption></figure>

### References

* Bivand R. S., Pebesma E. J., Gomez-Rubio V. (2008). *Applied spatial data analysis with R*. Springer.

* Chiles J.-P., Delfiner P. (2009). *Geostatistics: modeling spatial uncertainty*, Volume 497. John Wiley & Sons.
Goovaerts P. (2001). Geostatistical modelling of uncertainty in soil science. *Geoderma* 103(1), 3–26.

Griffith D., Wong D., Chun Y. (2015). Uncertainty-related research issues in spatial analysis. *Uncertainty Modelling and Quality Control for Spatial Data*, 3.

Kabacoff R. (2011). *R in Action*. Manning Publications Co.

Lloyd C. D. (2010). *geoENV VII – Geostatistics for Environmental Applications*, Chapter Multivariate Interpolation of Monthly Precipitation Amount in the United Kingdom, pp. 27–39. Dordrecht: Springer Netherlands.

Oliver M., Webster R. (2014). A tutorial guide to geostatistics: Computing and modelling variograms and kriging. *Catena* 113, 56–69.

Pebesma E. J. (2004). Multivariable geostatistics in s: the gstat package. *Computers & Geosciences* 30, 683–691.

R Development Core Team (2013). *R: A language and environment for statistical computing*. Vienna, Austria: R Foundation for Statistical Computing. http://www.R-project.org.

Ribeiro Jr. P., Diggle P. (2001). geoR: a package for geostatistical analysis. *R-NEWS* 1(2), 15–18.
# A Hybrid approach for land cover mapping based on the combination of soft classifiers outputs and uncertainty information

Luisa M. S. Gonçalves¹³⁴*, Cidália C. Fonte<fnref n="2" />,<fnref n="4" />

*Corresponding author: luisa.goncalves@ipleiria.pt

## Abstract

In this article, the authors present a hybrid approach to produce more accurate land cover maps of diverse landscapes representing Mediterranean environments. The innovation of the proposed methodology is to use, in the combination of soft classifiers outputs, information about the classification uncertainty associated with each pixel and its neighbourhood.

The hybrid combined classification method developed includes the following steps: 1) definition of the training areas; 2) pixel-based classification using several soft classifiers; 3) computation of the classification uncertainty; 4) application of rules to combine the outputs of the pixel-based soft classifications and the uncertainty information obtained with uncertainty measures, on a pixel and neighbourhood based approach; 5) image segmentation; and 6) object classification based on decision rules that include the results of the combined soft pixel-based classification and its uncertainty. The proposed methodology was applied to an IKONOS and a SPOT-4 multispectral image with, respectively, 4m and 20m spatial resolution.

The overall accuracy of the hybrid classification obtained with the proposed methodology was higher than the one obtained for the individual pixel-based classifications, which shows that this approach may increase classification accuracy. The approach showed to be particularly powerful to obtain Land Cover Map for landscapes representing Mediterranean environments.

## Keywords

Soft classifiers, uncertainty, land cover, IKONOS, SPOT

## I INTRODUCTION

The use of a hybrid classification approach, which combines pixels and objects, has been shown to be suitable for the identification of Landscape Units that contain a variety of land cover objects using Very High Spatial Resolution images. With the combination of a set of classifiers outputs it is possible to obtain a classification that is often more accurate than the individual classifications (Gonçalves et al., 2010; Gonçalves, 2011). Although several approaches have been proposed for combining hard classifications, the development of methods to combine soft classifications and their integration in a hybrid pixel and object based approach is still a field of research. This study tests whether the combination of the outputs of a set of soft classifiers in a hybrid classification approach using uncertainty and contextual information can improve the accuracy of the results. In the proposed hybrid classification

<footnote n="1">¹ Polytechnic Institute of Leiria, School of Technology and Management, Portugal</footnote>

<footnote n="2">² Department of Mathematics, University of Coimbra, Portugal</footnote>

<footnote n="3">³ NOVA IMS - NOVA Information Management School, Universidade Nova de Lisboa, Portugal</footnote>

<footnote n="4">⁴ Institute for Systems and Computers Engineering at Coimbra, Portugal</footnote>
method, the uncertainty information was used in two phases. First, for combining the outputs of three pixel-based soft classifications. This was done through the development of rules that incorporate the information provided by the pixel-based soft classifications and the results given by the application of an uncertainty measure. Second, in the classification of the obtained segmented objects, which represent the Land Units. These are classified through decision rules which include the results of the combined soft pixel-based classification and its uncertainty. The main objective of integrating uncertainty in the classification process is to avoid the use of misclassified pixels in the classification. The first step of the method was applied to different kinds of images, namely to two IKONOS multispectral images and a SPOT-4 multispectral image, with 4m and 20m of spatial resolution respectively, with different Mediterranean landscapes characteristics, to evaluate the outputs reliability. The second step was only applied to an IKONOS multispectral image to obtain a Land Unit Map (LUM).

## II DATA

The study was conducted using three multispectral satellite images with High Spatial Resolution (HSR) from two regions of Portugal (Figure 1), two located in the Alentejo region and one in the Central region. An image of the Alentejo region near Alcácer was used, obtained by the IKONOS sensor, having 4 spectral bands namely (1) blue (0.45-0.52 µm), (2) green (0.52-0.60 µm), (3) red (0.63-0.69 µm) and (4) near infrared (0.76-0.90 µm), with a spatial resolution of 4m and a dimension of 11.8 km by 8.7 Km. The two other images were obtained by the SPOT 4 sensor, having 4 spectral bands, namely (1) green (0.50-0.59 µm), (2) red (0.61-0.68 µm), (3) near infrared (0.78-0.89 µm) and (4) short-wavelength infrared (1.58-1.75 µm), with a spatial resolution of 20m and a dimension of 16.4 km by 15.1 Km.

<figure><img src="image_60.png" /><figcaption>Figure 1: Multispectral images (false color) : a) SPOT-4 image of the Alentejo region, located in the municipal district of Setubal; b) SPOT-4 image of the Central region located in the municipal district of Marinha Grande; c) IKONOS image of the Alentejo region located in the municipal district of Alcácer.</figcaption></figure>

The locations of the images chosen to carry on the study are representative of the main forest species and Mediterranean landscapes in Portugal. The SPOT-4 image is located in the municipal district of Marinha Grande, in the Central region, near the coast, and includes areas with different characteristics, such as: built-up areas; agricultural fields and forest. The dominant forestry species in the region are eucalyptus and coniferous trees.

The SPOT-4 image is located in the municipal district of Setubal, in the Alentejo region, is occupied mainly by agriculture, pastures, forest and agro-forestry areas. The dominant forest species are eucalyptus, coniferous and cork trees.

The IKONOS image is located in the municipal district of Alcácer, in the Alentejo region. It is also occupied mainly by agriculture, pastures, forest and agro-forestry areas, where the dominant forest species are eucalyptus, coniferous and cork trees.

## III METHODOLOGY

The main goal of the classification is to obtain a LUM, where the objects that represent Landscape Units (LU) have a mean area of 0.5 ha and are classified using a hybrid
classification approach that integrates the combination of the outputs of pixel-based soft classifications and their uncertainty. The combination of the pixel-based soft classification is made through the application of rules that incorporate the information provided by the pixel-based classifications and the results given by the uncertainty measures on a pixel and neighbourhood based approach. This approach considers not only each pixel-based classification and its uncertainty information but also the classification and uncertainty of its neighbourhood. Three soft classifiers were used in this application: 1) the neural network Multi-Layer Perceptron Classifier (MLPC); 2) a pixel-based supervised fuzzy soft classifier based on the underlying logic of Minimum-Distance-to-Means (FC); and 3) Maximum Likelihood Classifier (MLC). The classifiers were trained using the same sampling protocol that included 100 pixels per-class. The classes used in this study to identify Surface Elements (SE), which are the basic elements of landscape, are: Eucalyptus Trees (ET), Cork Trees (CKT), Coniferous Trees (CFT), Shadows (S), Shallow Water (SW), Deep Water (DW), Herbaceous (H), Sparse Herbaceous (SH), Non Vegetation Area (NVA), Irrigated Herbaceous (IH) and Non Irrigated Herbaceous (NIH). These classification methods assign to each pixel and to the each class under consideration, in the case of MLPC, different degrees of assignment, in the case of FC, different degrees of possibility and in the case of MLC, different degrees of probabilities. This extra data provides additional information at the pixel level that allows the assessment of the classification uncertainty. To analyze if the use of an output combination of a set of soft classifiers in a hybrid approach classification can improve the accuracy of the results, a similar method, where the combination of the three classifiers and the classification uncertainty was not considered, is also presented. The hybrid classification method that was developed includes the following steps: 1) definition of the training areas; 2) pixel-based soft classification using several soft classifiers; 3) computation of the classification uncertainty; 4) development and application of rules to combine the soft classifications, incorporating the information provided by the pixel-based classifications the results given by the uncertainty measures, on a pixel and neighbourhood based approach; 5) image segmentation; and 6) object classification based on decision rules that include the results of the combined soft pixel-based classification and its uncertainty. The second method that does not take into consideration the output combination of the pixel-based classifications and their uncertainty includes three steps: 1) pixel-based classification of the image; 2) image segmentation and 3) object classification based on decision rules. This article focuses only on the first methodology since the second one was already presented in (Gonçalves et al., 2009a).

## Combination of pixel-based soft classifications

The first phase of the algorithm developed to combine pixel-based soft classifications checks whether the same class is assigned to each pixel by all classifiers. If this condition is satisfied, the class is accepted. If the output classes for each individual pixel differed, the uncertainty information is compared and the class assigned with the lower value of uncertainty is chosen to be the one assigned to the pixel. In this approach the uncertainty measure E, developed by (Chow, 1970), was used to quantify the uncertainty at each spatial unit. This measure is given by

$$E=1-p(x_1) \quad (1)$$

where $p(x_1)$ is the largest degree of possibility or probability of the possibility distributions or probability distributions assigned to the pixel corresponding to the several classes. This measure is also called ambiguity measure.

If the classifiers have different results for a certain pixel but the ambiguity values are equal the results of the neighborhood classification is used to make a judgment.
To evaluate if the combined classification improves the results, the accuracy assessment was made with the same protocol used with the single classifiers and the results were compared.

One of the classifiers used was the MLP neural network, which is a non-parametric method and is the most commonly used in remote sensing. Details of the MLP can be found in Atkinson and Tatnall (1997) and in Brown et al. (2009). The MLP provides an activation level for every output class of each pixel, and for hard classifications each pixel is allocated to the class with the largest activation level. A soft classification may be derived from this classifier by considering the activation levels of the network output units for each pixel. These activation levels range from 0 to 1, and may be used as indicators of the uncertainty associated with the pixel allocation to the classes. The other classifier used was a pixel-based supervised fuzzy classifier based on the underlying logic of the Minimum-Distance-to-Means classifier. Details of the fuzzy classifier can be found in (Gonçalves et al., 2009b). The third was a supervised Bayesian classifier similar to the maximum likelihood classifier. The maximum likelihood classifier is based on the estimation of the multivariate Gaussian probability density function of each class using the classes statistics (mean, variance and covariance) estimated from the sample pixels of the training set. Details can be found in (Gonçalves et al., 2009b).

To evaluate the classification accuracy of the individual soft classifications and the combined results, a stratified random sampling with about 100 pixels per class was selected considering the entire image scene, which also included mixed pixels. The number of pixels was chosen to obtain a standard error of 0.05 for the estimation of the accuracy indexes of each class. Each land cover class was sampled independently and the accuracy assessment was made with an error matrix.

## Hybrid Classification to obtain a Landscape Unit Map

The LUM was built using the combined output of the pixel-based classification, its ambiguity information and the objects obtained with the segmentation algorithm. In the segmentation stage the whole image was partitioned into a series of closed objects, corresponding to the spatial patterns. The extraction of the objects was driven using the “Fractal Net Evolution Approach” (FNEA) segmentation method, implemented in eCognition software, which can be described as a region merging technique (Baatz and Schape, 2000).

In this study only one segmentation level was considered, chosen from a series of experiments done with different parameters, the results of which were visually analyzed. The criterium that led to their choice was the identification of meaningful image-objects i.e., groups of pixels that represented the LU existing in the study area, with a mean area of 0.5 ha. The next step was the development of rules that incorporate the information provided by the combined pixel-based classification within each object and the results given by the ambiguity measure E. The classification of the LU is similar to a decision tree which, for geographical objects, is a hierarchical structure consisting of several levels. At each level a test is applied to one or more attribute values. The application of a rule results either in a leaf, allocating an object to a class, or a new decision node, specifying a further decision rule. In this study eight LU classes were used: Water Bodies (WB), Agriculture and Pasture Areas (A), Non-Vegetated Areas (NVA), Broad-Leaved Forest (BF), Coniferous Forest (CFF), Cork Forest (CF), Agro-Forestry Areas (AFA) and Mixed Forest (MF). Table 1 shows the classification rules, based on the SE identified in the previous phase. The structure of the rules is based on the ones used on the study performed by (Gonçalves *et al.*, 2009a).

For the accuracy assessment, the sampling unit to assess the accuracy of the LUM was a fixed-area square plot sampling unit with an area of 0.5 ha. A stratified random sampling of about 50 samples per class was chosen, which guarantees a standard error of 0.07 for the Conditional
Probability of the Map (CPM) and Conditional Probability of Reference (CPR) estimates for each class, assuming that the classification accuracy is superior to 50%, which is acceptable because the construction of the LUM already involved a prior pixel-based classification and an analysis of the terrain. The accuracy assessment was made with an error matrix, where the pᵢⱼ entry is the proportion of area that is class i in the map and class j in the reference within the square areas with 0.5 ha. The CPM and CPR accuracy parameters were then derived from the error matrix.

<table><caption>Table 1. Object's classification rules.</caption><thead><tr><th>Rules</th><th>Test</th><th>Class if true</th></tr></thead><tbody><tr><td rowspan="2">Rule 1</td><td>Objects which have more than 10% of SE classified as tree crowns, regardless of species, with ambiguity of less than 0.5</td><td>Forest</td></tr><tr><td>Objects which do not satisfy the previous test</td><td>Non-Forest</td></tr><tr><td rowspan="3">Rule 2</td><td>The mode of the SE, inside the object, with ambiguity of less than 0.5 is Deep Water or Shallow Water</td><td>Water Bodies</td></tr><tr><td>The mode of the SE, inside the objects, with ambiguity of less than 0.5 is Herbaceous Vegetation or Sparse Herbaceous Vegetation</td><td>Agriculture</td></tr><tr><td>The mode of the SE, inside the objects, with ambiguity of less than 0.5 is Non-Vegetated Area or Shadow</td><td>Non-Vegetated Areas</td></tr><tr><td rowspan="4">Rule 3</td><td>Eucalyptus Trees represent more than 75% of the trees or objects that have only Broad- Leaved Trees inside</td><td>Broad-Leaved Forest</td></tr><tr><td>Coniferous Trees represent more than 75% of the trees or objects that have only Coniferous Trees inside</td><td>Coniferous Forest</td></tr><tr><td>Cork Trees represent more than 75% of the trees or objects that have only Cork Trees inside and the percentage of Herbaceous or Sparse Herbaceous is inferior to Cork Trees</td><td>Cork Forest</td></tr><tr><td>Objects that do not satisfy the previous test</td><td>Mixed or Non-</td></tr><tr><td rowspan="2">Rule 4</td><td>The percentage of trees is less than 50%; the percentage of Herbaceous or Sparse Herbaceous is superior to Cork Trees and 80% of trees are Cork Trees with ambiguity of less than 0.5</td><td>Agro-Forestry Areas</td></tr><tr><td>Objects that do not satisfy the previous test</td><td>Mixed Forest</td></tr></tbody></table>

## IV RESULTS

### Combined Pixel Classifications

The accuracy assessment of the combined classification was made with the same testing datasets used to evaluate the individual classifications. Table 2 shows the Global Accuracy (GA) results.

<table><caption>Table 2. Global Accuracy of the classifications.</caption><thead><tr><th>Classifiers</th><th>Marinha Grande (SPOT-4)</th><th>Setubal (SPOT-4)</th><th>Alcacer (IKONOS)</th></tr></thead><tbody><tr><td>FC</td><td>73.8%</td><td>70.4%</td><td>65.5%</td></tr><tr><td>MLPC</td><td>75.6%</td><td>72.4%</td><td>64.5%</td></tr><tr><td>MLC</td><td>86.9%</td><td>73.7%</td><td>66.9%</td></tr><tr><td>Combined</td><td>89.2%</td><td>82.7%</td><td>72%</td></tr></tbody></table>
The GA of the combined classification of the Marinha Grande, Setubal and Alcácer regions was 89%, 82% and 72% respectively which represents a higher value than that of the most accurate individual classification. Figure 2, Figure 3 and Figure 4 show the classification results when each pixel is assigned to the class with a higher degree of probability for the MLC classifier, higher degree of possibility for the FC classifier, and with the largest activation level for the MLPC classifier for the three areas classified. The spatial distribution of the ambiguity measure E is also presented in the same figures, below the classification results. The regions with the larger ambiguity (darker zones) are the ones where the assignment degrees were lower.

<figure><img src="image_61.png" /><figcaption>Figure 2: Hard version of the classification results of the Setubal area (above) and spatial distribution of ambiguity (below) with: a) MLC classifier b) MLPC classifier and c) FC classifier. Graph d) shows the mean ambiguity per class.</figcaption></figure>

<figure><img src="image_62.png" /><figcaption>Figure 3: Hard version of the classification results of the Marinha Grande area (above) and spatial distribution of ambiguity (below) with: a) MLC classifier b) MLPC classifier and c) FC classifier. Graph d) shows the mean ambiguity per class.</figcaption></figure>