# III Methodology

## Spectral Indices

In this study, three spectral indices (NDVI, NBR and BAI) were calculated using the following equations (Table 1) to distinguish between burned and unburned surfaces for Landsat ETM+ image. The indices were calculated using the bands wavelengths: NIR (0,75 to 0.90 µm), R (0.63 to 0.69 µm), and SWIR (1.55 to 1.75 µm) for each spectral region.

<table><caption>Table 1. Spectral indices used in this study (NIR: Near Infrared, R: Red, SWIR: Short Wave Infrared).</caption><thead><tr><td>Spectral Index</td><td>Abbreviation</td><td>Formula</td></tr></thead><tbody><tr><td>Normalized Difference Vegetation Index</td><td>NDVI</td><td>NDVI = <sup>NIR-R</sup>&frasl;<sub>NIR+R</sub></td></tr><tr><td>Normalized Burn Ratio</td><td>NBR</td><td>NBR = <sup>NIR-SWIR</sup>&frasl;<sub>NIR+SWIR</sub></td></tr><tr><td>Burned Area Index</td><td>BAI</td><td>BAI = <sup>1</sup>&frasl;<sub>(0.1-R)<sup>2</sup> + (0.06-NIR)<sup>2</sup></sub></td></tr></tbody></table>

Veraverbeke et al. (2011) state that enhancing the NBR with Landsat's thermal band provides better separability between burned and unburned lands. The latter index, the BAI, has a high discrimination ability for burned areas in the R-NIR spectral domain (Chuvieco et al., 2002). Furthermore, the NDVI has been widely used in monitoring fire affected areas (García-Haro et al., 2001; van Leeuwen et al., 2010).

## Image Segmentation and Classification

Segmentation, the first step in the OBIA, splits an image into random object primitives to build up homogenous (segments or image objects) regions. Multiresolution segmentation algorithm (Baatz and Schäpe, 2000), which is a bottom-up region merging technique was utilized in this paper. The segmentation of Landsat ETM+ image was constructed by adjusting the scale, band weights and color-shape parameters.

The scale parameter is an unitless abstract which assess the maximum allowed heterogeneity of image objects. Optionally, user can define the scale parameter according to required level of detail in image. Higher scale parameter values construct larger objects; smaller scale parameter values generate smaller objects. Some studies emphasized the relation between scale parameter and number of objects (Addink et al., 2007). Moreover, several studies have underlined the variation in classification accuracy due to the user-defined scale setting (Kim et al., 2009; Kavzoglu and Yildiz, 2014).

In this study, bottom-up approach was performed to define land use/land cover objects. Based on a “trial and error” procedure, scale parameter was set to 50 and 15 respectively for level 1 and 2 in segmentation process. Membership function classifier was implemented in the commercial software eCognition (Developer v. 9.1). Membership functions allow users to define the relationship between feature values and the degree of membership to a class using fuzzy logic. Minimum and maximum value set the upper and lower limits of the membership function (Trimble, 2011). The fuzzy sets were defined by membership functions to identify feature attributes according to spectral and contextual information. Classifications were performed on two segmentation levels: “water” class at Level 1 (i.e., scale 50); ‘burned area’, ‘slightly-burned area’ and ‘non-burned’ class at Level 2 (i.e., scale 15).
## IV Results and Discussion

Multiresolution segmentation and fuzzy membership function classifier were applied by using three different combinations of the selected indices (NDVI, BAI-NBR, NDVI-NBR) to discriminate water, burned, slightly-burned and non-burned areas from each other. In iterative steps, a two-level network of image objects was advanced. Membership functions were utilized by the thresholds of objects regarding to their spectral characteristics.

The Level 1 objects were initially categorized. SWIR band and NDVI values were utilized on membership functions of spectral information to distinguish water class, which was delineated similarly on the combinations. Level 2, main level of classification, was created smaller image objects. Three classes representing ‘burned area’, ‘slightly-burned area’ and ‘non-burned area’ were generated at this level.

In the first combination which is based on only NDVI, classification results were not sufficient to detect burn scars. It was observed that rocky land, vegetation and urban areas were mixed with burned area (Figure 2a). Mean brightness and NDVI values were used to define membership functions for burned area and slightly-burned area. NDVI based threshold values were selected as “0.089 / 0.20”, “0.20 / 0.34” and “-0.52 / 0.85” for burned, slightly burned and non-burned class, respectively. At the end of the first step, the overall accuracy of the classification was estimated as 94.59%. Secondly, BAI-NBR combination was implemented to delineate the burned areas (Figure 2b). Image segments were classified as burned area class using NBR based threshold values “-0.8 / 0.1” and excluding the BAI based threshold values “-1.0 / 1.8”. Similarly, slightly-burned area was determined using NBR based threshold values “-0.03 / 0.41” and excluding the BAI-based threshold values “-0.8 / 2.33”. The related BAI-based threshold values excluded to eliminate confusion between burned and non-burned pixels. In non-burned area class determination, NBR based threshold values were selected as “-0.39 / 0.99”. Although this approach produced good results for burned area discrimination, it was observed that several pixels were misclassified in coastal regions. In BAI-NBR combination, overall classification accuracy was estimated 97.23%. In third approach, NDVI-NBR combination was applied to distinguish burned areas from non-burned areas (Figure 2c). NBR index was utilized to differentiate burned areas. NBR based threshold value was selected as “-0.38 / -0.08” for burned area class. Slightly-burned area was determined by using NDVI and NBR indices together. Thresholds of NBR and NDVI for slightly burned area class selected as “-0.06 / 0.42” and “0.19 / 0.42”, respectively. In non-burned area discrimination, NBR based threshold value was chosen as “-0.3 / 0.79”. Finally, overall classification accuracy was calculated 98.37% for NDVI-NBR combination.

Since there is no information about post-fire perimeter map estimated by the Turkish Directorate of Forestry, the accuracy assessment was evaluated by using authors’ previous research (Kaya et al., 2014). The related study was investigated to predict burn severity with using the differenced Normalized Burn Ratio (dNBR) algorithm in same region. Training and test datasets were collected according to mentioned above study’s results. The overall classification accuracies of three combinations for burned area mapping were presented in Table 2.

Results presented in Table 2 showed that the NDVI-NBR approach increased the overall accuracy by about 4% compared to NDVI, while this increment was approximately 3% for the burned area class. It was also observed that producer’s accuracy calculated for each class was higher for BAI-NBR combination than NDVI-NBR combination.
<table><caption>Table 2. Accuracy assessment results for NDVI, BAI-NBR and NDVI-NBR analysis.</caption><thead><tr><th rowspan="2">Combination</th><th rowspan="2"></th><th colspan="4">Producer's Class Accuracy (%)</th><th rowspan="2">Overall Accuracy(%)</th></tr><tr><th>Burned</th><th>Slightly Burned</th><th>Non-Burned</th><th>Water</th></tr></thead><tbody><tr><td>1</td><td>NDVI</td><td>94.06</td><td>92.42</td><td>94.64</td><td>100</td><td>94.59</td></tr><tr><td>2</td><td>BAI-NBR</td><td>98.03</td><td>94.62</td><td>94.56</td><td>100</td><td>97.23</td></tr><tr><td>3</td><td>NDVI-NBR</td><td>96.66</td><td>97.48</td><td>98.31</td><td>100</td><td>98.37</td></tr></tbody></table>

It is apparent from Figure 2 that the classified objects of burned area (especially in NDVI-NBR and BAI-NBR combinations) could strongly be distinguished from other classes. Additionally, users and fire agencies can extract burned area borders by using burned area mask for the purpose of fire perimeter definition.

<figure><img src="image_113.png" /><figcaption>Figure 2: Classification results of three different combinations: (a) NDVI, (b) BAI-NBR, (c) NDVI-NBR.</figcaption></figure>
## V Conclusion

In this study, effectiveness of object-based classification was assessed in the determination of burned area mapping of Antalya-Taşağıl forest fire in Turkey. For this purpose, spectral indices of NDVI, NBR and BAI derived from Landsat ETM+ image were computed in analyses. Multiresolution segmentation and fuzzy membership function classifier were implemented to determine burned area.

The findings from this study make several contributions to the current literature. The combination of different spectral indices with OBIA approach showed very promising results. Furthermore, the use of membership functions contributed to reduction of misclassified pixels. The highest accuracy (98.37%) was obtained from NDVI-NBR combination in OBIA process. The following conclusions can be drawn from the present study, object-based classification can be used as a tool for rapid operational burned area mapping of Mediterranean forest fires. Unfortunately, no information about post-fire perimeter map was estimated by the Turkish Directorate of Forestry. The results obtained from these findings could provide valuable information for government agencies, planners and decision makers to construct management of forest fire strategies.

## References

* Addink E., Jong S. de, Pebesma E. (2007). The importance of scale in object-based mapping of vegetation parameters with hyperspectral imagery. *Photogrammetric Engineering & Remote Sensing*, 73(8), 905-912.

* Baatz M., Schäpe A. (2000). Multiresolution Segmentation: an optimization approach for high quality multi-scale image segmentation. In: Strobl, J., Blaschke, T., Griesebner, G. (Eds.), *Angewandte Geographische Informations-Verarbeitung XII*. Wichmann Verlag, Heidelberg, pp. 12-23.

* Chen X., Vogelmann J. E., Rollins M., Ohlen D., Key C. H., Yang L., Huang C., Shi H. (2011). Detecting post-fire burn severity and vegetation recovery using multitemporal remote sensing spectral indices and field-collected composite burn index data in a ponderosa pine forest. *International Journal of Remote Sensing* 32(23), 7905-7927.

* Chuvieco E., Martín M. P., Palacios A. (2002). Assessment of different spectral indices in the red-near-infrared spectral domain for burned land discrimination. *International Journal of Remote Sensing* 23(23), 5103-5110.

* García-Haro F.J., Gilabert M.A., Meliá J. (2001). Monitoring fire-affected areas using Thematic Mapper data. *International Journal of Remote Sensing* 22 (4), 533-549.
Kavzoglu T., Kaya S., Tonbul H. (2014, October 14-17). Detecting Burn Area and Burn Severity Using MODIS Satellite Image Via Spatial Autocorrelation Techniques. In *5th Remote sensing and Geographic Information Systems Symposium (UZALCBS'2014)*, Istanbul/Turkey (in Turkish).

* Kavzoglu T., Yildiz M. (2014, September 29-October 2). Parameter-based performance analysis of object-based image analysis using aerial and Quikbird-2 images. In: *ISPRS Annals of Photogrammetry, Remote Sensing and Spatial Information Sciences*, II-7, pp.31-37.

* Kaya S., Kavzoglu T., Tonbul H. (2014, December 15-19). Detection of Burn Area and Severity with MODIS Satellite Images and Spatial Autocorrelation Techniques. In *AGU Fall Meeting*, San Francisco/USA.

* Kim M., Madden M., Warner T. (2009). Forest type mapping using object-specific texture measures from multispectral IKONOS imagery: segmentation quality and image classification issues. *Photogrammetric Engineering & Remote Sensing*, 75(7), 819-829.

* Mitri G., Gitas I. (2002, November 18-23). The development of an object-oriented classification model for operational burned area mapping on the Mediterranean island of Thasos using Landsat-TM images. In *Proceedings of IV International Conference on Forest Fire Research / Wildland Fire Safety Summit*, Coimbra/Portugal.

* Mitri G., Gitas I. (2004). A performance evaluation of a burned area object-based classification model when applied to topographically and non-topographically corrected TM imagery. *International Journal of Remote Sensing* 25, 2863-2870.

* Pleniou M., Koutsias N., (2013). Sensitivity of spectral reflectance values to different burn and vegetation ratios: A multi-scale approach applied in a fire affected area. *ISPRS Journal of Photogrammetry and Remote Sensing* 79, 199-210.
OGM. (2014). Turkish General Directorate of Forestry Reports; Available from: http://www.ogm.gov.tr/ekutuphane/FaaliyetRaporu/Forms/AllItems.aspx. Accessed: 2016.04.05.

Trimble (2011). *eCognition developer 8.64.0* reference book. München: Germany Trimble Documentation.

van Leeuwen W., Casady G., Neary D., Bautista S., Alloza J., Carmel J., Wittenberg L., Malkinson D., Orr B., (2010). Monitoring post-wildfire vegetation response with remotely sensed time series data in Spain, USA and Israel. *International Journal of Wildland Fire* 19 (1), 75-93.

Veraverbeke S., Harris S., Hook S. (2011). Evaluating spectral indices for burned area discrimination using MODIS/ASTER (MASTER) airborne simulator data. *Remote Sensing of Environment*. 115(10), 2702-2709.
# Irregularly Sampled Data in Space and Time: Using Poisson Kriging to Reduce the Influence of Uncertain Observations in Assessing the Risk of Aflatoxin Contamination of Corn in Southern Georgia, USA

Ruth Kerry<fnref n="1" />, Brenda Ortiz<fnref n="2" />, Ben Ingram<fnref n="3" />, Brian Scully<fnref n="4" />, EunHye Yoo<fnref n="5" />

*Corresponding author: ruth_kerry@byu.edu

## I. INTRODUCTION

Aflatoxin is a mycotoxin produced by fungi (A. flavus or A. parasiticus) which can contaminate several staple crops such as corn (Payne, 1992) and can cause liver cancer in humans and animals (Barrett, 2005). The Food and Drug administration office of the USA have set a limit of 20 ppb, total Aflatoxin, for interstate commerce of food and feed (FDA, 2015). Infection of corn with A. flavus or A. parasiticus is driven by high temperature and drought conditions (Guo et al. 2008) associated with particular climatic areas (Abbas et al. 2007), agro-ecological zones (Setamou et al. 1997) and soil types (Palumbo et al. 2010). In the Southeast of the U.S.A, corn planted as a summer crop is highly susceptible to Aflatoxin contamination. Rainfall variability and high temperatures in this region during summer, along with light textured soils can induce water stress and contamination. Also, lack of irrigation infrastructure in some areas further aggravates the situation (Brenneman et al. 1993). Salvacion et al. (2011) found that June maximum temperatures and precipitation were key predictors of Aflatoxin levels in southern Georgia (GA), USA.

There are several methods for accurate Aflatoxin measurement, but most are time-consuming and expensive (Papadoyanis, 1990) and conducted at harvest which does not allow implementation of in-season adaptation strategies to reduce the risk. Given the expense of Aflatoxin measurement, an important goal for agricultural extension services is to identify those years and counties most at risk of contamination to reduce unnecessary expense on testing in years and areas when there is little risk of contamination. Identifying such years and counties would also allow adaptation of management strategies in season to reduce contamination risk and identify areas where more resistant varieties (Chen et al. 2002) of corn should be planted or irrigation infrastructure improved.

This paper aims to use Poisson kriging to give a space-time summary of Aflatoxin contamination data collected irregularly, both spatially and temporally, over a 27 year period (1977-2004) in 53 counties in southern GA. This will then be used to evaluate a risk factors approach to identify counties with the greatest risk of Aflatoxin contamination based on weather, soil and crop variables. Due to the irregular sampling in space and time, and highly skewed data approaching a Poisson distribution typical geostatistical methods are unstable for analysing spatial and temporal patterns in this data.

<footnote n="1">¹Brigham Young University, USA</footnote>

<footnote n="2">²Auburn University, USA</footnote>

<footnote n="3">³Universidad de Talca, Chile</footnote>

<footnote n="4">⁴USDA-ARS, USA</footnote>

<footnote n="5">⁵University at Buffalo, SUNY, USA</footnote>
## II. METHODS

### Data Collection

From 1977-2004, corn samples were collected at harvest using a grab sampling technique with an average of 3 replications in 53 counties in southern GA (Fig. 1a) for Aflatoxin assessment by the Plant Pathology Department, University of Georgia at Tifton, GA. Data was collected for 23-45 counties each year but exact locations within counties was not recorded. For all years combined there was a total of 705 observations and these data approached a Poisson distribution (Fig. 1b).

<figure><img src="image_114.png" /><figcaption>Figure 1. (a) Counties in Southern Georgia where Aflatoxin was measured, the proportion of years measurements were made and the location of weather stations (b) Frequency distribution of Aflatoxin values measured in southern Georgia 1977-2004.</figcaption></figure>

Based on the findings of Salvacion et al. (2011), monthly maximum temperatures for June (June Tmax, °C) and June rainfall data (June RF, mm) were obtained for each year 1977-2004 from the Georgia Weather Network (http://georgiaweather.net). Fig. 1a shows the location of weather stations within GA. As all counties do not have a weather station, some have more than one, the weather stations are not located at the center of the county and have different installation dates, the weather variables were ordinary kriged (OK) to county centroids and a 1 km grid.

The area planted with corn per county was determined using The CropScape - Crop Land data layer produced by the National Agricultural Statistics Service (NASS, http://nassgeodata.gmu.edu/CropScape/). Unfortunately this information was only available from 2008-2009 onwards and so does not coincide with the period of Aflatoxin collection. So this and the soils data derived using it are quite uncertain and may not give an exact reflection of the proportion of each county under corn production because it assumes that the same rotations are used and that the proportion of agricultural land per county has not greatly changed over time.

A geo-corrected 1:250,000 map of soil associations (NRCS, 2006) was simplified and used to generate a map with 3 drainage classes: excessively, well and poorly drained soil. The percentage of
land areas with soil in each drainage class in the corn growing area (as identified using Cropscape above) was calculated for each county.

## Geostatistical Analysis

Making sense of the Aflatoxin data with geostatistical methods was difficult. If a map of Aflatoxin contamination for each year for all counties in southern GA is to be produced then a variogram must be computed with an average of 37 data points and as few as 23. Variograms for individual years were unreliable and showed little spatial structure. This is typical of highly skewed (Kerry and Oliver, 2007a,b) and sparse data (Webster and Oliver, 1992). Poisson kriging (Monestiez et al. 2006) is ideal for data with a Poisson distribution and that have been irregularly observed in space or time. The proportion of years a county had Aflatoxin levels > 20ppb and > 100 ppb were Poisson kriged to county centroids and also to a 1 km grid. As proportions have a numerator and a denominator, the numerator was the number of years Aflatoxin levels were above one of the thresholds in a given county and the denominator was the number of years Aflatoxin was measured in that county. The influence of proportions for counties with fewer observations is down-weighted during variogram computation using the following weighted estimator

$$ \hat{\gamma}_{Rv}(\mathbf{h}) = \frac{1}{2 \sum_{\alpha, \beta}^{N(h)} \frac{d(v_{\alpha})d(v_{\beta})}{d(v_{\alpha}) + d(v_{\beta})}} \sum_{\alpha, \beta}^{N(h)} \left\{ \frac{d(v_{\alpha})d(v_{\beta})}{d(v_{\alpha}) + d(v_{\beta})} [r(v_{\alpha}) - r(v_{\beta})]^2 - m^* \right\} \quad (1) $$

where $N(h)$ is the number of pairs of counties $(v_{\alpha}, v_{\beta})$ whose denominator weighted centroids are separated by the vector $\mathbf{h}$, and $m^*$ is the denominator-weighted mean of the $N$ area ratios. The usual squared differences, $[r(v_{\alpha}) - r(v_{\beta})]^2$, are weighted by a function of their respective denominator sizes, $d(v_{\alpha})d(v_{\beta})/[d(v_{\alpha})+d(v_{\beta})]$, which gives more importance to more reliable data pairs based on larger denominators. Poisson kriging is a form of kriging with non-systematic errors and is parametric, modelling the noise attached to each observation with a Poisson distribution. Observations with small denominators receive less weight in kriging, by adding an error variance term to the diagonal of the kriging system (Monestiez et al. 2006, Goovaerts, 2005). Poisson kriging was done in SpaceStat (BioMedware, 2013).

## Risk Factors Approach

By applying OK, risk factor data (June TMax, June RF, % Corn and Soil Type) were generated for each county and the nodes of a 1 km grid. The OK data was used to determine if risk factors exceeded a certain threshold converting it into indicator data (0/1). Table 1 shows the thresholds chosen for each variable. The thresholds for June TMax and June RF were chosen with respect to 30-year normals in southern GA and values receiving a (1) show hotter or drier than normal years. The indicator thresholds for other risk factors were determined based on natural marked breaks in the frequency distribution or were based on values associated with the tails of a normal distribution.

Once the number of risk factors above the specified threshold for each county in each year was determined, the relationship between these data and the Poisson kriged Aflatoxin data was assessed. This suggested broad groupings of years and counties with different levels of Aflatoxin contamination risk. These broad groupings were used to define grouping variables for Mann-Whitney U and Kruskal-Wallis H comparison tests to determine if there were significant differences in Aflatoxin levels based on the thresholds identified by the risk factors approach.
<table><caption>Table 1. Threshold values used for Risk Factor Indicators</caption><thead><tr><td>Risk Factor</td><td>Threshold for Indicator (1/0)</td></tr></thead><tbody><tr><td>June monthly maximum temperature (°C)*</td><td>>33°C</td></tr><tr><td>June monthly Rainfall (mm)*</td><td><50 mm</td></tr><tr><td>Percent of county area growing corn (%)</td><td>>1.75%</td></tr><tr><td>Percent of county with well-drained soils (classes 1-4) (%)</td><td>>90%</td></tr><tr><td>Percent of county with excessively drained soils (classes 1-2.5)</td><td>>2.5%</td></tr><tr><td>Percent of years with 2 weather risk factors</td><td>>30%</td></tr></tbody></table>

### III. RESULTS AND DISCUSSION

In trying to summarize the risk of Aflatoxin contamination in space and time and to verify if the risk factor approach proposed in this work is useful, an obvious starting point is to calculate summary statistics for each year and county. The means (not shown) were influenced by the maximum value and this was particularly pronounced in counties where a smaller proportion of observations had been made. This suggests that there can be great spatial and temporal variability in the Aflatoxin levels and that the small number problem is an issue with data for counties where few observations have been made giving very high or low Aflatoxin contamination levels.

Temporal patterns in Aflatoxin risk are mostly associated with weather variables, specifically June Tmax and RF (Salvacion et al. 2011) which change from year to year and can help identify the specific years at greatest risk of contamination. The temporal results are not shown here but will be shown in the associated presentation to illustrate this point.

#### Analysis of Risk Data - Spatial Patterns

Analysing spatial patterns in risk of Aflatoxin contamination to identify the counties at greatest risk is more complicated than temporal analysis as it requires analysis of factors that are relatively stable in time (% corn and soil type) as well as weather variables. Fig. 2 a, c, e show the patterns in percent corn, well- and excessively-drained soil and Fig. 2 b, d, f show the indicators produced from these data using the thresholds in Table 1. Based on % corn, there is greatest risk (1) in south western GA (Fig. 2b) for well-drained soil risk is greatest in the west (Fig. 2d) and for excessively-drained soil the north is the highest risk area (Fig. 2f). Fig. 3a shows a map where the two weather risk factors (June Tmax and June RF) have been combined and the proportion of years with two weather risk factors has been kriged to the 1 km grid. This gives a temporal summary of the areas most prone to drought in southern GA. When converted to an indicator based on the thresholds in Table 1, Fig. 3b shows the areas with a high (1, >30% of years) risk of drought. Fig. 3c shows the percentages of years with three or more risk factors for Aflatoxin contamination where June Tmax and RF are considered separate risk factors. The map has similarities to Figs. 2f and 3b suggesting that weather and excessively drained soils are the greatest risk factors for Aflatoxin contamination. Nevertheless, crucial importance of both weather factors is shown by the striking similarity in the patterns shown in Fig. 3a and b which show drought summary and Fig. 3e and f which show the Poisson kriged summary of Aflatoxin levels in all years.

<footnote>*Thresholds for June Tmax and June RF were chosen with respect to 30-year normal Tmax and RF in the area to show hotter and drier years than normal</footnote>
<figure><img src="image_115.png" /><figcaption>Figure 2. Risk factors and their associated indicators plotted kriged to a 1 km grid (a-b) percent of area growing corn, (c-d) percent of area with well-drained soils, (e-f) percent of area with excessively drained soils</figcaption></figure>
<figure><img src="image_116.png" /><figcaption>Figure 3. (a) Percent of years with 2 weather risk factors and (b) associated indicator, (c) percent of years with ≥ 3 risk factors, (d) Relationship between percent of years with ≥ 3 risk factors and percent of years and >20 ppb Aflatoxin, (e) percent of years with Aflatoxin > 20 ppb and (f) percent of years with Aflatoxin > 100 ppb.</figcaption></figure>
Fig. 3d shows the relationship between the proportion of years with ≥3 risk factors and proportion of years with Aflatoxin levels >20ppb for each county. The correlation coefficient for this relationship was 0.59 which is significant at *p*<0.001. Distinct groupings of counties are visible in the plot and have been circled. Such grouping were not very well defined when just weather factors were considered and the correlation coefficient was lower, suggesting that the other risk factors (% corn and soil drainage types) help to distinguish the spatial differences between counties. The higher risk counties, circled in red and green in Fig. 3d, are the northern most counties in southern GA as well as those in the central area of the southern half of the state.

Comparison tests were performed to determine whether there were significant differences in Aflatoxin levels defined in particular ways based on risk levels identified from the scatter-graph produced (Fig. 3d) with the risk factor approach. Counties at 'high risk' of Aflatoxin contamination (with > 3 risk factors) were compared with those not at high risk (< 3 risk factors) by Mann-Whitney U tests. When average Aflatoxin levels were used for comparison there was no significant difference (*p*=0.569) but when Poisson kriged % years with > 20ppb and 100 ppb Aflatoxin data were used for comparison there were significantly higher, *p*=0.002 and *p*=0.012, proportions of years with Aflatoxin levels exceeding thresholds for 'high' risk years. Kruskal Wallis H tests comparing all risk classes identified in Figure 3d also showed more significant differences (*p*<0.001) and the expected order of classes when Poisson kriged data were used rather than average Aflatoxin data. This shows that average Aflatoxin levels do not give a good summary of the counties most and least at risk of Aflatoxin contamination while the Poisson kriged data which down-weight the influence of proportions based on low numbers of observations, do give a good summary.

## IV. CONCLUSIONS

This study showed that when data have been irregularly sampled in space and time and have a Poisson distribution, Poisson kriging is a reliable way to generate a temporal summary of spatial patterns. Simple averages were shown to be unreliable where fewer observations were made and standard geostatistical methods do not work well when data have a Poisson distribution or have few data for individual years. Comparison tests showed that counties and years identified as having the greatest risk levels using the risk factors approach did have significantly higher Aflatoxin levels. This was not the case however, for average Aflatoxin due to unreliable averages because of irregular sampling. Identifying the weather conditions and counties associated with the highest contamination risk will allow for in-season adaptation strategies such as irrigation to avoid drought as temperatures and rainfall in June are carefully monitored with respect to 30 year normals. Also, testing can be focused in the highest risk counties and very little expensive Aflatoxin testing will be needed in low risk years and the highest risk counties could plant more resistant varieties of corn.

Future work should investigate including new variables in the risk factors approach, fine tuning of the thresholds of existing variables and quantifying the relative uncertainty of predictions associated with each risk factor. Using OK risk factor data means that uncertainty would obviously be greatest where distance from sampling points is greatest and for data like % corn where data from a different time period to sampling had to be used. The data kriged to a 1 km grid shows the potential for identifying high risk areas at the sub-county level. This should be more reliable when sub-county data on in-season corn NDVI values are used as an indicator of in-season drought stress in the risk factors approach. An online interactive Aflatoxin risk assessment tool that uses the risk factors approach outlined here is currently being developed and will include NDVI data. There is the potential that such an online tool could be adapted to other crops, states and even farms so that Aflatoxin levels may be better managed.
## References

Abbas, H.K., Shier, W.T and Cartwright, R.D. (2007) Effect of temperature, rainfall and planting date on aflatoxin and fumonisin contamination, in commercial Bt and non-Bt-corn hybrids in Arkansas, *Phytoprotection*, 88, 41-50.

Barrett, J.R. (2005) Liver Cancer and Aflatoxin: New Information from the Kenyan Outbreak. *Environmental Health Perspectives*. 113, A837–A838.

BioMedware, Inc. (2013) *SpaceStat User Manual version 3.6*, BioMedware, Inc. Available at: http://www.biomedware.com/files/documentation/spacestat/default.htm

Brenneman, T.B., Wilson, D.M. and Beaver, R.W. (1993) Effects of diniconazole on aspergillus populations and aflatoxin formation in peanut under irrigated and non-irrigated conditions, *Plant Disease*, 77, 608-612.

Chen, Z.Y., Brown, R.L., Damann, K.E. and Cleveland, T.E. (2002) Identification of unique or elevated levels of kernel proteins in aflatoxin-resistant maize genotypes through proteome analysis, *Phytopathology*, 92, 1084-1094.

FDA (2015) http://www.fda.gov/ICECI/ComplianceManuals/CompliancePolicyGuidanceManual/ucm074703.htm accessed November 2015.

Goovaerts, P. (2005) Geostatistical Analysis of Disease Data: Estimation of Cancer Mortality Risk from Empirical Frequencies Using Poisson Kriging. *International Journal of Health Geographics*. 4, 31.

Guo, B., Chen, Z., Lee, R.D and Scully, B.T. (2008) Drought stress and preharvest aflatoxin contamination in agricultural commodity: Genetics, genomics and proteomics, Journal of *Integrative Plant Biology*, 50, 1281-1291.

Kerry, R. and Oliver, M.A. (2007) Determining the Effect of Skewed Data on the Variogram. I. Underlying Asymmetry. *Computers & Geosciences*, 33, 1212-1232.

Kerry, R. and Oliver, M.A. (2007) Determining the Effect of Skewed Data on the Variogram. II. Outliers. *Computers & Geosciences*, 33, 1233-1260.

Monestiez, P., L. Dubroca, E. Bonnin, J. P. Durbec, and C. Guinet. (2006). Geostatistical Modelling of Spatial Distribution of Balaenoptera Physalus in the Northwestern Mediterranean Sea from Sparse Count Data and Heterogeneous Observation Efforts. *Ecological Modelling*, 193, 615–28.

NRCS (2006) Digital General Soil Map of U.S., U.S. Department of Agriculture, Natural Resources Conservation Service, Fort Worth, Texas. Online_Linkage: URL:http://SoilDataMart.nrcs.usda.gov/

Palumbo, J.D., O'Keeffe, T.L., Kattan, A., Abbas, H. K. and Johnson, B.J. (2010) Inhibition of Aspergillus flavus in Soil by Antagonistic Pseudomonas Strains Reduces the Potential for Airborne Spore Dispersal, *Phytopathology*, 100, 532-538.

Payne, G.A. (1992) Aflatoxin in maize, *Critical Reviews in Plant Sciences*, 10, 423-440.

Salvacion, A. R., Ortiz, B. V., Scully, B. T., Wilson, D. M., Fraisse, C. W., Hoogenboom, G., Lee, R. D. (2011) Modeling Probability of Corn Aflatoxin Contamination Using Drought Index in South Georgia. *South East Climate Consortium Poster*

Setamou, M., Cardwell, K.F., Schulthes, F. and Hell, K. (1997) Aspergillus flavus infection and Aflatoxin contamination of preharvest maize in Benin, *Plant Disease*. 81, 1323-1327.

Webster, R., Oliver, M.A. (1992) Sample adequately to estimate variograms of soil properties. *Journal of Soil Science*. 43, 177–192.
# A comparison of optimal map classification methods incorporating
uncertainty information

Yongwan Chun*¹, Hyeongmo Koo¹, Daniel A. Griffith¹

¹ University of Texas at Dallas, USA

*Corresponding author: ywchun@utdallas.edu

## Abstract

Uncertainty in spatial data attributes can produce unreliable spatial patterns in choropleth maps, but only a few studies have considered uncertainty in map classification processes. Unfortunately, a less desirable classification result often is generated by existing methods. For example, most observations are assigned to a single class while the remaining classes have a very small number of observations allocated to them. Also, selection of proper criteria for an optimal map classification is difficult. The purpose of this paper is to expand the discussion about incorporating data uncertainty for map classification by extending optimal map classification strategies with Bhattacharyya distance. The proposed method is illustrated with an application of soil lead contamination measurements in the City of Syracuse.

## Keywords

Uncertainty, Map classification, Choropleth map

# I INTRODUCTION

Most spatial data attributes inevitably contain uncertainty due to sampling and/or measurement error, among other sources of error (e.g., specification). The preferred form of visualization for these data is uncertainty as well as attribute values, which is limited because of existing technical and conceptual deficiencies. One common approach is to utilize a bivariate mapping technique that simultaneously represents estimates with their corresponding uncertainty information using additional visual variables. For example, using Bertin's (1983) graphic variables, such as color value and texture, with a general choropleth map (e.g., MacEachren et al., 1998, 2005; Xiao et al., 2007).

Sun et al. (2014) propose a class separability measure to incorporate uncertainty information in a map classification to produce more reliable spatial patterns. This approach produces a more elaborated classification result. But it is heavily affected by outliers, and, at worst, results in a number of classes with a single observation, allocating most observations to a single class. Sun et al. (2016) further propose a heuristic approach to overcome this imbalance issue of the class separability classification results; but untrained users still might find selecting values of criteria for achieving an optimal map classification to be difficult.

The purpose of this paper is to expand the discussion about incorporating data uncertainty for choropleth mapping. More specifically, this study incorporates optimal map classification strategies with Bhattacharyya distance.
## II METHOD

Bhattacharyya distance is effective to quantify dissimilarities between two probability distributions (Bhalerao and Rajpoot, 2003), resulting in it commonly being used in feature selection and extraction (e.g., Choi and Lee, 2003; Reyes-Aldasoro and Bhalerao, 2006). When two observations conform to a normal distribution, Bhattacharyya distance between the two observations can be calculated with the following equation (Coleman and Andrews, 1979):

$$D_B(i, j) = \frac{1}{4} \ln \left( \frac{1}{4} \left( \frac{SE_i^2}{SE_j^2} + \frac{SE_j^2}{SE_i^2} + 2 \right) \right) + \frac{1}{4} \left( \frac{(\bar{x}_i - \bar{x}_j)^2}{SE_i^2 + SE_j^2} \right) \quad (1)$$

where $\bar{x}_i$ and $\bar{x}_j$ are the estimates and $SE_i$ and $SE_j$ are corresponding standard errors for observations i and j. Bhattacharyya distance can index the dissimilarity between two observations having the same estimate with the first term of equation (1). In general, an optimal classification can be achieved by minimizing within class Bhattacharyya distance. A within class distance may be determined in two different ways. The first is the sum of pairwise distances within a class (SB), and the other is the maximum pairwise distance within a class (MB). Although the first method calculates a measure from all pairwise distances, the second method intends to control for the worst case to improve homogeneity within a class. Class breaks can be determined in order to minimize the sum of costs, which are within class distances here.

This study utilizes optimal classification methods incorporating uncertainty information by extending the optimal classification method developed by Cromley (1996). He shows that an optimal univariate map classification can be modelled as an analogy to a constrained shortest path problem, when an acyclic network is constructed by the rank ordered estimates of observations. In this network, a node represents an observation, and lines correspond to assigning observations to a class. Here, a cost (or impedance) of a line is defined by within class Bhattacharya distance. An optimal classification can be constructed by identifying a combination of lines minimizing a total cost.

Extending Cromley (1996), optimal map classification can be formulated as follows:

$$\text{Minimize:} \quad \sum_{i,j} c_{i,j} d_{i,j} \qquad \forall i, j \in N, i \le j \qquad (2)$$

$$\text{Subjected to:} \quad \sum_{i} d_{i,k} = \sum_{j} d_{k,j} \qquad \forall i \in In_k; j \in Out_k \qquad (3)$$

$$\sum_{j} d_{s,j} = 1 \qquad \forall j \in Out_s \qquad (4)$$

$$\sum_{i} d_{i,e} = 1 \qquad \forall i \in In_e \qquad (5)$$

$$\sum_{i,j} d_{i,j} = T \qquad \forall i, j \in N, i \le j \qquad (6)$$

$$d_{i,j} \in \{0,1\} \qquad \forall i, j \in N, i \le j \qquad (7)$$

where $d_{i,j}$ is a binary decision variable, $c_{i,j}$ is the cost of a line from i to j, $In_k$ is a set of lines that terminate at node k, $Out_k$ is a set of all lines that originate at node k, and nodes s and e respectively are the minimum and maximum values. Equation (2) is the objective function, and equation (3) ensures that all observations are assigned to a class. Equations (4) and (5) ensure that the first node is assigned to the first class, and the last node to the last class. Equation (6) constrains the number of classes, and equation (7) is a binary integer restriction.
### III APPLICATION

The proposed method is illustrated with lead (Pb) contamination measurements collected from soil samples across the City of Syracuse (Griffith, 2008). These measurements are in milligrams per kilogram of soil (ppm), and are log-transformed so that their frequency distribution better conforms to a bell-shaped curve. A Pb surface was krigged with a Bessel function semivariogram model (26,402 pixel values were interpolated), and then the krigged surface was aggregated to 264 grid cell polygon values. Predicted measurements and prediction standard errors for these grid cells are utilized for map classification purposes. The proposed methods have been implemented as an extension of ArcGIS 10.1 using C# in the Microsoft .Net Framework 4. The optimization problems are solved by a branch-and-bound algorithm using Gurobi optimizer 6.5.0.

<figure><img src="image_117.png" /><figcaption>Figure 1: Map classification results for an aggregated Pb krigged surface across Syracuse</figcaption></figure>

The classification results are compared with those from the class separability method by Sun et al. (2014), and from the standard natural breaks classification method. Figure 1 displays map
classification results for Pb levels. The class separability method produces a spatial pattern that is remarkably different from the other results. This method produces classes that contain only outliers on both sides of its frequency distribution. The other three map classifications show a relatively similar spatial pattern, although they have a noticeable difference for the last class. The last class by the natural breaks method has the greatest number of observations among the map classification results. In contrast, class 6 in MB has only three observations.

Table 1 presents the numbers of observations for output classes. The results of the class separability method display considerable variation in terms of the number of observation counts across classes. More specifically, most of the observations are assigned to class 3. However, classes 1 and 5 have only a single observation. In contrast, proposed classification methods have a more balanced number of observations across classes. The natural breaks classification also produces well balanced results, but this method considers only variances among estimates without a consideration of uncertainty (Jenks, 1977).

<table><caption>Table 1. The number of observations in classification result classes</caption><thead><tr><th></th><th>Class 1</th><th>Class 2</th><th>Class 3</th><th>Class 4</th><th>Class 5</th><th>Class 6</th></tr></thead><tbody><tr><th>SB</th><td>43</td><td>56</td><td>46</td><td>48</td><td>44</td><td>27</td></tr><tr><th>MB</th><td>32</td><td>59</td><td>55</td><td>76</td><td>39</td><td>3</td></tr><tr><th>Class separability</th><td>1</td><td>3</td><td>239</td><td>18</td><td>1</td><td>2</td></tr><tr><th>Natural breaks</th><td>38</td><td>50</td><td>42</td><td>43</td><td>50</td><td>41</td></tr></tbody></table>

## IV DISCUSSION

This study extends map classification by integrating an optimal map classification (Cromley, 1996) and uncertainty information with various costs (e.g., distance here) and objective functions. Generally, the proposed methods produce homogenous classes by their respective criteria, and also achieve visually balanced classification results. Some limitations should be investigated in future studies. First, a performance evaluation for the proposed methods needs to be investigated. Because widely used evaluation methods examine only estimates (e.g. Jenks and Caspall, 1971; Xiao et al., 2007), these methods are not effective for evaluating the performance of the proposed map classification methods. Second, the proposed map classification methods are solved using a branch-and-bound algorithm. Thus, as the number of nodes increases, the amount of time needed to find an optimal solution tends to increase often at roughly an exponential rate.

## V 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

* Bertin, J. (1983). *Semiology Of Graphics: Diagrams, Networks, Maps*. Madison: University of Wisconsin press.

* Bhalerao, A. H., Rajpoot, N. M. (2003). Discriminant feature selection for texture classification. *In Proceedings Of The British Machine Vision Conference 2003*. United Kingdom.

* Choi, E., Lee, C. (2003). Feature extraction based on the Bhattacharyya distance. *Pattern Recognition*, 36(8), 1703–1709.

* Coleman, G. B., Andrews, H. C. (1979). Image Segmentation by clustering. *In Proceedings of IEEE*, 67, 773–788.
Cromley, R. G. (1996). A comparison of optimal classification strategies for choropleth displays of spatially aggregated data. *International Journal of Geographical Information Systems*, 10(4), 405–424.

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.

Jenks, G. F. (1977). *Optimal Data Classification for Choropleth Maps*. Occasional Paper No. 2, Department of Geography, University of Kansas.

Jenks, G. F., Caspall, F. C. (1971). Error on choropleth maps: definition, measurement, reduction. *Annals of the Association of American Geographers*, 61(2), 217–244.

MacEachren, A. M. A., Robinson, A., Hopper, S., Gardner, S., Murray, R., Gahegan, M., Hetzler, E. (2005). Visualizing geospatial information uncertainty: What we know and what we need to know. *Cartography and Geographic Information Science*, 32(3), 139–160.

MacEachren, A. M., Brewer, C. A., Pickle, L. W. (1998). Visualizing georeferenced data: representing reliability of health statistics. *Environment and Planning A*, 30, 1547–1561.

Reyes-Aldasoro, C. C., Bhalerao, A. (2006). The Bhattacharyya space for feature selection and its application to texture segmentation. *Pattern Recognition*, 39(5), 812–826.

Sun, M., Wong, D. W., Kronenfeld, B. J. (2014). A classification method for choropleth maps incorporating data reliability information. *The Professional Geographer*, 67(1), 72–83.

Sun, M., Wong, D., Kronenfeld, B. J. (2016). A heuristic multi-criteria classification approach incorporating data quality information for choropleth mapping. *Cartography and Geographic Information Science*, 1–13.

Xiao, N., Calder, C. A., & Armstrong, M. P. (2007). Assessing the effect of attribute uncertainty on the robustness of choropleth map classification. *International Journal of Geographical Information Science*, 21(2), 121–144.
<figure><img src="image_118.png" /></figure>
<figure><img src="image_119.png" /></figure>

# Scale effects in spatial accuracy
<figure><img src="image_120.png" /></figure>
# Getting the right spatial mix: optimising the size, type and location of renewable energy facilities

Alexis Comber*<fnref n="1" />, Jen Dickiie<fnref n="2" />

**Abstract**

Supply and demand modelling for facilities that require land related resources need to take into account the spatial distribution of the resource. Using a Scottish case study, this short paper presents a method extension that to the p-median problem that identifies optimal locations and combinations of different types and sizes of land based biomass renewable energy facilities. Whilst there are many decision models and tools for siting other types of renewable energy (solar, wind, hydro), supply and demand tools decisions around land based biomass renewable energy do not exist. The p-median extension and optimisation allows the trade-offs between different decisions around land, for example related to agriculture, food security, biodiversity, flood risks, to be evaluated. In the context of Scotland, this methodology supports the many policy agendas around community initiatives, agri-renewables, circular economy, supply chains, local food agendas, carbon sequestration and green infrastructures.

**Keywords**

Supply and demand; renewable energy; biomass; land use

# I INTRODUCTION

In much facility location analysis, the resources needed by a facility at potential locations are assumed to exist or that their supply is trivial. Much greater focus is given to identifying potential facility locations that meet the spatial distribution of demand, and where facility resources are considered these are typically done so from the perspective of transportation overheads, infrastructure and road networks.

Land based and biomass Renewable Energy (RE) facilities typically require consideration of the spatiality of supply as well as demand. The facilities require biomass feedstocks (inputs) that are usually derived from animal manures, crop residues, forestry / timber production. In the context of RE, as well as climate change mitigation, the objectives should also include net energy gains and minimal carbon impacts. Thus it is important to consider the resource catchments needed to supply RE facilities, rather than subsume financial transportation costs into site evaluation.

Recent work has extended the p-median problem to be able to handle the resource catchments needed to support facilities at potential locations in location allocation analyses (Comber et al., 2015) as well as demand. It has been further extended to locate multiple types of facility and to determine optima locations for pre-defined groups of different sized facilities which have

<footnote n="1">¹ University of Leeds, Leeds, LS2 9JT, UK</footnote>

<footnote n="2">² University of Stirling, Stirling, FK9 4LA, UK</footnote>

<footnote>*Corresponding author: a.comber@leeds.ac.uk</footnote>
essentially imposed an optimal packing dimension onto the p-median problem (Comber et al, 2016).

This paper further extends these developments to the p-median problem to identify the optimal mix of facilities. The problem remains a packing problem but the extension seeks to identify the optimal mix as well as spatial arrangement of facilitates. The objective function remains the same (to minimise distances to demand and to forbid resource catchment overlap), but the search space is extended from sets of predefined facilities (e.g. 3 large, 10 medium and 20 small) to identifying the best combination of facilities for the study area being considered. The context for this extension is the need to accommodate competing land use demands, for example food security, flood defence or biodiversity.

The spatial accuracy context of this work relates to the sensitivity of location-allocation algorithms and their associated evaluation functions to spatial consideration of resource catchments.

## II INITIAL RESULTS

Consider the location, catchments and distribution of potential biomass RE facilities in Figure 1a along with the underlying spatial distribution of biomass from forestry. This shows an optimal arrangement for 3 Combined Heat and Power (CHP) units, with different capacities (energy outputs) and different input requirements for a case study in north east Scotland. The total potential amount of annual biomass from forestry is 1,895,677 t yr⁻¹ and the mix of 3 20MW CHPs, 10 4MW CHPs and 30 1MW CHPs equates to a biomass feedstocks of 650,000 t yr⁻¹ or 34% of the available supply and the production of 130MW.

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

a)

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

b)
<figure><img src="image_123.png" /><figcaption>Figure 1: Optimal spatial arrangements for sets of biomass CHP facilities with the same overall generating capacity a) for a predefined mix of capacities, b) for a different mix of capacities, and c) the optimal, optimal mix</figcaption></figure>

The question addressed by this paper is whether a different spatial arrangement of different mix of facilities sizes, such as suggested in Figure 1b, would achieve the same energy generating capacity output but use the biomass resource more efficiently, ie with less energy expended on transportation.

To achieve this, the algorithm has been further extended to search through *all possible combinations* of the 3 types of facility to achieve the same capacity. The optimal spatial configuration is shown in Figure 1c and suggests that the following combination of CHPs to achieve 120MW and minimising resource transportation costs is as follows: 2 x 20MW, 17 x 4MW and 12 x 1MW.

## III DISCUSSION

There are a number of critical methodological issues and spatial accuracy considerations related to this work.

In terms of method development, this paper describes a further extension to the algorithm, first described in Comber et al., (2015), which accommodates the spatial distribution of the resources needed to supply the RE facility and satisfy the energy demand at any given location. It extended the p-median problem to prevent facility locations with overlapping resource
catchments from being selected. The first extension, described in Comber et al (2016), supported the allocation of a predefined number of different sizes of facility (3 small, 2 medium and 5 large, for example). However, this was based on a packing problem heuristic, which identified the location for $n_1$ largest facilities first, then the $n_2$ medium sized facilities and finally the $n_3$ largest facilities.

The extension reported here identifies the optimal mix of different sized facilities and considers all sizes of facility together. In each case the evaluation function was to minimise population distance to the resource supply. In this case this was done through a deterministic search. The resource catchment for each size of facility at each potential location was pre-computed and then the ‘best’ combination of $n_1$ 20MW $n_2$ 4MW and $n_3$ 1MW facilities was determined. This was possible because of the relatively small number of potential sets of each {$n_1, n_2, n_3$} combination. However for a larger study area, the number of potential combinations of facilities may preclude a deterministic search suggesting the need for search heuristics such as genetic algorithms, or perhaps more pertinently to this type of study grouping genetic algorithms (Falkenauer, 1998) which have been shown to more effectively identify optimal groups of facilities than standard genetic algorithm (Comber et al., 2011).

In terms of the accuracy of the results generated by this method, they are subject to the usual considerations in any location-allocation application, namely the extent to which the evaluation function matches the problem specification and the reliability of the input data.

Here, optimality (the evaluation function) was determined by the degree to which potential facility locations minimised distance to population centres (ie they were not weighted by the population at those centres). There are two obvious areas to refine this function: first, to include population, for example, to generate demand weighted distances, reflecting differences in the spatial distribution of demand. However, this study used Output Area centroids as the demand locations. These are the smallest areal unit over which UK census data are reported (Martin, 1998) and were constructed such that they contain broadly an equal number of people (~350). Perhaps more interestingly data describing more nuanced measures of demand for RE could be used, such as the Call Credit Green and Ethical geo-demographic classification of the UK. Second, the object of any planning and policy in relation to renewable energy should be to minimise net carbon costs and maximise net carbon gains. This is rarely the case in the RE literature, which much effort seems to be spent on modelling transport costs in RE facility location decisions (eg Sultana and Kumar, 2012; Panichelli and Gnansounou, 2008; Sliz-Szkliniarza and Vogt, 2012). This is plainly daft as it renders the result of any analysis meaningless if there are radical changes in the price of oil. The optimal selection of suitable sites for RE facilities using these kinds of paradigms is critical if land based biomass resources are to be efficiently and maximally used to support a diverse range of objectives including food and energy security as well as environmental protection.

In summary the current method allows locations for multiple sizes of facilities to be evaluated. The algorithm could be applied to select optimal sites for multiple types as well as sizes of renewable energy facilities: CHPs, anaerobic digesters, gasification units etc. This would support truly holistic, strategic regional planning and well as community level energy initiatives. The latter are increasingly being supported in Scotland (eg http://www.localenergyscotland.org/cares).

The next steps in this work are multiple. First, to consider network distances to resources rather Euclidean distances, to explore how asymmetric, amorphous catchments may be incorporated, allowing them to fill the available space between already selected sites and to consider different combinations of feedstocks – domestic, forest and agricultural and household wastes – would could be used to refine the results of the analysis. Second, to develop methods for searching more efficiently, for example through heuristics, the use of *net*
*energy gain* as an evaluation function, the spatial distribution of demand that is receptive to RE and whether network distances to resources improve the outcomes.

This paper presents a framework for doing this: compute all possible combinations of RE supply, identify the resource catchments needed for each of individual facility at each potential location, and then develop a search through this very highly dimensional decision space. Future work will consider the development and application of evaluation functions that minimise the distances (and net energy and carbon costs) that resources have to travel to supply the facility as well as more nuanced measures of potential demand.

This research was funded by the EPSRC SECURE project, Reference EP/M008347/1 (http://www.gla.ac.uk/research/az/secure/).

## References

Comber A, Dickie J, Elston D and Miller D (2016). What to put where? Extending the p-median problem to consider multiple facilities, multiple sizes and associated resource needs. Short paper AGILE 2016, Conference Helsinki, June 2016.

Comber A, Dickie J, Jarvis C, Phillips M and Tansey K. Locating bioenergy facilities using a modified GIS-based location-allocation-algorithm: considering the spatial distribution of resource supply. *Applied Energy*, 154: 309-316, 2015.

Comber, A.J., Sasaki, S., Suzuki, H. and Brunsdon, C., (2011). A modified grouping genetic algorithm to select ambulance site locations. *International Journal of Geographical Information Science*, 25(5): 807–823

Falkenauer, E., 1998. *Genetic algorithms and grouping problems*. London: John Wiley and Sons.

Martin, D., 1998. 2001 Census output areas: from concept to prototype. *Population Trends* 94, 19–24.

Panichelli L and Gnansounou E, (2008). GIS-based approach for defining bioenergy facilities location: A case study in Northern Spain based on marginal delivery costs and resources competition between facilities. *Biomass and Bioenergy*, 32: 289-300.

Sliz-Szkliniarz B and Vogt J, (2012). A GIS-based approach for evaluating the potential of biogas production from livestock manure and crops at a regional scale: A case study for the Kujawsko-Pomorskie Voivodeship. *Renewable and Sustainable Energy Reviews* 16: 752–763.

Sultana A and Kumar A, (2012). Optimal siting and size of bioenergy facilities using geographic information system, *Applied Energy*, 94: 192-201
# Impact of compositional and configurational data loss on downscaling accuracy

Amy E. Frazier<fnref n="1" />* and Peter Kedron<fnref n="2" />

*Corresponding author: amy.e.frazier@okstate.edu

## Abstract

Datasets collected at widely varying spatial scales are often merged to address questions related to global environmental change. Integrating data requires aggregation, which reduces data quality and introduces statistical biases collectively known as the Modifiable Areal Unit Problem (MAUP). These biases result from different forms of compositional and configurational data loss that occur during aggregation, but little is known about the relationship between data loss and MAUP biases for downscaling. This study uses the well-established process of landscape and surface metric scaling to examine how uncertainties related to the composition and configuration of land cover patterns propagate across scales when data are aggregated and ultimately impact downscaling results. Results suggest a link between compositional data loss and downscaling accuracy, particularly in the patch-based landscape paradigm. Further work is needed to determine if relationships exist between compositional and configurational data loss measures and downscaling error in the surface paradigm.

**Keywords:** scale and scaling; landscape ecology, spatial pattern metrics, heterogeneity, remote sensing

## I. Introduction

In the age of data-driven science, diverse datasets collected at widely varying spatial scales are increasingly being merged to address questions related to global environmental change. However, integrating data collected at different spatial scales requires aggregation, which reduces data quality and introduces statistical bias. These biases, collectively known as the Modifiable Areal Unit Problem (MAUP), result from different forms of compositional and configurational data loss that occur during aggregation. MAUP biases can be minimized when aggregated data have a lower degree of heterogeneity (Holt et al. 1996; Steel and Holt 1996), but beyond this recognition, studies on the relationship between forms of data loss and MAUP biases remain limited.

Building on recent theoretical and methodological advancements in spatial science, remote sensing, and landscape ecology, we investigate how MAUP biases propagate across resolutions and whether the spread of bias, in the form of compositional and configurational data loss, can be used to forecast downscaling accuracy. Spatial pattern scaling provides a useful platform from which to examine these issues because research has established that several land cover pattern metrics exhibit consistent and robust scaling relationships across resolutions (Wu 2004). Yet, when these scaling relationships are extrapolated (i.e., downscaled) to predict metrics at a finer resolution, large errors typically result (Frazier 2014), likely from MAUP-driven aggregation biases (Frazier 2014, 2015a).

Until recently, examination of this hypothesis was limited to hard-classified landscape images, but the emergence of sub-pixel remote sensing classification techniques that preserve greater heterogeneity than their traditional, pixel-based counterparts have improved our ability to quantify data loss and preserve landscape heterogeneity, thereby opening the door for renewed investigations. Simultaneously, the proliferation of surface metrics in landscape ecology provides a means from which to compare this hypothesis across two different landscape paradigms: patch-based and surface. We use the well-established process of landscape and surface metric scaling (Turner et al. 1989; Wu 2004; Frazier 2015b) to examine how uncertainties related to the composition and configuration of land cover patterns

<footnote n="1">¹ Oklahoma State University, Stillwater, Oklahoma, USA</footnote>

<footnote n="2">² Ryerson University, Toronto, Ontario, Canada</footnote>
propagate across scales when data are aggregated and ultimately impact downscaling results. Predetermination of the impact of these biases on downscaling may eventually allow assessment of whether a landscape is a satisfactory candidate for downscaling.

## II. Methods

Data were collected through a geographically stratified sampling of forest land cover in four forested ecoregions in the eastern United States (Omernik 1987). Within each ecoregion, we randomly sampled 125 20x20km plots from a continuous grid (Fig. 1a,b). We removed grid squares comprising urbanized areas greater than 500,000 people. We then clipped the national land cover map (NLCD), aggregated to parent classes including a ‘forest’ class, and the tree canopy cover (TCC) product, both 30m resolution, to sample boundaries. Each NLCD and TCC plot was then aggregated to 60, and the 60m raster aggregated to 120, 180, 240, 420, and 480m using majority rules (NLCD) or mean (TCC) aggregation (Fig. 1c).

We computed analogous patch-based and surface metrics to measure downscaling accuracy (Table 1) and a suite of metrics from each paradigm that measure landscape composition and configuration (McGarigal et al. 2009). Patch-based metrics were computed using Fragstats (McGarigal et al. 2012) and surface metrics using SPIP software (Image Metrology).

<figure><img src="image_124.png" /><figcaption>Figure 1. (a) Eastern U.S. forest ecoregions, (b) grid sampling scheme, and (c) NLCD aggregation.</figcaption></figure>

Reserving the 30m raster, a scalogram was plotted for PD and Sds for the seven coarse resolutions for each plot, and various scaling functions fit to the scalogram (i.e., power law, first-, second-, and third-order polynomials). Wu (2004) and Frazier (2015b) demonstrate these curve types accurately model metric behavior across scale. To generate a measure of MAUP-induced error, those curves were downscaled to predict each metric value at 30m and those predicted values compared to true 30m PD and Sds values through a measure of relative error (Eq. 1), where $M_p$ is the predicted metric from the scaling function, and $M_t$ is the true metric value. Low $E_{rel}$ values indicate the scaling function accurately predicts the metric at a finer resolution. The curve type producing the lowest average $E_{rel}$ for each metric was selected as the candidate model for further examination.
<table><caption>Table 1. Analogous metrics for measuring downscaling accuracy and data loss in the two paradigms.</caption><thead><tr><th>Patch-Based</th><th>Surface</th></tr></thead><tbody><tr><td><em>Downscaling metrics</em></td><td></td></tr><tr><td>Patch Density (PD)</td><td>Peak Density (Sds)</td></tr><tr><td><em>Compositional data loss metrics</em></td><td></td></tr><tr><td>Total area of forest (AREA)</td><td>--</td></tr><tr><td>Largest patch index (LPI)</td><td>Maximum peak height (Sp)</td></tr><tr><td><em>Configurational data loss metrics</em></td><td></td></tr><tr><td>Edge Density (ED)</td><td>Peak-Peak range height (Sy)</td></tr><tr><td>Percent of Like Adjacencies (PLADJ)</td><td>Moran's I index (Smi)</td></tr><tr><td>Mean Fractal Dimension (FRAC)</td><td>Surface fractal dimension (Sfd)</td></tr></tbody></table>

To capture compositional data loss within the patch-based paradigm, we computed the linear rate of change of the data loss metrics (Table 1) across the seven coarse resolutions. Finally, we used conventional, ordinary least squares regression to measure the influence of data loss on downscaling accuracy. Expressed as Equation 2, $\hat{y}_i$ is the estimated value of the dependent variable $E_{rel}$ for plot $i$, $\beta_0$ represents the intercept, and $\beta_{config(k)}$ and $\beta_{comp(j)}$ are coefficients for the independent variable $x_{ik}$ and $x_{ij}$, our measures of either compositional or configurational data loss, and $\varepsilon_i$ represents the error term. Independent regressions were completed for each ecoregion for each paradigm and a combined regression for all regions for each paradigm.

$$E_{rel}(\%) = |(M_p - M_t)/M_t| * 100 \quad (1)$$

$$\hat{y}_i = \beta_0 + \sum_k \beta_{config(k)} x_{ik} + \sum_k \beta_{comp(j)} x_{ij} + \varepsilon_i \quad (2)$$

### III. Results and Discussion

Across the four ecoregions, a third-order polynomial model performed best for downscaling PD, and power law performed best for Sds, which is consistent with recent findings (Frazier 2015a). PD mean $R^2$ values were consistently 0.97. Model fit improved with the surface paradigm, and $R^2$ values were >0.99 for Sds (Table 2). Relative error statistics (Table 2) show average downscaling accuracies ranged from 28.3 to 39.8 for PD, and 27.7 to 39.3 for Sds. Thus, there is little correlation between model fit ($R^2$) and downscaling accuracy ($E_{rel}$), but results confirm prior findings that even when scaling relationships are strong and model fits are high, metric downscaling of aggregated data is not particularly accurate (Frazier 2014).

<table><caption>Table 2. Scaling model fit ($R^2$) and downscaling relative error ($E_{rel}$) across ecoregions and paradigms.</caption><thead><tr><th rowspan="2">Ecoregion</th><th colspan="4">PD</th><th colspan="4">Sds</th></tr><tr><th>8.1</th><th>8.3</th><th>8.4</th><th>8.5</th><th>8.1</th><th>8.3</th><th>8.4</th><th>8.5</th></tr></thead><tbody><tr><td>Mean R<sup>2</sup></td><td>0.97</td><td>0.97</td><td>0.97</td><td>0.97</td><td>0.998</td><td>0.996</td><td>0.997</td><td>0.998</td></tr><tr><td>Std. Dev. R<sup>2</sup></td><td>0.01</td><td>0.01</td><td>0.02</td><td>0.01</td><td>0.002</td><td>0.004</td><td>0.005</td><td>0.002</td></tr><tr><td>Mean E<sub>rel</sub></td><td>28.3</td><td>31.8</td><td>39.8</td><td>28.8</td><td>27.7</td><td>39.3</td><td>34.4</td><td>28.1</td></tr><tr><td>Std. Dev. E<sub>rel</sub></td><td>10.7</td><td>9.85</td><td>14.9</td><td>7.94</td><td>8.7</td><td>13.1</td><td>16.1</td><td>10.3</td></tr><tr><td>Sample (n)</td><td>124</td><td>125</td><td>125</td><td>125</td><td>124</td><td>125</td><td>125</td><td>125</td></tr></tbody></table>

Tables 3 and 4 summarize the results of OLS regression for the two paradigms. For the patch paradigm, the rate of AREA loss was consistently a significant predictor of relative downscaling error. ED was also significant for most models. Across all ecoregions, PLADJ was also significant. Model fit values ranged from 0.41 to 0.596. For the surface paradigm, models were poorly fit. No measures of compositional or configurational data loss were significant across all regions. Variance inflation factors indicated possible
collinearity of variables. These findings suggest further examination of alternative metrics is necessary. Of the two models with reasonable fit, the rate of loss of maximum peak height (Sp), a measure of compositional loss, was the best predictor of downscaling accuracy.

<table><caption>Table 3. Impacts of data loss on downscaling accuracy of PD (patch paradigm).</caption><thead><tr><th>Ecoregion</th><th>8.1</th><th>8.3</th><th>8.4</th><th>8.5</th><th>Total</th></tr></thead><tbody><tr><td>Intercept</td><td>18.174***</td><td>34.948***</td><td>40.101***</td><td>36.102***</td><td>40.802***</td></tr><tr><td>AREA</td><td>0.666*</td><td>1.250***</td><td>2.243**</td><td>1.914***</td><td>1.361***</td></tr><tr><td>LPI</td><td>-0.243</td><td>0.678**</td><td>-0.644</td><td>-0.602</td><td>0.029</td></tr><tr><td>ED</td><td>-2.566***</td><td>1.351***</td><td>0.984</td><td>0.934**</td><td>1.609***</td></tr><tr><td>PLADJ</td><td>0.625</td><td>1.173***</td><td>2.745</td><td>0.558</td><td>1.184***</td></tr><tr><td>Diagnostics</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>R<sup>2</sup></td><td>0.478</td><td>0.596</td><td>0.467</td><td>0.410</td><td>0.454</td></tr><tr><td>Sample (n)</td><td>124</td><td>125</td><td>125</td><td>125</td><td>499</td></tr></tbody></table>

<table><caption>Table 4. Impacts of data loss on downscaling accuracy of Sds (surface paradigm).</caption><thead><tr><th>Ecoregion</th><th>8.1</th><th>8.3</th><th>8.4</th><th>8.5</th><th>Total</th></tr></thead><tbody><tr><td>Intercept</td><td>34.375***</td><td>52.849***</td><td>41.33***</td><td>30.780***</td><td>39.288***</td></tr><tr><td>Sp</td><td>1.083</td><td>14.762***</td><td>10.790***</td><td>1.251</td><td>5.371***</td></tr><tr><td>Sy</td><td>-0.281</td><td>-6.769</td><td>-0.335</td><td>-0.114</td><td>-1.093</td></tr><tr><td>Sfd</td><td>-1.414***</td><td>-2.969***</td><td>0.361</td><td>-0.716</td><td>-1.140***</td></tr><tr><td>Diagnostics</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>R<sup>2</sup></td><td>0.064</td><td>0.246</td><td>0.264</td><td>0.024</td><td>0.112</td></tr><tr><td>Sample (n)</td><td>124</td><td>125</td><td>125</td><td>125</td><td>499</td></tr></tbody></table>

## IV. Conclusions

Results suggest a link between compositional data loss and downscaling accuracy, particularly in the patch-based paradigm—the greater the rate of AREA loss during aggregation, the more difficult it becomes to predict metric values at a finer resolution. Further work is needed to determine if any relationships exist between compositional and configurational data loss measures and downscaling error in the surface paradigm. An initial step would be to further establish correspondences between patch and surface metrics within these and other ecoregions. Future work should also address whether ecological characteristics impact the success these measures.

**Acknowledgements:** This work is supported by a grant to the authors from the U.S. National Science Foundation (#1561021) for “Data Complexity and Spatial Scaling: Prediction Accuracy and Implications for Emerging Landscape Paradigms”.

## References:

* Frazier, A. E. (2014) A new data aggregation technique to improve landscape metric downscaling. *Landscape Ecology*, 29(7), 1261-1276

* Frazier, A.E. (2015a) Landscape heterogeneity and scale considerations for super-resolution mapping. *International Journal of Remote Sensing*,

* Frazier, A. E. (2015b) Surface metrics: scaling relationships and downscaling behavior. *Landscape Ecology*, 1-13

* Holt, D., Steel, D.G., Trammer, M., & Wrigley, N. (1996) Aggregation and Ecological Effects in Geographically Based Data. *Geographical Analysis* 28(3):244-261

* McGarigal, K., Tagil, S., & Cushman, S. A. (2009). Surface metrics: an alternative to patch metrics for the quantification of landscape structure. *Landscape Ecology*, 24(3), 433-450

<footnote>*p < 0.10, **p<0.05, ***p<0.01; FRAC was not significant in any models</footnote>

<footnote>*p < 0.10, **p<0.05, ***p<0.01; SMI was not significant in any models</footnote>
McGarigal, K., SA Cushman, and E Ene. 2012. FRAGSTATS v4: Spatial Pattern Analysis Program for Categorical and Continuous Maps. Computer software program produced by the authors at the University of Massachusetts, Amherst. Available at the following web site: http://www.umass.edu/landeco/research/fragstats/fragstats.html

Omernik, J.M. (1987) Ecoregions of the conterminous United States. *Annals of the Association of American Geographers* 77(1), 118-125.

Turner, M. G., O'Neill, R. V., Gardner, R. H., & Milne, B. T. (1989). Effects of changing spatial scale on the analysis of landscape pattern. *Landscape ecology*, *3*(3-4), 153-162

Steel, D. G., & Holt, D. (1996) Analysing and adjusting aggregation effects: the ecological fallacy revisited. *International Statistical Review/Revue Internationale de Statistique*, 39-60

Wu, J. (2004) Effects of changing scale on landscape pattern analysis: scaling relationships. *Landscape Ecology*, 19:125-138.
# A hierarchical scale setting strategy for improved segmentation performance using very high resolution images

Taskin Kavzoglu¹*, Merve Yildiz Erdemir¹

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

(kavzoglu, m.yildiz)@gtu.edu.tr

## Abstract

Land use/land cover (LULC) classification is a specific implementation to define terrain features to the closest real world object. Object-based image analysis (OBIA) has been proved to improve classification accuracy, particularly for very high resolution remotely sensed images. In this study, multiresolution segmentation algorithm was utilized in the image segmentation process using a pan-sharped Quickbird-2 image. Segmentation scales were determined by widely-used estimation of segmentation parameter (ESP-1) tool that produces rate of change graph (LV-RoC) in terms of local variance of the image. In this study, the LV-RoC graph of the image was evaluated to determine optimal scale values ranging from fine to coarse levels. An attempt was made to estimate optimal scale parameter for an image considering not only single-scales but also multi-scales for an image using a hierarchical scale setting strategy. Nearest neighbour classifier was used on single-scale segmented images and fuzzy classifier employing membership functions was used on multi-scale segmented image. Equal numbers of pixels for each class were randomly selected to estimate accuracy metrics (i.e. overall accuracy and kappa coefficient). The differences in classifier performances (~ 6%) were statistically significant according to McNemar's test. It was found that the proposed strategy has a great potential for LULC classification using very high resolution imagery.

**Keywords:** Object-based classification, Segmentation, Scale parameter, ESP, Quickbird-2

## I Introduction

Parallel with the increasing use of very high spatial resolution images, object based classification has become more advanced (Kim et al., 2011). As stated by Blaschke et al. (2014), the object based classification technique that uses image objects instead of pixels has been widely used in last decade primarily applied to very high resolution images. Object-based image analysis (OBIA) offers several unique advantages in comparison to pixel-based classification, one of which is the removal of the so-called salt-and-pepper effect. In OBIA, shape, context and textural information of image objects are considered instead of individual pixels. At the same time, each object is considered both by its spectral, shape or texture features, and by its unique neighbours, its sub- and super-objects (Benz et al., 2004). The recent availability of OBIA provides various opportunities for sensitive and detailed LULC classification.

OBIA is generally applied in three stages: image segmentation, classification and accuracy assessment. Image segmentation process is the first and crucial step to define image objects. Although there exists several techniques, multiresolution segmentation introduced by Baatz and Schäpe (2000) has been one of the most commonly used image segmentation algorithm. Multiresolution segmentation produces highly homogeneous image objects in arbitrary resolution on different types of data and is applied to many problems in the field of remote
sensing (Baatz and Schäpe, 2000). It includes the setting of three major parameters, namely scale, shape and compactness. Setting of these parameters is of critical importance for the performance of a subsequent classification process (Witharana and Civco, 2014). Of these parameters, it is agreed that scale parameter is the most important one as it determines the objects size and thereby affects largely the resulting thematic map accuracy (Kim et al., 2011; Myint et al., 2011). The optimal value for scale parameter varies depending on the structure of the study area, land cover types and scale of the imagery (Myint et al., 2011).

Results of classification accuracy are affected by quality of segmentation process (Marpu et al., 2010). For this reason, many studies have already been carried out to determine optimal scale parameter. In a previous research, Kavzoglu and Yildiz (2014) endeavoured to prepare guidelines for the selection of scale parameter. Besides selection of the scale parameter without relying on any method, evaluation methods of segmentation quality can be divided into two categories as supervised (Clinton et al., 2010) and unsupervised (Johnson et al., 2011; Gao et al. 2011).

In recent years, researches mainly focused on the use of statistical methods for determination of optimal scale parameter in segmentation (Witharana and Civco, 2014). However, only a few studies have led to automatically determining results that produce fast and applicable solutions. Woodcock and Strahler (1987) utilized local variance (LV) graphs to ascertain the spatial structure of each type of image. Then, the estimation of scale parameter (ESP) tool was developed by Drăguț et al. (2010) for automated detection of scale parameter. The tool has been recently extended to multi-scale analysis based on single layers (Drăguț et al., 2014). With this tool, LV was generated for each scale value from an image. Furthermore, a graph is generated using LV of image and rate of change (RoC) values of LV. Thus, LV-RoC graph is obtained on each scale step in an image considering a single layer.

The purpose of this study is to determine optimal scale parameter(s) of the image using the LV-RoC graphs produced by the ESP tool and perform usage of single- and multi-scale analysis on classification process. Nearest neighbour classifier was employed in single-scale classification process and membership functions were used in multi-scale classification process. The differences in classifier performances were analysed with McNemar's test, which is a statistical test for differences.

## II Study area and Dataset

The study area chosen for this research covers approximately 80 ha area located in the Yomra district of Trabzon province, Turkey. A multi-spectral pan-sharpened Quickbird-2 satellite image having four spectral bands at 0.6 m spatial resolution acquired in May 2008 was used for the study site. It has high mountainous terrain covered by forest, sea and urban classes. For the purpose of evaluation, image was clipped to a subset of 1500x1500 pixels (Figure 1). eCognition Developer (v9.1), a widely-used object based image processing software, was used for all segmentation and classification experiments in this study.

The study area is mainly covered by eleven land use/cover features, namely bare ground, concrete surface, forest, asphalt road, gravel, pasture, shadow, water, red, blue and white roof.

## III Image segmentation

Image segmentation is the first and crucially important process of object based classification, since the initial objects are created on this stage. The most widely used algorithm in image segmentation is multiresolution segmentation introduced by Baatz and Schäpe (2000). It is a bottom up region merging technique based on local homogeneity criteria. It merges neighbour and similar spectral pixels into image objects. Size of image objects is adjusted in conjunction
with segmentation parameters including scale, shape, compactness and band weights. Scale parameter sets the image size and considering most crucial parameter. In general, higher values for the scale parameter result in larger image objects, smaller values in smaller ones (Yildiz et al., 2012).

<figure><img src="image_125.png" /><figcaption>Figure 1. Location of the study area, Trabzon province of Turkey.</figcaption></figure>

Estimation of scale parameter (ESP) tool was developed for automated detection of scale parameter by Drăguț et al. (2010) who states that “the ESP tool iteratively generates image objects at multiple scale levels in a bottom-up approach and calculates the local variance (LV) for each scale image”. Rate of change (RoC) values of LV is determined for each scale levels using Eq. 1. The RoC is calculated as:

$$RoC = \left[ \frac{L-(L-1)}{L-1} \right] * 100 \qquad (1)$$

where L is local variance at the target level and L-1 is local variance at next lower level. The peaks in the LV-RoC graph designates the object levels at which the image can be segmented in the most appropriate way, according to data properties at the scene (Drăguț et al., 2010).

## IV Classification

Classification process which is the next step after the image segmentation uses segmented images. Two types of classifiers were applied in this study. Nearest neighbour classifier uses a set of samples of different classes in an attempt to assign class values to a segmented object. The procedure contains two major steps: teaching the system by giving it certain image objects as samples and classifying image objects in the image object domain based on their nearest sample neighbours. Membership functions allow describing the relationship between feature values and the degree of membership to a class using fuzzy logic. It can be defined by the degree of membership; for example, any value between one and zero (Definiens, 2008). In other words, Myint et al. (2011) states that “the membership function describes intervals of feature characteristics that determine whether the objects belong to a particular class or not”.
## V Results

This study aims to investigate optimal scale parameter(s) for single- and multi-scale segmentation. A pan-sharpened Quickbird image includes many land use/cover features such as urban, sea, forest and gravel was used in this study. For detecting optimal scale parameter(s), local variances for each of scale parameter were estimated using ESP tool in the image. As can be seen in Figure 2, the prominent peaks are apparent in the LV-RoC graph created using Eq. 1. Peaks at 26, 41, 53 and 77 on the graph appeared to be suitable values. The four scale levels were denoted S_ESP1, S_ESP2, S_ESP3 and S_ESP4, where S_ESP1 represented the finest object scale and S_ESP4 the coarsest scale. These scale parameter values specified by LV-RoC graph was individually applied to multiresolution segmentation for single-scale segmentation. Furthermore, all scale values were used in multi-scale segmentation using a hierarchical strategy (M_ESP). It should be noted that shape and compactness parameters were kept constant as 0.1 and 0.5, respectively.

<figure><img src="image_126.png" /><figcaption>Figure 2. LV-RoC graph of the image. Scale values of 26, 41, 53 and 77 were selected as potential optimal scale values.</figcaption></figure>

For single-scale segmentation, segmented images obtained from the estimated scale parameters were classified using nearest neighbour classifier. Membership function classifier was applied on multi-scale segmentation and different segmentation scale levels were applied to detect feature classes. Each scale level was generated within a hierarchy, where scale parameters of 77, 53, 41 and 26 were set as scale levels of 1, 2, 3 and 4, respectively. Land use/cover classes were classified at mentioned levels. To be more specific, while water and gravel classes were extracted using scale level 1, blue roof and forest was extracted with scale level 2, asphalt road and white roof with scale level 3, and other classes with scale level 4.

Classifications were achieved using test datasets on the basis of contingency matrices. It should be noted that test datasets were prepared using random pixel selection. Equal numbers of samples for each class (700 pixels) were selected. For the assessment of classification results, overall classification accuracy (OA) and Kappa statistics were computed from the contingency matrices (Table 1).
<table><caption>Table 1. Overall accuracies of single- and multi-scale segmentations.</caption><thead><tr><th></th><th>Scale Unit</th><th>Overall Accuracy (%)</th><th>Kappa</th></tr></thead><tbody><tr><td rowspan="4">Single-scale</td><td>S_ESP1</td><td>86.10</td><td>0.847</td></tr><tr><td>S_ESP2</td><td>84.65</td><td>0.831</td></tr><tr><td>S_ESP3</td><td>81.18</td><td>0.793</td></tr><tr><td>S_ESP4</td><td>80.39</td><td>0.784</td></tr><tr><td>Multi-scale</td><td>M_ESP</td><td>91.56</td><td>0.907</td></tr></tbody></table>

The overall classification accuracies for single-scale parameters of 26, 41, 53 and 77 were estimated as 86.10%, 84.65%, 81.18% and 80.39%, respectively, suggesting superior performance when fine level scale values were used. For multi-scale segmentation applied through a particular strategy with fuzzy approach, a classification accuracy of 91.56% was achieved. As can be seen from Table 1, classification accuracies on single-scale segmentation range between 80% and 86% in terms of overall accuracy. All classified image results are shown in Figure 3.

<figure><img src="image_127.png" /><figcaption>Figure 3. Single-scale classification results for (a) S_ESP1, (b) S_ESP2, (c) S_ESP3, (d) S_ESP4 and multi-scale classification result with hierarchical scale (e) M_ESP.</figcaption></figure>

Although some classes (e.g., water and forest) were delineated easily, there were misclassification problems between some classes (e.g., gravel, bare soil, asphalt road and red building roofs). It should also be noted that low accuracies, compared with the accuracies produced for fine scale selection (i.e. 26 of scale value), were achieved in the classification of coarser scale selection (i.e. 53 and 77 of scale values). The reason could be related to training objects containing two or more classes in same object. Overall, the highest classification result
was estimated by multi-scale segmentation using membership functions. One reason for this can be that scale levels separated to identify objects of different classes and each land use feature was represented according to its spatial features. Another reason may be related to using membership functions in classification process. Thus, some classes were easily defined by object features such as spectral bands, NDVI, band ratio, brightness etc.

McNemar's test based on $\chi^2$ distribution is a popular non-parametric test that is generally performed to compare the classification errors (Kavzoglu et al., 2015). In this study, McNemar's test was implemented to determine statistical significance of the differences between classifications according to scale selections (Table 2). It was found that there is statistical significance among all five classification results at 95% confidence interval ($\chi^2_{0.05} = 3.84$).

<table><caption>Table 2. McNemar's test result for single- and multi-scale segmentation.</caption><thead><tr><th></th><th>S_ESP1</th><th>S_ESP2</th><th>S_ESP3</th><th>S_ESP4</th><th>M_ESP</th></tr></thead><tbody><tr><th>S_ESP1</th><td>-</td><td>210.72</td><td>463.59</td><td>539.13</td><td>930.44</td></tr><tr><th>S_ESP2</th><td></td><td>-</td><td>60.84</td><td>101.82</td><td>324.38</td></tr><tr><th>S_ESP3</th><td></td><td></td><td>-</td><td>17.19</td><td>132.64</td></tr><tr><th>S_ESP4</th><td></td><td></td><td></td><td>-</td><td>79.78</td></tr><tr><th>M_ESP</th><td></td><td></td><td></td><td></td><td>-</td></tr></tbody></table>

## VI Conclusion

Determination of the optimal scale parameter of an image is currently an important research for object based image analysis. The more representative segments are introduced to a segmentation process, the more accurate classification results can be produced. Land use/cover features of used imagery are usually composed of a complex combination of buildings, roads, water area, trees, and pasture. Scale parameter depends on certain structural features such as spatial resolution of image, land use/cover characteristics and study area.

In this study, LV-RoC graph derived from ESP tool was utilized for determine optimal scale parameter. Thus, each scale parameter was separately applied as single-scale segmentation and hierarchical usage of scale parameters was used in multi-scale segmentation. Nearest neighbour and membership function classifiers were applied on single- and multi-scale segmentation to perform the land cover classification, respectively.

The scale parameter of 26 was selected as an optimum scale parameter of single-scale segmentation after estimated segmentation levels, determined through LV-RoC graph analysis, were individually tested for multiresolution segmentation process. However, scale parameter of 77 produced the lowest accuracy of classification. The reason may be related to the fact that the large image objects produced for the coarse scale, thus several land cover properties were represented in a single object. Particularly, about 6% classification accuracy improvement was achieved when considering the optimal scale parameters in multi-scale segmentation process compared to single-scale segmentation. Consequently, overall accuracies of the classification based on multi-scale segmentation were higher than single-scale segmentation. As a result, it can be concluded that multi-scale segmentation shows better performance compared to single-scale segmentation. Main reason for poor performance of single-scale may be related to construction of objects including several land cover features, especially in coarse scale values. In summary, the analyses conducted in this study clearly showed that the multi-scale segmentation based on hierarchical scale setting strategy is an effective approach for classifying very high resolution imagery.
## References

Baatz M., Schäpe A. (2000). Multiresolution Segmentation: an optimization approach for high quality multi-scale image segmentation. In: Strobl, J., Blaschke, T., Griesebner, G. (Eds.), *Angewandte Geographische Informations-Verarbeitung XII*. Wichmann Verlag, Heidelberg, pp. 12-23.

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

Blaschke T., Hay G.J., Kelly M., Lang S., Hofmann P., Addink E., Queiroz Feitosa R., van der Meer F., van der Werff H., van Coillie F., Tiede D. (2014). Geographic Object-Based Image Analysis – Towards a new paradigm. *ISPRS Journal of Photogrammetry and Remote Sensing*. 87, 180-191.

Clinton N., Holt A., Scarborough J., Yan L., Gong P. (2010). Accuracy assessment measures for object-based image segmentation goodness. *Photogrammetric Engineering and Remote Sensing*, 76(3), 289-299.

Definiens 2008. eCognition Developer 7.0 User Guide. Definiens AG, Munich, Germany.

Drăguț L., Csillik O., Eisank C., Tiede D. (2014). Automated parameterisation for multi-scale image segmentation on multiple layers. *ISPRS Journal of Photogrammetry and Remote Sensing*. 88, 119-127.

Drăguț L., Tiede D., Levick S. R. (2010). ESP: a tool to estimate scale parameter for multiresolution image segmentation of remotely sensed data. *International Journal of Geographical Information Science*, 24(6), 859-871.

Gao Y., Mas J.F., Kerle N., Navarrete Pacheco J.A. (2011). Optimal region growing segmentation and its effect on classification accuracy. *International Journal of Remote Sensing*, 32(13), 3747-3763.

Johnson B., Xie Z. (2011). Unsupervised image segmentation evaluation and refinement using a multi-scale approach. *ISPRS Journal of Photogrammetry and Remote Sensing*, 66(4), 473-483.

Kavzoglu T., Yildiz M. (2014, Sept.29- Oct. 2). Parameter-based performance analysis of object-based image analysis using aerial and Quikbird-2 images. In: *ISPRS Annals of Photogrammetry, Remote Sensing and Spatial Information Sciences*, II-7, pp.31-37.

Kavzoglu T., Colkesen I., Yomralioglu T. (2015). Object based classification with rotation forest ensemble learning algorithm using very high resolution WorldView-2 image. *Remote Sensing Letters*, 6(11), 834-843.

Kim M., Warner T. A., Madden M., Atkinson D. S. (2011). Multi-scale GEOBIA with very high spatial resolution digital aerial imagery: scale, texture and image objects. *International Journal of Remote Sensing*, 32(10), 2825-2850.

Marpu P. R., Neubert M., Herold H., Niemeyer I. (2010). Enhanced evaluation of image segmentation results. *Journal of Spatial Science*, 55(1), 55-68.

Myint S.W., Gober P., Brazel A., Grossman-Clarke S., Weng Q. (2011). Per-pixel vs. object-based classification of urban land cover extraction using high spatial resolution imagery. *Remote Sensing Environment*, 115(5), 1145-1161.

Witharana C., Civco D. L. (2014). Optimizing multi-resolution segmentation scale using empirical methods: Exploring the sensitivity of the supervised discrepancy measure Euclidean distance 2 (ED2). *ISPRS Journal of Photogrammetry and Remote Sensing*. 87, 108-121.

Woodcock C.E., Strahler A.H. (1987). The factor of scale in remote sensing. *Remote Sensing Environment*. 21, 311-332.

Yildiz M., Kavzoglu T., Colkesen I., Sahin E.K. (2012, July 10-13). An assessment of the effectiveness of segmentation methods on classification performance, In: *10th International Symposium on Spatial Accuracy Assessment in Natural Resources and Environmental Sciences*, Florianópolis, SC, Brazil, pp. 133-138.
# Exploring Accurate Spatial Downscaling using Optimization

Michael Poss*3 and Didier Josselin¹,²

¹UMR ESPACE 7300, CNRS, Université d'Avignon, France

²LIA, Université d'Avignon, France

³UMR CNRS 5506 LIRMM, Université de Montpellier, France

*Corresponding author: michael.poss@lirmm.fr

## Abstract

This paper proposes to combine several downscaling processes based on expert knowledge and objective functions from Operations Research to fill a table of social data about illiteracy in the USA during the Thirties.

## Keywords

Downscaling, social data, Operations Research, optimization, data about illiteracy, USA, 1930

# I THE PROBLEM OF SPATIAL DOWNSCALING ON CENSUS DATA

This paper deals with the well-known downscaling problem (Bierkens et al. (2000)) related to the ecological inference fallacy (King (1997)), change of support problem (King et al. (2004)) or aggregation problem (Josselin et al. (2008)) in general. Many applications and research fields (environment, social science, geography) are concerned and this problem notably impairs census data analysis (Gehlke and Biehl (1934); Robinson (1950); Tranmer and Steel (1998)).

## 1.1 Example

Let us consider that we have information (attribute values) for a complete statistical population of a territory. We want to disaggregate this information into a set of smaller sub-areas that compose the whole territory, at a finer scale. To do so, we need to make assumptions on the way we distribute the values into the spatial partition elements. Fortunately, we also have a few relevant information about the complete population and about individuals aggregates, for all the categories of the variable, and also in each sub-area. But we do not know the exact distribution of the category-specific values in each area. We have to make assumptions to estimate it. This is illustrated in the Table 1 with an example about housing.

<table><caption>Table 1: 14 millions of people from 3 French regions to distribute in 9 cells of the table according to different types of housing.</caption><thead><tr><th>Region</th><th>House owner</th><th>Tenant</th><th>Free accommodation</th><th>Inhabitants (millions)</th></tr></thead><tbody><tr><td>PACA</td><td>?</td><td>?</td><td>?</td><td>5</td></tr><tr><td>Rhône-Alpes</td><td>?</td><td>?</td><td>?</td><td>6.5</td></tr><tr><td>Centre</td><td>?</td><td>?</td><td>?</td><td>2.5</td></tr><tr><td>Quantity</td><td>8</td><td>5.6</td><td>0.4</td><td>14</td></tr></tbody></table>

Depending on the number of sub-areas (e.g. regions) and classes of the attribute (e.g. housing types), there exist many ways to fill the partition of the Table 1. At this point, those are all
equivalent solutions to the problem since they respect the constraints of the (known) summed values in the row and column margins. For instance, the most simple table (that would correspond to a freedom degree of only one in a contingency table) with only 20 individuals and 4 cells, i.e. 2 categories for both variables, has already 5 possible solutions.

To find the solution that fits the reality the best way –in the unreachable case of knowing the complete table content, i.e. how many people really live in a given region with a given type of housing–, we sometimes need to use statistical tools or to fix complementary assumptions related to complementary knowledge.

However, we are not completely deprived of means. What is interesting here is that we can firstly set simple constraints based on aggregated information. Indeed, we know that:

* the sum of the values in every line or column should be equal to the corresponding margin values;

* the sum of the margin values should be equal to the total number of individuals.

To find what is the "good" solution, then we have to fix an objective to maximize or minimize, according to assumptions from experts. This is the main idea of this paper. On the one hand, getting reliable information using census is very costly and cannot be made very frequently. On the other hand, there are experts in social science who may know a lot about socio-spatial characteristics of people. It can strongly help in making assumptions to drive disaggregation process, although we know that we definitely cannot find the real and exact distribution of the data at a given time. However, it seems anyway better than reading in a crystal ball.

In this paper, we propose an approach mixing statistics and optimization to find an optimal data distribution in a table, according to aggregated constraints and expert knowledge. Our objectives are the following:

* getting a better accuracy in statistical data due to a downscaling process driven by a statistical criterion optimized using a mathematical solver;

* finding the solution the most related to a given statistical criterion or to expert knowledge translated to a new matrix of expected values;

* assessing the quality of the solution when taking into account uncertainty in downscaling estimations, providing *min* and *max* values in each cell of the contingency table.

This method is tested on data studied by Robinson (1950), for which we know the complete matrix. Thence it enables to compare the results of our method to the real observed data.

## 1.2 Our data set: USA census data about illiteracy in 1930 used by Robinson (1950)

We use the data about the illiteracy in the population of USA in 1930, which are reference data in the field of sociology <fnref n="1" />. They include the number and the percentage of illiterate people in each State and in regions, by population color and nativity. We study four groups of population:

* *native white with native parentage* people;

* *native white with foreign or mixed parentage* people;

* *foreign born white* people;

* *black* people.

In this paper, optimization methods are used for the whole matrix *M* and results are provided and only discussed about black people in USA in 1930. Indeed, we test a complementary hypothesis on the role of history (slavery in confederate states before civil war) in black population illiteracy (for instance, see map provided in the Figure 1).

<footnote n="1">¹http://www.ru.nl/sociology/mt/rob/downloads/</footnote>
Part (%) of black illiterate people in the population of USA in 1930

<figure><img src="image_128.png" /><figcaption>Figure 1: Illiterate black people in USA in 1930: a location related to history (unified vs confederate states in East-Southern).</figcaption></figure>

Using these data, Robinson (1950) noticed that the correlation between illiteracy rate of black people and the different regions was close to 1. But if we consider the statistics at another more accurate scale (the States), it goes down to 0.2. This lead Robinson to define the *ecological fallacy* in social science and more generally on census data. On a practical point of view, it seems indeed very difficult to apply a knowledge learnt using the same statistics provided at a certain aggregated level, on a more disaggregated and finer level. That is the purpose of our proposition: trying to improve local estimation accuracy and downscaling process by modifying the data according to hypothesis and new constraints from aggregated statistics and from expert knowledge. Because we know the complete information about illiteracy in each type of population and in all the States in USA in 1930, we can compare and discuss the results of our method to the real data of the census. In this first stage, we fix the objective function criteria and the expert assumptions by ourselves to test the method.

## II A METHOD TO DOWNSCALE CENSUS DATA

### 2.1 Downscaling in areal units

In geography, it is very common to disaggregate spatial data into more accurate layers. Usually, a hypothesis of proportionality is applied. The assumption is the following: space is somehow isotropic and densities are the same all over geographical space or inside a given type of characterized area. This corresponds to the *H0* hypothesis in *Chi²* statistics. For instance, as in a contingency table, the expected number of individuals $\hat{x}_{ij}$ in a cell of the table is estimated by the product of the margin values of the corresponding column $x_i$ and row $x_j$, divided by the total number of individuals $x_\cdot$. Here we consider that illiteracy and regions where individuals are located are independent statistical variables. If there are significant differences in the number of individuals between regions or between classes of the variable, this can lead to very various values in the cells of the table.
Another way to compute downscaling may be to find the solutions minimizing a certain statistical criterion. For instance, variance of the values in the matrix $M$ of the table can be minimized. In this case, experts consider a kind of homogeneity in downscaling, i.e. that cell values should not differ too much from each others whatever the margin value weights. This approach is slightly different from the previous one, because a global statistical criterion is now considered. Moreover, it is possible to minimize many different functions $f(\theta)$ on such a table, according to several $L_p$-norms as we shall see in the next sections.

In an analog way, experts can use information about sub-areas peculiarities to provide probabilities or ranges (probability $min$ and $max$) in each cell. Then the problem leads to find the solution that is the closest to the probability matrix $M$ given by the expert. Practically, experts can provide bounds in which (s)he expects the values to be, in each cell of $M$. Whatever the function to minimize and levels of uncertainty (intervals), the optimal solution can be reached using an optimizer, as we propose in the following section.

## 2.2 Variability criteria to optimize

We denote by $x_{ij}$ the value of row $i$ and column $j$ in the resulting matrix $M$ and $\bar{x}$ the mean of all $x_{ij}$.

$$ \theta_{ij} = |x_{ij} - \bar{x}| \qquad (1) $$

We aim at minimizing 4 different objective criteria $f(\theta)$, each of them representing one approach to weight the variability among the different components of $x$. Then, given each specific objective criteria, we analyze the resulting optimal distributions of $\hat{x}_{ij}$ in $M$.

A first criterion is based on the $L_\infty$-norm, which illustrates an equity process: outliers have an important weight in this solution.

$$ f(\theta) = \max_{i,j}(\theta_{ij}) \qquad (2) $$

A second is similar but takes into account squared residuals $\theta_{ij}$. Outliers are even more considered.

$$ f(\theta) = \max_{i,j}(\theta_{ij}^2) \qquad (3) $$

The following case corresponds to the least absolute deviation case ($L_1$-norm) which depicts an efficiency purpose.

$$ f(\theta) = \sum_{i,j} \theta_{ij} \qquad (4) $$

This last criterion is the variance and corresponds to the $L_2$-norm (equality).

$$ f(\theta) = \sum_{i,j} \theta_{ij}^2 \qquad (5) $$

The drawback of the objective functions (2) and (4) is that many optimal solutions may exist, and in particular, optimal solutions that contain many components of $x$ set to 0 which does not provide interesting insights. To avoid that situation, we smoothed the objective functions by replacing $\theta$ with $\theta'$ defined as

$$ \theta'_{ij} = \theta_{ij} + 0.001 \theta_{ij}^2 \qquad (6) $$
Smoothing ensures unicity of the optimal solutions and, more importantly, increase the number of positive coefficients.

These functions are tested with different data relating to different assumptions and in two cases:

* No information is provided about the range of the value in each cell (no information known or expressed about uncertainty of $\hat{x}_{ij}$);

* In each cell, the value is bounded and the solution must take into account this constraint to be included in the interval (uncertainty modeling for $\hat{x}_{ij}$).

## 2.3 Mathematical model: optimization under constraints

Let $C_j$ be the given value for the sum of all elements in column $j \in \{1, \dots, m\}$ and $R_i$ be the sum of all elements in row $i \in \{1, \dots, n\}$. We denote by $x_{ij}$ the value of row $i$ and column $j$ in the resulting matrix $M$. The problem of finding the optimal values for $x$ can be stated as the following optimization problem

$$ (P) \equiv \left\{ \begin{array}{ll} \min & f(\theta) \\ \text{s.t.} & \displaystyle\sum_{i=1}^{n} x_i = C_j \quad j \in \{1, \dots, m\} \\ & \displaystyle\sum_{j=1}^{n} x_j = R_i \quad i \in \{1, \dots, n\} \\ & x \ge 0 \\ & \theta_{ij} \ge x_{ij} - \bar{x} \quad i \in \{1, \dots, n\}, j \in \{1, \dots, m\} \\ & \theta_{ij} \ge -x_{ij} + \bar{x} \quad i \in \{1, \dots, n\}, j \in \{1, \dots, m\} \end{array} \right. $$

where the function $f$ models the preference used in the construction of the matrix $M$ as explained in the previous section, and the last two groups of inequalities are a linearization for the definition of $\theta$ from (1). Notice that the objective functions involving $\max$ can be linearized by introducing the artificial optimization variable $\Theta$ linked to $\theta_{ij}$ through

$$ \Theta \ge \theta_{ij} \quad i \in \{1, \dots, n\}, j \in \{1, \dots, m\}, $$

and optimizing the objective function $f(\Theta) = \Theta$. Finally, the model can be completed by adding bounds on variables $x_{ij}$ as explained in the previous section.

Since the four cases of functions $f$ presented in the previous section are convex quadratic functions, $(P)$ turns to a convex quadratic and linearly constrained optimization problem, which can be solved very quickly with state-of-the-art optimization solvers (e.g. CPLEX, Gurobi).

## 2.4 Several cases and hypothesis to test

For all the studied cases, the process is similar: we set an hypothesis that will change the estimation of each value $x_{ij}$ in $M$. Then we recompute the sums of individuals in rows and columns and we adjust the constraints accordingly (e.g. bounds). Several types of solutions are searched, based on margin values and according to different hypothesis:

* Models without optimization:

- [Variable independence $H0$]. As explained in the introduction, this model is based on the $H0$ table of contingency; so $\hat{x}_{ij} = x_i \cdot x_j / x_{ii}$;

- [Global Illiteracy Rate $GIR$]. In the second model, we apply the global rate (%) of illiterate people observed in the whole USA on each type of population in every State;
- [Population Illiteracy Rate PIR]. In this model, we take into account the part of illiteracy in each type of population (%);

- [Effect of Old Confederacy EOF]. Here we add a historical hypothesis based on whether or not a given State belong (or was close) to the Confederation area during the civil war in USA (1861-1865); indeed, this could partly explain the capacity of a State to potentially integrate black or foreign population including language learning facilities. The probability of illiteracy is multiplied by 2 for the Confederate states and by 1.5 for a few bordering States which used to allow slavery and to belong to the Union however.

* Models with optimization but without any local information about value range (no information about uncertainty, in all of them, we only use recalculated margin values in rows and columns):

- [Raw Data RW]. We try to find out an optimal solution;

- [Global Illiteracy Rate GIR]. We look for an optimal solution using GIR;

* Models with local bounded values (the solver finds a solution given new local constraints of uncertainty in each matrix cell); we use recalculated margin values in rows and columns and we fix min and max bounds in each cell of M using two different formula 7 and 8; are concerned by this procedure:

- [Global Illiteracy Rate GIR] data (for the whole USA);

- [Population Illiteracy Rate PIR] data (by type of population);

- [Population Illiteracy Rate PIR] data (by type of population) updated by [Effect of Old Confederacy EOF] data.

For the models with local bounded values, we apply a function to find the minimum inf($\hat{x}_{ij}$) and maximum sup($\hat{x}_{ij}$) values of the interval centered on the estimated $\hat{x}_{ij}$:

$$ \sup(\hat{x}_{ij}) = \text{round}(\hat{x}_{ij} \cdot k^{1/\text{length}(\hat{x}_{ij})}) \quad (7) $$

and

$$ \inf(\hat{x}_{ij}) = \text{round}(2\hat{x}_{ij} - \sup(\hat{x}_{ij})) \quad (8) $$

These functions (equations 7 and 8) allow to have symmetric intervals, larger for low values of $\hat{x}_{ij}$ (indeed length counts the digits of $\hat{x}_{ij}$) depending on k.

For models dealing with uncertainty, we consider two cases:

* k = 2; for instance [109 - 185] and [6023 - 8835]

* k = 4; resp. [0 - 294] and [2363 - 12495]

### III RESULTS

We focus on the relation between illiteracy and black people in USA in 1930. To do so, we compute a LS regression model between observed data and estimations we obtained. We study the coefficient of determination r² and the slope of the regression. Closer to 1, better the quality solution, for both indicators.

In the Figure 2, we can notice that using optimization increases the estimation quality when expert assumption is weak (case with H0 and Global Illiteracy Rate GIR). For other hypothesis (Population Illiteracy Rate PIR and Effect of Old Confederacy EOF), optimization has no effect compared to estimations based on accurate hypothesis from experts.

Globally, Figures 3 and 4 show that:
* The results are better for bounded optimization methods;

* The difference between objective functions is not marked; indeed they all somehow compute a sort of variability minimization;

* Generally, solutions are better with $k = 2$ (more narrow intervals);

* Most of the time, functions including the *max* operator seem to get better estimates;

* The best optimization reaches a coefficient of determination $r^2$ of 0.86 and a slope of 0.76, both quite close to 1.

<figure><img src="image_129.png" /><figcaption>Figure 2: Comparison of different estimation methods according to the values of the coefficient of determination $r^2$ and of the slope of regression line (on left, 4 solutions without optimization; on right, 4 best solutions using optimization). Only the black population in USA in 1930 is considered.</figcaption></figure>

## IV CONCLUSION

In this paper, we explore how we can mix optimization methods and expert knowledge translated in statistical method to improve data estimation in contingency table downscaling. Although the results do not strongly demonstrate that the use of optimization is determinant on the result quality compared to other current methods from experts of the domain, it shows that, in certain conditions (weak assumptions or bounds knowledge from experts for instance), estimations can be improved, especially prediction quality (slope of the regression line between estimations and observed data).

Moreover, using optimization allows to find the most optimal solution among a large set of possible solutions, according to a given objective. These solutions can be then compared and participate in the data exploration.

Other experiments will be made, on the whole set of data whatever the type of population nativity, also on other types of data and with other compared methods (multiple regressions for example) to assess in which conditions, such a mixed approach may be useful for predicting data in contingency tables.
<figure><img src="image_130.png" /><figcaption>Figure 3: Comparison of different estimation methods according to the $r^2$ of the matrix *M* and to 4 different objective functions of $\theta$. Only the black population in USA in 1930 is considered.</figcaption></figure>

<figure><img src="image_131.png" /><figcaption>Figure 4: Comparison of different methods to estimate the slope of the regression line according to 4 different objective functions of $\theta$. Only the black population in USA in 1930 is considered.</figcaption></figure>

## References

* Bierkens M. F., A. F. P., de Willigen P. (2000). *Upscaling and Downscaling Methods for Environmental Research*. Springer, Series Developments in Plant and Soil Sciences.

* Gehlke C., Biehl H. (1934). Certain effects of grouping upon the size of the correlation coefficient in census tract material. *Journal of the American Statistical Association supplement* 29, 169–170.

* Josselin D., Mahfoud I., Fady B. (2008). Impact of a change of support on the assessment of biodiversity with
shannon entropy. In *Spatial Data Handling, SDH'2008", Montpellier, June, 23-25, pp. 109-131.

King G. (1997). *A solution to the ecological inference problem*. *Reconstructing individual behaviour from aggregate data*. Princeton University Press.

King G., Rosen O., Tanner A. M. (Eds.) (2004). *Ecological Inference*. New Methodological Strategies. Cambridge University Press.

Robinson W. (1950). *Ecological correlations and the behaviour of individuals*. *American Sociological Review* 15, 351-357.

Tranmer M., Steel D. (1998). Using census data to investigate the causes of the ecological fallacy. *Environment and Planning A* 30, 817-831.
# Downscaling of soil moisture with area-to-point geographically weighted regression kriging and uncertainty analysis

Yan Jin<fnref n="1" />,<fnref n="2" />, Yong Ge*<fnref n="1" />, Jianghao Wang<fnref n="1" />, Yuehong Chen<fnref n="1" />,<fnref n="2" />, Xudong Guan<fnref n="1" />,<fnref n="2" />

*Corresponding author: gey@lreis.ac.cn

## Abstract

Combine the visible/infrared remote sensing and passive microwave remote sensing has demonstrated its potential for capturing the high spatial resolution of the near-surface soil moisture. A geostatistical approach by integrating geographically weighted regression (GWR) and area-to-point kriging (ATPK) was proposed to downscale soil moisture. The proposed area-to-point geographically weighted regression kriging (GWATPRK) was verified on the AMSR-2 soil moisture product at 25 km resolution to obtain the downscaled soil moisture of 1 km spatial resolution. The results showed the GWATPRK can enhance the accuracy of the downscaled estimates. The best downscaled estimates owned the root mean square error value of 0.072m³/m³ and the mean absolute error value of 0.056m³/m³. The different auxiliary variables would affect the accuracy of downscaled estimates. Comparison of the downscaled estimates using different covariates with the *in situ* observations, the topography-corrected land surface temperature can improve the downscaled performance in the study area.

**Keywords:** downscaling, soil moisture, area-to-point kriging, geographically weighted regression, AMSR-2, MODIS

## I. INTRODUCTION

Soil moisture (SM) is one of fundamental land surface parameters for modelling eco-hydrological processes and understanding impacts and conservation of ecosystems (Srivastava et al., 2013). Although the ground measurement methods (Bogena et al., 2010) can obtain the SM information with depth monitoring and high accuracy, the expanded applications at large scales included regional or global scales are impractical. The remote sensing of surface SM is considered to fill this gap, including visible/infrared remote sensing and microwave remote sensing. Microwave remote sensing technique has been widely applied to monitor SM owing to high SM sensitivity and ignorance of atmospheric conditions. Most of these satellites can provide the SM products at a fine temporal resolution (1 to 3 days) and a coarse spatial resolution (about tens of kilometers). The coarse spatial resolution couldn't satisfy the application requirements (Ge et al., 2015a). One the other hand, it is difficult to match the large-scale SM to the fine-scale ground observations in the remote sensing validation. Therefore, it is necessary to improve the spatial resolution of remote sensing products of SM.

A number of methods have been developed to downscale the passive microwave remote sensing products (Kim and Hogue, 2012; Yang et al., 2016). Combine the visible/infrared remote sensing and passive microwave remote sensing that has become more general, because it benefits the fine resolution of the former and takes advantage of the latter

<footnote n="1">¹ State Key Laboratory of Resources and Environmental Information Systems, Institute of Geographic Sciences & Natural Resources Research, Chinese Academy of Sciences, Beijing 100101, China</footnote>

<footnote n="2">² University of Chinese Academy of Sciences, Beijing 100049, China</footnote>
which is high sensitive to the SM. Many SM downscaling algorithms focused on the relationship between the SM and the relative variables (Piles et al., 2011). Meanwhile, the geostatistical methods which consider the spatial correlation have received extensive concern in downscaling field (Atkinson et al., 2008) and growing applications in SM. The geostatistics perform well in spatial prediction and uncertainty analysis (Chiles and Delfiner, 1999). And the area-to-point regression kriging (ATPRK) (Kerry et al., 2012; Wang et al., 2015) is a flexible and well-performing downscaling method which incorporates regression kriging (Hengl, et al., 2007; Ge et al., 2015b) and ATPK (Kyriakidis, 2004). It can consider both the correlated variables and the change of support problem during the downscaling process. ATPRK has been shown the potential in downscaling the irregular geographical units (Zhang et al., 2014) but still less in remote sensing images and never be used in SM.

In this paper, a new geostatistical approach by integrating geographically weighted regression (GWR) (Fotheringham et al., 1997) and area-to-point kriging (ATPK) was proposed to downscale SM with high resolution visible/infrared variables (i.e., the Moderate Resolution Imaging Spectroradiometer (MODIS) products). The proposed area-to-point geographically weighted regression kriging (GWATPRK) as a downscaling strategy can improve SM estimates, because it retains the advantages of ATPRK and considers spatial nonstationarity relationship of variables. The approach was illustrated by an application for downscaling the SM remote sensing product of the Advanced Microwave Scanning Radiometer 2 (AMSR-2) (Imaoka et al., 2010) and a validation with in situ observations in the Heihe Water Allied Telemetry Experimental Research experiment (HiWATER) (Li et al., 2013).

## II. STUDY AREA AND DATA DESCRIPTION

### A. Study area

The study area is the upper reaches of the Heihe River Basin (HRB) in northwestern China, within latitudes 37.66 to 39.16 N, and longitudes 98.50 to 101.25 E. It includes the Babao River basin which is one of foci experimental areas of the comprehensive eco-hydrological program HiWATER. The district covers approximately 2452 km² with an elevation ranging from 2640 to 5000 m and is in a typical landscape of cold region with natural grassland as main vegetation. The topography of the study area with ground-based site locations and the AMSR-2 grid pixels (SM product of one day) was shown in Figure 1.

<figure><img src="image_132.png" /><figcaption>Figure 1: (a) The study area; (b) the ground-based site locations and AMSR-2 grid pixels.</figcaption></figure>

### B. In situ observations

In the upstream area of HRB, seven ordinary automatic meteorological stations (AMSs) and a superstation at A'rou have been installed to capture the atmospheric states. The Ecological and Hydrological Wireless Sensor Network (ENWSN) has been designed and deployed to monitor
the surface SM (Li et al., 2013). There are forty nodes of ENWSN to obtain the variation dynamically at depths of 4cm, 10cm and 20cm. In this study, one month (July 20-August 21, 2014) of in situ SM observations were collected. For validating the downscaled SM estimates from the AMSR-2 SM product, the average observations with 4 cm depth during satellite passes were considered, i.e., from 9:30 p.m. to 2:30 a.m. for the descending. Specifically, one period was adopted for validation because of the availability of the auxiliary MODIS images, i.e., 21 August, 2014 which was the day of Year DOY233.

### C. Remote sensing data

The Level 3 descending surface SM product of AMSR-2 (version 001) at a 25 km spatial resolution and the Version 5 MODIS night products at a 1 km spatial resolution were employed. The daily SM product data were resampled into 25km × 25km regular grids using the nearest neighbor resampling technique for completely overlapping the MODIS pixels. The daily LST was extracted from the MYD11A1 and MOD11A1 products. The NDVI was acquired from MYD13A2 product at a temporal resolution of 16 days. All of the MODIS data were projected and extracted consistently with large scale AMSR-2 data and the aggregated coarse MODIS images achieved full coverage. The 16-day NDVI product is cloud free. However, the single LST product owned a bad cloud free condition within the study area, it needed to combine the Terra and Aqua products. The data coverage fraction was 86% on DOY 233 and it was already high in available dates. Then the spline interpolation was used to interpolate the pixel values in uncovered region for obtaining the fully downscaled estimates.

## III. DOWNSCALING STRATEGY

Assume that vector $Z(S_i)$ represents coarse SM at n spatial observation pixels $S_i = (U_i, V_i)$ ($i=1,2,\dots,n$) and $Z(s_j)$ represents fine SM at $nF^2$ spatial observation pixels $s_j = (u_j, v_j)$ ($j=1,2,\dots,nF^2$). F is the ratio between the coarse resolution and the fine resolution. In which, $(U_i, V_i)$ and $(u_j, v_j)$ state the geographic positions (latitudes, longitudes) of the centers of the corresponding pixels, respectively. The proposed downscaling strategy contains two phases. It first establishes GWR model among the SM and the visible/infrared variables to predict the spatial trend of the SM at fine resolution; then applies ATPK to downscale the regression residuals to obtain the fine resolution predictions. The general form of GWATPRK estimate is:

$$Z(s_j) = m(s_j) + R(s_j) \quad (1)$$

where $m(s_j)$ is the deterministic part which can be estimated by GWR and $R(s_j)$ is the ATPK parts for regression residuals. The detailed descriptions were displayed in the following sections.

The GWR model is a local regression that generates local coefficients to each pixel and has been successfully used in SM (Yu et al., 2008). Many Researches have indicated the feasibility of applying the universal relationship at different scale resolutions (Zhang et al., 2014). Let $x_k(S_i)$ ($k=1,\dots,p$) be the value of the auxiliary variables at pixel $S_i$ and $\beta_l(U_i, V_i)$ ($l=0,1,\dots,p$) be regression coefficients with geographic position $(U_i, V_i)$. The corresponding parameters of the fine resolution are $x_k(s_j)$ at pixel $s_j$ and $\beta_l(u_j, v_j)$ with position $(u_j, v_j)$. The changing regression coefficients $\beta_l(U_i, V_i)$ at different position can be estimated by using the least square method in basic GWR, in which one of the key steps is to choose a distance-decay function for representing the strength of the connectivity between pixels. The Gaussian function which is one of the most common distance-decay functions was selected in the following experiments. Furthermore, the GWR with a locally compensated
ridge term (Lu et al., 2014) was employed to address local collinearity problems in basic GWR. The relationship among the SM and the auxiliary variables at coarse resolution was given:

$$Z(S_i) = \beta_0(U_i, V_i) + \sum_{k=1}^{p} \beta_k(U_i, V_i) \cdot x_k(S_i) \quad (2)$$

Meanwhile, the regression residual was computed using the following equation:

$$R(S_i) = Z(S_i) - [\beta_0(U_i, V_i) + \sum_{k=1}^{p} \beta_k(U_i, V_i) \cdot x_k(S_i)] \quad (3)$$

Based on the hypothesis of the scale invariant, the relationship in Eq. 2 can be applied at fine spatial resolution:

$$m(s_j) = \beta_0(u_j, v_j) + \sum_{k=1}^{p} \beta_k(u_j, v_j) \cdot x_k(s_j) \quad (4)$$

where $\beta_t(u_j, v_j)$ equals $\beta_t(U_i, V_i)$ when the pixel $s_j$ is covered by the large one $S_i$.

If the regression model was perfect, there would be no bias between the fitting result and the observations at coarse resolution, meaning that the regression part (Eq. 4) can attain the aim for downscaling. Because of the impracticality to search such ideal regression model, the regression residual (Eq. 3) should not be neglected. Consider that adding the residual into the regression prediction at fine resolution directly would ignore the change of supports, it was no downscaling step actually. The ATPK method was employed to downscale the regression residual to acquire the prediction at fine resolution by summing the weighted residuals at coarse resolutions. If a given pixel $s_j$ was estimated as a combination of the $R(S_h)$ ($h=1,2,\dots,m$) in $m$ large neighboring pixels, the ATPK model can be written as follows:

$$R(s_j) = \sum_{h=1}^{m} \lambda_h \cdot R(S_h) \quad (5)$$

where $\lambda_h$ are the weight coefficients for the prediction of the fine resolution and satisfy $\sum_{h=1}^{m} \lambda_h = 1$. The weight coefficients can be estimated by minimizing the prediction error variance, in which the main step for the implementation of ATPK is to obtain the point support semivariogram. A presented deconvolution procedure (Goovaerts, 2008) was utilized here.

Combine the regression predictions and downscaled fine residuals from the Eq. 4 and Eq. 5, the general form of GWATPRK estimate (Eq. 1) can be further written as:

$$Z(s_j) = \beta_0(u_j, v_j) + \sum_{k=1}^{p} \beta_k(u_j, v_j) \cdot x_k(s_j) + \sum_{h=1}^{m} \lambda_h \cdot \{Z(S_h) - [\beta_0(U_h, V_h) + \sum_{k=1}^{p} \beta_k(U_h, V_h) \cdot x_k(S_h)]\} \quad (6)$$

where the symbols of parameters stay the same as stated above. The coarse SM can be downscaled to the fine one through the Eq. 6.

# IV. RESULTS AND DISCUSSION

## A. Downscaled soil moisture maps

There were three groups of auxiliary variables as different independent variable(s) in downscaling. For comparison, a mentioned commonly used downscaling strategy known as ATPRK and a quadratic regression model (QRM) (Piles et al., 2011) which is one of the SM
downscaling models were adopted. The ATPRK was employed an ordinary linear regression.

<figure><img src="image_133.png" /><figcaption>Figure 2: (a) The vegetation index-surface temperature trapezoidal feature space and (b) the calculated VTCI at 1km resolution on DOY 233</figcaption></figure>

<figure><img src="image_134.png" /><figcaption>Figure 3: The images at 1km resolution on DOY 233. (a) The auxiliary variables; (b) the downscaled SM with LST; (c) the downscaled SM with corrected LST; (d) the downscaled SM with VTCI.</figcaption></figure>

Firstly, the LST and NDVI were employed as the auxiliary variables which have already been closely linked to SM (Chauhan et al., 2003). The correlation analysis figured out that the LST and NDVI have correlation coefficients of 0.34 and 0.41, respectively. The variance inflation factors of LST and NDVI were both reported less than 10 which revealed the reasonability for employing these two independent variables without considering the multicollinearity. Secondly, considering the topographic effects in hilly surrounding area, a topography-corrected land surface temperature (Merlin et al., 2013) was used to replace the LST in the first group. Thirdly, we chose the vegetation temperature condition index (VTCI) which has a positive correlation with SM. It has been demonstrated to increase downscaling performance with VTCI as the only input (Peng et al., 2016). The normalized LST and the
vegetation fraction cover were employed to describe the feature space. Figure 2a and 2b showed the physically meaningful trapezoidal feature space and the calculated VTCI at 1km resolution, respectively.

The three downscaling approaches were realized on the three variables group. Figure 3a showed images corresponding to the auxiliary data at 1km resolution, meanwhile Figure 3b, 3c and 3d showed downscaled images for different approaches with varying variables. The downscaled results have a consistent tendency. The white parts in the images owned the highest SM values which were located in the low elevation district of the upstream area. All top right corners displayed the lowest SM values in accordance with the conditions of the highest temperature and less vegetation coverage. The GWATPRK and ATPRK could capture more similar spatial information of SM with the coarse image. And the downscaled images of GWATPRK represented stronger variability. The GWATPRK could obtain less minimum values of downscaled estimates which were closer to the minimum *in situ* observation. Furthermore, the values greater than around 0.7 m³/m³ of downscaled SM estimates may not accord with the natural cognition of the cold region. It caused by the large original SM values of the AMSR-2 data. And to focus on the downscaling methods, we skipped correcting for the AMSR-2 SM product here. In the following comparison, the GWATPRK with changed regression coefficients at each coarse pixel has obtained the good quality of which downscaled values were closer to the ground observations expectably.

## B. Validation and comparison

To assess the accuracy of the downscaling strategy, the *in situ* observations were compared with the downscaled SM, including the mean absolute error (MAE) (m³/m³), the standard deviation of the error (MSD) (m³/m³), the root mean square error (RMSE) (m³/m³) and the correlation coefficient (R). It desires the downscaled estimates to be more approximate to the *in situ* measurements. The less MAE, less RMSE, and higher R would be expected. Considering one day of most growing season days was chosen and the corresponding AMSR-2 SM values were all greater than the *in situ* measurements. The mean of these differences between them would be removed from the downscaled estimates for a more suitable comparison with ground observations. The downscaled SM estimates were plotted against the *in situ* measurements and the accuracy of SM results of different downscaling approaches described in Figure 4.

<figure><img src="image_135.png" /><figcaption>Figure 4: Downscaled soil moisture estimates versus *in situ* measurements and the summary of the comparison results. (a) The downscaled soil moisture with uncorrected LST; (b) the downscaled soil moisture with corrected LST; (c) the downscaled soil moisture with VTCI.</figcaption></figure>

Although the downscaled estimates of the QRM had the higher *R*, the *MAE* and the *RMSE* were also higher. And its estimates were almost all over or below the 1:1 line. The results demonstrated that the GWATPRK downscaling strategy achieved better downscaled estimates
than other two methods in every variables group. The corrected LST can improve the estimates of SM in each method. The VTCI just reduced the *MAE* of QRM and increased the *R* of ATPRK, albeit the downscaled result of GWATPRK was better than other two. It was also clear that with the least *MAE* value of 0.056 m³/m³, a less *MSD* value of 0.047 m³/m³, the least *RMSE* value of 0.072 m³/m³ and a higher *R* value of 0.896, the GWATPRK downscaling method with corrected LST resulted in the best performance.

The downscaled error mainly comes from the input images, the models and the scale effect during the downscaling process (Zhao et al., 2015). The latter two are always complex and complicated, which are still affected by all the participating variables. Because that the study area is a hilly surrounding area, the topography-corrected LST is actually a less error variable for impacting SM than the original LST by eliminating the terrain influence. It worked well. In addition, the VTCI as a synthetic index didn't improve the accuracy expectably, in which the LST has not been corrected. The VTCI of ignoring the topography may not adequately capture the interactions with SM in this mountainous area. Meanwhile, it is best to apply the VTCI with wide range required vegetation coverage. It also reminded us again to use the temperature-vegetation indexes with more analyses of the area conditions. To select the appropriate variable will conduce to enhance the downscaled accuracy.

## V. CONCLUSION

Downscaling passive microwave remote sensing SM product is an essential problem in SM monitor. To improve SM spatial resolution, GWATPRK was proposed to downscale SM by integrating GWR and ATPK. We evaluated different downscaling methods on the AMSR-2 SM product in the upper reaches of the HRB. The results showed the GWATPRK can produce downscaled images at a fine spatial resolution with greater qualities for heterogeneous region. And it would improve the accuracy of the downscaled estimates by using the topography-corrected LST in this study area. The input images as a significant part need considerations in downscaling process which would affect the accuracy of the downscaled results. Moreover, for an operational downscaled SM approach, further validation work will be carried out in the future studies.

## ACKNOWLEDGMENTS

This work was supported by the Key Program of the National Natural Science Foundation of China (No. 41531174). The authors would like to thank all the data producers, and the anonymous reviews for their comments and suggestions to improve the quality of the paper.

## References

* Atkinson P M, Pardo-Iguzquiza E, Chica-Olmo M. (2008). Downscaling cokriging for super-resolution mapping of continua in remotely sensed images. *Geoscience and Remote Sensing, IEEE Transactions on* 46(2): 573–580.

* Bogena H R, Herbst M, Huisman J A, et al. (2010). Potential of wireless sensor networks for measuring soil water content variability. *Vadose Zone Journal*, 9(4): 1002–1013.

* Chauhan N S, Miller S, Ardanuy P. (2003). Spaceborne soil moisture estimation at high resolution: a microwave-optical/IR synergistic approach. *International Journal of Remote Sensing*, 24(22): 4599–4622.

* Chiles, J. P., and Delfiner, A. (1999). Geostatistics: Modelling spatial uncertainty: Wiley Interscience. *New York*.

* Fotheringham, A. S., Brunsdon, C., & Charlton, M. (2003). Geographically weighted regression: the analysis of spatially varying relationships. *John Wiley & Sons*.

* Ge, Y., Wang, J. H., Heuvelink, G. B. M., Jin, R., Li, X., and Wang, J. F. (2015a). Sampling design optimization of a wireless sensor network for monitoring ecohydrological processes in the Babao River basin, China.
*International Journal of Geographical Information Science*, 29(1), 92–110.

Ge, Y., Liang, Y., Wang, J., Zhao, Q., and Liu, S. (2015b). Upscaling sensible heat fluxes with area-to-area regression kriging. *Geoscience and Remote Sensing Letters*, IEEE, 12(3), 656–660.

Goovaerts P. (2008). Kriging and semivariogram deconvolution in the presence of irregular geographical units. *Mathematical Geosciences*, 40(1): 101–128.

Hengl, T., Heuvelinkb, G. B. M., & Rossiter, D. G. (2007). About regression-kriging: From equations to case studies. *Computers & Geosciences*, 33,1301–1315.

Imaoka K, Kachi M, Kasahara M, et al. (2010). Instrument performance and calibration of AMSR-E and AMSR2. *International Archives of the Photogrammetry, Remote Sensing and Spatial Information Science*, 38(8): 13–18.

Kerry R, Goovaerts P, Rawlins B G, et al. (2012). Disaggregation of legacy soil data using area to point kriging for mapping soil organic carbon at the regional scale. *Geoderma*, 170: 347–358.

Kim, J., Hogue, T. S. (2012). Improving spatial soil moisture representation through integration of AMSR-E and MODIS products. *Geoscience and Remote Sensing, IEEE Transactions on*, 50(2), 446–460.

Kyriakidis P C. (2004). A Geostatistical Framework for Area-to-Point Spatial Interpolation. *Geographical Analysis*, 36(3):259–289.

Li X, et al. (2013). Heihe watershed allied telemetry experimental research (HiWATER): Scientific objectives and experimental design. *Bulletin of the American Meteorological Society*, 94(8): 1145–1160.

Lu, B., Harris, P., Charlton, M., & Brunsdon, C. (2014). The GWmodel R package: further topics for exploring spatial heterogeneity using geographically weighted models. *Geo-spatial Information Science*, 17(2), 85-101.

Merlin O, Escorihuela M J, Mayoral M A, et al. (2013). Self-calibrated evaporation-based disaggregation of SMOS soil moisture: An evaluation study at 3km and 100m resolution in Catalunya, Spain. *Remote sensing of environment*, 130: 25–38.

Peng J, Loew A, Zhang S, et al. (2016). Spatial downscaling of satellite soil moisture data using a vegetation temperature condition index. *Geoscience and Remote Sensing, IEEE Transactions on* 54(1): 558–566.

Piles M, Camps A, Vall-Llossera M, et al. (2011). Downscaling SMOS-derived soil moisture using MODIS visible/infrared data. *Geoscience and Remote Sensing, IEEE Transactions on* 49(9): 3156–3166.

Srivastava P K, Han D, Ramirez M R, et al. (2013). Machine learning techniques for downscaling SMOS satellite soil moisture using MODIS land surface temperature for hydrological application. *Water resources management*, 27(8), 3127-3144.

Wang Q, Shi W, Atkinson P M, et al. (2015). Downscaling MODIS images with area-to-point regression kriging. *Remote Sensing of Environment*, 166: 191–204.

Yang K, Zhu L, Chen Y, et al. (2016). Land surface model calibration through microwave data assimilation for improving soil moisture simulations. *Journal of Hydrology*, 533: 266-276.

Yu G, Di L, Yang W. (2008). Downscaling of global soil moisture using auxiliary data//Geoscience and Remote Sensing Symposium. IGARSS 2008. IEEE International. IEEE, 2008, 3: III-230-III-233.

Zhang J, Atkinson P, Goodchild M F. (2014). Scale in spatial information and analysis. CRC Press.

Zhao W, Li A. (2015). A comparison study on empirical microwave soil moisture downscaling methods based on the integration of microwave-optical/IR data on the Tibetan Plateau. *International Journal of Remote Sensing*, 36(19-20): 4986–5002.
<figure><img src="image_136.png" /></figure>

SPATIAL
ACCURACY

2016

5-8 JULY 2016

MONTPELLIER, FRANCE

# Uncertainty in 3D geodata
<figure><img src="image_137.png" /></figure>
# A proposal of a statistical test to control positional accuracy by means of 2 tolerances simultaneously

José Rodríguez-Avi*¹, Francisco Javier Ariza-López¹

¹University of Jaén, Jaén, Spain

*Corresponding author: jravi@ujaen.es

## Abstract

We propose a new method for positional accuracy quality control of spatial data. This method is valid for 1D, 2D, 3D and nD dimensional data. Normality of errors is not required. The method is an exact statistical hypothesis testing based on multinomial distribution. The proportions of the multinomial distribution are defined by means of several metric tolerances. The proposed statistical test is exact, so the p-value can be derived by exploring a space of solutions and summing up the probabilities of each isolated case of this space. The performance of the test has been analyzed by means of a simulation procedure. The validity and the power of the contrast seem to be good enough. Some application examples are presented for the case of working with two tolerances.

## Keywords

Multinomial distribution, Spatial Data, positional accuracy quality control

# I INTRODUCTION

Spatial data are data referring to features that have a position in space, and for this reason positional quality is one of the most desirable characteristics of spatial data sets (SDS). Positional accuracy of SDS has traditionally been evaluated using control points. Following this idea there are a great many statistical Positional Accuracy Assessment Methodologies (PAAM) like: National Map Accuracy Standard (NMAS) (USBB, 1947), Engineering Map Accuracy Standard (EMAS) (ASCE, 1983), National Standard for Spatial Data Accuracy (NSSDA) (FGDC, 1998), STANAG 2215 (STANAG, 2002), ASPRS Accuracy Standards for Large-Scale Maps (ASLSM) (ASPRS, 1990), the new ASPRS Positional Accuracy Standards for Digital Geospatial Data (ASPRS, 2015), and many others.

The majority of PAAMs take as an underlying hypothesis the Gaussian distribution (Normal Distribution) of positional errors, and here are some examples taken from the existing PAAM. The EMAS shows a control (acceptance/rejection) perspective and proposes a set of statistical hypothesis control test based on the normality of data, the first one for testing the presence of bias (a t-Student test) and the second one for comparing variational behaviour of data (a Chi squared test). In another way, the NSSDA shows an estimation perspective, determining an accuracy value scaled to a confidence value (95%) and leaving to the user the final decision about the appropriateness of the data. In the NSSDA the planimetric case is computed using a Chi2, which can only be used under the assumption that both X and Y error distributions follow a Normal Distribution where homogeneity of variances is needed, and the altimetric case is computed using the value 1.96 as a scaling factor, which is only valid under the assumption that Z error distribution follows a Normal Distribution Function. Another PAAM that requires the strict normality of errors is the procedure stated in ISO 3951-1 (ISO 2013b), being applied in some European countries. The ISO 3951-1 method is a sequence of tests allowing the acceptance/rejection of lots that can be applied to spatial data based on their positional behaviour.
It is also important to notice here that prior to any statistical analysis of the data assuming normality, the outliers must be eliminated. The elimination of outliers in parametric distributions is a usual issue, and for the case of the normal distribution there are many methods, such as those presented in ISO 16269-4 (ISO 2010) (for more detail on outliers elimination methods see Barnett and Lewis 1994).

But several studies indicate that this hypothesis (normality of errors) is not true. There are some proposals, for instance, in the case of Global Navigation Satellite Systems error distribution, the Rayleigh distribution (Logsdon, 1995); and for the case of geocoding errors, a log-normal distribution (Cayo and Talbot, 2003). For the case of vertical errors in digital elevation models there are many references (e.g. Bonin and Rousseaux 2005, Oksanen and Sarjakoski 2006) indicating that error distribution is not Normal, and for this reason it is proposed to perform and express the results of quality control checks by means of percentiles (Maune 2007) of the observed distribution. The latter was introduced by the ASPRS into the Guidelines for reporting vertical accuracy for Lidar data (ASPRS, 2004) and the most recent standard from the ASPRS (ASPRS, 2015) follows this method by informing separately the vegetated vertical accuracy (supposed to be non-Gaussian distributed) and the non-vegetated vertical accuracy (supposed to be normally distributed). This situation is somewhat confusing for non-expert users.

The Normal Distribution Function and other distribution functions mentioned above are parametric distribution functions, which means that they can be determined by a set of parameters, for instance the mean (μ) and the variance (σ²) in the case of the Normal Distribution Function. Parametric models summarize the behaviour of a population by means of a mathematical formula and parameters which control its shape (e.g. for the Normal Distribution Function the mean gives the position and the deviation the more or less flattened shape of the curve). Parametric distribution functions are very convenient when dealing with scarce data and computing capabilities are limited, as in the past. But now we live surrounded by very large data sets, and by very large computing capabilities, so we can use the observed distribution functions (Free-Distribution Functions or Parameters-Free-Distribution Functions) without major problems.

The PAAMs based on Free-Distribution Functions are scarce in positional accuracy assessment. For instance, the use of percentiles is proposed for dealing with height data of DEM when capture (e.g. by Lidar) is performed on a vegetal cover. However, the use of percentiles is merely descriptive, and no method for quality control is proposed. Recently, Ariza-López & Rodríguez-Avi (2014a) proposed a method that can be applied to data belonging to free-distribution functions. A step further would require a statistical method allowing control in the same way (mean value, variation, control of outliers) as performed when data come from a Normal Distribution Function.

Our aim is to propose a general positional accuracy control method for dealing with data following any kind of distribution function. It is a method based on the observed distribution function of the data and, in this way, can be applied to 1D, 2D or 3D data without the limitations of homogeneity of variances stated in traditional methods (e.g. EMAS, NSSDA, STANAG 2215, etc.). To achieve this aim we propose a method based on proportions of a multinomial distribution function in order to establish a strict control over data coming from any distribution function. The control is multiple and can test jointly proportions corresponding to tolerances related with, for instance, median values, extreme values (e.g. 95%), or the amount of outliers existing in the data set. The proposed control is based on an exact test, in the same way as the Fisher's exact test (Fisher 1922, Freeman and Halton 1951, Müller 2001).

## II PROPOSED METHOD

In order to analyse a pass/fail model applied to positional accuracy control, an approximation is given in terms of the Binomial distribution (Ariza-López & Rodríguez-Avi, 2014). For each sampling point taken in a k-dimensional space, the value:
$$E_i = \sqrt{\sum_{i=1}^{k} (x_i - x_i^T)^2}$$

is calculated (where $(x_1^T, ..., x_k^T)$ is the “exact value).

Given a tolerance $T$, we count the number of sampling points where $E_i$ is greater than $T$. This can be seen as a realization of a Binomial distribution with parameters $N$ and $p$. In this case we consider a hypothesis test where the null hypothesis is: the proportion of defectives points is less or equal than $p_0$, against $H_1$: the proportion of defectives is greater than $p_0$. In consequence, a right hand unilateral test is proposed and $H_0$ is rejected if $p = P[X > x|X \sim B(N, p_0)] < \alpha$ where $x$ is the number of defectives (i.e. number of points whose error is greater than the previously specified tolerance) found in the sample.

This approximation is a pass/fail procedure and implies that we are not able to distinguish the “degree of bad quality” of a defective point (in position). In several circumstances is interesting to propose an ordination in the degree of fail. We propose a gradation in the degree of positional defectiveness, splitting the interval, the excess of error greater than the tolerance, into two parts. In this way we consider two metric tolerances defining three intervals for the positional error values measured in the sample. Each interval defines a quality level, and we can determine the minimum percentage of “correct” points in the first level (best level), and the maximum percentage of points in each of the other two levels (worse levels). For instance, and in order to control outlier points, we can stablish a percentage of at least a $p_1$% of points with an error less than a tolerance $T_1$ and at most a $p_2$% of points with an error between tolerances $T_1$ and $T_2$ ($T_2 > T_1$) and, in consequence, at most a percentage of $p_3$% = 100 - $p_1$% - $p_2$% of points with error greater than tolerance $T_3$. So, we classify the sample into three (instead of two, as before) categories, and consequently, the multinomial approximation instead of the binomial approximation is adequate.

Let us recall that the multinomial distribution is a multivariate extension of a binomial distribution when we classify the result of an experiment into more than two categories. So if we realize $N$ independent experiments where we classify the results for exactly one of 3 categories, with probabilities $\pi_1, \pi_2, \pi_3$, and $\pi_1 + \pi_2 + \pi_3 = 1$, then the mass probability function of a such multinomial $M(N, \pi_1, \pi_2, \pi_3)$ is given by:

$$P[(X_1 = N_1, X_2 = N_2, X_3 = N_3) = \binom{N!}{N_1! N_2! N_3!} \pi_1^{N_1} \pi_2^{N_2} \pi_3^{N_3}$$

Where $N_i$ is the number of points that belongs to the category $i$, which has a probability $\pi_i$.

We use the multinomial distribution to propose an exact test to decide whether a specified hypothesis is false or not. Firstly we decide two tolerance levels, $T_1$ and $T_2$, $T_2 > T_1$ to decide if a point is adequate ($E_i < T_1$), slightly inadequate ($T_1 < E_i < T_2$) or roughly inadequate ($E_i > T_2$).

The null hypothesis implies that the proportion of points in each category is known. In consequence:

* $H_0$: The sampling statistics has a multinomial distribution with parameters $N, \pi_1^0, \pi_2^0, \pi_3^0$ where $\pi_i^0$ is the hypothesized probability for category $i, i = 1, 2, 3$

We propose an unilateral test (we will only rejected the null hypothesis if the true error distribution is worse, and this situation takes place when the proportion of elements with tolerance less than $T_1$ is less than $P_1$, or when the other two proportions account for more than $P_2$ or $P_3$, because this implies a worsening in tails.
A sampling of size N is dropped, and we count the number of point that belongs in each category. The sampling statistics is S = (N₁, N₂, N₃), so that N₁ + N₂ + N₃ = N, where Nᵢ is the number of points that, in a sample of size N belongs to category i, i = 1, 2, 3

Under the null hypothesis we can calculate the exact probability of obtaining a such point S and the probability of every point worse than S. As a result, to calculate the p-value we sum the probability, under the null hypothesis, of elements (M₁, M₂, M₃) that verify:

* M₁ < N₁

* M₁ = N₁ y M₂ < N₂

And we reject the null hypothesis if the p-value obtained is lesser than α.

To prove the validity of our proposal we show two complementary simulations:

* The first one is in order to prove the validity of the contrast with respect to the Type I error.

* The second one is in order to prove the power of the contrast, that it is to say, the probability of rejecting the null hypothesis when this hypothesis is false.

<table><caption>Table 1 - Validity of the exact contrast with respect to the Type I error</caption><thead><tr><th rowspan="2">Case<br>π<sub>0</sub> = (0.8; 0.15; 0.05)</th><th rowspan="2">Critical<br>value</th><th colspan="7">Sampling size (n)</th></tr><tr><th>20</th><th>50</th><th>75</th><th>100</th><th>150</th><th>200</th><th>300</th></tr></thead><tbody><tr><td></td><td>0.5</td><td>0.5190</td><td>0.5205</td><td>0.5195</td><td>0.4980</td><td>0.4950</td><td>0.4885</td><td>0.4725</td></tr><tr><td></td><td>0.3</td><td>0.3345</td><td>0.3015</td><td>0.3090</td><td>0.2880</td><td>0.2875</td><td>0.2905</td><td>0.2815</td></tr><tr><td></td><td>0.2</td><td>0.2205</td><td>0.2000</td><td>0.2030</td><td>0.2025</td><td>0.1925</td><td>0.1960</td><td>0.1960</td></tr><tr><td></td><td>0.1</td><td>0.1095</td><td>0.1020</td><td>0.0970</td><td>0.1050</td><td>0.0975</td><td>0.0890</td><td>0.1010</td></tr><tr><td></td><td>0.05</td><td>0.0655</td><td>0.0565</td><td>0.0555</td><td>0.0470</td><td>0.0520</td><td>0.0450</td><td>0.0505</td></tr><tr><td></td><td>0.01</td><td>0.0115</td><td>0.0130</td><td>0.0120</td><td>0.0095</td><td>0.0085</td><td>0.0090</td><td>0.0095</td></tr></tbody></table>

Table 1 shows the proportion of simulated points obtained when the null hypothesis is, in fact, the true population distribution (in this case, π₀ = 0.80, 0.15, 0.05) and Table 2 shows the power of the test when the true population distribution is worse than the null hypothesis, and in consequence, the test would be rejected.

<table><caption>Table 2. - Proportion of times where the null hypothesis is rejected</caption><thead><tr><th rowspan="3">Case<br>π<sub>0</sub> = (0.7; 0.25; 0.05)</th><th rowspan="3">Critical<br>value</th><th colspan="7">Sampling size (n)</th></tr><tr><th colspan="7">True population probability: π<sup>1</sup> = (0.69; 0.01; 0.30)</th></tr><tr><th>20</th><th>50</th><th>75</th><th>100</th><th>150</th><th>200</th><th>300</th></tr></thead><tbody><tr><td></td><td>0.5</td><td>0.6235</td><td>0.6235</td><td>0.6930</td><td>0.6255</td><td>0.6510</td><td>0.6380</td><td>0.6705</td></tr><tr><td></td><td>0.3</td><td>0.4240</td><td>0.3865</td><td>0.3925</td><td>0.4455</td><td>0.4295</td><td>0.4660</td><td>0.4865</td></tr><tr><td></td><td>0.2</td><td>0.2560</td><td>0.2810</td><td>0.2900</td><td>0.2870</td><td>0.3045</td><td>0.3435</td><td>0.3340</td></tr><tr><td></td><td>0.1</td><td>0.1395</td><td>0.1890</td><td>0.1450</td><td>0.1755</td><td>0.1885</td><td>0.1955</td><td>0.2050</td></tr><tr><td></td><td>0.05</td><td>0.1395</td><td>0.1160</td><td>0.0935</td><td>0.0820</td><td>0.1040</td><td>0.0960</td><td>0.1180</td></tr><tr><td></td><td>0.01</td><td>0.0230</td><td>0.0180</td><td>0.0205</td><td>0.0220</td><td>0.0260</td><td>0.0285</td><td>0.0295</td></tr></tbody></table>

Recall that in this last case, a binomial approach will not be rejected the null hypothesis because the probability of adequate point under it is similar to the true population distribution.

### III APPLICATION EXAMPLES

We illustrate this procedure working with two different examples. In the first example $H_0$ is true, and we are going to see what happens if a sample confirms this hypothesis or not and in the second example the true model is worse, in the sense of generating a higher number of errors than the desired model.
We work in a three-dimensional space, and we suppose that a spatial data product establishes that errors in X, Y and Z are distributed according to three Normal and independent distributions with μ=0 m and σ=1m. This is the null hypothesis for both cases. Under this hypothesis, the quadratic error in each element (e.g. point, line or whatever kind) is:

$$QE_i = x_i^2 + y_i^2 + z_i^2$$

which is distributed according to a chi-square distribution with 3 degrees of freedom. In this case, the probability that an element has a $E \le 2$ m ($QE_i < 4 m^2$) is approximately 0.74, and the probability that an element has a $E \le 2.5$ m ($QE_i < 6.25 m^2$) is 0.90. In consequence, if we take a sample of control elements of size N, calculate the errors and count the number of elements whose $QE \le 4 m^2$; the number of elements with $4 m^2 \le QE \le 6.25 m^2$ and the number of elements whose $QE > 6.25 m^2$, these quantities will follow a multinomial distribution with parameters (N, 0.74, 0.16, 0.10).

Let suppose a sampling of size 30 obtained for such population.

<table><caption>Table 3 Example with a hypothesized case (Null population = True population)</caption><thead><tr><th>Element</th><th>e<sub>X</sub> [m]</th><th>e<sub>Y</sub> [m]</th><th>e<sub>Z</sub> [m]</th><th>QE [m<sup>2</sup>]</th><th>Category</th></tr></thead><tbody><tr><td>1</td><td>0.378</td><td>-1.044</td><td>-0.104</td><td>1.244</td><td>I</td></tr><tr><td>2</td><td>-0.607</td><td>-0.528</td><td>1.127</td><td>1.917</td><td>I</td></tr><tr><td>3</td><td>0.069</td><td>0.003</td><td>-0.711</td><td>0.511</td><td>I</td></tr><tr><td>4</td><td>-1.211</td><td>1.644</td><td>-1.073</td><td>5.322</td><td>II</td></tr><tr><td>5</td><td>-0.600</td><td>-0.744</td><td>-0.887</td><td>1.701</td><td>I</td></tr><tr><td>6</td><td>0.238</td><td>0.806</td><td>-0.278</td><td>0.784</td><td>I</td></tr><tr><td>7</td><td>-0.431</td><td>-0.187</td><td>-0.397</td><td>0.378</td><td>I</td></tr><tr><td>8</td><td>0.798</td><td>-0.059</td><td>2.211</td><td>5.529</td><td>II</td></tr><tr><td>9</td><td>1.518</td><td>0.004</td><td>0.266</td><td>2.376</td><td>I</td></tr><tr><td>10</td><td>1.302</td><td>1.442</td><td>0.217</td><td>3.822</td><td>I</td></tr><tr><td>11</td><td>-0.880</td><td>0.662</td><td>-0.063</td><td>1.218</td><td>I</td></tr><tr><td>12</td><td>-1.594</td><td>-0.549</td><td>-1.194</td><td>4.269</td><td>II</td></tr><tr><td>13</td><td>-1.274</td><td>-0.699</td><td>1.272</td><td>3.732</td><td>I</td></tr><tr><td>14</td><td>-1.633</td><td>0.855</td><td>-0.914</td><td>4.234</td><td>II</td></tr><tr><td>15</td><td>-0.884</td><td>-0.592</td><td>1.395</td><td>3.077</td><td>I</td></tr><tr><td>16</td><td>2.109</td><td>0.210</td><td>-0.667</td><td>4.938</td><td>II</td></tr><tr><td>17</td><td>0.497</td><td>-0.203</td><td>-0.498</td><td>0.536</td><td>I</td></tr><tr><td>18</td><td>0.296</td><td>-0.719</td><td>-0.441</td><td>0.800</td><td>I</td></tr><tr><td>19</td><td>1.429</td><td>1.264</td><td>0.098</td><td>3.649</td><td>I</td></tr><tr><td>20</td><td>-0.594</td><td>0.704</td><td>1.139</td><td>2.145</td><td>I</td></tr><tr><td>21</td><td>1.286</td><td>-1.149</td><td>-0.082</td><td>2.982</td><td>I</td></tr><tr><td>22</td><td>-0.717</td><td>-1.082</td><td>0.630</td><td>2.083</td><td>I</td></tr><tr><td>23</td><td>1.088</td><td>-0.793</td><td>-1.565</td><td>4.261</td><td>II</td></tr><tr><td>24</td><td>-0.546</td><td>-1.149</td><td>-0.561</td><td>1.935</td><td>I</td></tr><tr><td>25</td><td>-0.543</td><td>-1.332</td><td>0.220</td><td>2.117</td><td>I</td></tr><tr><td>26</td><td>0.533</td><td>-1.214</td><td>1.457</td><td>3.881</td><td>I</td></tr><tr><td>27</td><td>2.795</td><td>0.639</td><td>0.653</td><td>8.648</td><td>III</td></tr><tr><td>28</td><td>-0.196</td><td>1.201</td><td>1.514</td><td>3.775</td><td>I</td></tr><tr><td>29</td><td>0.970</td><td>-0.605</td><td>-0.190</td><td>1.343</td><td>I</td></tr><tr><td>30</td><td>-1.016</td><td>0.220</td><td>1.261</td><td>2.671</td><td>I</td></tr></tbody></table>

In this case, the sampling statistics is $S = (23, 6, 1)$, To obtain the p-value, we calculate the probability in a multinomial (30; 0.74; 0.16;0.10) of the point S and all points where the first value is less than 23 or if the first point is 23 the second point is lesser than 6. Finally, all these probabilities are summed and the p-value is 0.5362. In consequence, we do not reject the null hypothesis.
For the second example, we want to test if the same null hypothesis as before is true, but in a second population. In this population, errors are truly distributed according to three normal standard distributions, as before, but now, errors in X and Y are highly correlated (correlation coefficient equals to 0.6). We obtain a sample of size 30 for this new population, and we obtain the sampling that appears in Table 4.

For this sample we classify each point into its corresponding class according to its QE value, and the test estimator value is S = (18, 9, 3). As before, to obtain the p-value, we calculate the probability in a multinomial (30; 0.74; 0.16;0.10) of the point S and all points where the first value is less than 18 or if the first point is 18 the second point is lesser than 9. Finally, all these probabilities are summed and the p-value is 0.02950. In consequence, we can reject the null hypothesis, to an alpha level of 5%.

<table><caption>Table 4 Example with a hypothesized case (Null population worse than True population)</caption><thead><tr><td>Element</td><td>e<sub>X</sub> [m]</td><td>e<sub>Y</sub> [m]</td><td>e<sub>Z</sub> [m]</td><td>QE [m<sup>2</sup>]</td><td>Category</td></tr></thead><tbody><tr><td>1</td><td>-0.3869</td><td>0.6570</td><td>1.7462</td><td>3.6306</td><td>I</td></tr><tr><td>2</td><td>-0.7862</td><td>-0.5534</td><td>0.2419</td><td>0.9828</td><td>I</td></tr><tr><td>3</td><td>1.4546</td><td>-0.3234</td><td>0.8716</td><td>2.9802</td><td>I</td></tr><tr><td>4</td><td>-0.0815</td><td>0.7061</td><td>-0.2738</td><td>0.5803</td><td>I</td></tr><tr><td>5</td><td>0.6979</td><td>0.6140</td><td>0.7460</td><td>1.4207</td><td>I</td></tr><tr><td>6</td><td>-0.0322</td><td>0.0424</td><td>0.2619</td><td>0.0714</td><td>I</td></tr><tr><td>7</td><td>1.0974</td><td>1.7072</td><td>1.0563</td><td>5.2346</td><td>II</td></tr><tr><td>8</td><td>-1.6764</td><td>0.3839</td><td>-1.2768</td><td>4.5880</td><td>II</td></tr><tr><td>9</td><td>-1.3978</td><td>-1.6323</td><td>-1.0792</td><td>5.7830</td><td>II</td></tr><tr><td>10</td><td>-1.0880</td><td>-0.6952</td><td>0.2615</td><td>1.7353</td><td>I</td></tr><tr><td>11</td><td>0.8909</td><td>-0.9406</td><td>1.4338</td><td>3.7342</td><td>I</td></tr><tr><td>12</td><td>0.6538</td><td>1.3700</td><td>-1.1415</td><td>3.6072</td><td>I</td></tr><tr><td>13</td><td>1.3760</td><td>1.1411</td><td>-1.2979</td><td>4.8800</td><td>II</td></tr><tr><td>14</td><td>-1.2167</td><td>-1.5128</td><td>0.0603</td><td>3.7727</td><td>I</td></tr><tr><td>15</td><td>1.5711</td><td>1.7410</td><td>-1.4285</td><td>7.5398</td><td>III</td></tr><tr><td>16</td><td>1.3679</td><td>1.5336</td><td>0.0297</td><td>4.2240</td><td>II</td></tr><tr><td>17</td><td>1.0336</td><td>-0.1400</td><td>0.5406</td><td>1.3802</td><td>I</td></tr><tr><td>18</td><td>1.5528</td><td>0.6278</td><td>-1.1301</td><td>4.0823</td><td>II</td></tr><tr><td>19</td><td>-1.9135</td><td>-2.0818</td><td>-0.4297</td><td>8.1803</td><td>III</td></tr><tr><td>20</td><td>1.4945</td><td>1.2784</td><td>-1.0525</td><td>4.9755</td><td>II</td></tr><tr><td>21</td><td>-1.6507</td><td>-1.2067</td><td>-0.2056</td><td>4.2232</td><td>II</td></tr><tr><td>22</td><td>-0.2347</td><td>-1.6093</td><td>1.2372</td><td>4.1755</td><td>II</td></tr><tr><td>23</td><td>1.1281</td><td>0.9067</td><td>1.0737</td><td>3.2474</td><td>I</td></tr><tr><td>24</td><td>1.5816</td><td>0.8454</td><td>-0.8331</td><td>3.9101</td><td>I</td></tr><tr><td>25</td><td>1.3497</td><td>0.6396</td><td>-0.9129</td><td>3.0641</td><td>I</td></tr><tr><td>26</td><td>1.3805</td><td>-0.0765</td><td>-0.2301</td><td>1.9645</td><td>I</td></tr><tr><td>27</td><td>1.6766</td><td>0.0924</td><td>2.7859</td><td>10.5810</td><td>III</td></tr><tr><td>28</td><td>0.9029</td><td>1.2889</td><td>1.0359</td><td>3.5496</td><td>I</td></tr><tr><td>29</td><td>-1.6151</td><td>-1.1481</td><td>0.1718</td><td>3.9560</td><td>I</td></tr><tr><td>30</td><td>-0.1513</td><td>0.4965</td><td>-0.5257</td><td>0.5458</td><td>I</td></tr></tbody></table>
# IV CONCLUSIONS

We propose a general positional accuracy control method, when we are interested in splitting points into three categories in respect with measurement errors. This method has the advantage that may be employed without any previous hypothesis about the undelaying distribution of errors, and it can be applied to 1D, 2D or 3D error data.

The method is an exact statistical hypothesis testing based on multinomial distribution with three parameters, N, $π_1, π_2$. In the case of error, we proceed to discretize data, through the classification into three categories given for two tolerances previously defined by the user. So, probabilities of the multinomial distribution are related to the probability of each category under the null hypothesis. The proposed statistical test is exact, so the p-value can be derived by exploring a space of solutions and summing up the probabilities of each isolated case of this space.

To probe the viability of this contrast a simulation procedure has been carried out and two examples are presented.

This contrast has some advantages, such as its easy realization and implementation, for instance in R program. Albeit we have applied it in a continuous underlying situation (measurement errors), it is also useful for discrete and categorical situations.

# ACKNOWLEDGMENTS

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

# REFERENCES

* Ariza-López, F.J., Rodríguez-Avi, J. (2014). A Statistical Model Inspired by the National Map Accuracy Standard. *Photogrammetric Engineering & Remote Sensing* 80(3), 271-281.

* ASCE (1983). *Map Uses, Scales and Accuracies for Engineering and Associated Purposes*. American Society of Civil Engineers, Committee on Cartographic Surveying, Surveying and Mapping Division, New York.

* ASPRS (1990). ASPRS accuracy standards for large scale maps. *Photogrammetric Engineering and Remote Sensing*, 56(7):1068-1070.

* ASPRS (2004). *ASPRS Guidelines: Vertical Accuracy Reporting for Lidar Data*. American Society for Photogrammetry and Remote Sensing. Bethesda, Maryland.

* ASPRS (2015). ASPRS Positional Accuracy standards for digital Geospatial Data, *Photogrammetric Engineering and Remote Sensing*, 81(3):53-75.

* Barnett, V., Lewis, T. (1994). *Outliers in Statistical Data* (3rd edition). John Wiley & Sons, Ltd, Chichester.

* Bonin O., Rousseaux, F. (2005). Digital terrain model computation from contour lines: How to derive quality information from artifact analysis. *GeoInformatica*, 9:253-268.

* Cayo M.R., Talbot, T.o. (2003). Positional error in automated geocoding of residential addresses. *International Journal of Health Geographics*, 2:10.

* FGDC (1998). *FGDC-STD-007: Geospatial Positioning Accuracy Standards, Part 3. National Standard for Spatial Data Accuracy*. Federal Geographic Data Committee, Reston, USA.

* Fisher, R.A. (1922). On the interpretation of $\chi^2$ from contingency tables, and the calculation of P. *Journal of the Royal Statistical Society* 85(1): 87-94.

* Freeman G.H., Halton J.H. (1951). Note on an exact test treatment of contingency, goodness of fit and other problems of significance. *Biometrika* 38:141-149.

* ISO (2010). *ISO 16269-4:2010. Statistical interpretation of data -- Part 4: Detection and treatment of outliers*. International Organization for Standardization, Geneva.
ISO (2013b). *ISO 3951-1:2013. Sampling procedures for inspection by variables -- Part 1: Specification for single sampling plans indexed by acceptance quality limit (AQL) for lot-by-lot inspection for a single quality characteristic and a single AQL.* International Organization for Standardization, Geneve.

Logsdon T. (1995). *Understanding the Navstar: GPS, GIS, IVHS*. Springer US.

Maune D.F. (ed) (2007). *Digital Elevation Model Technologies and Applications: The Dem Users Manual*. American Society for Photogrammetry and Remote Sensing. Bethesda, Maryland.

Müller M.J. (2001). Exact Tests for Small Sample 3x3 Contingency Tables With Embedded Fourfold Tables: Rationale and Application. *German Journal of Psychiatry*, 4(1):57-62.

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

STANAG (2002). *Standardization Agreement 2215: Evaluation of land maps, aeronautical charts and digital topographic data*. North Atlantic Treaty Organization, Bruxelles.

USBB (1947). *United States National Map Accuracy Standards*. Bureau of the Budget, Washington DC.
# Practical global elevation data error simulation

Ashton Shortridge, Joe Messina, Xue Li, and Nick Ronnei

Department of Geography, Michigan State University, USA

*Corresponding author: ashton@msu.edu

Digital elevation models (DEMs) contain error, and that this error results in uncertainty in scientific and applied applications that use DEMs. Although the fitness-for-use paradigm has been widely recognized in the geospatial data production community for over two decades, today there are no publicly-available datasets that actually meet this criterion. We report on two research threads to efficiently generate well-parameterized error models for medium-resolution global DEM products: the first thread concerns the development of error models to account for relationships between SRTM error and local environmental characteristics, and the second concerns the efficient generation of spatially autocorrelated error realizations.

## I GLOBAL TERRAIN DATA AND ERROR

It is widely recognized that DEMs contain error. This is a particular problem for terrain data products of global extent, as it is difficult to conduct rigorous “ground-truth” exercises to assess data quality in remote portions of the Earth. The impact of spatial data error on applications is typically not reported, and the burden of dealing with the potential impact of data error falls squarely on users of spatial data (Fisher & Tate 2006). Although the fitness-for-use paradigm has been widely recognized in the geospatial data production community for more than 15 years, today there are no publicly-available datasets that actually meet this criterion. We think this problem has three main components: the difficulty in translating the many DEM accuracy assessments into error models, the computational complexity of implementing Big Data stochastic simulation, and the difficulty in using realizations in scholarly and applied research with terrain data.

This paper reports on our team's progress towards addressing the first two of these challenges. We describe a global SRTM error model that captures important aspects of error using secondary, globally available data, and we introduce an efficient simulation implementation based on process convolution, enabling the rapid generation of very large error realizations. The paper concludes with comments on the method's limitations and potential extensions.

## II ELEVATION ERROR MODELING

In a continental-scale accuracy assessment, Shortridge & Messina (2011) linked SRTM error with globally available ecoregion, forest cover percentage, and topographic data. While a regression model with these variables (Table 6 in Shortridge & Messina 2011) accounted for substantial variance in SRTM error, the residuals were spatially autocorrelated, as Figure 2d demonstrates for a sample of locations in the United States. We extend this approach by proposing a basic regression kriging model (Hengl et al. 2007) to characterize SRTM or GDEM error as a linear function (parameters estimated via generalized least squares) of these variables along with a separate spatially autocorrelated error term, (calculated via simple kriging):

$$ \hat{y}(s_0)=\hat{m}(s_0)+\hat{e}(s_0) \qquad (1) $$
Where y-hat is the DEM error estimate at location $s_0$, m-hat is a linear model accounting for covariates measured at $s_0$, and e-hat is an estimate of the error in the linear model at $s_0$, with the estimate derived from simple kriging. This model can be developed in data-rich areas where reference elevations are available at locations $s_0$ and then efficiently deployed in areas without reference elevation data to generate realizations. The e-hat component is estimated via unconditional simulation, and individual elevation estimates are found as:

$$ \hat{z}(s_0) = DEM(s_0) - \hat{y}(s_0), (s_0) \in DEM \qquad (2) $$

where z-hat is an estimate of the unknown actual elevation at $s_0$ and DEM is the SRTM or GDEM elevation at $s_0$.

<figure><img src="image_138.png" /><figcaption>Figure 1: a) 3" SRTM DEM near Trafalgar, Indiana, USA; b) Resampled 3" NED DEM; c) 2000 MODIS forest cover %; d) SRTM errors.</figcaption></figure>

We illustrate this model with a case study on rolling terrain in south-central Indiana, USA near the town of Trafalgar. The region extent is 8,670 m x 7,700 m, comprising 10,043 pixels of 3" SRTM data (Figure 1a). For reference data we use USGS National Elevation Dataset (NED) data, resampled to 3" to match SRTM (Figure 1b); mismatches are apparent between these elevation surfaces. MODIS forest cover percentage for the year 2000 (Figure 1c) is clearly linked to larger positive SRTM errors (Figure 1d). Implementing the continental error model described above, parameters are multiplied by local SRTM slope, aspect, and MODIS forest cover percentage values to produce a predicted error surface (Figure 2a), which is subtracted from the original SRTM to produce a modeled DEM (Figure 2b). This model is not perfect; residuals from (SRTM - Predicted Error) - GDEM are presented in Figure 2c. RMSE of the original SRTM is 5.66 m, while RMSE for the model is 4.43 m.
Further, errors are spatially autocorrelated. Using a USA-wide sample of prediction errors, an average variogram model was formulated; this nationally-derived model is plotted against the empirical variogram of errors for the study site on Figure 2d.

<figure><img src="image_139.png" /><figcaption>Figure 2: a) Predicted SRTM error using Shortridge & Messina (2011) regression model; b) Model DEM (SRTM - Predicted Error); c) Residuals (Model - GDEM; d) Empirical variogram (dots) for the study area; model (line) derived from USA-wide sample of locations.</figcaption></figure>

### III EFFICIENT ERROR REALIZATION PRODUCTION

A significant challenge in global error propagation modeling via stochastic simulation is the computational complexity of implemented geostatistical simulation models. Generating dozens of error realizations for large (greater than one million cell) DEMs presents practical challenges for many algorithms. We implement an approach of Ver Hoef & Barry (1998), further developed by Higdon (2002) to construct realizations of valid spatial covariance models by passing specially constructed moving averages over gridded white noise. Convolution filters are commonly used by geospatial analysts for feature detection and classification in image processing, but their utility in error modeling is less well-developed (Oksanen & Sarjakoski 2005 is the notable exception). As explained in Higdon, particular covariance models (e.g., Gaussian with a specified range) may be convolved into specific smoothing filter. The result is an unconditional realization with spatial structure corresponding to a valid covariance model.

We have implemented this approach in R. Figure 3a demonstrates that realizations comprising millions of cells can be generated rapidly with this algorithm, making global deployment possible. The filter realization methodology can readily be adapted for specific areas. Figure 3b shows off a single error realization for the Trafalgar study site using a Gaussian covariance
model with parameters defined from a USA-wide sample: nugget 4.45, partial sill 7.41, and effective range 416 meters. One hundred similar realizations from this model for this study area were simulated in less than five seconds on the first author's older desktop PC.

<figure><img src="image_140.png" /><figcaption>Figure 3: a) Simulation times for spatially-autocorrelated fields using process convolution; b) Example error realization for the Trafalgar study site (UTM zone 16 projection).</figcaption></figure>

## IV DISCUSSION AND CONCLUSIONS

This brief paper demonstrates how DEM accuracy research may be developed into DEM uncertainty model research. Regression kriging offers an appealing framework for capturing important error properties of global DEM products and efficiently reproducing observed pattern in new areas. Further, the use of convolution filters to produce the spatially-structured random component is very promising. There are limitations to both approaches: RK can not work on non-colocated data, and it uses a very simple measure of spatial dependence. Process convolution models do not appear capable of conditioning realizations to locally available reference data, and several papers have identified approaches that capture pattern with greater fidelity (e.g., Bolin & Lindgren 2013). A significant limitation of the current implementation is that a single covariance model is employed. However, these models are computationally cheap and implementable on a large scale, and clear potential exists for nonstationary variogram modeling. Ultimately, this is part of a larger effort to enable efficient deployment of DEM uncertainty propagation; we expect to have a pilot system running very soon.

## References

* Bolin, D. Lindgren, F. (2013) A comparison between Markov approximations and other methods for large spatial data sets. *Computational Statistics & Data Analysis*, 61: 7-21.

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

* Hengl, T., Heuvelink, G. B., Rossiter, D. G. (2007) About regression-kriging: from equations to case studies. *Computers & Geosciences*, 33(10): 1301-1315.

* Higdon, D. M. (2002) Space and space-time modeling using process convolutions. In Anderson, C. W, Barnett, V, Chatwin, P. C and El-Shaarawi, A. H., (Eds) *Quantitative Methods for Current Environmental Issues*. London: Springer, 37-54.

* Oksanen, J. Sarjakowski, T. (2005) Error propagation of DEM-based surface derivatives. *Computers & Geosciences*, 31(8): 1015-1027.

* Shortridge A., Messina J. (2011) Spatial structure and landscape associations of SRTM error. *Remote Sensing of Environment* 115(6): 1576-1587.

* Ver Hoef, J. M. and Barry, R. P. (1998) Constructing and fitting models for cokriging and multivariable spatial prediction. *Journal of Statistical Planning & Inference*, 69: 275-294.
# How far is far enough? Towards an adaptive and “site-centric” modelling integrating co-visibility constraints for optimal land use

Valerio Signorelli*<fnref n="1" />, Thomas Leduc*<fnref n="1" />, Guillaume Chauvat<fnref n="2" />

*Corresponding author: valerio.signorelli@crenau.archi.fr

In this article, we propose a renewed site-centric solution that allows us to characterize a specific region of interest by defining the extent of the surroundings that influence sunlight exposure. The proposed method is a mix of an adaptive refinement and a visual-based clipping technique. This method has been implemented in the SketchUp context and applied to three sites located next to the French historical thermal town of Aix-les-Bains.

## I INTRODUCTION

Urban planning policies aim at defining the conditions of human settlements. Coherence and convergence of public action rely on the spatial continuity of its field of application for common issues. That is why French planning regulations tend to follow a concerted process led on several municipalities. They operate through planning regulations which are the expression of a political answer to issues emerging from territorial analyses. These analyses reveal the complexity of the territory by describing the local implications of the relations between cultural, physical and anthropological phenomena. Terrain features, settlement patterns, vegetation and infrastructures do not only influence environmental parameters, such as the amount of solar energy potential and daylight availability, in existing and planned urban fabric, but they also define the way in which inhabitants perceive their environment.

For this reason, urban planning practices should consider at the same time, and alongside urban regulations, perceptual, environmental and climate features in available and potential construction sites. Within this context solar exposure takes on great importance in terms of energy efficiency, quality of public and private spaces, and physical-perceptual enhancement of the local eco-system.

The availability of high resolution terrain and buildings models, the improvement in computation capability, and the development of 2.5D modelling simulation tools, based on image processing computation have provided, in the last decades, essential means for obtaining effective solar analysis from territorial to city scale (Prévot *et al.*, 2011; Morello *et al.*, 2010; Morello & Ratti, 2009; Floriani & Magillo, 2003; Tandy, 1967). However, simulations conducted on large and detailed raster grid models, the format commonly used by public agencies to deliver these basic data, are still costly operations, due to the amount of data needed, and tools based on image processing approaches do not permit yet an effective integration with vector-based solutions used in urban planning practices.

<footnote n="1">¹UMR AAU - CRENAU, Ecole Nationale Supérieure d'Architecture de Nantes, France</footnote>

<footnote n="2">²Cittànova, Nantes, France</footnote>
In large and complex topographic areas the influences of the various territorial features, even far from the actual position of the chosen sites, should be carefully considered. That is to say, by considering just the immediate surrounding of a site we can over- or under- estimate the solar contribution in terms of energy and sunlight. But which are the features to consider? Which is the level of detail they should provide? And how can we ensure an accurate simulation model able to provide reliable results in a reasonable amount of time? In other words, how far is far enough?

The aim of this explorative research is to propose a site-centric simplification method, based on 3D visibility analysis, in order to obtain a vector-based terrain model easier to handle, faster to compute. The entire process is integrated, as a series of extensions, for a well-known CAAD system, Trimble SketchUp. Three case studies, illustrating different levels of topographic constraints, will be used. Potentialities and limitations of the method will be highlighted in the discussion and further investigation will be proposed.

## II METHOD

### Pre-processing

The spatial datasets are provided by the IGN<fnref n="1" />, a national French institute in charge of the management and updating of geodesic and leveling networks, aerial photographs, and geospatial data. More precisely the aforementioned datasets are part of the French Large Scale Reference system (RGE): Digital Elevation Model or DEM (RGE® ALTI 5 m) in raster format for the representation of the landform (supposedly free of vegetation, buildings, etc.), and 3D vector models of significant spatial features such as footprints of individual buildings, forest cover, etc. (BD TOPO® 3D).

The needed pre-processing operations has been made using the Geospatial Data Abstraction Library GDAL/OGR (GDAL, 2016). Specifically the merging of the geospatial data, obtained through “gdal_merge.py” command line tool in the context of raster-based tiles or “ogr2ogr -update -append” command line tool in the context of vector-based layers; the clipping of data sources to some specified bounding box by using “gdalwarp” in the context of raster-based tiles or “ogr2ogr -clipsrc” in the context of vector-based layers.

To generate 3D vector contour files from the input raster DEM, the “gdal_contour.py” command line has been used and the resulting polylines have been simplified using the “ogr2ogr -simplify” tool. The simplified contour polylines are then reused to build the various Terrain Models presented hereafter.

### A two-step process: adaptive refinement and visual-based clipping techniques

The objective, after the conversion of the raster-grid model in a vector-based model, is to refine the virtual model of the terrain with an acceptable trade-off between the amount of data, and thus computation time, and data accuracy. Several terrain configurations have been produced and compared. However, only three of them will be detailed in this presentation.

<footnote n="1">¹ The Institut national de l'information géographique et forestière (National Institute of Geographic and Forestry Information, IGN), is a French public state establishment to produce and maintain geographical information for France.</footnote>
The “reference” Terrain model (M1) has the same planimetric resolution (5 meters) for the whole region, as provided by the original IGN dataset. In this first model there is no difference between studied parcels and the surrounding landscape, and the entire complexity of the terrain is represented. A side effect of this model is that even the areas that do not influence the chosen sites are precisely modelled. The “mixed” Terrain model (M2), combines the high resolution model of the three selected sites, with the low resolution model of the rest of the whole region. The third model is the “local” Terrain model (M3), where just the selected site and their immediate surroundings are considered. The mixed Terrain model (M2) is therefore some sort of intermediary between two extreme solutions. On the one hand, the "local" Terrain model (M3) does not take into account far distances' masks. On the other hand, the "reference" Terrain model (M1) is unnecessarily precise all over the wide region.

The three models have been developed using the contour lines obtained during the preprocessing phase, and then imported as SHP file (through a tool developed by one of the authors) in SketchUp. Through the existing tool “Sandbox From Contours”, the contour lines have been converted into a Terrain Model (Triangulated Irregular Network or TIN). Lastly, the building footprints have been imported as SHP file, drape on the terrain surface and extruded using the elevation values given as attributes.

The adaptive refinement of the Terrain models consists in the spatial union of two different sets of contours lines layers, with distinct spatial resolutions. Connections between this two datasets are automatically handled by the SketchUp “Sandbox from Contours” tool. The three images presented in Fig. 1 show a) the boundaries of the immediate surroundings of the studied site (the red polygon was obtained as a 300 meter radius buffer) b) the spatial union of the two contours lines datasets, with different resolutions, and c) a 3D sketch of the result in SketchUp.

<figure><img src="image_141.png" /><figcaption>Figure 1: The refinement technique used to adapt the resolution of the Terrain model. The contour lines outside the red polygon of the studied site, have been replaced by lower resolution contour lines.</figcaption></figure>

In order to further reduce the amount of data to be processed, a visibility-based clipping technique is used, considering that the hidden parts of landscape, from a given position, will not influence the studied area in terms of solar exposure.

A 3-step method has been implemented in SketchUp. First of all, we placed, over the three selected areas, a point grid with a fixed step equal to the resolution of the terrain obtained by IGN. For each of the sampling nodes, a 3D viewshed is computed using our extension based on the native ray casting engine of SketchUp. Finally, the spatial union of all these viewsheds is assessed and, to avoid any interpolation effect in all concavities during the TIN building phase, the convex hull of the resulting spatial union is delineated (see Fig. 2). The convex hull is then reused to clip the coarse-grained contour lines, and therefore (potentially) divide by two the area of the region to be taken into account.
<figure><img src="image_142.png" /><figcaption>Figure 2: The visual-based clipping technique shows the amount of contour lines, and therefore of terrain, that will not be taken into account in the simulation phase.</figcaption></figure>

## Post-processing

The three terrain models developed (M1, M2, and M3) have been used in the solar simulation and the obtained outcomes compared. Two indicators have been considered in order to test the reliability of our simplification method: the beam (direct) solar irradiation values (Wh/m²) and the daylight duration (min.). Both indicators consider a standard clear sky model and do not take into account sky or model reflections. We decided to conduct the simulation on December, 21st as the day with the lowest sun angles of the year, thus lowest amount of irradiation and daylight.

The measure of irradiation per unit area, depends obviously on the Terrain model itself. More precisely, it depends on its own direction towards the various sun positions, insofar as irradiation understood as the sum of instantaneous density of solar radiation incident on the surface over the given time period is the scalar product of the normal to the face by the sunlight direction).

## III USE CASE

The region of interest, of about 386.5km² (a 20 km-width square), is located in the French Alps, on the shores of the wide *Lac du Bourget*. In its south part, it embeds a 8-km wide valley oriented north-south which spreads between two mountain ranges (*Le Mont du Chat* and *Les Bauges*) peaking at around 1.5 km over the sea level. In the northern part, on the contrary, the valley broadens and the surrounding ridges sink down (see Fig. 3).
<figure><img src="image_143.png" /><figcaption>Figure 3: The region of interest with the 3 studied sites (Ontex on the NW side, Pugny-Chatenod on the East side, and Tresserve close to the Lac du Bourget).</figcaption></figure>

This region covers a specific inter-communality (the Communauté d'Agglomération du Lac du Bourget: Grand Lac) consisting of 17 municipalities. Planning regulations mostly tend to concentrate urban development around the administrative center of municipalities. This aims to allow planning optimization, spaces preservation for natural and agricultural purposes. On this constrained territory, land pressure requires a global project which can be embodied by many strategic locations. We propose to evaluate the sun exposure on those strategic parcels, in the existing fabric or on its boundaries.

In order to select relevant parcels, we assessed the expected impact of terrain on municipalities' center by evaluating the highest aspect ratio (H/W) to the closest relevant ridges (see Table 1).

<table><caption>Table 1. Characterization of the impact of closest ridges for each municipality center.</caption><thead><tr><td>Municipality</td><td>Aspect ratio</td><td>Municipality</td><td>Aspect ratio</td></tr></thead><tbody><tr><td>Tresserve</td><td>0.865</td><td>Bourdeau</td><td>0.274</td></tr><tr><td>Grésy-sur-Aix</td><td>0.579</td><td>Viviers-du-Lac</td><td>0.252</td></tr><tr><td>Méry</td><td>0.445</td><td>Aix-les-Bains</td><td>0.228</td></tr><tr><td>Drumettaz-Clarafond</td><td>0.433</td><td>Voglans</td><td>0.220</td></tr><tr><td>Le Montcel</td><td>0.430</td><td>Saint-Offenge</td><td>0.208</td></tr><tr><td>Mouxy</td><td>0.316</td><td>Brison-Saint-Innocent</td><td>0.184</td></tr><tr><td>Pugny-Chatenod</td><td>0.305</td><td>Ontex</td><td>0.147</td></tr><tr><td>Trévignin</td><td>0.295</td><td>La Chapelle-du-Mont-du-Chat</td><td>0.088</td></tr><tr><td>Le Bourget-du-Lac</td><td>0.284</td><td></td><td></td></tr></tbody></table>

Instead of *La Chapelle-du-Mont-du-Chat*, whose administrative center is constrained by the topography, we chose parcels in *Ontex*, on the slope of *Le Mont du Chat* (see Table 2, Fig. 4). We also chose locations on *Tresserve's* hill and in hillside *Pugny-Chatenod*. Those different locations all embody development opportunities inside the existing urban fabric.
<table><caption>Table 2. Various nomenclatures (models names) in use.</caption><thead><tr><th></th><th>Ontex site</th><th>Pugny-Chatenod site</th><th>Tresserve site</th></tr></thead><tbody><tr><th>Reference TIN</th><td>M11</td><td>M12</td><td>M13</td></tr><tr><th>Mixed TIN</th><td>M21</td><td>M22</td><td>M23</td></tr><tr><th>Local TIN</th><td>M31</td><td>M32</td><td>M33</td></tr></tbody></table>

<figure><img src="image_144.png" /><figcaption>Figure 4: Zoom in the three sites (from left to right: *Ontex*, *Pugny-Chatenod*, and *Tresserve*). The close horizon limits are represented by red circles.</figcaption></figure>

## IV DISCUSSION

Predictably, values obtained from M1 are the lowest and indices' values provided by M3 are the greatest (see Fig. 5). Indeed, a fine modeling of the mountainous Terrain model add new masks to the mock-up and therefore decreases the solar potential of the Terrain patches.

<figure><img src="image_145.png" /><figcaption>Figure 5: Comparison of the two indicators (duration of sunlight exposure and irradiance, Dec. 21st) for the three sites.</figcaption></figure>
The respective average elevations of the three sites vary from 310m (*Tresserve*), to 604m (*Pugny-Chatenod*), and 711m (*Ontex*). The corresponding standard deviation in each site is included between 2 and 3.8m (with a maximum value in the site of *Tresserve*). Height range between lowest and highest points within a 300m buffer varies along the three sites. This range is 115m (from 645m to 760m) high at *Ontex*, 125m (from 555m to 680m) at *Pugny-Chatenod* and 75m (from 250m to 325m) at *Tresserve*. On that last site, the closest obstacles do not alter direct irradiance, since they are located in a northerly area. The *Pugny-Chatenod* site is located on a West-East slope, which alters early morning exposure (lower energy input). Eventually, the parcels chosen in *Ontex* are affected by southerly masks (see Fig. 6), that reduce mid-day exposure (higher energy input).

<figure><img src="image_146.png" /><figcaption>Figure 6: The southern edges of the “cumulative viewshed” of the *Ontex* site (NW) are close from measurement locations (from 180m to 670m with a peak of about 1150m high located just south at about 2.5 km).</figcaption></figure>

As one can notice, for the site located on the top of a small hill close to the lake (*Tresserve*), the coarse-grained model M2 provides simulations results that perfectly match those obtained using the reference Terrain model M1. When the Terrain model gets hillier in the immediate surroundings, simulation outcomes are not as clear-cut. Thus, in the specific case of *Pugny-Chatenod*, where the site is located on a west-facing slope, the coarse-grained model M2 is accurate enough for the irradiation indicator (because it impacts mainly the early morning exposure with lower energy input). However, the mixed TIN M2 is obviously not precise enough to assess the precise daylight duration (there is indeed an average difference of about 30-minutes with the reference values). At last, in the particular case of *Ontex* site, where the Terrain model is particularly hilly nearby the measurement locations, the mixed TIN M2 is undoubtedly inaccurate.