<figure><img src="image_63.png" /><figcaption>Figure 4: Hard version of the classification results of the Alcacer 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>

The comparison of the mean ambiguity per class shows that for the Setubal area, almost all the classes were assigned with lower ambiguity with the FC classifier. The only exception was the Non-Vegetation Area which was assigned to the pixels with lower ambiguity by the MLC classifier. For the Marinha Grande area, SH was the class assigned to the pixels with more uncertainty by all classifiers, almost all the classes were assigned with lower ambiguity by the MLC classifier. The only exception was Herbaceous (H) that was assigned with lower mean ambiguity by the FC classifier.

For the Alcacer area, the comparison of the mean ambiguity per class shows that forest species, such as CKT and CFT, were assigned to the pixels with lower ambiguity with MLPC classifier. The class DW, SW and NVA classes were assigned to the pixels with lower ambiguity with MLC. The FC classifier presents the higher values of ambiguity. The mean ambiguity presented by the classifiers depends essentially on the characteristics of the image.

## Land Unit Map

The Global Probability (GP) classification accuracy for the LUM, obtained with the hybrid approach which classifies the LU (segmented objects) using the combined classification and the ambiguity information, was 72%. The best result for the GP classification with the hybrid approach, using an individual classification (the ones that had the best GA results) and without the ambiguity information, was 59%. This shows that the accuracy increased significantly with the combined classification and the inclusion of the ambiguity information.

Figure 5 allows the comparison between the results of the CPR and CPM accuracies for the LUM obtained with both hybrid classification approaches. These show that the classification results obtained with the method using uncertainty are considerably better for almost all LU classes and this improvement is more evident for the forest classes. Figure 5 shows the final results of the classification with the proposed hybrid classification method.

<figure><img src="image_64.png" /><figcaption>Figure 5: Conditional probability of reference (CPR) and Conditional probability of the map (CPM) obtained with the hybrid approach (LUM) with and without combined classification and uncertainty.</figcaption></figure>
<figure><img src="image_65.png" /><figcaption>Figure 6: Land Unit Map.</figcaption></figure>

## V CONCLUSIONS

The goal of this study was to evaluate if the proposed hybrid combined classification method, using uncertainty information associated with the pixel-based soft classification, produces more accurate land cover maps of diverse landscapes representing Mediterranean environments. Another goal was to investigate if the proposed method that combines the outputs of different soft classifications, using the uncertainty information to choose the best class to assign to each pixel, when applied to different multispectral satellite images always produces more accurate Land Units. The method was applied to IKONOS and SOPT-4 multispectral satellite images with 4m and 20m spatial resolution and the results have shown that with the combined method the accuracy results in both cases were higher. With the proposed hybrid pixel-object classification, the global accuracy of the classification of Mediterranean Landscape Unit classes increased by 11% when compared to a similar classification method that does not use combined classification and ambiguity information. These results have confirmed that the information provided by the uncertainty measure was useful for the combined classification process and in the construction of the Landscape Unit Map because it allowed the determination of the best class to assign to the pixels.

## Acknowledgements

This work has been supported by the Fundação para a Ciência e a Tecnologia (FCT) under project grant UID/MULTI/00308/2013.

## References

* Atkinson P. M., Tatnall A. R. L. (1997). Neural networks in remote sensing. *International Journal of Remote Sensing* 18, 699-709.

* Baatz M., Schape A. (2000). Multiresolution segmentation – an optimization approach for high quality multi-scale image segmentation. In *Proceedings of Angewandte Geographische Informationsverarbeitung XII*. Beitrage Zum AGIT – Symposium Salzburg 2000, Karlsruhe: Herbert Wichmann Verlag, pp. 12-23

* Brown K. M, Foody G. M., Atkinson P. M. (2009). Estimating per-pixel thematic uncertainty in remote sensing classifications. *International Journal of Remote Sensing* 30, 209-229.

* Chow C. K. (1970). On optimum error and reject tradeoff. *IEEE Transactions on Information Theory*, 16, 41-46

* Gonçalves L. M., Fonte C. C., Júlio E. N. B. S., Caetano M. (2009a). A method to incorporate uncertainty in the classification of remote sensing images. *International Journal of Remote Sensing* 30, 5489-5503.

* Gonçalves L M, C C Fonte, E N B S Júlio & M Caetano (2009b). Spatial Data Quality From Process to Decisions, In *Proceedings of 6th International Symposium on Spatial Data Quality 2009*, Boca Raton: CRC Press, pp. 163-177

* Gonçalves L. M., Fonte C. C., Caetano M. (2010). Using uncertainty information to combine soft classifications. In *Proceedings of the 13th International Conference on Information Processing and Management of Uncertainty in Knowledge-Based Systems*, IPMU 2010, Berlin Heidelberg: Springer-Verlag, LNAI 6178, pp. 455-463.

* Gonçalves L. M. S. (2011). A hybrid approach classification of remote sensing images. In *Proceedings of the 31st EARSEL Symposium of the European Association of Remote Sensing Laboratories*, H. Halounová editors Czech Technical University, pp. 328-335.
# Sampling to validate a global cropland map.

Javier Gallego*<fnref n="1" />, Anne Schucknecht<fnref n="1" />, François Waldner<fnref n="2" />

*Corresponding author: javier.gallego@jrc.ec.europa.eu

**Abstract**

A global cropland map is being produced with a 300 m resolution in the SIGMA Project of the EU (Stimulating Innovation for Global Monitoring of Agriculture). The map, based on PROBA-V 300 m images of 2015, is not yet available but the collection of reference data is ongoing on a two-tier sample of plots photo-interpreted on publicly available Very High Resolution (VHR) images with the help of two versions of a Geo-wiki tool: a large sample (initially 40000 plots) photo-interpreted by volunteers (crowdsourcing) and a sub-sample of 4000 plots photo-interpreted by experts. A stratification based on a cropland probability map has been used to select a higher rate in the most difficult areas where the cropland probability is between 25% and 75%. Stratified systematic sampling with a variable number of replicates has been used. The sampling scheme has been assessed using as pseudo-truth the cropland probability map in the European Union. The conclusions suggest that the variance efficiency of the sampling scheme due to the improved spatial homogeneity is moderate. A more important gain appears in terms of traceability of the sample.

**Keywords**

Global cropland map; stratified systematic sample; crowdsourcing for validation.

# I INTRODUCTION

The SIGMA project (Stimulating Innovation for Global Monitoring of Agriculture, http://www.geoglam-sigma.info/Pages/default.aspx) is elaborating a global cropland map based on PROBA-V 300 m images (Dierkx et al., 2014). The project is funded by the research framework program FP7 of the European Union. This paper presents the sampling scheme that has been elaborated for the validation of this map. The aim is ensuring at the same time 1) flexibility, 2) homogeneous spatial layout and 3) improved traceability. The solution we chose was a systematic sample with ranked replicates based on 1° latitude-longitude cells.

The concept of cropland for this map corresponds rather to “annual crops in the current year” rather than the standard concept of cropland in agricultural statistics (FAO, 1996). The chosen definition includes a minimum size threshold of 0.25 ha and a minimum width of 30 m. The SIGMA global cropland map is a binary map (cropland/non cropland) with the reference year 2015. The map will be delivered in latitude-longitude coordinates. Therefore the area of each pixel changes with the latitude. For the sample selection the same latitude-longitude coordinates have been used rather than a cartographic projection.

<footnote n="1">¹ Joint Research Centre. Ispra, Italy</footnote>

<footnote n="2">² Université Catholique de Louvain, Belgium</footnote>
A two-tier sample is applied for the validation: a large sample of around 40,000 units is being photo-interpreted on very high resolution (VHR) images by volunteers that label as crop/non-crop a grid of 9 points inside each pixel. Experts photo-interpret the whole pixel after automatic segmentation for a smaller subsample of around 4000 units. Validators will interpret sub-pixels or polygons as cropland or non cropland based on a very high resolution imagery such as Google Earth and medium resolution time-series. Additionally, there will be the possibility to opt for “unknown/undetermined”, if the validator does not know whether the sample is cropland or not.

## II SAMPLING SCHEME

The collection of reference data for validation is being carried out before the cropland map to be validated is available. For this reason the map itself is not used for stratification as it is usually recommended (Olofsson et al., 2014). However a stratification has been built on the basis of the IIASA cropland probability map that has been produced by comparing different thematic maps (Fritz et al., 2015). Areas with very high or very low probability are considered easier to classify and need a lower sampling rate. Pixels with an attributed crop probability between 25% and 75% are sampled with a higher rate (table 1). The fact that the sampling plan does not use the SIGMA cropland map for the stratification might make the reference data set more easily usable to validate other cropland maps, as long as they have a compatible geometry.

A Global Agro-Environmental Stratification (GAES) provided by Alterra and FAO (Mücher et al., 2016), which defines strata based on agrosystem characteristics is applied. These zones will be used as post-strata for the computation, but were not used at the sampling stage.

<table><caption>Table 1. Strata definition from the IIASA cropland probability map</caption><thead><tr><td>Stratum</td><td>Cropland probability</td></tr></thead><tbody><tr><td>1</td><td>0</td></tr><tr><td>2</td><td>1-25%</td></tr><tr><td>3</td><td>25-75%</td></tr><tr><td>4</td><td>&gt;75%</td></tr></tbody></table>

The sampling unit is the pixel of the map (300 m). The option of using 3 x 3 pixels windows to limit the impact of location inaccuracy has been discussed and discarded to reduce the photo-interpretation effort. Comparing the accuracy computed from the centre and the peripheral part of the 300 m unit should help quantifying the impact of location errors.

The usual approach to sample for the validation of a land cover map is independent random sampling in each of the classes of the map, considered as strata. When field work is necessary, sampling units are grouped into clusters to optimize logistics, but this has little interest when the reference data are obtained by photo-interpretation on publicly available images. An alternative approach is systematic sampling with a random starting point. Systematic sampling provides unbiased estimators of the parameter being assessed and is always more efficient than random sampling if the spatial correlation of the parameter of interest decreases with the distance (Bellhouse, 1988) and in particular for land cover data obtained from remote sensing (Dunn and Harrison, 1993). However systematic sampling has some drawbacks:

* It lacks flexibility if the sample size needs to be modified for practical reasons, adding more units in certain strata or reducing the initially foreseen sample (Stehman, 2009).

* Stratified systematic sampling may result in a loss of the spatial homogeneity of the sample distribution that is the basis of its good performance. This happens in particular when the strata are not defined by a limited number of large compact patches, but by scattered small areas.
* There is no unbiased estimator for the variance. The usual estimator for the variance may strongly overestimate it, although some alternatives have been proposed to reduce the overestimation (Wolter, 1984). We use a bi-dimensional adaptation of one of the formulas Wolter had assessed for unidimensional systematic sampling:

$$ \hat{V} = \frac{\sum_{j \neq j'} w_{jj'} \delta_{jj'} (y_j - y_{j'})^2}{2n \sum_{j \neq j'} w_{jj'} \delta_{jj'}} \quad (1) $$

where the weight $w_{jj'}$ is an average of the weights $w_j$ and $w_{j'}$; $\delta_{jj'}$ is a decreasing function of the distance between j and j', usually with zero value beyond a moderate distance.

An important advantage of systematic sampling appears when the estimates are politically sensitive and their acceptance is problematic or if there is a conflict of interest. For example if the validation is performed by the same team that has produced the map under validation, they can be tempted to eliminate difficult points, in particular if an economic interest is involved. It is more difficult to proof that a random sample has not been manipulated than it is with a systematic sample.

Figure 1 illustrates the irregular spatial layout of a random sample compared with a systematic sample and with a stratified systematic sample with the strata defined above based on the IIASA cropland probability map. In the stratified systematic sample the higher sampling rate for stratum 3 (higher uncertainty of cropland presence) has been applied by applying a grid step of 40 km instead of 70 km used for strata 2 and 4. A simple visual inspection suggests that the regular spatial distribution that gives some advantage to the systematic sampling is lost to a large extent with strongly mixed-up strata.

<figure><img src="image_66.png" /><figcaption>Figure 1: Random and systematic samples in central Europe in an equal-area projection.</figcaption></figure>

An additional issue with systematic sampling appears if we work with global products, such as the one we are considering, and we want to use geographic coordinates (latitude-longitude) in the sampling process. A straight forward systematic sampling in latitude-longitude is strongly unbalanced because the area of each cell changes with the latitude (111 km x 111 km x cos(lat) assuming a spherical shape).

One possible way to deal with this issue is keeping the systematic sample in latitude-longitude coordinates. This means we have an unequal sampling probability $\pi$ and we have to apply $1/\pi$ as weight to the sample elements (Horvitz-Thompson estimator, Cochran, 1977). However this solution is not very efficient because the higher sampling rate far from the equator is not justified.

Figure 2 illustrates the approach we have used: Figure 2A represents a systematic sample built with three replicates in a 1° x 1° cell. Each replicate is the set of locations with the same relative position in each cell of the grid. Even if the area covered by the illustration is not very large, the higher density of the sample in northern latitudes is visible. This uneven density can be
rebalanced by modifying replicates in such a way that the number of points per replicate for a given latitude band (ring) is proportional to the land area in the corresponding ring (Figure 2B). Since the parallel at latitude $α$ has approximately a length $cos(α)$ compared to the equator, a proportion $1-cos(α)$ of points belonging to replicate 1 is downgraded to replicate 2, a proportion $2*(1-cos(α))$ is downgraded from replicate 2 to 3 and so on. Very quickly we have a proportion $rep*(1-cos(α))>1$ and points may be downgraded for example from replicate 5 to replicate 8. In the current scheme the selection of points to be downgraded is purely random.

Figure 2B represents a stratified sample based on two modified replicates: stratum 1 is disregarded, 1 replicate is kept for strata 2 and 4 and 2 replicates are kept for stratum 3 (stratified sample 1-2-1). We can see that the geographical homogeneity is broken with some empty rows of cells. Some of the corresponding points have been dropped in the random selection of points to be downgraded from one to another replicate.

<figure><img src="image_67.png" /><figcaption>Figure 2: Low density sampling schemes in latitude-longitude coordinates</figcaption></figure>

Figure 3 gives an example of stratified systematic sample in latitude-longitude coordinates with 5 modified replicates for strata 2 and 4 and 10 for stratum 3 (stratified sample 5-10-5). It is similar to the stratified 1-2-1 above, but the irregularities are much smaller because the random subsampling for downgrading replicates happens only in the last replicate.

<figure><img src="image_68.png" /><figcaption>Figure 3: A higher density stratified sample with a systematic of replicates.</figcaption></figure>
### III SOME SIMULATIONS

At the moment of writing this paper, the SIGMA cropland map is not yet available and the reference data collection is still ongoing. For this reason we limit the scope of this paper to assessing the behaviour of the sampling schemes described above with two sets of data: the IIASA cropland probability map that we have used for stratification in SIGMA is used as pseudo-truth to be estimated. Considering a target on which we have an exhaustive knowledge allows us to make as many simulations as we wish and to distinguish the real variance of a sampling approach from the estimated variance computed with an empirical formula. We are aware however that the spatial structure of this map is not necessarily similar to the spatial structure of the inaccuracies of the SIGMA cropland map that are the target of our study. For the stratified versions we have chosen CORINE Land Cover (CLC) as data source (EEA 2007). For this reason, we have limited the simulation to the European Union (EU). The CLC land cover classes have been grouped into 3 strata: arable land, other agricultural classes and non-agricultural, that have been sampled with rates proportional to 3-2-1.

Table 2 summarises some of the results obtained. Simple random sampling (srs) is the benchmark with which the other standard errors are compared. Sampling in latitude-longitude coordinates, as if it were an equal area projection, introduces a bias of 7%. Therefore we do not even compute variances or relative efficiency, since such a bias is a sufficient reason to reject the method whatever the variance of the estimates. The same bias appears both with random and with systematic sampling: it has nothing to do with the sampling approach itself, but with the use of a reference system that is not equal-area. The other sampling schemes confirmed in the simulation to be unbiased. Simple systematic sampling has a relative efficiency of 1.22 (i.e. the variance is divided by 1.22 with the same sample size). This standard error could be computed in our example only because it was a simulation and we had an exhaustive knowledge of the sampled population. We can notice that using for systematic sampling the usual variance estimator for random sampling has a positive bias (overestimation) of 23%, leading to an “apparent efficiency” of 0.99. The estimator (Equation 1) based on local variance has a slightly positive bias and reports an estimated efficiency of 1.19. It should be noticed that the sample size in some approaches is not the same as the reference sample size of 1000 units we have chosen for the comparison. This highlights the rigidity of the usual systematic sampling towards the sample size.

<table><caption>Table 2. Comparison of mean and std. error of cropland probability for some sampling schemes.</caption><thead><tr><th>Sampling method</th><th>Mean</th><th>n</th><th>Std error</th><th>relative efficiency</th></tr></thead><tbody><tr><td>Simple random Sample (srs) equal area</td><td>53.1</td><td>1000</td><td>1.306</td><td></td></tr><tr><td>Srs lat-long</td><td>49.45</td><td>****</td><td>****</td><td>****</td></tr><tr><td>Systematic latlong</td><td>49.45</td><td>****</td><td>****</td><td>****</td></tr><tr><td>Systematic equal-area</td><td>53.1</td><td>889</td><td>1.256</td><td>1.22</td></tr><tr><td>Syst. Eq-area usual variance estimator</td><td>53.1</td><td>889</td><td>1.391</td><td>0.99</td></tr><tr><td>Syst. Eq-area local variance</td><td>53.1</td><td>889</td><td>1.268</td><td>1.19</td></tr><tr><td>Systematic latlong balanced replicates</td><td>53.1</td><td>1000</td><td>1.213</td><td>1.16</td></tr><tr><td>Systematic latlong balanced replicates</td><td>53.1</td><td>100</td><td>4.07</td><td>1.03</td></tr><tr><td>Stratif random</td><td>53.1</td><td>1000</td><td>1.103</td><td>1.40</td></tr><tr><td>Stratified syst. equal area independent</td><td>53.1</td><td>959</td><td>1.056</td><td>1.59</td></tr><tr><td>Stratified systematic latlong balanced</td><td>53.1</td><td>1000</td><td>1.049</td><td>1.55</td></tr></tbody></table>
The systematic sampling schemes that we have labelled in the table as “balanced replicates” refers to the process described above to reduce the number of units in each replicate by a factor cos(latitude) (Figure 2B and 2C). If we consider a relatively large sample (n=1000 corresponding to a bit less than 3 replicates), we get a reasonably good efficiency of 1.16. In exchange a small sample of 100 units (a portion of one replicate) does not significantly improve the efficiency of a random sample of the same size. This happens because the subsample inside the replicate is random. At the same time the traceability of systematic sampling is lost to a large extent.

The stratification defined from CLC provides a more substantial improvement (efficiency=1.40). Combining stratification with systematic sampling roughly collects the advantages of both aspects in terms of efficiency. Independent systematic samples in different strata with different sampling steps in each stratum (figure 1C) gives a good efficiency in spite of the visually unpleasant pairs of points in different strata that are too close to each other. The type of scheme chosen for SIGMA (stratified systematic samplingflat-long balanced, illustrated in figure 3) behaves slightly worse, but still significantly better than stratified random.

## IV SOME PROVISIONAL CONCLUSIONS.

The better performance of systematic sampling compared to random sampling is confirmed even if the relative efficiency is not particularly high. A relative efficiency above 1.5 for systematic sampling is usually not realistic. The rigidity of systematic sampling in terms of sample size can be removed by using a system based on a pattern of replicates. This type of system slightly reduces the efficiency compared to the pure systematic sampling. The usual variance formula for random sampling give a pessimistic bias if applied to systematic sampling. Better estimations are computed on the basis of local variance. The main advantage of systematic sampling is the traceability, i.e. the possibility to proof that no irregularities have been introduced in the sampling process.

Spatial sampling in a geographic reference that is not an equal-area projection, such as latitude-longitude can introduce a significant bias. This can be avoided with different adjustments, both in random and systematic sampling, but such adjustments in systematic sampling can strongly degrade its advantages, both in terms of variance and in terms of traceability.

The conclusions reported in this paper are based on an example of pseudo-truth, even if they are consistent with the conclusions of the wide literature on systematic sampling (see Bellhouse, 1988, for some references). Simulations with a pseudo-truth behaving in a more similar way to the errors of a land cover map would be useful, as would be, of course, the calculation of accuracy indicators on the SIGMA cropland map.

## References

* Bellhouse D.R., (1988). Systematic sampling, *Handbook of Statistics*, vol. 6, ed. P.R. Krisnaiah, C.R. Rao, pp. 125-146, North-Holland, Amsterdam

* Cochran W., (1977) *Sampling Techniques*. New York: John Wiley and Sons

* Dierckx, W., Sterckx, S., Benhadj, I., Livens, S., et al. (2014). PROBA-V mission for global vegetation monitoring: Standard products and image quality. *International Journal of Remote Sensing, 35*(7),

* Dunn R., Harrison A.R. (1993). Two-dimensional systematic sampling of land use. *Journal of the Royal Statistical Society Series C: Applied Statistics*, vol. 42 n. 4, pp. 585-601.

* EEA (2007) CLC2006 Technical guidelines; EEA technical report 17/2007. http://www.eea.europa.eu/publications/technical_report_2007_17

* FAO (1996). *Conducting Agricultural Censuses and Surveys*. FAO Statistical Development Series n. 6. FAO Publication, Rome.
Fritz, S., See, L., Mccallum, I., You, L., Bun, A., Moltchanova, E., Obersteiner, M. (2015). Mapping global cropland and field size. *Global Change Biology*, 21(5), 1980-1992.

Mücher, C.A., de Simone, L., Kramer, H., de Wit, A., Roupioz, L., Hazeu, G., Boogaard, H., Schuiling, R., Fritz, S., Latham, (2016). Global Agro-Environmental Stratification (GAES), SIGMA report D31.1.

Olofsson, P., Foody, G.M., Herold, M., Stehman, S.V., Woodcock, C.E., Wulder, M.A., (2014). Good practices for estimating area and assessing accuracy of land change. *Remote Sensing of Environment* 148, 42 - 57

Stehman, S. V. (2009). Sampling designs for accuracy assessment of land cover. *International Journal of Remote Sensing*, 30(20), 5243-5272

Wagner, J. E., & Stehman, S. V. (2015). Optimizing sample size allocation to strata for estimating area and map accuracy. *Remote Sensing of Environment*, 168, 126-133.

Wolter K.M., (1984). An investigation of some estimators of variance for systematic sampling. *Journal of the American Statistical Association*, Vol. 79 No 388, pp. 781-790.
# Exploring the Usefulness of Land Parcel Data for Evaluating Multi-Temporal Built-up Land Layers

Johannes H. Uhl*<fnref n="1" />, Stefan Leyk<fnref n="1" />, Aneta J. Florczyk<fnref n="2" />, Martino Pesaresi<fnref n="2" />

*Corresponding author: johannes.uhl@colorado.edu

The increasing amount of free and open remote sensing data suggests that a number of multi-temporal land use and built-up land datasets derived from remotely sensed imagery will be made available soon. However, little research has been done regarding the approaches to evaluate spatiotemporal uncertainty of such products. Employing publicly available cadastral data with information on construction date may be useful for this purpose but requires developing a proper validation protocol. In this work we present preliminary results from a feasibility study that examines the potential use of land parcel data as reference data for spatiotemporal evaluation of such built-up land data, exemplified by the novel Global Human Settlement Layer (GHSL). The results indicate that alternative strategies shall be considered, as the use of parcel data tends to bias the evaluation results due to inherent mismatches between the two data sources, especially in rural areas.

## I INTRODUCTION

The Global Human Settlement Layer (GHSL) is a methodology developed to identify and map built-up areas from Landsat satellite imagery (Pesaresi et al. 2013) and create a new global information baseline describing the spatial evolution of the human settlements in the past 40 years (Pesaresi et al. 2016). This global dataset is available at high spatial resolution (38m) and for various periods of time (around 1975, 1990, 2000, and 2014). GHSL data may provide new opportunities for population projections (Freire et al. 2016), disaster management and risk assessment (Freire et al. 2015), as well as for analysing and modelling urban dynamics and land use change.

Before making such novel data products available to the research community, an extensive quality assessment is needed to demonstrate their usability. However, such assessments are difficult and rarely done due to the lack of reliable reference data, particularly for earlier time periods and in less developed regions. We carried out a first experiment to evaluate multi-temporal spatial data on built-up land such as GHSL or developed land cover classes in a typical land cover database using publicly available tax parcel (cadastral) data. How meaningful reference layers based on parcel data can be, represents an open question. This study is meant to shed light on the feasibility of such evaluations in order to establish a thorough protocol for validation studies in the near future.

## II DATA AND METHOD

Open data policy makes cadastral and tax assessment data publicly available – often as GIS-compatible format – for many regions in the U.S. Often these parcel data contain rich attribute

<footnote n="1">¹University of Colorado Boulder, Department of Geography, Boulder, Colorado, USA</footnote>

<footnote n="2">²European Commission – Joint Research Centre (JRC), Institute for the Protection and Security of the Citizen (IPSC), Global Security and Crisis Management Unit, Ispra, Italy</footnote>
information related to the type of land use, characteristics of the structure and the year when a
structure in a parcel has been established (built year). To create spatiotemporal reference layers,
the built year attribute is used to reconstruct snapshots of parcel-derived built-up areas that
correspond to the GHSL time spans of built-up land (before 1975=class 6, 1975-1990=class 5,
1990-2000=class 4, 2000-2014=class 3) and its non-built-up land (class 2).

For each time span, GHSL built-up land is compared with parcel-derived built-up areas. This
has been implemented in a raster-based approach where parcel polygons are rasterized to match
the spatial resolution of GHSL (approx. 38m) using the time-span class as raster value (values
2-6). In order to comply with the GHSL-based built-up land definition, a pixel is classified as
built-up when it overlaps with a built-up entity (Pesaresi et al. 2016). In this first experiment the
parcel area is used as a proxy for built-up land despite well-known limitations (see below). This
reference raster dataset is overlain with the original GHSL data and pixel-based confusion
matrices are built to derive various accuracy metrics (Fielding and Bell 1997) for each time
period. These metrics can be derived to quantify the classification accuracy of GHSL for each
time span (cumulatively e.g., class 4 labels built-up land before 2000 including classes 5 and 6)
and provide rich material for assessing its quality across space and time. Figure 1 shows GHSL
built-up labels (a) and created reference data (b and c) for a subset of Boulder County, CO.

<figure><img src="image_69.png" /><figcaption>Figure 1: GHSL built-up labels and corresponding reference data for North Boulder, Colorado.</figcaption></figure>

Parcel data including built year information is publicly available for different administrative
regions in the U.S. (states, counties, and cities) but not a standard product. In order to test parcel
data for evaluation of GHSL in different settings with regards to development intensity, we
included two rather rural (Figure 2) and two urban counties (Figure 3).
<figure><img src="image_70.png" /></figure>

<figure><img src="image_71.png" /><figcaption>Figure 2: Rural study areas: Baltimore County (excludes the city of Baltimore), Maryland (left), and Boulder County, Colorado (right). Basemap: Open Street Map.</figcaption></figure>

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

<figure><img src="image_73.png" /><figcaption>Figure 3: Urban study areas: City of Madison, Wisconsin (left), and Milwaukee County, Wisconsin (right). Basemap: Open Street Map.</figcaption></figure>

III RESULTS

The results are based on the comparison of GHSL with the parcel-derived reference surface and intended to illustrate the feasibility of using parcel data to evaluate multi-temporal built-up land data. Confusion matrices are created for all time spans in both settings (Figures 4 and 5) and overall agreement metrics are shown in Table 1 for each study area. There are noticeable differences in agreement between rural and urban settings (e.g., Producer's Accuracy or Kappa).

<table><caption>Table 1: Overall accuracy measures for the four study areas.</caption><thead><tr><th></th><th>Average<br>Producer's<br>Accuracy</th><th>Average<br>User's<br>Accuracy</th><th>Kappa</th><th>Normalized<br>Mutual<br>Information</th><th>Overall<br>Agreement</th><th>Average<br>Omission<br>Error</th><th>Average<br>Commission<br>Error</th></tr></thead><tbody><tr><td>Madison (urb.)</td><td>0.497</td><td>0.428</td><td>0.284</td><td>0.153</td><td>0.568</td><td>0.503</td><td>0.572</td></tr><tr><td>Milwaukee (urb.)</td><td>0.369</td><td>0.304</td><td>0.220</td><td>0.095</td><td>0.538</td><td>0.631</td><td>0.696</td></tr><tr><td>Boulder (rur.)</td><td>0.275</td><td>0.636</td><td>0.155</td><td>0.062</td><td>0.673</td><td>0.725</td><td>0.364</td></tr><tr><td>Baltimore (rur.)</td><td>0.284</td><td>0.387</td><td>0.140</td><td>0.035</td><td>0.462</td><td>0.716</td><td>0.613</td></tr></tbody></table>
Confusion matrices for urban areas (Figure 4) show that while there is some variation among counties, in general the agreement for class 2 (non-built-up) and class 6 (built before 1975) is high. Yet, large portions of non-built-up areas (class 2) in parcel data have been classified as built-up before 1975 (class 6) in GHSL.

For rural areas very high agreement can be seen in non-built-up land. However, higher portions of built-up pixels in the parcel surface in all epochs (particularly class 6) have been classified as non-built-up pixels in GHSL (Figure 5) decreasing Producer's Accuracy. A possible explanation could be that portions of large rural residential parcels, which are mostly of agricultural use, are falsely assumed (and thus overestimate) built-up land in the reference data. These portions are often correctly classified as not built-up land (class 2) in GHSL resulting in this kind of disagreement.

<figure><img src="image_74.png" /><figcaption>Figure 4: Overall confusion matrices for the urban study areas Madison, Wisconsin (left) and Milwaukee County, Wisconsin (right).</figcaption></figure>

<figure><img src="image_75.png" /><figcaption>Figure 5: Overall confusion matrices for the rural study areas Boulder County, Colorado, and Baltimore County, Maryland.</figcaption></figure>

Binary confusion matrices are used to derive agreement metrics comparing GHSL and parcel data for non-built-up labels (class 2) and built-up land in each epoch (classes 3-6). These metrics (Figures 6 and 7) also illustrate some interesting differences between rural and urban areas.

User's Accuracy (i.e., a pixel labelled class 4 in GHSL is also labelled class 4 in the parcel reference data) in rural areas increases over time and is higher for the built-up classes when compared with urban areas. Pixels of higher development intensity are more reliably detected as built-up in GHSL, and these locations are most likely within residential parcels in rural settings.

Producer's Accuracy (a pixel of class 4 in parcel data is classified as class 4 in GHSL) increases slightly in both settings over time but is much lower in rural areas due to the overestimated built-
up land in parcel data. Overall agreement (PCC) decreases towards more recent epochs due to an accumulation effect of disagreement in earlier epochs. Interestingly, Kappa, which accounts for chance agreement and is more conservative than PCC, increases over time in rural settings but decreases in urban areas. Possibly, this reflects higher detection rates in rural areas with improved technology, and at the same time a cumulative effect in urban areas similar to PCC.

<figure><img src="image_76.png" /><figcaption>Figure 6: Trends in agreement metrics for urban and rural areas: a) PCC, and b) User's Accuracy.</figcaption></figure>

<figure><img src="image_77.png" /><figcaption>Figure 7: Trends in agreement metrics for urban and rural areas: a) Producer's Accuracy, b) Kappa.</figcaption></figure>

IV DISCUSSION

Measuring the agreement in built-up land labels between publicly available parcel data and GHSL over time allows to investigate the usefulness of parcel data for evaluation of GHSL classes in rural and urban settings. While the built-year attribute provides a unique key feature for spatio-temporal evaluation purposes, the results show that the use of parcel data may bias the evaluation results due to inherent mismatches between the two data sources. Higher degrees of disagreement become evident in rural areas, where the residential parcels can be very large. Thus in order to establish robust protocols for thorough evaluation of GHSL the spatial integration of
building data, which are becoming increasingly available (e.g., LiDAR derived), will be tested to identify built-up areas within parcels more precisely thus creating more reliable reference datasets. Such an evaluation will also need to address uncertainty due to spatial offsets between datasets and the inherent uncertainty in the temporal information (e.g., (i) due to tear-downs and rebuilt structures, not reflected by the built-year attribute, (ii) the exact timestamp of an image and the assigned nominal class). Such an evaluation procedure will develop a broader understanding of the uncertainty across space and time in GHSL built-up area labels and inform the future user community about fundamental quality aspects if GHSL is applied to similar settings in other countries where no reference data are available. In future steps, the study area will be further extended, alternative criteria to differentiate urban vs. rural settings will be examined, and the appropriateness of the accuracy metrics employed will be critically reviewed.

## References

Fielding, A. H., Bell, J. F. (1997). A review of methods for the assessment of prediction errors in conservation presence/absence models. *Environmental Conservation*, 24(01), 38-49.

Freire, S., Florczyk, A., Ehrlich, D. Pesaresi, M. (2015). Remote sensing derived continental high resolution built-up and population geoinformation for crisis management. In *Geoscience and Remote Sensing Symposium (IGARSS), 2015 IEEE International*, 2677-2679.

Freire, S., MacManus, K., Pesaresi, M., Doxsey-Whitfield, E., Mills J., (2016). Development of new open and free multi-temporal global population grids at 250 m resolution. Proceedings of the 19th AGILE Conference on Geographic Information Science. Helsinki, Finland, June 14-17, 2016.

Pesaresi M., Ehrlich D., Ferri S., Florczyk A., Carneiro Freire S. M., Halkia S., Julea A. M., Kemper T., Soille P., Syrris V. (2016). Operating procedure for the production of the Global Human Settlement Layer from Landsat data of the epochs 1975, 1990, 2000, and 2014. JRC Technical Report EUR 27741 EN; doi:10.2788/253582 (online)

Pesaresi, M., Huadong, G., Blaes, X., Ehrlich, D., Ferri, S., Gueguen, L., Halkia, M., Kauffmann, M., Kemper, T., Lu, L., Marin-Herrera, M. A., Ouzounis, G. K., Scavazzon, M., Soille, P., Syrris, V., Zanchetta, L. (2013). A global human settlement layer from optical HR/VHR RS data: concept and first results. In *IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing*, 6(5), 2102-2131.
# Urban objects recognition feasibilities by airborne hyperspectral and multispectral remote sensing

Sébastien Gadal*¹, Walid Ouerghemmi*¹

¹Aix-Marseille University, CNRS ESPACE UMR 7300

*Corresponding authors: sebastien.gadal@univ-amu.fr, walid.ouerg@gmail.com

## Abstract

This paper explores the recognition uncertainty of urban objects by multiband imagery. The purpose is to recognize the urban objects by their spectral signature, using an external spectral library. Two Vis-NIR images were used for the study: a four bands Kompsat-2 multispectral image and a 16 bands Ricola's airborne hyperspectral image, two supervised classifiers were tested; a spectral based classifier, called the Spectral Angle Mapper (SAM), coupled to an external spectral library and a machine learning based classifier called the Support Vector Machine (SVM), in a second step the classification results obtained by the two classifiers were merged, the goal was to take advantage of both techniques, to optimize the classification result. The classifiers performance and the objects recognition feasibility were discussed for both images.

## Keywords

Objects recognition uncertainty, spectral library, multiband imagery, urban environment.

## I. INTRODUCTION

The cities space structure and territories result in a heterogeneous spatial and spectral mosaics of objects, which are combination of different manmade materials, urban vegetation species, and street networks. The recognition and characterization of urban objects becomes a difficult task in remote sensing image processing, requiring the uses of powerful data processing methods (e.g. Baillard et al, 1998). The use of multiband imagery is one of the approaches for identifying objects and urban space structures (e.g. Benediktsson et al, 2005), since it offer a precise characterization of the objects thanks to their spectral signatures. Two methods were selected for this study to deal with the complex context of urban objects recognition and classification. The first method called the Spectral Angle Mapper (SAM) (Kruse et al, 2003) is a simple and powerful classifier which is based on metric calculation, which is the angle formed between a spectral reference and a given pixel, each pixel is then classified depending on the value of the measured metric. Usually, the SAM method is applied using training samples acquired from the image, nevertheless, the sampling process could be complicated due to the presence of mixed pixels or to a lack of ground truth knowledge (e.g., Nidamanuri & Zbell, 2011). For this study, we used an external spectral library of urban materials in addition to the image-based spectral sampling. The use of such libraries in the exploitation of hyperspectral data is a growing and promising approach (Hueni et al, 2009; Zomer et al, 2009; Kotthaus et al, 2014). The second method called the Support Vector Machine (SVM) (e.g. Wu et al, 2004; Hsu et al, 2010) is a powerful classification method based on a machine learning algorithm, the method use training samples acquired from the image, to build optimal surfaces called the hyperplanes, and assign the
image pixels to these hyperplanes. Both classifiers were tested on the multispectral and hyperspectral images. Many studies have tested and compared the efficiency of these methods for urban land cover detection, burned areas detection, agricultural crops recognition, and geology mapping (e.g., J. Murphy et al, 2012; Hegde et al, 2014; Anggraeni & Lin, 2014). However, the fusion between these methods wasn't made. The second part of this work concerns the fusion of the SAM and SVM classification results, to enhance the overall classification accuracy and to take advantage of both methods.

## II. STUDY ZONE AND MATERIALS

The study zone concerns the city of Kaunas (Lithuania) which is characterized by a large variety of urban materials including roofing materials, roofing's painting, roads, pavements, and urban vegetation. For the study we used a 16 bands airborne hyperspectral image acquired over the city centre of Kaunas, the camera includes a Finnish Vis-NIR sensor (RICOLA LTD), with a spectral range of 500-900nm, and a spatial resolution of 1m². In addition to the hyperspectral image we used a 4 bands multispectral KOMPSAT-2 image, with a spectral range of 500-900nm, and a spatial resolution of 16m². Both images were originally of radiance. They were converted to reflectance using a MODTRAN 4 radiative transfer model (Matthew et al, 2000).

To characterize the urban materials of the city, a spectral library was generated over Kaunas in July 2015. The spectral measurements were done in a black room using the superspectral Themis-Vision VNIR400H imaging camera (Fig 1.a). The sensor was developed by the NASA at the John Stennis Space Center. The camera is able to measure the electromagnetic radiation from 400 nm to 1000 nm using 1000 narrow bands; the spectral resolution is equal to 0.6 nm. The collected spectra were converted to relative reflectance values; using a standard white reference panel (99% of the Avian Technologies LLC white reference panel). The relative reflectance is equal to the brightness ratio of the sample grid to the white reference panel grid, multiplied by the standard spectral reflectance values (99%). The library includes thirty common urban materials, including roofing's material, roads, and vegetation types (Fig 1.b).

<figure><img src="image_78.png" /><figcaption>Figure 1: a) Themis Vision superspectral camera, b) Common urban objects in Kaunas city center (Lithuania) measured in laboratory.</figcaption></figure>
In addition to the spectral library, we used a Geographic Information system (GIS) over the same study zone including many information's about the urban materials (i.e. photo, material type, color, date of acquisition, state) collected over Kaunas from December 2015 to March 2016. The GIS was used to validate the classification results and as a decision tool to extract training samples -if needed- for the classifications process.

### III. IMAGE PROCESSING METHODS

#### a. PRE-TREATMENTS

Multiband images and spectral library were both converted to reflectance, concerning the images, an atmospheric correction using the Flaash module (Exelis Visual Information Solutions) were conducted, the correction is based on a modified version of the MODTRAN 4 radiative transfer code. Mid-latitude summer atmospheric model, 820 nm water feature, and no predefined aerosol parameters were used. The aerosol amount in air was estimated by visibility and was set to 40 km.

Concerning the spectral library, the relative reflectance was calculated with the Themis-Vision camera software, using a white reference panel that have a reflection coefficient of 99%, the relative reflectance value for a given pixel is then expressed as the brightness ratio of the sample grid to the white reference panel grid, multiplied by the standard spectral reflectance values (99%) (see Eq 1).

$$ \rho_{\lambda} = \frac{NV_{sample}}{NV_{white reference}} \cdot 99\% \quad (1) $$

#### b. CLASSIFICATION

Two classification methods were used for this study, which are the SAM and SVM, concerning the SAM method, the algorithm is based on a spectral classification scheme, and uses metric calculation. This is a powerful tool that permits to distinguish the classes by their spectral signatures. It decreases the shade influence to enhance the targets detection, and it didn't require any statistics about the distribution of the data. The metric used is the spectral angle $\theta$ between a reference spectrum y and an unknown spectrum x following this equation (see Eq 2):

$$ \theta(x, y) = \arccos\left(\frac{\langle x, y \rangle}{||x||||y||}\right) \quad 0 \le \theta \le \frac{\pi}{2} \quad (2) $$

Where $\langle \dots \rangle$ is the dot product and $||\dots||$ is the 2-norm operator, smaller $\theta$ value indicate better match between x and y.

Concerning the SVM method, the algorithm is based on machine learning theory, below we introduce briefly the problem formulation, introduced by Vapnik (Vapnik, 1995) (see Eqs 3 and 4). Considering a two class problem in n dimensional space $R^n$, let's assume l training samples $x_i \in R^n$ with their corresponding labels $y_i \in \{-1, +1\}$, the method consists in finding the optimal hyperplane (i.e. surface) that separates the classes and maximizes the margin between them (i.e., distance to the closest training data points for both classes). The closest training data points to the hyperplane are called support vectors. We used the Radial Basis Function as a Kernel, the penalty parameter that controls the margin rigidity and misclassification was set to 100, no pyramidal
levels were fixed. The image was processed at full resolution without rescaling. The hyperplane was derived following Lagrangian problem below:

$$ \text{maximize } f(\alpha) = \sum_{i=1}^{l} \alpha_i - \frac{1}{2} \sum_{i,j=1}^{l} \alpha_i \alpha_j y_i y_j K(x_i, x_j) \quad (3) $$

With the constraint: $0 \le \alpha_i \le C \quad \forall i \in [1, l]$ and $\sum_{i=1}^{l} \alpha_i y_i = 0$, the normal vector of the hyperplane which give the maximal margin is given by:

$$ w = \sum_{i=1}^{l} \alpha_i y_i x_i \quad (4) $$

Where C is the penalty parameter, K is the Kernel function, and $\alpha_i$ is non-zero vectors called the support vectors.

### c. METHODS IMPLEMENTATION

In this study we present two classification strategies: a simple classification based on the separate use of SAM and SVM, and a classification based on the merging of SAM and SVM methods (Fig 2). The first step concerns the image pre-treatments. The radiance images are converted to reflectance using the MODTRAN 4 model; the multiband images generated by the laboratory camera Themis Vision were also converted to reflectance using the camera software (Hyper Visual Imaging Systems, Themis Vision Systems). The second step concerns the application of SAM and SVM methods. The SAM is coupled to an external spectral library in addition to an image-based library. The SVM is coupled to GIS database for a precise acquisition of training samples. The classification is followed by an artifact elimination to enhance the results interpretation and to delete some misclassified pixels. The last step concerns an original fusion strategy between the SAM and SVM which is guided by statistical accuracy indicators and ground truth knowledge.

<figure><img src="image_79.png" /><figcaption>Figure 2: Method implementation</figcaption></figure>
# IV. URBAN MATERIALS CLASSIFICATION

## a. CLASSIFICATION OF HYPERSPECTRAL RICOLA IMAGE

Concerning the SAM method, the results (Fig 3.a) reveal good identification of tile and red painted metal roofing's. The recognition is carried using the spectral library, nevertheless old tile with altered color was hardly recognizable by SAM, the incorporation of spectral samples of old tile to our library should enhance the detection of tile material. White painted metal and dark painted roofing's are recognized by spectral training samples extracted from the image; they were almost well recognized both. The trees are well recognized but sometimes confounded with grass. The pavements are weakly detected due to their high correlation with asphalt. Asbestos roofing's are not detectable due to their high correlation with roads. The results are evaluated by the Kappa coefficient and overall accuracy, which are respectively of 0.47 and 54%. The classification performance is moderately accurate due to the weak detection of some materials (i.e. vegetation and pavements mainly).

Concerning the SVM method, the results (Fig 3.b) reveal a good identification of both recent and old tile, nevertheless, an overestimation is observed, and the red steel roofing's are often misclassified with tile. The recognition is altered compared to the SAM results for this material. White steel roofing's, dark steel roofing's, and pavement are better recognized with SVM, the trees (coniferous) and grass are better separable, deciduous trees couldn't be detected due to the lack of dedicated validation points. Theoretically the SVM performs better than SAM, with a Kappa coefficient and an overall accuracy of respectively 0.79 and 82%. Nevertheless the algorithm tends to overestimate the pixels compared to the SAM. In the next step, the results of SAM and SVM will be merged to take advantage of both techniques.

<figure><img src="image_80.png" /><figcaption>Figure 3: Urban materials classification using a) SAM, b) SVM and c) SAM and SVM fusion.</figcaption></figure>

The merging of the SAM and SVM (Fig 3.c) is done following the observations noticed previously. The SVM tends to overestimate the pixels classes for materials that are not
homogeneous in color or not highly reflective, which is the case of tile and red steel roofing's (Fig 3.b). These materials were detected with SAM method (Fig 3.a). The SVM estimates correctly the highly reflective materials (e.g. white roofing's, pavements) and the materials which are homogeneous in color (e.g. vegetation, some of the dark roofing's) (Fig 3.b). The white steel roofing's, the dark steel roofing's, the pavements and the vegetation were classified by SVM. The classification of coniferous trees was only possible with SAM, thanks to the spectral library (Fig. 3.a). The combination of SAM and SVM results (Fig 3.c) offered and interesting and precise pattern of urban materials over Kaunas city, with reasonable performances in terms of estimation indicators; the Kappa coefficient and overall accuracy were respectively of 0.71 and 75%.

## b. CLASSIFICATION OF MULTISPECTRAL KOMPSAT-2 IMAGE

Concerning the SAM method, the results (Fig 4.a) reveal a correct identification of tile and red painted metal roofing's. The recognition is carried using the spectral library, an over estimation of the clay is noticeable. The detection of white steel roofing's is also possible. The detection of dark steel roofing's and asbestos roofing's is not possible using this 4 bands image. The pavements are partially detectable. Trees are detected partially and sometimes confounded with grass. The results are evaluated by the Kappa coefficient and overall accuracy, which are respectively of 0.48 and 57%, the classification performance is moderately accurate.

Concerning the SVM method, the results (Fig. 4.b) reveal a good identification of tile with moderate correlation with other materials. The red steel is weakly detected, in the other hand; the white steel roofing's and pavements are almost correctly detected, as for the SAM method, dark steel roofing's and asbestos roofing's couldn't be detected with this image. Concerning the vegetation; grass and trees are well distinguished, with an overestimation that could be noticed for the grass in comparison to the SAM method. Theoretically, the SVM method performs better than SAM, with a Kappa coefficient and an overall accuracy of respectively 0.8 and 84%. In the next step the SAM and SVM methods will be merged together.

<figure><img src="image_81.png" /><figcaption>Figure 4: Urban materials classification using a) SAM, b) SVM and c) SAM and SVM fusion.</figcaption></figure>
The SVM shows a better identification of tile, pavements and vegetation than SAM (Fig 4.b). The red steel roofing's is better detected with SAM (Fig 4.a), the white steel roofing's detection is correctly done with both methods. So both class of white steel roofing's are merged. The fusion of SAM and SVM (Figure 4.c) seems homogeneous in terms of class repartition with an accurate theoretical accuracy; the Kappa coefficient and overall accuracy were respectively of 0.81 and 85%.

## V. CONCLUSION AND PERSPECTIVES

This study explored urban objects recognition using multiband imagery with two classification methods tested: the SAM and the SVM. The results obtained by both methods are merged to enhance the classification accuracy in a second step of processing. The SAM and SVM showed encouraging results in terms of detection with Ricola's 16 bands hyperp spectral image. Concerning the SAM, the materials were precisely recognized, nevertheless the external library couldn't cover all urban materials at the moment, therefore, image-based spectra have been used in addition to the spectral library to initialize the SAM. The overall accuracy of SAM was moderate due to non-sufficient detection of vegetation and pavements, the method accuracy is also sensible to the maximum angle affected to each class, which is tunable. Concerning the SVM, the method requires longer processing time than SAM, nevertheless it presented a better accuracy in terms of overall accuracy. The method inclines to overestimate certain classes, especially tile and vegetation. The results must be carefully examined and compared to ground truth knowledge before validation. In a second step we merged the SAM and SVM results. This fusion presents an encouraging opportunity for hyperspectral classification, and permits to overcome the limitations of both methods, while preserving an acceptable overall accuracy.

In addition to the hyperspectral image, we have tested our methods on a 4 bands multispectral image. The SVM seemed more suitable and powerful to process multispectral images than SAM. The moderate spectral resolution of the image is not suited for metric based classification. The fusion of SAM and SVM enhanced however -slightly- the overall accuracy and enhanced the identification of some materials (e.g. red steel roofing's). The enrichment of the spectral library will probably enhance the SAM classification performances.

## Acknowledgments

This research was supported by the French National Research Agency (ANR) through the Hyep project (ANR-14-CE22-0016). We are thankful to Dr. Gintautas Mozgeris, Aleksandras Stulginskis University (Lithuania), for providing us with the hyperspectral image of 2015 and for allow us to perform the laboratory measurements with the Themis-Vision Camera.

## References

Anggraeni, A., and Lin C. (2011). Application of SAM and SVM Techniques to Burned Area Detection for Landsat TM Images in Forests of South Sumatra. *International Conference on Environmental Science and Technology*, Singapore: pp. 160-164.

Baillard, C., Dissard, O., Jamet, O., and Maître, H. (1998). Extraction and textural characterization of above-ground areas from aerial stereo pairs: a quality assessment. *ISPRS Journal of Photogrammetry and Remote Sensing* 53, 130-141.

Benediktsson, J.A., Palmason, J.A., and Sveinsson, J.R. (2005). Classification of hyperspectral data from urban areas based on extended morphological profiles. In *IEEE Transactions on Geoscience and Remote Sensing*, vol. 43, no. 3, pp. 480-491.
Hegde, G., Ahamed, J. M., Hebbar, R., and Raj, U. (2014). Urban land cover classification using hyperspectral data. *ISPRS-International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences*, Volume XL-8, 751-754.

Hsu, C.-W., Chang, C.-C., and Lin, C.-J. (2010). *A practical guide to support vector classification*. National Taiwan University. http://ntu.csie.org/~cjlin/papers/guide/guide.pdf.

Hueni, A., Nieke, J., Schopfer, J., Kneubühler, M., Itten, K.I. (2009). The spectral database SPECCHIO for improved long-term usability and data sharing. *Computers & Geosciences* 35, 557-565.

Kotthaus, S., Smith, T.E.L., Wooster, M.J., Grimmond, C.S.B. (2014). Derivation of an urban materials spectral library through emittance and reflectance spectroscopy. *ISPRS Journal of Photogrammetry and Remote Sensing* 94, 194-212.

Kruse, F. A., Lefkoff, A. B., Boardman, J. B., Heidebrecht, K. B., Shapiro, A. T., Barloon, P. J., and Goetz, A. F. H. (1993). The Spectral Image Processing System (SIPS) - Interactive Visualization and Analysis of Imaging spectrometer Data. *Remote Sensing of Environment* 44, 145-163.

Matthew, M. W., Adler-Golden, S. M., Berk, A., Richtsmeier, S. C., Levine, R. Y., Bernstein, L. S., Acharya, P. K., Anderson, G. P., Felde, G. W., Hoke, M. P., Ratkowski, A., Burke, H.-H., Kaiser, R. D., and Miller, D. P. (2000). Status of Atmospheric Correction Using a MODTRAN4-based Algorithm. *SPIE Proceedings, Algorithms for Multispectral, Hyperspectral, and Ultraspectral Imagery VI*. Vol. 4049, pp. 199-207.

Murphy, R. J., Monteiro, S. T., and Schneider, S. (2012). Evaluating Classification Techniques for Mapping Vertical Geology Using Field-Based Hyperspectral Sensors. In *IEEE Transactions on Geoscience and Remote Sensing*, vol. 50, no. 8, pp. 3066-3080.

Nidamanuri, R. R., & Zbell, B. (2011). Transferring spectral libraries of canopy reflectance for crop classification using hyperspectral remote sensing. *Biosystem Engineering* 110(3), 231-246.

Plaza, A., Benediktsson, J.A., Boardman, J.W., Brazile, J., Bruzzone, L., Camps-Valls, G., Chanussot, J., Fauvel, M., Gamba, P., Gualtieri, A., et al. (2009). Recent advances in techniques for hyperspectral image processing. *Remote Sensing of Environment* 113, Supplement 1, S110-S122.

Samsudin, S. H., Shafri, H. Z. M., Hamedianfar, A., and Mansor, S. (2014). Spectral feature selection and classification of roofing materials using field spectroscopy data. *Journal of Applied Remote Sensing* 9(1), 095079.

Villa, A., Chanussot, J., Benediktsson, J.A., and Jutten, C. (2011). Spectral Unmixing for the Classification of Hyperspectral Images at a Finer Spatial Resolution. In *IEEE Journal of Selected Topics in Signal Processing*, vol. 5, no. 3, pp. 521-533.

Vapnik, V. (1995). *The Nature of Statistical Learning Theory*, New York: Springer-Verlag.

Wu, T.-F., C.-J. Lin, and R. C. Weng. (2004). Probability estimates for multi-class classification by pairwise coupling. *Journal of Machine Learning Research* 5, 975-1005.

Zomer, R.J., Trabucco, A., Ustin, S.L. (2009). Building spectral libraries for wetlands land cover classification and hyperspectral remote sensing. *Journal of Environmental Management* 90, 2170-2177.
International Conference on Spatial Accuracy (2016):
Space-Time Kriging of Temperature over Van-Turkey

Pınar ASLANTAS¹, Okan YELER²

¹Yuzuncu Yil University, Faculty of Agriculture, Landscape Architecture Department,
TURKEY

²Yuzuncu Yil University, Muradiye Vocational School, TURKEY

*Corresponding author: pinaraslantas@yyu.edu.tr

Predictions of the variables at points with no measurements are obtained by interpolation techniques. Space-time interpolation techniques that consider variation both in space and time provided a new research area. Temperature is an important climatic parameter varying both in space and time. Like other meteorological, hydrologic and environmental variables, temperature is measured at specific locations. In order to obtain predictions for all grid locations, kriging methods have been applied for a long time. In here, space-time Ordinary kriging (ST-OK) and space-time Universal kriging (ST-UK) have been used in annual temperature estimation over Lake Van Basin using at 13 meteorological station observations for 2001-2011. Elevation, land cover, distance to Van Lake are used as secondary information in ST-UK. Elevation at 500 m resolution are obtained by Nearest Neighbour resampling of 3 arc second Shuttle Radar Topography Mission (SRTM). MOD12Q1 land cover data set has been downloaded from USGS organization. This dataset has 500 m spatial resolution and 17 subclasses. Distance to nearest coast variable is obtained by calculating the Euclidean distances of each SRTM pixel to the nearest boundary of the Van Lake coast vector. Annual temperature values are analysed and predicted at 500m*500 m resolution for 11 year-period. One-fold cross-validation is used to assess accuracy performance of both methods. R-square and Root Mean Square Error (RMSE) are calculated and evaluated for each technique. Comparison of kriging methods and inclusion of secondary information is assessed.

Key words: Space-time kriging, Van Lake Basin, Temperature

**Introduction**

Spatial kriging methods have been used for many years to predict variables at unmeasured locations in many disciplines. The first geostatistics and spatial kriging applications started in mining and geology. The variables used in these sciences can often be assumed constant in time. After understanding the usefulness and reliability of kriging in these disciplines, it was also introduced to many other disciplines within the earth and environmental sciences, such as meteorology, climatology, agronomy, soil science, hydrology, etc. Generally variables in these sciences vary both in time and space. Therefore the requirement of kriging methods for space-time interpolation is raised (Heuvelink and Griffith, 2010). If the data have been measured in different time and space locations, then more data may be used for prediction, and this allows obtaining more accurate predictions, helps to parameter estimation and helps to define spatial and/or temporal auto-correlation in measurements (Gething et al., 2007). In case of space-time kriging, to predict the value of the variable of interest at a specific location and time, past and
future measurements are used to predict on the specified time. This may add more complexity
to the kriging procedure but may help to gain more accurate results.

In this study space-time Universal kriging (ST-UK) method is applied to average annual temperature values measured from 9 meteorological stations from 2001 to 2011 over the Lake Van basin of Turkey. The aim is to discuss applicability of space-time kriging methods on annual temperature values by using limited number of meteorological station.

Study Area and Data

The study area is Lake Van Basin that is located at the Far East part of Turkey (Figure 1). The area of basin is about 16.000 km². Lake Van basin has a high topography. The high mountains are located at the northern and southern parts of the basin. The mean elevation of basin is about 2200-2400 m., minimum elevation is about 1500 m. and maximum elevation is approximately 4000 m (Figure 2). Lake Van which is the biggest lake of country is located at the basin (Figure 2). The Lake is a depression state in the middle of high mountains. Lake has a surface of 3574 km², length of shoreline is 505 km, and a volume of 607 km³. The lake stands at 1650 m. above sea level. The Lake is a closed lake without any significant outflow. With a maximum depth of 451 m and a volume of 607 km³, it ranks fourth in water content among all the closed lakes of the world (Degens et al., 1984).

<figure><img src="image_82.png" /><figcaption>Figure 1: Location of Lake Van Basin on Turkey</figcaption></figure>
<figure><img src="image_83.png" /><figcaption>Figure 2: Location of Lake Van Basin on Turkey Lake Van, SRTM (90m) of basin, and distribution of meteorological stations over basin.</figcaption></figure>

## Data

The temperature data used in this study were obtained from the Turkish State Meteorological Service. The primary dependent data source was monthly temperature measured at 9 meteorological stations between 2001 and 2011. The spatial distribution of stations is not fairly uniform over the basin; when looking the overall distribution condensed placement can be seen near the Lake boundary (Figure 2). The highest average monthly temperature is between 24-25 °C and is measured generally at the June. The lowest average monthly temperature is about -17 °C and is measured at the January and Far East part of basin. It is observed in many studies that secondary information can often improve the spatial interpolation of environmental variables (Bostan et al., 2012; Lloyd, 2005; Hofierka, 2002; Boer, 2001). Therefore, three secondary variables correlated with temperature are selected to improve accuracy of predictions. As independent data source, an elevation map with 500 m spatial resolution was used (Figure 2). It was obtained by resampling the 3 arc second SRTM (the Shuttle Radar Topography Mission) (approximately 90 m spatial resolution) to 500 m spatial resolution. MOD12Q1 land cover data set has been downloaded from USGS organization. This dataset has 500 m spatial resolution and 17 subclasses. Distance to nearest coast variable is obtained by calculating the Euclidean distances of each SRTM pixel to the nearest boundary of the Van Lake coast vector.

## Methodology

Annual temperature predictions are made on spatio-temporal framework. Space-time universal kriging (ST-UK) method is used to obtain predictions over the basin. Elevation, distance to nearest Lake Coast and land cover data sets are used as secondary variable. Separable, productSum, Metric, sumMetric and simpleSumMetric variogram models are applied to sample variogram (Figure 3 and 4).

## Space-time Kriging

Consider a variable z which varies in the spatial (s) and time (t) domain. Let z be observed at n space-time points (sᵢ, tᵢ), i=1, ..., n. These measurements constitute a space-time network of observations. However it is practically impossible to measure data point z at each spatial and temporal point. In order to obtain a complete space-time coverage, interpolation of z is required. The aim of space-time interpolation is to predict z(s₀, t₀) at an unmeasured point (s₀,
$t_0$), which is a node of a space-time grid. To predict $z$ at these nodes, it is assumed to be a realization of a random function $Z$ which has a known space-time dependence structure. Next $Z(s_0, t_0)$ is predicted from the observations and using the assumed space-time model (Heuvelink and Griffith, 2010).

The random function $Z$ can be defined with a deterministic trend $m$ and a zero-mean stochastic residual $V$ as follows (1):

$$Z(s,t) = m(s,t) + V(s,t) \qquad (1)$$

The deterministic trend $m$ represents large-scale variations whereas the stochastic component $V$ represents small-scale variations (Heuvelink and Griffith, 2010).

## Results and Discussion

Space-time Universal kriging is performed to dependent variable: annual temperature. Figure 3 and 4 represent the plots of sample and modelled variograms to data. Optimum parameter value of the modelled variograms are used to select the best variogram model (Table 1). In addition, visual interpretation of modelled variograms are considered while selecting the appropriate model. According to these values productSum and sumMetric models resulted minimum and very similar optimum values. Therefore productSum model is selected for kriging operation (Figure 5).

<figure><img src="image_84.png" /><figcaption>Figure 3: Plot of sample and modelled variogram models</figcaption></figure>
<figure><img src="image_85.png" /><figcaption>Figure 4: Wireframes of sample and modelled variogram models</figcaption></figure>

<table><caption>Table 1: Optimum Parameter Values of Variograms</caption><thead><tr><td>Separable</td><td>ProductSum</td><td>Metric</td><td>SumMetric</td><td>SimpleSumMetric</td></tr></thead><tbody><tr><td>0,62</td><td>0,37</td><td>0,58</td><td>0,36</td><td>0,88</td></tr></tbody></table>

Prediction maps of some selected years are represented at Figure 5. According to these maps, it is clearly seen that around the lake, temperature values are higher than the other regions. Around the Lake the area is more flat, after some distances elevation gets higher. So it can be said that distance to lake and elevation are highly correlated with temperature variable.
<figure><img src="image_86.png" /><figcaption>Figure 5: Prediction maps for some observation years</figcaption></figure>

## Conclusions

In this study space-time kriging method was applied to predict annual temperature of the Lake Van Basin, Turkey. Measurements obtained from nine meteorological stations were used for 2001-2011 time period. Secondary variables that vary in space but are static in time (elevation and distance to lake) and variable changes in space and time (land cover) were used by the space-time Universal kriging method. ST-UK method resulted with reasonable prediction values at space and time; however prediction values for each time scale are similar to each other. Limited number of observations were used at space-time kriging and gives satisfactory results with regard to spatial and temporal framework. However it is thought that using more dense observations would give more accurate results.

## References

* Heuvelink, G. B. M., Griffith, D. A., “Space-Time Geostatistics for Geography: A Case Study of Radiation Monitoring Across Parts of Germany”, *Geographical Analysis* 42 (2010) 161–179.

* Gething, P.W., Atkinson, P.M., Noor, A.M., Gikandi, P.W., Hay, S.I., Nixon, M.S., “A local space-time kriging approach applied to a national outpatient malaria data set”, *Computers & Geoscienses* 33, (2007), pp 1337–1350.

* Degens, E.T., Wong, H.K., Kempe, S., Kurtman, F., “A Geological study of Lake Van, Easten Turkey”, *Geologische Rundschau* 73, 2, pp 701–734, (1984).

* Bostan, P.A., Heuvelink, G.B.M., Akyurek, S.Z., “Comparison of Regression and Kriging Techniques for Mapping the Average Annual Precipitation of Turkey”, *International Journal of Applied Earth Observation and Geoinformation* 19, pp 115–126, 2012.

* Lloyd, C.D., ‘Assessing the effect of integrating elevation data into the estimation of monthly precipitation in Great Britain’, *Journal of Hydrology* 308, pp 128–150, 2005.

* Hofierka, J., Parajka, J., Mitasova, H., Mitas, L., ‘Multivariate interpolation of precipitation using regularized spline with tension’, *Transactions in GIS* 6 (2), pp 135–150, 2002.

* Boer, E. P. J., Beurs, K. M., Hartkamp, A. D., ‘Kriging and thin plate splines for mapping climate variables’, *International Journal of Applied Earth Observation and Geoinformation* 3 (2), pp 146–154, 2001.
<figure><img src="image_87.png" /></figure>

# Spatial inference under uncertainty
<figure><img src="image_88.png" /></figure>
# Where (we think) the wild things are: comparing citizen generated and formal measures of wilderness

Alexis Comber¹*, Steve Carver¹

¹School of Geography, University of Leeds, UK

*Corresponding author: a.comber@leeds.ac.uk

## Abstract

This paper compares mapped data describing the degree of wilderness from a formal data analysis – the Wilderness Quality Index (WQI) – against crowdsourced estimates of wilderness collected by the Geo-Wiki initiative. The analysis examines the quality of the crowdsourced wilderness reporting, and how that quality varies spatially, using a geographically weighted model. The results indicate globally a positive relationship between the 2 datasets, but with wilderness being under-estimated by the crowdsourced data. However, when local patterns are examined it is evident that in more wild areas (the Alps, eastern Norway, etc) the crowdsourced data have a stronger relationship with the formal data. This suggests that, while citizens may be good at identifying locales that at the extremes of a continuum of wilderness (*wild* and *not-wild*), they may not be so reliable in describing the wilderness in areas with a wilderness character between these extremes.

## Keywords

Geographically weighted model; citizen science; wildness; Geo-Wiki

## I BACKGROUND

This paper compares citizen measures of human impact from Geo-wiki with formal wilderness data to determine the degree to which citizens perceive wilderness. Geo-Wiki contributors record the degree of human impact [0-100%] and the land cover they perceive to be present. Different approaches have been used to map wilderness and activities at continental / global scales use a continuum / environmental modification spectrum (Kuiters et al. (2011). Other approaches such as the 2012 Google data that applied a 10km buffer to its road dataset under the assumption that roadless areas represent natural ecosystems that are undisturbed. The aim of this work was to compare two different datasets describing wilderness and wild land, one collected under the Geo-Wiki initiative (Perger et al., 2012; Fritz et al., 2012) and the other as part of a pan-European wildness project described in Kuiters et al. (2011).

Geo-Wiki uses Google Earth imagery and asks contributors to label features at randomly selected points. There have been a number of campaigns and interface developments. While much work has examined the utility and quality of the land cover information that can be derived from Geo-Wiki land cover data (See et al 2013; Foody et al 2014, Comber et al 2013), where users had to allocate scenes in Google Earth imagery to one of a predefined number of classes, as yet little work has considered the uncertainties, accuracies and errors associated with user contributed measures within a continuum, such as wildness. Specifically this work seeks to determine whether the effort involved in developing traditional maps of wilderness (as a suitability layer) results in significantly higher quality data products. Or whether, more
informal approaches, such as are facilitated by crowdsourced data, are able to support the development of data of quality to be scientifically useful.

## II DATA & ANALYSIS

Figure 1 shows ~13,000 Geo-Wiki data points and formal data the wilderness quality index (WQI, Kuiters et al. (2011). It shows that, broadly, both datasets are describing similar distributions of wildness. Iceland, the Alps and North western Scotland has all have high wildness values in these areas.

The WQI data reports wildness measure in the range [0, 1] at 1km scale. The Geo-Wiki data has a number of attributes including land cover class (from a pre-defined set of 10 classes), a Human Impact score which was treated as the inverse of a wildness measure and subtracted from 1, and a self-reported measure of Geo-Wiki volunteer confidence for each point that was labelled.

The analysis used a standard OLS regression to model the relationship between Geo-Wiki volunteer measures of Human Impact, under the assumption that this describes the inverse of wildness, with the formal measure of wilderness quality index. For each Geo-Wiki data point, the coincident wilderness measure was extracted from the WQI data. The model describing the degree to which Geo-Wiki Human Impact predicts WQI data was constructed using the following:

$$y = \beta + \beta_1 x_1 \quad (\text{Eqn 1})$$

where y is the WQI at each data point and $x_1$ is the Geo-Wiki wildness value derived from the Human Impact value, and $\beta_1$ is the coefficient estimate for Geo-Wiki wildness. This identified constructed a global model of the relationship between these 2 variables.

<figure><img src="image_89.png" /><figcaption>Figure 1: a) the Wildness Quality Index data, darker, bluer areas are more wild, and b) the wildness measures from the Geo-Wiki platform, with a small transparency term</figcaption></figure>

Next a Geographically Weighted model was constructed to examine whether and how the coefficient estimates vary. The idea here is that global statistical models implicitly make an
assumption that the relationship between the dependent and independent variables is constant over the study area. In reality this is frequently not the case. In such situations, where spatial non-stationarity is thought to occur, local models can provide a more detailed and informative analysis. The Geographically Weighted (GW) framework (Brunsdon et al., 1996) supports a suite of local models in order to account for local variation in statistical relationships between variables (Gollini et al, 2015). A GW regression (GWR) was used to examine how the coefficient estimates of the degree to which Geo-Wiki measure predict WQI measures from of a standard OLS regression vary spatially. A GWR version of the OLS (Equation 1) is defined as follows

$$y(u_i, v_i) = \beta_0(u_i, v_i) + \beta_1 x_1(u_i, v_i) \quad (\text{Eqn 2})$$

where $(u_i, v_i)$ is a vector of two dimensional co-ordinates describing the location of i over which the coefficient estimates are assumed to vary. The outputs provide spatially varying estimates of the coefficients and the aim of generating these was to better understand how informal, citizen generated measures of landscape wildness relate to formal measures, how these vary spatially and how measures of wildness interact with the self-reported confidence value attached by volunteers to their labelling.

### III RESULTS

Figure 2 shows the Geo-Wiki wildness scores (inverted from the Human Impact value) against the WQI data. The OLS model suggested a coefficient estimate $\beta_1$ of 0.356 for the Geo-Wiki wildness value, indicating that globally, each increase of 0.1 in the WQI data is associated with an increase of 0.036 in the Geo-Wiki data. The results of the all the regressions are summarised in Table 1. The distribution of the GWR coefficient estimates indicates that there is considerable local variation from the global model.

<figure><img src="image_90.png" /><figcaption>Figure 2: WQI against rescaled Geo-wiki wildness scores, with the slope of regression coefficient.</figcaption></figure>

<table><thead><tr><th>OLS</th><th>Estimate</th><th>Std. Error</th><th>t value</th><th>Pr(&gt;|t|)</th><th>Sig</th></tr></thead><tbody><tr><td>(Intercept)</td><td>0.3400</td><td>0.0031</td><td>108.5000</td><td>&lt;2e-16</td><td>***</td></tr><tr><td>Geo-Wiki wildness</td><td>0.3563</td><td>0.0051</td><td>69.9700</td><td>&lt;2e-16</td><td>***</td></tr></tbody></table>
<table><caption>Table 1. The results of the OLS regression and the GW regression</caption><thead><tr><th>GWR</th><th>Min.</th><th>1st Qu.</th><th>Median</th><th>3rd Qu.</th><th>Max</th></tr></thead><tbody><tr><td>(Intercept)</td><td>0.1602</td><td>0.3190</td><td>0.3797</td><td>0.4718</td><td>0.7934</td></tr><tr><td>Geo-Wiki wildness</td><td>-0.0643</td><td>0.2046</td><td>0.2813</td><td>0.3385</td><td>0.5746</td></tr></tbody></table>

Next the coefficient estimates were mapped to show how the coefficient estimates linking wildness score from Geo-Wiki volunteers relate to formal WQI wildness index scores. This is shown in Figure 3 and indicates distinct areas and gradients where volunteers found it difficult to accurately assign wilderness values. Clearly there is high agreement in the Alps and in eastern Norway but large differences between formal measures and volunteered ones in for example eastern France. Under the assumption that the formal WQI data are correct, this suggests that there are locales where people's perceptions of wilderness differ hugely from the way that it is formally recorded.

This paper presents the first analysis of cognitive aspects of how citizens perceive wilderness. Future work will seek to unpick the origins of these differences in relation to landscape context (for example the land cover and land use in these areas), user self-assigned confidence in their scores and their interaction. It will also consider the use of Geo-Wiki data to validate other datasets as has been done for many global land cover initiatives, and how to objectively quantify the accuracy of either dataset.

<figure><img src="image_91.png" /><figcaption>Figure 3: The spatial variation of degree to which informal Geo-Wiki measures of wildness predict formal measures from the WQI dataset.</figcaption></figure>

## References

* Brunsdon, C.F., Fotheringham, A.S. and Charlton M. (1996). Geographically Weighted Regression - A Method for Exploring Spatial Non-Stationarity, *Geogr Anal.*, 28, 281-298.
Comber, A., See, L., Fritz, S., Van der Velde, M., Perger, C., Foody, G.M. (2013). Using control data to determine the reliability of volunteered geographic information about land cover. *International Journal of Applied Earth Observation and Geoinformation*, 23: 37–48 http://dx.doi.org/10.1016/j.jag.2012.11.002

Foody, G.M., See, L., Fritz, S., Van der Velde, M., Perger, C., Schill, C., Boyd, D.S. and Comber, A., (2014). Accurate attribute mapping from volunteered geographic information: issues of volunteer quantity and quality. *The Cartographic Journal* DOI: http://dx.doi.org/10.1179/1743277413Y.0000000070

Fritz S, McCallum I, Schill C, Perger C, See L, Schepaschenko D, van der Velde M, Kraxner F and Obersteiner M 2012 Geo-Wiki: An online platform for improving global land cover. *Environmental Modelling and Software* 31: 110-123.

Kuiters, A. T., van Eupen, M., Carver, S., Fisher, M., Kun, Z., & Vancura, V. (2011). *Wilderness register and indicator for Europe*. Final Report - http://ec.europa.eu/environment/nature/natura2000/wilderness/pdf/Wilderness_registerzmz indicator.pdf

Perger C, Fritz S, See L, Schill C, Van der Velde M, McCallum I and Obersteiner M 2012 A campaign to collect volunteered geographic Information on land cover and human impact. In: Jekel T, Car A, Strobl J and Griesebner G (Eds.) *GI Forun 2012: Geovizualisation, Society and Learning*. Herbert Wichmann Verlag, VDE VERLAG GMBH, Berlin/Offenbach, pp.83-91.

See, L., Comber, A.J., Salk, C., Fritz, S., Van der Velde, M., Perger, C., Schill, C., McCallum, I., Kraxner, F. and Obersteiner M. (2013). Comparing the Quality of Crowdsourced Data Contributed by Expert and Non-Experts. *PLoS ONE* 8(7): e69958. doi:10.1371/journal.pone.0069958
# Uncertain clustering of social specialization in metropolitan areas

Giovanni Fusco¹*, Cristina Cao¹

¹ UMR ESPACE, CNRS / University of Nice Sophia-Antipolis, France

*Corresponding author: giovanni.fusco@unice.fr

**Abstract:** Instead of defining a priori target populations, sociodemographic variables and sector of residence are combined in order to identify geographically meaningful clusters of households on the French Riviera. A Bayesian classifier produces uncertainty-based clusters whose uncertain knowledge is represented through an interactive geo-dataviz solution. Cluster characteristics, sociodemographic content of geographic sectors, socio-spatial contrasts between neighboring sectors and proximity in variable space are explored taking into consideration the uncertain character of the available knowledge.

**Keywords:** Social Specialization, Uncertainty, Clustering, Bayesian Networks, Visualization, French Riviera

## I INTRODUCTION

Social specialization of residential space within an urban area is the concentration of households according to some characteristics like social status, demography, ethnicity, etc. in different urban subspaces. The metropolization process is often associated to increased social specialization in large urban areas (Lacour and Gaschet 2008). Taken to the extreme, social specialization of residential space can produce residential segregation, creating socioeconomic dividing lines within the metropolitan space and undermining social and territorial cohesion (Massey 1985). In France, for example, important policies are carried out in order to assure that municipalities within metropolitan areas have a minimum social mix in terms of household income (Blanc 2010). Nevertheless, knowledge of social specialization of space in metropolitan areas is still incomplete and general assumptions have to be confronted to empirical data in case studies. Understanding the logics of social specialization is of course crucial in order to define policies like urban planning or housing. More specifically, which groups of households are most opposed in residential space? And which socio-demographic factors contribute most to these oppositions in space? How segregated are they in space? If we recognize that the answers to these questions are affected by high levels of uncertainty, which uncertainty-based methods could be used in order to describe our uncertain knowledge of social specialization? Finally, how could we represent and communicate most effectively this uncertain geographic knowledge?

In order to answer these questions, a new clustering approach is proposed integrating both sociodemographic descriptors of households and geographic distribution of their place of residence. The case study of the analysis is the metropolitan area of the French Riviera, a coastal conurbation of more than 1 million inhabitants stretching over 60 km west of the French-Italian border, and including the coastal cities of Nice (348 000 inhabitants in 2013),
Antibes (76 000 inhabitants), Cannes (74 000 inhabitants) and the Principality of Monaco (38 000 inhabitants), which is an independent city-state within French territory. These cities form nowadays an urban continuum and extend their influence over their alpine hinterland, which absorbs a growing part of the metropolitan population. Being a traditional residential destination for affluent retirees, a more recent hub of high-tech development and the necessary home for large populations of low-skilled workers of the tourist and residential economy, the French Riviera is particularly concerned by social specialization of residential space (Centi 1993, Billard and Madoré 2009, Fusco and Scarella 2011).

## II METHODOLOGY: BAYESIAN CLUSTERING WITH SPATIAL CONSTRAINT AND INTERACTIVE GEO-DATAVIZ

Analyses of social specialization normally start by identifying target populations, whose segregation indicators are later calculated (Apparicio 2000). In our research, we opted for a bottom-up uncertainty-based approach: clusters of households were identified through data mining of selected sociodemographic variables within the 2008 Household Mobility Survey (which unfortunately does not cover the Principality of Monaco) combined with place of residence. Our starting hypothesis is that wealth differences are not the only factor beyond social specialization. A sample of 7539 households, representative of the population of the 94 sectors of the metropolitan area, has been analyzed through 16 variables describing social status, household structure, household demography and place in the workforce. Weights are attributed to variables in order to give approximately the same total weight to the four different dimensions of the analysis (Table 1), social status being slightly overweighted as its three variables cover more diverse issues. The *a priori* clustering of variables in four groups has been validated through Bayesian clustering algorithms based on mutual information among variables, given the empirical data. A 10-fold cross-validation procedure yields an average fit score of over 90% for this variable grouping.

<table><caption>Table 1 – The 16 indicators used to cluster households on the French Riviera.</caption><thead><tr><th>Thematic Area</th><th>Indicator</th><th>Weight</th></tr></thead><tbody><tr><td rowspan="4">Place in the Workforce</td><td>Occupation of the person of reference</td><td>1</td></tr><tr><td>Number of working active people</td><td>1</td></tr><tr><td>Number of unemployed people</td><td>1</td></tr><tr><td>Number of inactive people</td><td>1</td></tr><tr><td rowspan="3">Social Status</td><td>Maximum profession and socio-professional category among the spouses</td><td>2</td></tr><tr><td>Maximum qualification among the spouses</td><td>2</td></tr><tr><td>Occupancy status of the dwelling</td><td>2</td></tr><tr><td rowspan="4">Household Composition</td><td>Presence of spouse in the household</td><td>1</td></tr><tr><td>Number of household members</td><td>1</td></tr><tr><td>Number of children</td><td>1</td></tr><tr><td>Number of other members</td><td>1</td></tr><tr><td rowspan="4">Household Demography</td><td>Number of minors (0-17 years)</td><td>0.8</td></tr><tr><td>Number of young adults (18-29 years)</td><td>0.8</td></tr><tr><td>Number of adults 30-59 years</td><td>0.8</td></tr><tr><td>Number of seniors (60-75 years)</td><td>0.8</td></tr><tr><td>Total Weight 4</td><td>Number of elderly (more than 75 years)</td><td>0.8</td></tr></tbody></table>

Once the residence sector of the household is added as a further variable, different strategies of Bayesian clustering (Korb and Nicholson 2004) have been explored in order to produce an uncertainty-based socio-geographic clustering. More particularly, a naïve Bayesian classifier has been used with the variable weights of Table 1, with different constraints on maximal
number of clusters and minimal cluster content and with different weights for the new variable sector of residence. A minimal cluster content of 4% of the household population has finally been selected to avoid overfitting on the sample. As for the sector variable, weights beyond 3 forces the algorithm to consider it as the leading variable of the clustering: the likelihood of the resulting clustering is maximized by assigning the population of a given sector to one or two clusters only. This is clearly not the goal of our clustering, as we want to find dividing lines in the resident population that take into account place of residence together with and not instead of sociodemographic differences (the poor quality of clustering results based on place of residence only is also witnessed by the low contingency table fit on the 17 variables). The optimal compromise was found with a weight of 2 for the sector variable. In this case, the sector variable is only the fourth strongest variable in terms of mutual information with the clusters and the contingency table fit is 36,4% on all the variables.

This approach is different from previously developed research (Pallez et al. 2015): clustering is uncertain because household assignment to a given cluster is probabilistic and households can have several non-zero probabilities of being assigned to different clusters. Average probabilities in assigning a given household to a cluster range between a maximum of 0.97 (for cluster 11) to a minimum of 0.86 (for cluster 4). Individual households can have much lower probabilities, and sometimes even have similar probabilities of being assigned to two different clusters. Passing from the sample to the household population introduces additional uncertainties. Cluster labels are vague, too, in the sense that they are synthetic descriptions combining different variable values which are often (but not constantly) associated. The sociodemographic characteristics of clusters are thus described through Bayesian probabilities. Social specialization of metropolitan sectors with respect to these clusters is evaluated in terms of the classical dissimilarity index (Duncan and Duncan 1955), but different evaluations are proposed for different levels of uncertainty.

One of the main difficulties of the analyses was to convey the uncertainty associated with the results obtained. Several authors have already proposed approaches for graphical representation of uncertain information (MacEachren 1992, MacEachren and Howard 1993, Ehlschlaeger et al. 1997, Cedelnik and Rheingans 2000, Ward 2002). Within our research, we propose an interactive online geo-data visualization solution in order to explore the results of uncertainty-based analyses (Cao and Fusco 2015). Systems of dashboards for interactive visualization seem particularly useful in representing uncertain and complex phenomena like social specialization of space. First attempts of interactive representation of uncertain geographic data were proposed in the 1990s (Ehlschlaeger et al. 1997) even if they were not considered particularly effective in their applications (Evans 1997). Advances of the software interfaces have since been considerable and new applications seem to be both more user-friendly and scientifically sophisticated (Ban and Ahlqvist 2009, Kunz et al. 2011, Fusco et al. 2016). These applications link together interactive representations in the forms of maps, diagrams and text. In our case, different analyses are developed with respect to uncertain knowledge and represented in the geo-dataviz solution.
### III RESULTS FROM THE FRENCH RIVIERA

Results of the Bayesian clustering of households in the French Riviera will be presented and commented through four visualizations, taken from the interactive geo-dataviz solution (Cao and Fusco 2015).

<figure><img src="image_92.png" /><figcaption>Figure 1: Probabilistic description of household clusters.</figcaption></figure>

Figure 1 gives a probabilistic description of household clusters. The Bayesian classifier identified 12 clusters of households, which were later regrouped in 11 clusters after interpretation: five clusters concern families with children and differ in terms of social status, from highest (Cluster 9, families of skilled executives/professionals, very often owners of their dwelling) to lowest (Cluster 3, families of employees, with lower level of education and more frequently tenant than owner, sometimes social housing tenant); two clusters concern couples of retirees (of different social status); two clusters concern single adults (here social status is at least partially correlated with age); two clusters concern single retirees (they have been regrouped in one cluster only as their sole difference was the different age class 60-75 years or 75 years and more); a last cluster is specific to households of single parents with children, with difficult social situation (low education levels, employee not always with a full-time job and sometimes unemployed, tenant of social housing or in the private sector). The latter is, by no surprise, the most segregated cluster, with a dissimilarity index of 0.75. The exploration of the probabilistic values of the different variables for each cluster is of course richer than this simple summary of main features and lets the user better understand the vague content of the cluster labels and the subtle differences that sometimes exist between two relatively similar clusters.

Figure 2 represents through maps and diagrams the socio-demographic content of metropolitan sectors. Clusters of households are differently distributed in space. Every sector of the metropolitan area has a different probability profile in terms of cluster belonging of its resident population. Clusters can also be of particular importance for a given sector. In this visualization, like in the following ones, sectors can be linked to clusters according to different criteria: modal cluster, most overrepresented cluster (largest positive deviate from metropolitan
average), most underrepresented cluster (largest negative deviate from metropolitan average),
most characteristic cluster (highest location quotient compared to the metropolitan average).

<figure><img src="image_93.png" /><figcaption>Figure 2: Sociodemographic content of metropolitan sectors.</figcaption></figure>

Finally, the dissimilarity index measures which proportion of a sectors' population should
move in order to attain the same probability distribution as the whole metropolitan area.
Uncertain knowledge of population content of sectors can be explored in the visualization.
Uncertainties derive both by the uncertain cluster assignment of households and by the sample
process. Statistical significance of the observed deviates from metropolitan-wide values is
evaluated against the null hypothesis that the percentage of a given cluster in the household
population of the sectors is given by a binomial distribution whose expected value is the share
of the cluster in the whole metropolitan household population. Given the sample size of every
sector (always between 70 and 103), the binomial distribution can be approximated by a
normal distribution. In the data-viz, the significance level can be set by the user: every time
more certainty is required in knowledge, non-significant differences are omitted. The city-
center of Cannes, for example, with strong and extremely significant over-representation of
single retirees, is not necessarily the sector with the highest dissimilarity index, but becomes
such once only the most significant deviations from the metropolitan values are retained. The
user can thus interactively identify the deviates characterized by the lowest levels of
uncertainty, which contribute significantly to the social specialization or residential space.

Socio-spatial contrasts between neighboring metropolitan sectors are represented in Figure 3.
The diverse spatial repartition of clusters within the metropolitan area results in socio-
demographic differences among contiguous sectors, which can have important impacts on the
social functioning of the metropolitan area. The population content of sectors being known as
probability distributions, their differences are measured as distribution divergence. More
particularly, we use the Jensen-Shannon divergence (Lin 1991), which has the advantage of
being symmetrical and defined even when some clusters are absent in a given sector. Socio-
demographic distance between contiguous sectors is thus coherently measured with a
probabilistic divergence. Once again, the uncertainties in socio-demographic content of sectors
can be used to produce different measures of socio-spatial contrasts, according to significance level of differences from the metropolitan-wide values. The user can thus identify those contrasts which are identified with the lowest uncertainty. The socio-spatial contrasts in the city of Nice (at the center of the metro area) are thus determined with higher levels of certainty than those west of Nice, where some differences are relatively uncertain.

<figure><img src="image_94.png" /><figcaption>Figure 3: Socio-spatial contrasts between neighbouring metropolitan sectors.</figcaption></figure>

<figure><img src="image_95.png" /><figcaption>Figure 4: Proximity of metropolitan sectors in variable space.</figcaption></figure>
Figure 4 represents proximity of metropolitan sectors in variable space. Jensen-Shannon divergence values make up a complete distance matrix among the metropolitan sectors. A principal component analysis can be performed in order to visualize the relative proximity of sectors in terms of socio-demographic content, on the main factorial plane. Geographical space proximity can thus be evaluated in conjunction with sociodemographic proximity in a multi-dimensional scaling approach. Once again, reducing the level of uncertainty through minimum significance levels modifies the results of the analysis. Through this analysis we can highlight the extreme diversity of social specialization within the city of Nice, a phenomenon which is not observable in other main cities like Cannes and Antibes. At the same time, the city of Cannes presents a marked opposition between its westernmost sector, strongly marked by the presence of single parents with children, in difficult social situation, and the rest of the city, where all sectors show very high presence of single retirees, and low presence of couples with children. All these results are particularly robust to uncertainty levels and can be considered among the most salient characteristics of our study area.

## IV CONCLUSIONS AND PERSPECTIVES

In conclusion, with the example of the French Riviera metropolitan area, our research shows how knowledge of social specialization of residential space is uncertain. In this context, instead of defining a priori target populations, soft clustering techniques could be used to identify the most important sociodemographic divides in a given metropolitan area. For the French Riviera, position in the life-cycle as described by age and household composition seems even more important than social status in defining social specialization of space. Results of soft clustering benefit from appropriate geo-data-visualization solutions. In our example, a system of dashboards seems an appropriate way to describe complex phenomena like social specialization of space. Within this solution, knowledge uncertainty can be conveyed by interactive representations, whenever an appropriate calculus (in our case through the use of probabilities) can quantify the most relevant uncertainties.

The present work opens important perspectives in the search for uncertain geographically meaningful clusters. Concerning the analysis of social specialization of space, ethnicity is a crucial aspect which is missing in the French statistical information system and could only be integrated in the study through *ad hoc* surveys. As far as the methodology is concerned, the use of a naïve Bayesian classifier is a first solution and more sophisticated methods should produce better results. Multi-level clustering with hierarchical naïve Bayesian classifiers (Langseth and Nielsen 2006) could for example better exploit the structure of available information, where strong relations exist within groups of variables. The problem of identifying a correct weight for the geographic variable could be eliminated altogether by the use of algorithms of multi-objective optimization. We could thus optimize at the same time clustering likelihood based on the sociodemographic variables and dissimilarity index based on the geographic distribution of the clusters. The exploration of the Pareto front would then be instructive of role of the two criteria in the identification of socio-geographic clusters.
## References

Apparicio P. (2000). Residential segregation indices : a tool integrated into a geographical information system. *Cybergeo: European Journal of Geography*, 134, http://cybergeo.revues.org/12063

Ban H., Ahlqvist O. (2009). Representing and negotiating uncertain geospatial concepts - Where are the exurban areas?, *Computers, Environment and Urban Systems*, 33(4), pp. 233-246.

Billard G., Madoré F. (2009). Les Hauts du Vaugrenier: un exemple atypique de fermeture résidentielle en France, *Mappemonde*, 1-2009, http://mappemonde.mgm.fr/num21/lieux/lieux09101.html

Blanc M. (2010). The Impact of Social Mix Policies in France, *Housing Studies*, 25(2), pp. 257-272.

Cao C., Fusco G. (2015). *Representing Uncertain Clustering. The case of Social Specialization on the French Riviera*, https://public.tableau.com/profile/fusco#!/vizhome/RepresentingUncertainClustering/Story

Cedilnik, A., Reinghans, P. (2000). Procedural annotation of uncertain information, In *Proceedings of Visualization '00*, IEEE Computer and Society Press, pp. 77-84.

Centi C. (1993). Les enjeux du modèle niçois. L'approche localiste du développement en question. *Revue Economique*, 44(4), pp. 687-712.

Duncan, O., B. Duncan (1955). A methodological analysis of segregation indexes. *American Sociological Review*, 20(2), pp. 210-217.

Ehlschlaeger C., Shortridge A., Goodchild M. (1997). Visualizing Spatial Data Uncertainty Using Animation, *Computers & Geosciences*, 23(4), pp. 387-395.

Fusco G., Cao C., Dubois D., Prade H., Scarella F., Tettamanzi A. (2016). Social Polarization in the Metropolitan Area of Marseille. Modelling Uncertain Knowledge with Probabilistic and Possibilistic Networks, ECTQG 2015 Proceedings, *Plurimondi*, 8 p.

Fusco, G. et F. Scarella (2011). Métropolisation et ségrégation sociospatiale. Les flux des migrations résidentielles en PACA. *L'Espace Géographique*, 40(4), pp. 319-336.

Korb K., Nicholson A. (2004). *Bayesian Artificial Intelligence*, Chapman & Hall / CRC.

Kunz, M., Regamey-Gret, A., Hurni, L. (2011). Visualization of uncertainty in natural hazards assessment using an interactive cartographic information system, *Natural Hazards*, 59(3), pp. 1735-1751.

Lacour, C. et F. Gaschet (2008). *Métropolisation et ségrégation*, Bordeaux: PUB.

Langseth H., Nielsen T. (2006). Classification using Hierarchical Naïve Bayes Models, *Mach Learn*, 63, pp. 135-159.

Lin J. (1991). Divergence Measures Based on the Shannon Entropy, *IEEE Transactions on Information Theory*, 37(1), pp. 145-151.

MacEachren, A. (1992). Visualizing uncertain information. *Cartographic Perspective*, 13, pp. 10-19.

MacEachren, A. M., D. Howard, et al. (1993). Visualizing the health of Chesapeake Bay: An uncertain endeavor, In *GIS/LIS Proceedings*, vol. 1, American Society of Photogrammetry and Remote Sensing/American Congress on Survey and Mapping, Bethesda MD, pp. 449-458.

Massey D. (1985). Ethnic residential segregation: A theoretical synthesis and empirical review. *Sociology and Social Research*, 69(3), pp. 315-350.

Pallez D., Serrurier M., Da Costa Pereira C., Fusco G., Cao C. (2015). Social Specialization of Space: Clustering Households on the French Riviera, *GECCO Companion 15 - Proceedings of the Companion Publication of the 2015 Annual Conference on Genetic and Evolutionary Computation*, ACM, New York, pp. 1447-1448.

Ward, M.O. (2002). A taxonomy of glyph placement strategies for multidimensional data Visualization, *Information Visualization*, 1(3/4), pp. 194-210.
# Detection of outliers in crowdsourced GPS traces

Stefan S. Ivanović*<fnref n="1" />, Ana-Maria Olteanu Raimond<fnref n="1" />, Sébastien Mustière<fnref n="1" />, Thomas Devogele<fnref n="2" />

**Abstract**

Nowadays, crowdsourced GPS data are widely available in a huge amount. A number of people recording them has been increasing gradually, especially during sport and spare time activities. The traces are made openly available and popularized on social networks, blogs, sport and touristic associations' websites. However, their current use is limited to very basic metric analysis like total time of a trace, average speed, average elevation, etc. The main reasons for that are a high variation of spatial quality from a point to a point composing a trace and a need for referential data for evaluation of their quality. In this paper we present a novel approach for filtering and detection of outliers in crowdsourced GPS traces in order to assess their spatial quality intrinsically and make them more suitable for more advanced uses such as updating referential road network of French Mapping Agency – IGN. In addition, we propose a new definition of an outlier in GPS data, adapted to intrinsic assessment of spatial quality.

**Keywords**

Crowdsourced GPS trace, filtering, outlier detection, machine learning

# I INTRODUCTION

Intensive sports activities of professionals and amateurs are very frequently recorded by using GPS devices. The traces obtained are then made openly available to the community through social networks, blogs, sport and touristic communities' web-sites. Currently, they are mostly used for visualization purposes and some basic metric data analysis (e.g. total time of the trace, distance, speed). On the other hand, they have a huge potential to be used for more advanced purposes. (Bergman and Oksanen, 2016). The aim of our research is to use crowdsourced GPS traces as a potential data source for highlighting updates in a referential road network of the French National Mapping Agency (IGN).

In our work, secondary road network in mountainous area such as hiking, bicycle and tractor paths are in our focus. Even if not always necessary, this network is very important for production of touristic maps and for other different applications such as defence and sport activities. These paths are very challenging for continuous update due to their intermittent nature (e.g. they appear and disappear very often) and various landscape (e.g. forest, high mountains, seashore) which make the update process time consuming and the traditionally used methods such as stereo-restitution, insufficient.

The crowdsourced GPS traces are collected without any protocol, with low and heterogeneous frequency sampling and mostly by low class GPS devices. In addition, they are made available with few or inexistent metadata. Moreover, various errors are introduced by different external factors such as topography, canopy, etc. These errors can cause a significant bias and may limit

<footnote n="1">¹ Université Paris-Est, IGN, COGIT Lab., France</footnote>

<footnote n="2">² Université François Rabelais de Tours, Laboratoire d'informatique France</footnote>

<footnote>*Corresponding author: stefan.ivanovic@ign.fr</footnote>
the usability of the traces in different analysis. Thus, using GPS traces for updating authoritative data makes their quality issues very important. To assess the quality of GPS traces, a first step is to detect and filter outliers. This paper focuses on detection of outliers in crowdsourced GPS traces in order to improve their spatial quality.

## II OUTLIERS DETECTION IN GPS DATA

Many studies have analysed factors that influence the quality of GPS data. Environmental factors have been found most influencing. First, topography that reflects on positional error and a number of fixed positions (Lewis et al. 2007; Cain et al. 2005; DeCesare et al. 2005). Second, canopy cover, especially the density of its coverage, type of species (Klimànek 2010; Tucek and Ligos 2002) and height (Janeau et al. 2004) affect positional error. Third, obstacles degrade the quality of GPS signal by causing multipath effect (Tucek and Ligos 2002). Since they all confirmed the influences of the environment on GPS data quality, it is necessary to take these influences into account when dealing with such quality anomalies like outliers are.

We define an outlier as a GPS point whose metrics and geometrics characteristics differ significantly from the characteristics of other points composing a GPS trace. Overall, most of works on outlier's detection in GPS data has been considered GPS measurements' errors as outliers and treated them by various filtering methods (Ordonez 2011; Duran 2012; Knight 2009). In addition, few works defined outliers from geometric point of view, as points that differ from the rest of the traces following the same path (Etienne 2013; Gil de la Vega et al. 2015). In our work, we only consider a single trace which does not require presence of other traces. That is important in situations where only one trace or few traces exist for the same path not allowing to define a pattern, like in some challenging mountainous areas.

Red arrows in Figure 1 point examples of outliers whose positions cause values of speed, distance and angle between them and their consecutive points significantly different compared to other points of the same trace.

<figure><img src="image_96.png" /><figcaption>Figure 1: Examples of outlier points.</figcaption></figure>

## III TEST DATA

Crowdsourced GPS data recorded during sport activities such as running, cycling, walking and hiking in Vosges Mountains are collected from different websites and portals of French sport associations. We chose this area due to various landscapes with a plenty of secondary roads, as well as an important number of crowdsourced GPS traces.
In total, we selected 436 traces composed by 292337 points. GPS points are theoretically described by their position (WGS84), elevation and timestamps. Practically, some points have missing attributes, that is, missing elevation or timestamp or even both. Six indicators measuring the completeness of the attributes are computed. We found that 106206 points (36%) lacked timestamp, whereas the situation with elevation was far better – only 6580 points (2%) lacked elevation. Regarding the traces, we noticed that 157 (36%) had no timestamps at all, whereas 287 (66%) had at least one missing timestamp.

## IV METHODOLOGY AND RESULTS

Due to the heterogeneous nature of crowdsourced GPS data and more random than systematic influences of environment and other known factors, the detection of outliers is a complex task. Many intrinsic and extrinsic indicators may be computed to detect outliers. However, determining relevant criteria and thresholds as well as how to combine them is also a challenging task. Thus, to manage these difficulties, we propose an approach based on supervised machine learning techniques in order to generate generic rules and thresholds.

Our approach is composed by four steps: i) the first one consists in filtering noise such as redundant points, negatives values, etc.; ii) the second is the computation of different intrinsic and extrinsic indicators; iii) the third one consists in using machine learning, which involves a supervised sampling followed by applying a classification algorithm; iv) finally, the generated rules and thresholds are applied on non-classified points of the test area.

First, within the step of noise filtering, redundant points (e.g. overlapping consecutive points) are filtered by means of zero distance values in a small time window. Negative speed values are used to detect structural errors in GPX files such as two traces irregularly merged into one or errors of GPS clock. This is an important step since it was spotted that in some cases negative speed values caused geometric anomalies of the traces. In total 4594 points (2%) were filtered thanks to this step.

To formalize outliers’ detection, we propose intrinsic metrics and extrinsic indicators. The former are calculated from GPS measures only, such as distance, speed, direction and elevation between consecutive points. The later are calculated based on the analysis of the spatial context in which GPS points are recorded such as type and density of forest, slope and its curvature, proximity of obstacles (e.g. cliffs, buildings, forest) and other features (e.g. river, lake, building).

In total, 15 different indicators were computed. One representing direction change, three based on distance and speed respectively, taking into account their mean and variations between consecutive points, four derived from elevation (GPS and corresponding DTM) and four based on spatial context such as: proximity of obstacles (e.g. buffer of cliffs, buildings, forests), point in lake/river/building, point in the forest, type of forest. It is important to stress that due to missing attributes (e.g. lack of timestamps) for some points, some measures cannot be calculated (e.g. speed). The indicators are calculated for each point taking into account previous and next point in the trace.

In order to apply the machine learning techniques, first a sampling zone was chosen so that it represented faithfully the state and heterogeneities of the entire test area and the entire pattern of GPS points. Percentage of missing attributes ( timestamps and elevation) of sampled points compared to the total percentage of missing attributes differed for only 2%. In total 2342 points were sampled and manually classified as outliers (3%) or regular points (97%). The learning process was conducted in WEKA software package using JRIP algorithm proposed by Cohen (1995). Evaluation of results by means of 30 fold cross validation is presented in a confusion matrix in Figure 2.
<table><caption>Figure 2: Confusion matrix of 30 fold cross validation</caption><thead><tr><th rowspan="2"></th><th colspan="2">Predicted</th></tr><tr><th>Outlier</th><th>Not an outlier</th></tr></thead><tbody><tr><th rowspan="2">Actual class</th><th>Outlier</th><td>61</td><td>16</td></tr><tr><th>Not an outlier</th><td>16</td><td>2249</td></tr></tbody></table>

The number of correctly classified outliers is almost four times bigger than the number of misclassified, while overestimation and underestimation are balanced. In addition, a global precision of the approach calculated from the matrix is 98%. Both results confirm high performance of the approach and generated rules, particularly taking into account presence and random distribution of missing attributes.

In total, five rules as well as thresholds were generated and presented in Figure 3.

$$1. \text{AngleMean} \ge 87^{\circ}.54 \text{ and } \text{DistDiffMed} \ge 1.05 \Rightarrow \text{outlier}$$

$$2. \text{AngleMean} \ge 71^{\circ}.23 \text{ and } \text{SpeedRate} \ge 1.5 \Rightarrow \text{outlier}$$

$$3. \text{AngleMean} \ge 74^{\circ}.80 \text{ and } \text{DistDiffN} \le 0.21 \Rightarrow \text{outlier}$$

$$4. \text{AngleMean} \ge 83^{\circ}.15 \text{ and } \text{SpeedRate} \le 0.85 \Rightarrow \text{outlier}$$

$$5. \text{AngleMean} \ge 56^{\circ}.43 \text{ and } \text{DistMean} \ge 8847.31 \text{m} \Rightarrow \text{outlier}$$

<caption>Figure 3: Rules generated in machine learning</caption>

Where: i) AngleMean represents an average value of 3 direction changes; ii) DistDiffMed a relation between two consecutive distances (before and after a point that is evaluated) and a median distance of a trace; iii) DistDiffN a normalized value of two consecutive distances; iv) DistMean defines a mean distance of two consecutive distances and v) SpeedRate represents the velocity change rate being proposed by Winden et al (2016).

Finally, the generated rules were applied on the entire test area. As a result, 9309 points (3%) were recognized as outliers. In Figure 4, an example of a successfully detected outlier is illustrated.
<figure><img src="image_97.png" /><figcaption>Figure 4: An outlier detected within non-classified points</figcaption></figure>

Although the successfulness of the approach is high, there is a specific situation where it does not perform properly causing overestimations. This happens when treating traces with very high sinuosity and low spatial resolution such as illustrated in Figure 5.

<figure><img src="image_98.png" /><figcaption>Figure 5: Misclassified outliers</figcaption></figure>

All five points composing a part of the trace that fits referential road are identified as outliers. Such geometries of referential roads are not frequent (in our test area only 11), thus the overestimation produced are not numerous. However this should be treated by taking into account the sinuosity of referential roads.

Concerning the influence of external (environmental) factors, it is obvious that based on the generated rules, the impact and importance of the external factors on appearance of outliers were not discovered. It is likely that a link between them exists, however it is difficult to be detected and modelled in crowdsourced traces due to their extreme heterogeneity, and a lack of two very important information such as precision of GPS device and the quality of signal.
For example, the same obstacle produces a multipath and subsequently outlier while collecting data by low precision GPS, whereas it is not a case for high precision GPS under the same conditions. The same situation is with canopy cover. More precise device using better quality signal can produce a trace without outliers even in closed coniferous forest, while lower precision device supplied with low quality signal would cause outliers even in an open area.

## V CONCLUSION

GPS traces are widely spread as crowdsourced data due to the significant number of sportsmen's and amateurs recording them during sport or leisure activities. Our research is aiming to use them for highlighting updates in referential secondary road network. To do that, it is necessary to assess the traces quality, particularly to deal with significant anomalies and of their geometry that we defined as outliers. In this paper, we presented a novel approach for formalization and detection of the outliers despite the high heterogeneity of crowdsourced GPS traces and dominantly random influences of external factors on their quality. The outliers were modelled and detected using 15 intrinsic metric and spatial context indicators calculated for each point in a trace. The successfulness of the approach is generally high, except in the case where a sinuosity of a referential road is very high. By analysing generated rules and the nature of crowdsourced data, we conclude that causes of the outliers are obviously highly heterogeneous. Thus, the most efficient way to model and detect them is by means of their intrinsic indicators, particularly taking into account that links between spatial context and crowdsourced traces' quality proved to be weak.

## References

- Bergman C., Oksanen J. (2016). Conflation of OpenStreetMap and Mobile Sports Tracking Data for Automatic Bicycle Routing. *Transactions in GIS* 00(00): 00-00

- Cain III J.W., Krausman P.R., Jansen B.D., and Morgart J.R. (2005). Influence of topography and GPS fix interval on GPS collar performance. *Wildlife Society Bulletin* 33(3), 926-934

- Cohen, W. W. (1995). Fast Effective Rule Induction. In *Twelfth International Conference on Machine Learning*, Tahoe City, California

- DeCesare N.J., Squires J.R., and Kolbe Y.A. Effect of forest canopy on GPS-based movement data. *Wildlife Society Bulletin* 33(3), 935-941

- Duran, A., Earleywine, M. (2012). GPS Data Filtration Method for Drive Cycle Analysis Application. In *SAE 2012 World Congress*, Detroit, Michigan

- Etienne, L. (2011). Motifs spatio-temporels de trajectoires d'objets mobiles, de l'extraction àla détection de comportements inhabituels - Application au trafic maritime. PhD thesis. Institut de Recherche de l'Ecole Navale.

- Gil de la Vega, P., Ariza-Lopez, F.J., Mozas-Calvache, A.,T. (2015). Detection of outliers in sets of GNSS tracks from volunteered geographic information. In *Agile Conference 2015*, Lisbon

- Janeau G., Adrados C., Joachim J., Gendner J.P., Pépin D. 2005. Performance of differential GPS collars in temperate mountain forest. *C. R. Biologies* 327, 1143-1149

- Klimánek M. 2010. Analysis of the accuracy of GPS Trimble JUNO ST measurement in the conditions of forest canopy. *Journal of Forest Science* 56(2), 84-91

- Knight, N. L., Wang, J. (2009). A comparison of outlier detection procedures and robust estimation methods in GPS positioning. *J. Geodesy* 62(4), 699-709

- Lewis J.S., Rachlow J.L., Garton E.O. and Vierling L.A. (2007). Effects of habitat on GPS collar performance: using data screening to reduce location error. *Journal of Applied Ecology* 44, 663-671

- Ordoñez, C., Martínez, J., Rodríguez-Pérez, J., and Reyes, A. (2011). Detection of Outliers in GPS Measurements by Using Functional-Data Analysis. *Journal of Surveying Engineering*137(4), 150-155

- Tuček, J. Ligoš J. 2002. Forest canopy influence on the precision of location with GPS receivers. *Journal of Forest Science* 48(9), 399-407

- Van Winden, K., Biljecki, F. and van der Spek, S. (2016). Automatic Update of Road Attributes by Mining GPS Tracks. *Transactions in GIS* 00(00): 00-00
# Locational Error Impacts on Local Spatial Autocorrelation Indices: A
Syracuse Soil Sample Pb-level Data Case Study

Daniel A. Griffith*¹, Yongwan Chun¹, and Monghyeon Lee¹

¹University of Texas at Dallas, USA

*Corresponding author: dagriffith@utdallas.edu

**Abstract**

This paper focuses on propagation of location errors in spatial data analysis. Specifically, it investigates how location errors impacts local spatial autocorrelation indices that often are used to identify spatial clusters. Results of a simulation experiment using heavy metal soil sample points in Syracuse, NY, are summarized. In the simulation experiment, artificial location errors were introduced to perturb points, and then local Moran's *I* and Getis-Ord local statistics were calculated. The results show that location errors have an impact on the identification of spatial clusters. Some significant spatial clusters with no location error become insignificant ones with location errors, and some insignificant ones with no location error become significant ones with location errors. More severe deviations from the true results are observed with larger location error, as expected.

**Keywords**

Location error, local spatial autocorrelation, heavy metal soil sample, uncertainty

## I INTRODUCTION

Location error occurs when the spatial information of a geotagged observation deviates from its true locational information. If spatial data have locational error, the error may increase any uncertainty associated with a spatial data analysis. Although uncertainties may render only slightly incorrect modelling results, they also can be completely fatal to an analysis of georeferenced data, and undermine the outcome of the spatial data analysis (Fisher, 1999). Any spatial model output may have incorrect results that have been corrupted by uncertainty propagated to the output. Local indicators of spatial association (LISA; Anselin, 1995) and Getis-Ord statistics ($G_i^*$; Ord and Getis, 1995) are well-known indices measuring local spatial autocorrelation and identifying spatial clusters. This paper summarizes simulation experiment results demonstrating how location error affects these spatial autocorrelation indices.

## II DATA

Simulation experiments furnishing the basis of the analysis summarized in this paper utilized soil samples collected from across the City of Syracuse, NY over three years, mainly during the summers of 2003 and 2004. These samples were used to measure a suite of six heavy metals (Fe, Pb, Rb, Sr, Zn, and Zr), assayed in a chemistry laboratory using a NITON XL-700-series X-ray fluorescence (XRF) instrument (Griffith, 2008). These samples were assayed based on a 120-s testing time and NIST 2,711 standard reference materials (SRM); measurements are in milligrams of heavy metal per kilogram of soil, or ppm (Griffith, 2008). We focus on Pb because lead poisoning is one of the critical issues in public health today, and soil Pb levels tend to be strongly related to human Pb poisoning. The total number of collected
soil sample points is 3,574. We exclude five points that are outside of the city boundary, and averaged 284 duplicate sample assays by location. Consequently, this study utilized 3,290 distinct soil sample points (see Figure 1). The Pb levels were subjected to a logarithm transformation—LN(Pb/s<sub>Pb</sub> - 12)—so that their frequency distribution better conforms to a bell-shaped curve (Griffith, 2008).

The simulation experiments involved the following steps: (1) randomly adding location error to a subset of the soil sample location points, (2) aggregating the transformed Pb measures by census geography (i.e., census tracts and census block groups), and (3) calculating local Moran's I and G* statistics. The City of Syracuse has 57 census tracts and 147 block groups in its 2000 census geography maps. One census block group (ID: 29002) has been merged with its neighborhood block group (ID: 29001), which is in the same census tract as it, because this census block group does not contain a soil sample point.

<figure><img src="image_99.png" /><figcaption>Figure 1: Soil sample points, census tracts, and block group in Syracuse, NY.</figcaption></figure>

### III METHOD

After extensive data cleaning, we artificially and randomly added locational deviations to the soil sample points in the simulation experiments. The simulation experimental design may be described as follows:

1. Randomly sample a selected percentage (beginning with 10%) of the 3,290 soil sample points;

2. Assign a specified amount of location error (beginning with 10m) with a random direction to each sampled point, constraining the perturbed locations to remain within the city limits;
3) Aggregate transformed Pb measures in each census geography administrative unit (i.e., census tract or census block group), and then calculate local Moran's *I* and $G_i^*$ statistics, recording clusters (e.g., hotspots and coldspots) for each simulation replicate;

4) Repeat Steps 1) to 3) 1,000 times (the number of replicates);

5) Repeat Steps 2) to Step 4) for 25m, then 50m, then 75m, and finally 100m of location error; and,

6) Repeat Steps 1) to 4) for 20%, then 30%, then 40%, and finally 50% of the soil sample points.

This simulation experimental design generates 25 different sample sizes (five percentages by five levels of location error), each with 1,000 replicates (to exploit the Law of Large Numbers).

## IV RESULTS

This section summarizes aspects of the two extreme simulation cases: results for the minimum (a 10% sample size, and 10m of added location error) and maximum (a 50% sample size, and 100m of added location error) error levels. Intermediate error level results are between these minimum and maximum results. Statistical significant levels for evaluating the local Moran's *I* have been adjusted using a Bonferroni correction based upon an effective sample size that adjusts for latent spatial autocorrelation (Chun and Griffith, 2013).

Figure 2 represents the original aggregated Pb level maps for 2000 census tracts and census block groups. The eastern and southern areas of Syracuse have relatively low Pb levels, whereas northwest and downtown areas of the City have relatively higher Pb levels.

<figure><img src="image_100.png" /><figcaption>Figure 2: Pb level and its local Moran's *I* cluster maps for 2000 census tracts (a), and census block groups (b) in Syracuse, NY (Note: no $G_i$ clusters exist).</figcaption></figure>
In Figures 3 and 4, red and blue boundary regions represent original (i.e., true) statistically significant positive High-High (HH; a high local Moran's *I* value region is surrounded by high local Moran's *I* value regions) and Low-Low (LL; a low local Moran's *I* value region is surrounded by low local Moran's *I* value regions) local Moran's *I* geographic clusters. Meanwhile, the $G_i^*$ maps reveal no original geographic clusters. Polygons filled with green, rather than white, are geographic clusters; the numbers that overlay them represent their aggregate detection frequencies as significant geographic clusters across the simulation experiments. For example, if a region has a red boundary, is filled with dark green, and 1,000 overlays it, this region is an original HH cluster identified in every one of the 1,000 simulation replications; in other words, location error does not affect its detection. A blue boundary region filled with light green and overlaid with 200 is an original LL cluster identified in only 200 of the 1,000 simulation replications; location error obscures or hides this geographic cluster in eighty percent of the simulation replicates. Another possible scenario is that a region has a grey boundary, which means this region is not an original cluster, but filled with green and overlaid with a number; this is a false cluster that location error creates in a certain percentage of the simulation replicates. Finally, many regions have a grey boundary and are filled with white; these are neither original nor location error created clusters.

<figure><img src="image_101.png" /></figure>
<figure><img src="image_102.png" /><figcaption>Figure 3: Local Moran's *I* and $G_i^*$ location error simulation results for census tracts.</figcaption></figure>

Figure 3 portrays census tract simulation results. When the minimum level of location error is introduced (see Figure 3a), three original LL local Moran's *I* cluster regions are colored with dark green and overlaid with 1,000. However, one non-original cluster region has been detected 227 times out of 1,000 replicates. In this latter case, location error affects the final outcome. When we introduce the maximum level of location error (see Figure 3b), even though the three original cluster census tracts are not affected, one more non-original cluster region emerges five times as a cluster. Furthermore, the census tract that has been identified as a cluster 227 times with the minimum level of location error now has been detected 662 times with the maximum level of location error. One implication here is that as the amount of location error increases, uncertainty propagation also increases and impacts spatial analyses and output. Statistically significant original hotspots or coldspots in terms of $G_i^*$ are not identified at the census tract geographic resolution. However, $G_i^*$ identifies a census tract as a significant coldspot 20 times in the 1,000 simulation replicates (see Figure 3d).

**Minimum level of location error**
(10% sampled and 10m distance added)

**Maximum level of location error**
(50% sampled and 100m distance added)
<figure><img src="image_103.png" /><figcaption>Figure 4: Local Moran's *I* and $G_i^*$ location error simulation results for census block groups.</figcaption></figure>

Figure 4 portrays the census block group simulation results. Significant $G_i^*$ cluster regions do not exist at this geographic resolution, either. However, the local Moran's *I* simulation results clearly show location uncertainty propagation. In Figures 4a and 4b, the City has 10 significant original local Moran's *I* cluster regions: three of them are HH cluster regions, and seven of them are LL cluster regions. With the minimum level of location error, all of the original clusters are identified in every one of the 1,000 simulation replicates. In other words, the location error does not affect identification of the original clusters. However, two non-original clusters have been detected three and 70 times in the 1,000 replicates. They emerge because of location error. Uncertainty propagation tends to increase as the level of location error increases. With the largest level of location error, none of the original cluster regions is identified as a meaningful cluster in all 1,000 replicates. Nevertheless, detection of these clusters still occurs over 900 times; location error impacts cluster detection in this case. In terms of false positives, 16 census block groups are detected as significant geographic clusters even though they are not original clusters. In addition, the census block groups that have been
detected as significant cluster regions with the minimum level of location error are identified only 665 and 888 times out of 1,000 replicates. These findings imply that location error can generate severely biased outcomes. Similarly, $G_i^*$ results also have been affected by location error. One non-original cluster census block group has been detected 75 times with the minimum level of location error, whereas 11 originally nonsignificant census block groups have become significant geographic clusters when introduction of the maximum level of location error.

Table 1 summarizes the number of detected significant geographic clusters in the City, both the original ones and the location error created ones. This table includes both false positives and false negatives: location error creating geographic clusters that are not true clusters; and, location error obscuring or hiding geographic clusters that are true cluster. All red numbers in Table 1 denote that location error affects local spatial autocorrelation cluster detection results. Because the census block group geographic resolution tends to have units smaller than census tracts, in terms of area, location error impacts tend to be more severe at this geographic resolution. Census block groups have more soil sample points that potentially can be perturbed across unit boundaries with a certain magnitude of location error.

<table><caption>Table 1. A cross-tabulation of the number of significant soil Pb cluster regions in Syracuse, NY.</caption><thead><tr><th rowspan="3"></th><th rowspan="3"></th><th rowspan="3"></th><th rowspan="3"></th><th rowspan="3"></th><th colspan="4">Location error added to data (1,000 replicates)</th></tr><tr><th colspan="2">Minimum (10%; 10m)</th><th colspan="2">Maximum (50%; 100m)</th></tr><tr><th>Cluster</th><th>Not cluster</th><th>Cluster</th><th>Not cluster</th></tr></thead><tbody><tr><td rowspan="4">Census tract</td><td>local</td><td>Cluster</td><td>3</td><td>3</td><td>0</td><td>3</td><td>0</td></tr><tr><td>Moran's I</td><td>Not cluster</td><td>54</td><td>1</td><td>53</td><td>2</td><td>52</td></tr><tr><td rowspan="2">G<sub>i</sub><sup>*</sup></td><td>Cluster</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Not cluster</td><td>57</td><td>0</td><td>57</td><td>1</td><td>56</td></tr><tr><td rowspan="4">Block group</td><td>local</td><td>Cluster</td><td>10</td><td>10</td><td>0</td><td>0</td><td>10</td></tr><tr><td>Moran's I</td><td>Not cluster</td><td>136</td><td>2</td><td>134</td><td>16</td><td>120</td></tr><tr><td rowspan="2">G<sub>i</sub><sup>*</sup></td><td>Cluster</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Not cluster</td><td>146</td><td>1</td><td>145</td><td>11</td><td>135</td></tr></tbody></table>

## V FINDINGS

A majority of areal units remain unchanged in terms of geographic cluster detection in the presence of location error. Nevertheless, location error does tend to change a considerable number of true geographic clusters, and it also artificially creates geographic clusters that actually do not exist. The principal implication is that location error can propagate to the final output of a spatial analysis, compromising its quality and precision. Furthermore, location error seems to affected local Moran's I results more than $G_i^*$ results.

## VI ACKNOWLEDGEMENTS

This research was supported by the National Institutes of Health, grant 1R01HD076020-01A1; any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors, and do not necessarily reflect the views of the National Institutes of Health.

## References

* Anselin, L. (1995). Local indicators of spatial association—LISA. *Geographical analysis* 27, 93–115.
Chun, Y., Griffith, D.A. (2013). *Spatial Statistics and Geostatistics: Theory and Applications for Geographic Information Science and Technology*. Thousand Oaks: SAGE Publications.

Fisher, P.F. (1999). Models of uncertainty in spatial data. *Geographical information systems* 1, 191–205.

Griffith, D.A. (2008). Geographic sampling of urban soils for contaminant mapping: How many samples and from where. *Environmental Geochemistry and Health* 30, 495–509.

Ord, J.K., Getis, A. (1995). Local Spatial Autocorrelation Statistics: Distributional Issues and an Application. *Geographical Analysis* 27, 286–306.
# Modeling spatial risk of the Foot-Mouth-Disease epidemic in South Korea

EunHye Yoo<fnref n="1" /> and JiYoung Lee<fnref n="2" />

*Corresponding author: eunhye@buffalo.edu

**Abstract:** The 2010/2011 food and mouth (FMD) disease epidemic in South Korea spread nationwide and resulted in substantial economical damage. Its rapid transmission throughout the country despite the government's control policy argues for the need to improve our understanding of its spatial dimensions. The spatial point patterns of FMD incidences revealed that FMD incidences formed statistically significant spatial clusters during the early phase of epidemic but its significance and the spatial dimensions dynamically evolved over the course of epidemic. We developed a log-Gaussian Cox process (LGCP) model to quantify and to make inferences on the spatial distribution of FMD risk, while accounting for both known risk factors and unexplained variation. The FMD risk surface estimated from the LGCP model captured the spatial heterogeneity of FMD risk driven by environmental and landscape factors and unknown risk variations at both local and regional scales. This study demonstrated that point pattern analyses and statistical models enabled investigators to improve our understanding on the spatial dynamics of FMD, which may provide useful information for decision makers to make efficient and effective control strategies for future outbreaks.

**Keywords:** foot and mouth disease, spatial point pattern analysis, log-Gaussian Cox process

## I INTRODUCTION

Foot and mouth disease (FMD) is considered by far the most serious infection between livestock due to its rapid transmission (Keeling, 2005). In South Korea (S. Korea) the outbreaks of FMD during November 2010 and April 2011 spread nationwide and resulted in a total of 3.48 millions susceptible animals being culled (Park et al., 2013). The spatial pattern analysis of infectious diseases may reveal the shape of disease and guide how to control their future spread. However, most studies on the recent FMD incidences in S. Korea have overlooked the potential association between the observed FMD cases and spatial risk factors, and merely focused on descriptive statistics (Moutou and Durand, 1994; Alexandersen et al., 2003). In this study, we aim to fill the gap of FMD literature by examining the spatial patterns of infection risk and their dynamic changes. We will make an inference on the biological processes underlying the observed FMD cases using a log-Gaussian Cox process (LGCP) model.

## II MATERIAL AND METHODS

### 2.1 Data

The official public health records of the FMD epidemic in S. Korea were used, which include the location of farms on which FMD infection was confirmed, the reported date of infection, the type of livestock infected. The types of livestock infected include cattle, pig, deer, and goat, but we focused on cattle and pigs in this paper. Assuming that the spatial variation of FMD prevalence is affected by both the environmental factors and the spatial configurations of

<footnote n="1">¹Department of Geography, University of Buffalo (SUNY), USA</footnote>

<footnote n="2">²Department of Geoinformatics, Seoul University, S. Korea</footnote>
livestock farms, we incorporated 2010 national census data on the number of domestic animals and farms. Figure 1(A) shows the elevation, (B) and (C) are the spatial distribution of livestock density and the average temperature of December of 2010 and January of 2011. The driving distance between each grid cell at the same resolution to the index farm in Andong is illustrated in Figure 1(D). The index farm was considered as an initial seed of infection because there was hardly any FMD cases reported in S. Korea prior to this epidemic in 2010 (Park et al., 2013).

<figure><img src="image_104.png" /><figcaption>Figure 1: Spatial distributions of risk factors; (a) Elevation (in meter), (b) Density of livestock farm holdings per non-overlapping 5 × 5 km² areas, (c) Average temperature in winter months (in Fahrenheit), (d) Driving distance from the point source of infection (in kilometer).</figcaption></figure>

## 2.2 Spatial point pattern analysis

The K function describes the extent to which there is spatial dependence in the arrangement of georeferenced point locations of events in the study area (Gatrell et al., 1996). This function can be estimated from an observed event distribution and its behavior in a particular situation is assessed. One of the standard reference distributions is homogeneous poisson model. Given that FMD is infectious disease which forms a spatial cluster, homogeneous model is least sensible to the observed FMD cases. As an alternative, we estimated the inhomogeneous K-function $\hat{K}_{\text{inhom}}^{\text{iso}}(r)$ (Baddeley et al., 2000) where the spatially varying intensity is estimated from the observed FMD cases. To make an inference about observed point patterns under the consideration of non-stationarity, we simulated envelopes with the inhomogeneous K-function. An observed process is assumed to be different than the null process model, if the estimated K function at a specific distance lie outside the envelope (Moller and Waagepetersen, 2003).
## 2.3 A log-Gaussian Cox Process model

A special case of Cox process where log-intensity is a Gaussian process is a log-Gaussian Cox Process (LGCP) model (Møller et al., 1998; Diggle, 2013). We modeled the spread of FMD in S. Korea as a poisson process that was partially driven by observed risk factors. This conceptualization enabled us to link the observed global heterogeneity in the dispersion of FMD to the spatial patterns of environmental factors. On the other hand, we viewed the unexplained spatial variability of the point process as outcomes of infectious cases transmitting the disease to nearby susceptible. Assuming that FMD incidence locations {$s_i, i = 1, ..., n$} are independently distributed conditional on a Gaussian spatial random field $U(s)$, the LGCP model for FMD is written as

$$ \{s_i; i = 1, ..., n\} | \mathbf{U}(s) \sim \text{Poisson}(\Lambda(s)) \quad (1) $$

$$ \log \Lambda(s) = \mu + \mathbf{X}(s)\beta + U(s) \quad (2) $$

$$ \text{Cov}\{U(s), U(s')\} = \sigma^2 \rho[|s - s'|/\theta] \quad (3) $$

the vector of covariates $\mathbf{X}(s) = [X_1(s), ..., X_4(s)]$ consists the elevation $X_1$, the density of livestock holdings $X_2$, average temperature of winter months $X_3$ at the event location s, and $X_4$ denotes a driving distance from the index farm. We included $X_4$ as a risk factor for FMD risk, because the initial FMD case on November 28th, 2010 served as a source of infection in subsequent cases and the risk of infection decreases as susceptible is further away. Based on related studies (Ferguson et al., 2001; Grubman and Baxt, 2004), which showed that cold temperatures aided the virus to persist outside the host for a longer time period, we incorporated the spatial variability of winter temperature as one of known risk factors in our LGCP model. The unexplained spatial variation $U(s)$ in FMD risk is a stochastic component (the latent Gaussian process) whose spatial structure was modeled as a function of the separation vector between any two incidence locations with a correlation function specified with two parameters — a range $\theta$ and a sill $\sigma^2$. In geostatistical jargon, a range refers to the maximum distance at which spatial dependence exists and a sill is a maximum variance of the process.

The target of inference of the LGCP model in Equation 2 and 3 includes {$\beta, \theta, \sigma, U(s)$}. Møller et al. (1998) originally used the method of minimum contrast estimation, which is often computationally formidable (Diggle, 2013). We addressed this computational challenge by using two approximations: the Gaussian Markov Random Field (GMRF) approximation to the Matern correlation (Banerjee et al., 2004) and integrated nested Laplace approximation (INLA) for the marginal posterior distributions. More specifically, we used the lgcp function implemented in geostatsp library (Brown, 2015) where a computational grid that consisted of cells 5 × 5 km² in dimensions is superimposed over the study region A. The FMD risk surface $\Lambda(\cdot)$ was defined over a collection of cells $g_l, k = 1, ..., L$ that form a disjoint partition of A. Let $Y(g_l)$ be the case count within the lth cell, where the count follows a Poisson distribution with a FMD risk $\Lambda_l$. The risk is defined as product of an offset $\delta_l$ and the intensity function (or spatial random function) $\Lambda_l$, which is a function of known risk factors and a latent Gaussian process as

$$ Y_l \sim \text{Poisson}(\delta_l \Lambda_l) \quad (4) $$

$$ \log(\Lambda_l) = \mu + X_{1l}\beta_1 + X_{2l}\beta_2 + X_{1l}\beta_3 + X_{2l}\beta_4 + U_l \quad (5) $$

The spatial random effects $U(s)$ is approximated by $U_l$ using GMRF where $s \in g_l$ with Gaussian priors with a zero mean, $\theta \in [5, 20]$ km for range of the spatial random effects, and $\sigma \in [1, 4]$ for its standard deviation.
### III RESULTS

The FMD epidemic in 2010/2011 has shown that more than 90% of FMD cases occurred at five provinces throughout the epidemic and over 60% of the total number of livestock farms in S. Korea reside within these five provinces. We further examined their relationships in the LGCP model by taking the total number of the farms in each town and treated them as population at risk. The observed FMD epidemic curve in Figure 2(a) showed that the major peaks reached on the week 4 and week 7.

<figure><img src="image_105.png" /><figcaption>Figure 2: Temporal and spatial profile of FMD incidences in S. Korea</figcaption></figure>

We defined four time periods of the FMD epidemic to represent the natural phases of the epidemic: weeks 1-2 were the period of rapid spread followed by a sharp decrease (Period 1); in weeks 3-5 there was a medium- to long-distance spread combined with localized transmission from newly affected farms (Period 2); in weeks 6-8 the epidemic peaks (Period 3); in weeks 9-22 eradication was achieved (Period 4). The locations of FMD incidences are shown in the map of Figure 2(b) with different symbols per period. We characterized their spatial patterns using the inhomogeneous K-function, and the results are summarized in Figure 3. Spatial clustering is commonly found during the first three periods, although the spatial extents at which the clustering is observed are different with varying statistical significances. K-function for Period 1 in Figure 3(A), for example, shows statistically significant clustering across all distances, but the spatial cluster exists up to 20 km in Period 2. The K-function estimates for Period 3 (Figure 3(C)) shows statistically significant spatial clustering at medium- to long-distance. Lastly, the K-function in the last period in Figure 3(D) showed regularity without statistical significance.

Bayesian inference on the LGCP model parameters were conducted using the R-INLA package in the statistical computing language R (R Core Team, 2015). The posterior distribution of spatial intensity of FMD and spatial random effects are shown in Figures 4: The lower and upper bounds of the posterior mean of spatial intensity of FMD cases are shown in Figures 4 (A), (C), and (B), respectively. The corresponding posterior distribution of spatial random effects are shown in Figures 4(D), (F), and (E). The posterior distribution of parameters are also summarized in Table 1. The posterior mean of parameter estimates $\hat{\beta}$ are similar to $\hat{\beta}^*$
<figure><img src="image_106.png" /><figcaption>Figure 3: Inhomogeneous K function for four periods; (A) Period 1, (B) Period 2, (C) Period 3, and (D) Period 4. The solid and dotted lines, respectively, represent the inhomogeneous K function estimate and the K function under the null process, a theoretical expected value. The upper and lower simulation envelopes were shaded in a light grey colour.</figcaption></figure>

except the coefficient of elevation $\hat{\beta}_1$ which was not statistically significant. It is clear that the influence of temperature appears to be relatively substantial compared to other risk factors in LGCP model, too. Both the FMD risk map in Figure 4 and the non-zero estimate of intercept $\hat{\mu}$ indicate that the risk of FMD was not constant over the study region but instead was spatially varying. Lastly, the posterior mean of range parameter ($\hat{\theta}$) is approximately 14.49 km with lower and upper bounds of 9.94 km and 20.12 km, respectively.

<table><caption>Table 1: Posterior distribution of LGCP model parameters</caption><thead><tr><th></th><th>Mean</th><th>Std. Dev.</th><th>0.025 quantile</th><th>0.975 quantile</th></tr></thead><tbody><tr><td>(Intercept)</td><td>12.25</td><td>2.99</td><td>6.38</td><td>18.15</td></tr><tr><td>Elevation</td><td>-0.08</td><td>0.08</td><td>-0.23</td><td>0.07</td></tr><tr><td>Livestock holdings</td><td>0.30</td><td>0.08</td><td>0.14</td><td>0.47</td></tr><tr><td>Average temperature</td><td>-3.84</td><td>0.83</td><td>-5.48</td><td>-2.21</td></tr><tr><td>Driving distance</td><td>-1.10</td><td>0.23</td><td>-1.56</td><td>-0.66</td></tr><tr><td>Range</td><td>14,492.42</td><td>2,602.11</td><td>9,944.41</td><td>20,120.24</td></tr><tr><td>Std. dev.</td><td>1.12</td><td></td><td>0.94</td><td>1.32</td></tr></tbody></table>

## IV DISCUSSION AND CONCLUSIONS

We described the FMD transmission process in S. Korea during 2010-2011 using spatial point pattern analyses. The application of inhomogeneous K function to the farm level FMD cases
<figure><img src="image_107.png" /><figcaption>Figure 4: Posterior distribution of the spatial LGCP prediction; (A) 0.025 quantile, (B) Mean, and (C) 0.975 quantile of FMD risk (intensity) $E\{\Lambda|Y\}$; (D) 0.025 quantile, (E) Mean, and (F) 0.975 quantile of spatial random effects $E\{U|Y\}$</figcaption></figure>

across four time periods suggested the presence of dynamically evolving the spatial patterns of FMD transmission over the course of the epidemic. We found that clusters at short scales during Period 2 and 3 were not statistically significant, which contrasts to the common assumption that contagious diseases form a spatial cluster. The absence of short scale clusters (less than 5 km) in both periods might be due to the enforcement of the government's control policy including culling and the movement restriction (Park et al., 2013).

Inhomogeneous K function analysis enabled us to account for the non-stationarity of the disease spread, but the spatial intensity was purely based on the observed FMD cases. For a further investigation of FMD transmission with respect to the behaviors and appearance of intensity function, we developed a LGCP model where spatially varying intensity of the disease was estimated via a log-linear function of environmental and socioeconomic factors. We also took a Bayesian approach to explicitly take into account parameter uncertainty. We found that LGCP models for FMD can be used as a tool for improving our understanding of the spread of FMD and optimizing disease control. For instance, our analyses demonstrated that temperature plays a critical role in determining the spread of FMD — probably because the virus tends to persist in lower temperatures. As shown in the map of FMD intensity in Figure 4(A)-(C), the area with high intensity is centered at the source of infection where the onset of FMD epidemic started. On the other hand, the maps of unexplained risk in Figure 4(D)-(F) illustrate the areas with high to low risk whose shapes and sizes vary over the study area.

The Bayesian estimates of the LGCP model parameters and prediction of random effects enabled us to assess the spatial variability of FMD incidences risk and to identify the spatial clustering under the explicit consideration of population at risk. The non-zero estimate of spa-
tial variability confirmed that FMD incidences formed a surprising tendency to cluster together even after the spatially inhomogeneous livestock farm density was taken into account. Perhaps LGCP may not be the only means of testing the hypothesis of spatial independence or interactions between events as variant K-functions can provide similar outcomes. However, LGCP models provide statistically rigorous inference and estimation and allow investigators to estimate the second order properties of the underlying process, that is, the spatial dependence of events. Furthermore, the LGCP inference using INLA approximation has relatively tractable moment properties, which allowed a simple method of parameter estimation and a flexible specification of the space-time covariance structure. We simplified the mean structure as a set of spatial covariates but this approach is naturally extendable to incorporate additional covariate information.

The proposed analysis/modeling framework is general and can be applied to other infectious diseases or FMD cases in other countries. The proposed LGCP model can fully account for the temporal aspect of FMD incidences by incorporating relevant temporal covariates and specifying spatio-temporal covariance structure (Diggle, 2013). In doing so, for example, we might have been able to take into account the changes in the susceptible population due to the culling policy. We expect that the longer time series of FMD incidence will improve our understanding (given situation) and allows us to evaluate a wide range of control strategies using the model.

## References

Alexandersen S., Zhang Z., Donaldson A., Garland A. (2003). The pathogenesis and diagnosis of foot-and-mouth disease. *Journal of comparative pathology* 129(1), 1-36.

Baddeley A. J., Møller J., Waagepetersen R. (2000). Non-and semi-parametric estimation of interaction in inhomogeneous point patterns. *Statistica Neerlandica* 54(3), 329-350.

Banerjee S., Carlin B. P., Gelfand A. E. (2004). *Hierarchical Modeling and Analysis for Spatial Data*. Boca Raton, FL.: Chapman and Hall/CRC.

Brown P. E. (2015). Model-based geostatistics the easy way. *Journal of Statistical Software* 63(12), 1-24.

Diggle P. J. (2013). *Statistical Analysis of Spatial and Spatio-temporal Point Patterns*. Boca Raton, FL: CRC Press.

Ferguson N. M., Donnelly C. A., Anderson R. M. (2001). The foot-and-mouth epidemic in Great Britain: pattern of spread and impact of interventions. *Science* 292(5519), 1155-1160.

Gatrell A. C., Bailey T. C., Diggle P. J., Rowlingson B. S. (1996). Spatial point pattern analysis and its application in geographical epidemiology. *Transactions of the Institute of British Geographers*, 256-274.

Grubman M. J., Baxt B. (2004). Foot-and-mouth disease. *Clinical Microbiology Reviews* 17(2), 465-493.

Keeling M. J. (2005). Models of foot-and-mouth disease. *Proceedings of the Royal Society of London B: Biological Sciences* 272(1569), 1195-1202.

Møller J., Syversveen A. R., Waagepetersen R. P. (1998). Log Gaussian Cox processes. *Scandinavian Journal of Statistics* 25(3), 451-482.

Moller J., Waagepetersen R. P. (2003). *Statistical inference and Simulation for Spatial Point Processes*. Boca Raton, FL: CRC Press.

Moutou F., Durand B. (1994). Modelling the spread of foot-and-mouth disease virus. *Veterinary Research* 25(2), 279-284.

Park J.-H., Lee K.-N., Ko Y.-J., Kim S.-M., Lee H.-S., Shin Y.-K., Sohn H.-J., Park J.-Y., Yeh J.-Y., Lee Y.-H., et al. (2013). Control of foot-and-mouth disease during 2010-2011 epidemic, South Korea. *Emerging Infectious Diseases* 19(4), 655.

R Core Team (2015). *R: A Language and Environment for Statistical Computing*. Vienna, Austria: R Foundation for Statistical Computing.
# Sensitivity of DBSCAN in identifying activity zones

**using online footprints**

David W. S. Wong<fnref n="1" />*, Qunying Huang<fnref n="2" />

*Corresponding author: dwong2@gmu.edu

**Abstract**

In mining spatial(-temporal) data for trajectory and activity analyses, a common task is to determine spatial clusters, which may represent activity zones. DBSCAN is a popular clustering algorithm. How the two parameters of DBSCAN that control the clustering algorithm may affect the results spatially has not been thoroughly investigated. This paper reported an incremental effort in conducting a sensitivity analysis by changing the parameter values. Preliminary results show that the two parameters work against each other to a certain degree. Increasing one parameter value (*minpts*) will break up larger clusters into smaller ones, while the other parameter (*eps*) controlled the spatial scale of clusters that can be detected.

**Keywords**

Spatiotemporal trajectories, spatial clustering, twitter, scale of cluster

# I INTRODUCTION

Recent studies in GIScience and spatial analysis often exploited Twitter data, which may be considered as individual-level spatiotemporal data, to understand the mobility patterns of population. Many time-geography studies depict individual trajectories using space-time paths (sets of connected line segments) with data of similar nature to Twitter data. A fundamental issue is that this traditional representation depicts the trajectories of individuals with absolute certainty in space and time regardless if the data were gathered just for one day or several weeks, or derived from data like Twitter.

Twitter data may be regarded as data collected from semi-random spatiotemporal sampling of individual movements. Huang and Wong (2015) proposed a framework to use such data to depict the regular mobility patterns of individuals with certainty levels. The critical process is to identify zones that individuals visited regularly within the 24-hour period. Zones are formed by clustering visited locations within the same temporal window. These zones over different periods are connected with 3D cones to show the spatiotemporal (ST) trajectories. The variable sized ST cones depict the spatial variability of the trajectories and different colour hues on the cones indicate the levels of (un)certainty. A critical step in determining the zones regularly visited by the individuals employed a highly popular clustering method proposed by Ester et al. (1996), the density-based spatial clustering of applications with noise (DBSCAN). This clustering method has been adopted in many mining applications of spatial data (e.g., Birant and Kut, 2007; Huang and Wong, 2015; Zhou et al., 2004). How DBSCAN is executed in determining the activity zones will affect the sizes of the ST cones and thus the depiction of the ST trajectories, introducing another uncertainty dimension through the method, in addition to the uncertainty in the data.

DBSCAN requires two parameters: the minimum number of points that can form a cluster (*minpts*) and the maximum distance within which two points belong to a cluster (*eps*). Using different parameter values may obtain different results and thus the spatiotemporal trajectory

<footnote n="1">¹George Mason University, USA</footnote>

<footnote n="2">²University of Wisconsin-Madison, USA</footnote>
identification results are likely dependent on these parameter values. Although some studies have suggested appropriate values for *minpts*, these recommendations were data dependent and are not generally applicable. Although conducting a full-scale sensitivity analysis of DBSCAN is warranted to derive some general rules or guidelines in using this clustering method, the objective of this paper is more limited: evaluating the impacts of varying DBSCAN parameters on zone identification and thus ST trajectory variability depiction.

## II METHOD AND DATA

Following the recommendations provided by Ester *et al.* (1996), Huang and Wong (2015) set the *minpts* value to 4 and the *eps* value was determined using an iterative procedure. On the other hand, Zhou *et al.* (2004) used 20 meters as the *eps* value as this value approximated the positional error in GPS readings. As these recommended methods in determining the parameter values were based upon specific studies, they are not necessarily applicable to other research contexts, particularly in determining the regular activity zones of individuals based upon locations reported through social media data. Therefore, in this study, we will conduct a limited scope sensitivity analysis using a range of parameter values for DBSCAN to explore how the activity patterns and ST trajectories depiction will vary.

To assess the impacts of varying *minpts*, we use values range from 3 to 10, as Ester *et al.* (1996) suggested using 4. While Zhou *et al.* (2004) suggested an *eps* value of 20 meters, we test its impacts with a much larger range from 10 meters to 60 meters. The maximum *eps* value of 60 meters is probably sufficient to accommodate the locational uncertainty of human activity patterns at the intra-urban scale, but it may not be sufficient for activities conducted in the suburban or even rural areas. Therefore, in choosing our test data, we limited our choices to urban settings. For the sensitivity analysis, Twitter users in Washington, DC and Chicago, IL were used. From the tweets posted between January 1, 2014 and March 31, 2014 with users who used "Washington, DC" and "Chicago" in their profiles, we identified more than 9000 and 17,000 unique users in the two cities. After screening the data with other requirements (e.g., minimum number of geo-tagged tweets has to be 3 or more), 4,442 and 4,088 users from DC and Chicago, respectively, were used. Each of these users posted more than 3 geo-tagged tweets throughout the approximately two-year period (we retrieved the maximum number of tweets, 3,200, for each user).

In addition, we also selected one user likely with his/her residential locations in Washington, DC (we rarely can be definitely sure about the home locations of users) for detailed analysis. Specifically, we want to examine in detail how the clustering results are affected by different *minpts* and *eps* values. Using these data, we test different parameter values using DBSCAN to determine their activity clusters.

## III RESULTS OF SENSITIVITY ANALYSIS

The two concerned parameters for DBSCAN control two properties of clusters to be identified. The *minpts* parameter controls how dense point locations will constitute a cluster. The *eps* parameter controls the spatial extent of a cluster. With smaller *minpts* values, more but smaller clusters are expected to be identified. Therefore, using smaller *minpts* values will likely increase the probability of committing type I error (false positive), including random locations that may not be regular activity locations. Using larger *minpts* values will likely produce fewer but larger clusters, likely increasing the probability of type II error (false negative) and collapsing distinctive clusters. On the other hand, how *eps* value may affect clustering results has not been clear, although larger *eps* values may possibly produce smaller numbers of clusters but larger in size.
Figure 1 summarizes the results by varying *minpts* from 3 to 10, and *eps* from 10 to 60 for both Washington, DC and Chicago, IL. When *minpts* increases from 3 to 10, the averaged number of clusters decreases. This result is intuitively expected. When the clustering process requires a larger minimum number of points in order to form a cluster, given that the total number of points to be clustered is fixed, fewer clusters can be formed. This general pattern was found in both Washington, DC and Chicago, IL and is valid regardless of *eps* value. While the impact of changing *minpts* value on the number of clusters seems obvious, how the changes in *minpts* value affect the process is not apparent.

<figure><img src="image_108.png" /><figcaption>Figure 1: Averaged number of clusters with varying *minpts* and *eps* values (left: Washington, DC; right: Chicago, IL)</figcaption></figure>

Figure 1 also shows that the average number of clusters increases slightly with increasing *eps* values when *minpts* is kept constant. When *eps* is raised, even points farther away can be included to meet the *minpts* requirement to form a cluster. Then, clusters with larger *eps* values may likely be more spatially dispersed, or geographically larger in size. This effect will become more apparent when we examine the clustering results of locations from an individual Twitter user.

We selected a set of locations of tweets from a user to examine the changes in clustering detection when changing the DBSCAN parameter values. Figure 2 shows the clusters identified by DBSCAN using a subset of locations from geo-tagged tweets of a user in Washington, DC. The cluster boundaries were determined using the minimum bounding polygon. The *eps* value was kept at 20 meters, but *minpts* was altered from 5 to 9. In the aggregate analysis described above, we expected that when *minpts* increases, the number of cluster decreases because fewer clusters can be formed as more neighbouring points are required to be reachable within 20 meters for each point to be included in a cluster. Figure 2 shows that when *minpts* value increases, the process ignores smaller cluster and thus fewer
clusters are retained. Thus, fewer clusters were formed was not because clusters were competing for points to form larger but fewer clusters. However, the areal extents of clusters shrink when *minpts* increases, quite an unexpected outcome. In fact, larger clusters broke up into smaller clusters when *minpts* value increased (See the cluster(s) on the upper Figures d and e with *minpts* = 7 and *minpts* = 8).

<figure><img src="image_109.png" /><figcaption>Figure 2: Detected clusters using different *minpts* values and a fixed *eps* value of 20 meters on dispersed locations (The geo-tagged tweets are displayed as red dots and the boundary of a cluster is depicted in blue)</figcaption></figure>

Results reported above show that when value of *minpts* increases, smaller clusters will be removed, but larger clusters will be broken up into more fragmented clusters. This unexpected outcome can be traced back to the design of the DBSCAN algorithm and how *minpts* is used in the clustering process. For each point to be added to a cluster, the point has to reach *minpts* number of points given the *eps* value. So it is possible that a point originally in a large cluster, likely at the edge, with *minpts* value, say 7, could find 7 points, but can no longer find 8 or 9 points within the *eps* value when *minpts* increases to 8 or larger values.

We also examined the varying clustering results when *eps* value changes while keeping the *minpts* value to 4. Figure 3 reports part of the results. With increasing *eps* value from 20 to 70 meters, smaller clusters are merged to form larger but fewer clusters. To a large degree, this general pattern is expected. However, results described in Figure 3 illustrate a fundamental issue in clustering analysis – spatial scale. Cluster detection is a scale-dependent process (e.g., Donnelly, 1978). The *eps* value essentially controls the spatial scale in which the clusters are to be determined. Thus, when *eps* values are small, clusters are small in their spatial extents, and vice versa.
<figure><img src="image_110.png" /><figcaption>Figure 3: Detected clusters using different *eps* values and a fixed *minpts* value of 4 on dispersed locations</figcaption></figure>

IV CONCLUSIONS

In this short paper, we report only partially results from the sensitivity study. Using larger *minpts* values not only remove smaller clusters, lowering the probability of identifying false positives, but break up larger clusters into smaller clusters. Breaking up larger clusters into smaller ones may not be regarded as “concerning” although the spatial structure of clusters of the entire study area is likely affected. While the results from changing the *eps* value were not surprising, the results highlight the role of the *eps* parameter in determining the spatial scale of cluster detection. Conceptually, determining an “optimal” scale for clustering detection is challenging, not to mention the spatial heterogeneity that may be presented in the data. Taking all these together, we may not be able to qualify the reliability of clustering results using the current tool.

The reported results just focus on the impacts of changing *minpts* or *eps* values, but not simultaneously, or their interaction. On the surface, the two parameters seem to be working against each other: increasing value of *minpts* not only reduces the number of clusters but also areas that fall within the clusters, while increasing *eps* value produces extensive clusters. How the two parameters interact needs to be scrutinized in the future. In addition, the spatial distribution of points to be clustered should also be considered since it may affect the performance of DBSCAN or other clustering detection algorithms.

References

Birant, D., Kut, A. (2007). ST-DBSCAN: An algorithm for clustering spatial-temporal data. *Data & Knowledge Engineering* 60, 208-221.

Donnelly, K. (1978) Simulations to determine the variance and edge-effect of total nearest neighbor distance. In I Holder (ed) *Simulation Methods in Archaeology*, pp. 91-95.
Ester, M., Kriegel, H-P, Sander, J., Xu, X. (1996). A density-based algorithm for discovering clusters in large spatial databases with noise. In *Proceedings of the Second International Conference on Knowledge Discovery and Data Mining (KDD-96)*. AAAI Press. pp. 226–231.

Huang, Q., Wong, D. W. (2015). Modeling and visualizing regular human mobility patterns with uncertainty: an example using Twitter data. *Annals of the Association of American Geographers* 105(6), 1179-1197

Zhou, C., Frankowski, D., Ludford, P., Shekhar, S., Terveen, L. (2004). Discovering personal gazetteers: an interactive clustering approach. Paper presented at the *Proceedings of the 12th Annual ACM International Workshop on Geographic Information Systems*, ACM, pp. 266-273.
# On Reliability of Remote Sensing Data and Classification Methods for estimating transition rules of the land-use Cellular Automata

Yulia Grinblat*<fnref n="1" />,<fnref n="2" />, Michael Gilichinsky<fnref n="3" /> and Itzhak Benenson*<fnref n="1" />

*Corresponding author: juliagri@post.tau.ac.il, bennya@post.tau.ac.il

**Abstract**

Typically, the Cellular Automata (CA) models of Land-Use/Land-Cover (LULC) changes focus on estimating the rules of the LULC changes and analysis of the simulation results. However, the models put aside the uncertainty of the LULC maps that are used for establishing the transition rules. Our study questions the reliability of Remote Sensing (RS) data sources and classification methods applied for constructing these maps.

Based on four time intervals within a 36-year period, we construct LULC maps and estimate the transition probabilities between six LULC states. The LULC maps and transition probabilities matrices (TPM) are built based on the manual interpretation of high-resolution aerial photos and classification of multispectral Landsat images for the same years.

We consider the maps and TPM derived from the aerial photos as a reference data, and compare them to those constructed from the Landsat images classified with several methods: mean-shift segmentation algorithm followed by Random Forest classification method, and three pixel-based methods of classification: K-means, ISODATA, and maximum likelihood. Then, for each classification the TPM were compared to the referenced TPM.

The accuracy assessment of all maps obtained with the pixel-based methods is insufficient for estimating the LULC TPM. The LULC map obtained with the object-based classification method fit well to that based on the aerial photos, but the estimates of TMP are qualitatively different from those constructed from the aerial photos.

This article raises doubts regarding the adequacy of Landsat data and standard classification methods for establishing LULC CA model rules, and calls for the careful reexamination of the entire land-use CA framework.

**Keywords**

Cellular automata, Landsat images, land-use/land-cover changes, Markov transition probabilities matrices, validation of RS classification methods

# I INTRODUCTION

Conceptual simplicity and the ability of explicit representation of landscapes and their changes make Cellular Automata (CA) a standard tool for simulating urban and regional land-use dynamics (Clarke et al. 1996; White and Engelen 1997), which potential for modeling Land-Use/Land-Cover (LULC) dynamics is widely recognized (Wu and Webster 1998; Pijanowski et al. 2002; Verburg et al. 2002). The major source of data for the CA modeling is Remote Sensing (RS) multispectral imagery classified for establishing land uses and covers.

<footnote n="1">¹Department of Geography and Human Environment, Tel Aviv University, Tel-Aviv,</footnote>

<footnote n="2">²The Porter School of Environmental Studies, Tel Aviv University, Tel-Aviv</footnote>

<footnote n="3">³Elbit Systems, Rehovot</footnote>
It is often reported that the CA models are quite successful in predicting LULC changes, with the high overall fit (80-90%) between the real LULC data and model outputs. This is indeed true when the validation is based on comparing the *entire modeled area*. However, the period of time covered by the CA model is, usually, between one and few decades and the fraction of the modeled area that has been changed during such a period is, typically, few percent of the entire city area. As far as initial area is excluded from the comparison, the spatial fit between the predicted and real *changes* drops down (Hagen-Zanker et al. 2005; Pontius and Petrova 2010).

A hierarchy of reasons of limited capacity of the CA models for predicting LULC changes can be proposed: (1) CA framework as a whole is insufficient for predicting LULC dynamics, due to, say, an essential part of human bounded rational decisions in land planning and management; (2) The CA framework works, but wrong CA rules are chosen; (3) The CA framework works, the rules are properly established, but the data chosen for estimating parameters of the rules do not represent the real of the LULC changes. In this paper we deal with the latter and investigate the adequacy of the RS data for calibration and validation of the CA models.

## II Testing the adequacy of classifications methods

Strangely enough, the adequacy of the RS classification for representing LULC *changes* remains on the margin of the CA modeling studies. Despite a series of publications that regard the erroneous consequences of misclassification (van de Voorde et al. 2009; van der Kwast et al. 2009) and sensitivity of the CA dynamics to the parameters of the CA rules (Liu and Andersson 2004; Jantz and Goetz 2005; Dietzel and Clarke 2006), the majority of modeling studies carelessly exploit the simplest methods of the RS images classification, take their outputs for granted, and focus on model calibration. This may evidently result in inadequate transition rules regardless of the calibration methods. The standard source of data for the CA model calibration and validation is 30m resolution LANDSAT multispectral imagery. We thus investigated the adequacy of different methods of LANDSAT images classification for establishing CA model rules.

The background of the CA model is Transition Probability Matrix (TPM) {$p_{ij}$} - a set of probabilities, per time unit, of transition $S_i \rightarrow S_j$ between the states $S_i$ and $S_j$ of the LULC CA. Our study compares TPMs estimated based on the LANDSAT maps obtained by the different classification methods to the ground truth - the TPM that is estimated based on the manual interpretation of high-resolution aerial photos of the same area.

## III Study area and Data preparation

The experimental area is the 15x6 km transect that starts in the center of the city of Netanya, Israel, and extends to surrounding agriculture areas. The period of comparison 1972–2008 (36 years) is divided into 4 intervals of 6 - 11 years, depending on availability of the LANDSAT images and aerial photos. Based on the manual interpretation of the high-resolution aerial photos, we have constructed the maps of Netanya LULC dynamics of LULC states. Six LULC states are considered: built-up areas (BU), roads (RD), agricultural (AG) and vegetation (VG) areas, open spaces (OS) and water surfaces (WA). In this short paper we present the results aggregated into three states only - built-up (BU), agriculture (AG) and the rest states (RE) that aggregate the rest four LULC states.

Four popular in the CA studies pixel-based methods and one object-based method were applied for classifying the LANDSAT images. To remind, pixel-based classification methods consider pixels individually, while object-based methods recognize spatially continuous homogeneous domains of pixels, and then assign a land use to these segments (Lu and Weng 2007). All exploited pixel-based methods are traditional first choice of a CA modeler: K-means, ISODATA, Maximum Likelihood (ML) and hybrid classification. The object-based method we apply is two-staged: mean-shift clustering segmentation (Comaniciu and Meer 2002) is followed by a Random Forest classification (Breiman 2001).
## IV The results

As can be seen in Figure 1 the fit between the LANDSAT-based maps and the map that is based on manual classification varies depending on the method. We do not present here the results of numerical analysis of this fit and focus on comparison of transition probabilities matrices. For this comparison, the TPMs all estimated for the time periods of different length, were normalized to the 10-year period.

<figure><img src="image_111.png" /><figcaption>Figure 1. Land-use maps for the part of the study area obtained with the manual classification of the aerial photos and with several classification methods applied to the LANDSAT images</figcaption></figure>

Table 1 presents the TPMs obtained with two of the applied methods, Maximum Likelihood and Segmentation, together with the TPM obtained for the manually classified map. Due to limited space, the TPMs are presented for the LULC uses aggregated into three classes and for the last time interval of the investigated period - 2000-2007, only. The aggregated land uses classes are as follows: BU - built-up areas and roads; AG - agricultural areas; RE - the rest of land uses - open spaces, vegetated areas and water surfaces.

<table><caption>Table 1. TPM for transitions between three land uses based on the 2000-2007(LANDSAT) / 1999-2008(Aerial photos) data. Probabilities are normalized to the 10 year period</caption><thead><tr><th rowspan="2">Area of transition probability</th><th rowspan="2"></th><th colspan="4">Maximal Likelihood</th><th colspan="4">Segmentation</th><th colspan="4">Aerial photos</th></tr><tr><th>BU</th><th>AG</th><th>RE</th><th>Tot</th><th>BU</th><th>AG</th><th>RE</th><th>Tot</th><th>BU</th><th>AG</th><th>RE</th><th>Tot</th></tr></thead><tbody><tr><td rowspan="3">Transition probability</td><td>BU</td><td><strong>0.75</strong></td><td>0.05</td><td>0.20</td><td></td><td><strong>0.78</strong></td><td>0.13</td><td>0.09</td><td></td><td><strong>0.99</strong></td><td>0.00</td><td>0.01</td><td></td></tr><tr><td>AG</td><td>0.14</td><td><strong>0.64</strong></td><td>0.22</td><td></td><td>0.08</td><td><strong>0.86</strong></td><td>0.06</td><td></td><td>0.02</td><td><strong>0.95</strong></td><td>0.03</td><td></td></tr><tr><td>RE</td><td>0.32</td><td>0.09</td><td><strong>0.58</strong></td><td></td><td>0.13</td><td>0.13</td><td><strong>0.74</strong></td><td></td><td>0.07</td><td>0.10</td><td><strong>0.83</strong></td><td></td></tr><tr><td rowspan="3">Area of transition n, km²</td><td>BU</td><td><strong>7.10</strong></td><td>0.46</td><td>1.90</td><td>9.50</td><td><strong>16.9</strong></td><td>2.80</td><td>1.8</td><td>21.5</td><td><strong>19.2</strong></td><td>0.02</td><td>0.11</td><td>19.3</td></tr><tr><td>AG</td><td>5.50</td><td><strong>8.90</strong></td><td>25.9</td><td>40.3</td><td>3.7</td><td><strong>40.6</strong></td><td>3.0</td><td>47.3</td><td>0.87</td><td><strong>41.2</strong></td><td>1.30</td><td>43.4</td></tr><tr><td>RE</td><td>11.9</td><td>3.30</td><td><strong>21.4</strong></td><td>36.6</td><td>2.3</td><td>2.20</td><td><strong>13.2</strong></td><td>17.7</td><td>1.70</td><td>2.30</td><td><strong>19.9</strong></td><td>23.9</td></tr></tbody></table>

As can be seen from Table 1, for the presented period, the TPMs obtained with the ML and Segmentation methods are qualitatively and quantitatively different from the TPM estimated based on the aerial photos. Most important, in reality, LULC states are changing in time
essentially less frequently than it is obtained based on the RS images classified with the ML method; for example, in reality, the probability of the AG→AG transition per 10 years is 0.95, while according to the ML map this probability is 0.64 only. The same is true for the rest of the pixel-based methods (not presented here). The TPM obtained with the Segmentation method is much closer to the reference TPM that the TPM of the pixel-based methods, but yet essentially biased towards more changes than in reality. Similar differences are characteristic of all other periods, as well as for the TPMs constructed for the basic, non-aggregated, set of the LULC states.

## V Conclusions

We thus conclude that none of the maps obtained, based on the LANDSAT images, with the help of the popular pixel-based classification methods can be exploited for establishing CA transition rules. Object-based method provided better, but yet insufficiently precise estimates. We call for the revision of approach to the CA calibration and validation. An open depository of high-resolution, carefully validated, long-term series of the land-use/cover maps that reflect different types of LULC dynamics, and represent different types of land planning systems for different periods of population growth and economic development should be established. Instead of establishing a new database for every new CA model, one has to use these data series for calibration and validation of her/his new model. Only then, the model can be applied to the new dataset which, as we have demonstrated, must be constructed with the great care.

## References

* Breiman, L. 2001. Random Forests. Machine Learning 45 (1):5-32.

* Clarke, K. C., S. Hoppen, and L. Gaydos. 1996. A self-modifying cellular automata model of historical urbanization in the San Francisco Bay area. *Environment and Planning B: Planning and Design* 24:247 - 261.

* Comaniciu, D., and P. Meer. 2002. Mean shift: a robust approach toward feature space analysis. *Pattern Analysis and Machine Intelligence, IEEE Transactions on* 24 (5):603-619.

* Dietzel, C., and K. Clarke. 2006. The effect of disaggregating land use categories in cellular automata during model calibration and forecasting. *Computers, Environment and Urban Systems* 30 (1):78-101.

* Hagen-Zanker, A., J. van Loon, A. Maas, B. Straatman, T. de Nijs, and G. Engelen. 2005. Measuring performance of land use models: An evaluation framework for the calibration and validation of integrated land use models featuring Cellular Automata. Paper read at 14th European Colloquium on Quantitative Geography, September 9-13, at Tomar, Portugal.

* Jantz, C. A., and S. J. Goetz. 2005. Analysis of scale dependencies in an urban land-use-change model. *International Journal of Geographical Information Science* 9 (2):217-241.

* Liu, X., and C. Andersson. 2004. Assessing the impact of temporal dynamics on land-use change modeling. *Computers, Environment and Urban Systems* 28 (1-2):107-124.

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

* Pijanowski, B. C., D. G. Brown, B. A. Shellito, and G. A. Manik. 2002. Using neural networks and GIS to forecast land use changes: a Land Transformation Model. *Computers, Environment and Urban Systems* 26 (6):553-575.

* Pontius Jr, R. G., and S. H. Petrova. 2010. Assessing a predictive model of land change using uncertain data. *Environmental Modelling & Software* 25 (3):299-309.
van de Voorde, T., F. Canters, J. van der Kwast, G. Engelen, M. Binard and Y. Cornet. 2009. Quantifying intra-urban morphology of the Greater Dublin area with spatial metrics derived from medium resolution remote sensing data. *2009 Joint Urban Remote Sensing Event*, Shanghai: 1-7.

van der Kwast, J., I. Uljee, G. Enegelen, T. Van de Voorde, F. Canters, and C. Lavalle. 2009. Using Remote Sensing Derived Spatial Metrics for the Calibration of Land-Use Change Models. Paper read at 7th International Urban Remote Sensing Conference.

Verburg, P. H., W. Soepboer, A. Veldkamp, R. Limpiada, and V. Espaldon. 2002. Modeling the Spatial Dynamics of Regional Land Use: The CLUE-S Model. *Environmental Management* 30 (3):391 - 405.

White, R., G. Engelen, and I. Uljee. 1997. The use of constrained cellular automata for high-resolution modelling of urban land-use dynamics. *Environment and Planning B: Planning and Design* 24 (3):323-343.

Wu, F., and C. J. Webster. 1998. Simulation of land development through the integration of cellular automata and multicriteria evaluation. *Environment and Planning B: Planning and Design* 25 (1):103-126.
# Evaluating performances of spectral indices for burned area mapping using object-based image analysis

Taskin Kavzoglu¹*, Merve Yildiz Erdemir¹, Hasan Tonbul¹

¹Gebze Technical University, Department of Geomatics Engineering, 41400, Kocaeli, Turkey

kavzoglu, m.yildiz, tonbul@gtu.edu.tr

## Abstract

Determining post-fire information is crucial for post-fire management activities and rehabilitation treatments. The use of robust and advanced approaches is needed to determine fire severity and thoroughly analyze post-fire rehabilitation period. Object-based image analysis (OBIA) is a powerful approach that has been successfully applied in many research problems in remote sensing arena. However, its use in forest fire and related studies including fire severity and burned area estimation is quite limited. This study was carried out in Antalya's Taşağıl district (Turkey) where according to Directorate of Forestry reports one of the largest wildfires in the Turkey occurred in 2008. The objectives of the present work are (i) to investigate the performance of object based analysis for burned area mapping; (ii) to compare the performances of widely-used burned area related spectral indices in identifying burned, slightly burned, water and non-burned areas from each other, and (iii) to delineate the boundaries of burned area. In this context, spectral indices of Normalized Burn Ratio (NBR), Normalized Vegetation Index (NDVI), Burned Area Index (BAI) derived from the satellite image were employed in analyses. Multiresolution segmentation and fuzzy membership function classifier were applied to the combinations of the selected indices (NDVI, BAI-NBR, NDVI-NBR) to discriminate burned, slightly-burned and non-burned areas from each other. Results showed that all combinations constructed in this study produced satisfactory results in terms of classification accuracy. However, the highest accuracy (98.37%) was achieved by NDVI-NBR index combination whilst the lowest accuracy (94.59%) was achieved when only the NDVI index was employed in OBIA process. It is hoped that with this work a contribution will be made for the government agencies to delineate fire perimeter and determine risk of wildfire for post-fire damage management.

**Key words:** Forest fire, Wildfire, NBR, NDVI, BAI, Object-Based Image Analysis

## I Introduction

Wildfires are one of the most important natural disasters with respect to catastrophic consequences which cause serious social, economic and environmental problems. Particularly, the Mediterranean ecosystem of Turkey is suffered from increasing number of forest fires and severities due to human-induced activities or natural conditions. According to Turkish General of Directorate of Forestry report, 142,409 hectares of forest was burned between the years 2002 and 2014, most of which are located in the Aegean and Mediterranean regions of Turkey (OGM, 2014). In this context, burned area mapping is of crucial importance for fire management and post-fire damage estimation to determine fire behavior.

Several studies have investigated the use of remote sensing in burned area mapping on Mediterranean region (Mitri and Gitas, 2004; Kavzoglu et al., 2014; Pleniou and Koutsias, 2013). In addition; various spectral indices (e.g. NBR, NDVI, BAI) have been widely utilized
to monitor fire-induced vegetation changes, including burn severity and regeneration (Chen et al., 2011; Chuvieco et al., 2002; Veraverbeke et al., 2011).

OBIA considers spectral, textural and hierarchical information of objects, a group of neighboring pixels with similar characteristics. In contrast to pixel-based image analysis, it constructs segments from the objects and show high performances by producing more accurate thematic maps. OBIA, which deals with spectrally homogenous image objects instead of single pixels, may be more effective for burned area mapping (Mitri and Gitas, 2004). Mitri and Gitas (2002) states that “the combination of the object-oriented approach and the multispectral resolution data of Landsat TM showed very promising results in burned area mapping and in discriminating between burned and the other classes of confusion.”

The aim of this study was to examine the efficiency of OBIA together with various spectral indices for the Antalya-Taşağıl wildfire in Turkey by using Landsat ETM+ images. Main objectives of this study are twofold; i) to distinguish burned and non-burned areas precisely, (ii) to assess the effectiveness of object-based image classification for burned area mapping.

## II Study Area and Dataset

This study focuses on a large fire that occurred on 31th July 2008, in the Antalya-Taşağıl region (the central coordinates: 37° 03′N, 31° 10′E) of Turkey (Figure 1). The fire was human-induced and worsened by the prevailing winds. Pinus brutia is the dominant tree species at the region and other types of Mediterranean vegetation, (e.g. maquis) exist in the study area. The region is characterized by the Mediterranean climate, mean annual temperatures ranging between 10° and 18.5° C and with dry and hot summers. The mean annual precipitation ranges between ~450 and ~1020 mm.

<figure><img src="image_112.png" /><figcaption>Figure 1: Location of the study area.</figcaption></figure>

A cloud-free Landsat ETM+ image with an acquisition date of 12 September 2008 was obtained at no cost from the United States Geological Survey (USGS) archive (http://earthexplorer.usgs.gov/). The preprocessed (radiometrically corrected and geometrically registered) image was registered to UTM projection with WGS84 datum. A widely-used object based image processing software, eCognition Developer (v9.1), was used for segmentation and classification experiments in this study.