Taking account of uncertainties in digital land suitability assessment

Brief review of land suitability evaluation

For a contextual background we provide a brief overview of land and enterprise suitability evaluation. Rossiter (1996) and Mueller et al. (2010) have discussed the key concepts. Henceforth, we use the term land suitability assessment (LSA) as a collective term for its other connotations such as land evaluation, enterprise suitability evaluation or assessment or any such similar term that more-or-less conforms to the definition given by Van Diepen et al. (1991) as “all methods to explain or predict the use potential of land”. While there is a blurry concept of what ‘land’ is or means, we take its meaning from Triantafilis & McBratney (1993), paraphrasing from Brinkman & Smyth (1973) as “all reasonably stable, or predictably cyclic attributes of the biosphere above and below the earth surface, including those of the atmosphere, the soil and underlying geology, the hydrology, floral and faunal populations, and the result of previous human activity in which the effects significantly influence present or future land uses by man (sic)”.

Soil generally features heavily in a LSA with the reason possibly traceable to the rudimentary beginnings of soil science (and pedology) and subsequent grouping and classification of soils on the basis of their productivity (Brevik & Hartemink, 2010). Incorporation of climate information into LSA is also apparent. Climatic information can include variables such as temperature and precipitation (Zabel, Putzenlechner & Mauser, 2014), to other more complexly derived variables such as frequency of frost occurrence (e.g., Kidd et al., 2012). Less common is the incorporation of socioeconomic factors (D’Haeze et al., 2005). In general, greater understanding is needed in the area of coupled human and biophysical systems (Stuart et al., 2015). Intuitively, socioeconomic factors would feature heavily in a LSA because they compel activity despite biophysical constraints. Some factors include commodity price changes, institutional reforms and reforming of trade practices and agreements amongst governments and organisations. However, it is not altogether clear how these types of factors can be incorporated in the more common biophysical centric LSAs as they are often very localized to a particular area, and difficult to track.

LSA projects are practiced throughout the world due principally to the need for better management of natural resources which are finite—yet the pressures and demands emplaced upon them are increasing—for example, through human population increase. Underlying this is a fundamental need to determine whether a given site or locality would be productive or unproductive for a given land use, and then secondly evaluate its biophysical status to determine what the limitations are, and what the land may be best suited for. Broad global LSA programs include agro-ecological zoning (Fischer & Sun, 2001) and the Fertility Capability Classification (Buol et al., 1975). More recently a global high resolution LSA for a number of arable cropping enterprises together with future projections was presented by Zabel, Putzenlechner & Mauser (2014).

Mueller et al. (2010) describe numerous examples of LSA programs that are nationally rolled out throughout the world. There are also LSA projects that are quite locally situated. For example, LSA projects have been developed in both Northern Australia (Queensland) and Tasmania respectively as described in Harms et al. (2015) and Kidd et al. (2012) for enhancing agricultural pursuits in those areas.

While LSA may have had its origins in the early development of soil science and agricultural land capability classifications, the terminology and general framework were formalized with the land evaluation guidelines prepared by the Food and Agriculture Organization of the United Nations (FAO) in 1976 (FAO, 1976). These guidelines have strongly influenced and continue to guide LSA projects throughout the world. They have also informed the development of automated LSAs via such software as ALES (Automated land Evaluation System, Rossiter, 1990) and MicroLEIS (De la Rosa et al., 2009) as examples. At its core, the FAO framework is a crop specific LSA system with a 5-class ranking of suitability (FAO Land Suitability Classes) from 1: Highly Suitable to 5: Permanently Not Suitable. Given a suite of biophysical information from a site, each attribute is evaluated against some expert-defined thresholds for each suitability class. The final evaluation of suitability for the site is the one in which is most limiting. There are a number of variants to this approach such as the number of suitability classes, the weighting attributed to the inputs, and also the degree of complexity of the crop specific suitability criteria. This last variant is usually determined by availability of data, which is also related to spatial scale at which the LSA is applied. As we have seen recently however in Zabel, Putzenlechner & Mauser (2014) and with high resolution global soil mapping efforts such as the GlobalSoilMap project (Arrouays et al., 2014), the scale issue will not be as important into the future. A further generalisation of the FAO framework is evidenced by the incorporation of fuzzy threshold values for suitability parameters, as opposed to discrete threshold values (Lark & Bolam, 1997; Triantafilis & McBratney, 1993; Zabel, Putzenlechner & Mauser, 2014). Fuzzy cut-offs are intuitively sensible because one cannot easily justify (as an example) that a crop will or will not grow if a subsoil has pH of 5.4 as opposed to 5.5, given the complexity of soil and the interactions of other soil physical and chemical parameters in the system. These ideas are explored in Lark & Bolam (1997) with examples.

Other LSA frameworks include parametric grading systems that may be multiplicative or additive which ultimately rate the suitability of a land on a continuous scale rather than as discrete classes. The Storie index rating (Storie, 1932) of soils is an example of a multiplicative parametric LSA. A more data driven approach to the Storie Index for digital soil applications is described in O’Geen, Southard & Southard (2008). An additive parametric system was developed for NSW, Australia by Zhang (1989). The last LSA frameworks to be discussed are those associated with crop growth and ecosystem models that estimate local productivity of specific crops given soils and weather information (Mueller et al., 2010). These models are more sensitive to temporal variation, which make them a more dynamic LSA approach to those discussed already. While this may seem advantageous, their dynamic behavior can only be realized given appropriate input information which is often difficult to obtain. Mueller et al. (2010) describes a number of candidate models of this type that are applicable at different scales from global to local extents. The widespread application of these models in LSA projects is limited due to their sophistication and inflexibility about input data requirements as many are designed for a specific purpose and applicable only to their origin.

Other than fuzzy classification methods (Triantafilis & McBratney, 1993; Zabel, Putzenlechner & Mauser, 2014), the assumption of LSAs is that the input data and parameter thresholds are free of error. However, soil is complex and can vary quite erratically in the context of space and time (Webster, 2000), and subsequent model-based predictions of soil phenomena are anything but ‘error free’ (Brown & Heuvelink, 2005). In noting this, Harms et al. (2015) provided additional mapping of a confidence measure (based on Mahalanobis distance calculations) of the suitability classifications. They were able to state explicitly where suitability classifications were likely to be good and where they were likely to be uncertain. Indeed this is helpful and a good step forwards, but possibly a more direct way to derive a measure of land suitability is to incorporate the measures of prediction uncertainty into the final assessment.

Digital soil and climate modeling for land suitability assessment

An extensive soil sampling and climate modeling program was established in the Meander Valley in 2010 (Kidd et al., 2012; Kidd et al., 2015). In total 576 soil cores were extracted from various locations throughout the area and analysed for a number of physical and chemical soil properties both in the laboratory and with chemometric techniques dependent on mid-infrared soil spectral calibrations (Kidd et al., 2015). Climate and temperature monitoring throughout the Meander was carried out over 2010 and 2011 using a network of distributed temperature sensors (Webb et al., 2014), and calibrated to long-term climate data (Webb et al., 2015). Rainfall information was sourced from Australian Bureau of Meteorology monitoring sites located within and surrounding the study area (Webb et al., 2014). Soil sampling and climate monitoring network were conducted and installed respectively by field officers of the Tasmanian Government. For many sites that were situated on public or Government administered land, no permission was required for sampling or installation of temperature sensors. For sites that were situated on private landholdings, Government officers sought permission for access. In the rare situation where permission was not granted, the sampling or monitoring site was moved to an alternative and agreeable location.

The common workflow for all digital soil and climate variable mapping entailed:

1. The randomized splitting of observational data into calibration and validation datasets. Here a 75% and 25% split was used respectively for calibration and validation datasets. For consistency, the same calibration and validation datasets (soil, climate, rainfall) were used for all target variables, before removal of missing values.

2. Environmental covariates were sourced from DPIPWE and other government repositories which included derivatives from a digital elevation model (STRM DEM (Gallant et al., 2011)), gamma radiometric information (Minty et al., 2009), and spectral indices derived from Landsat 7 ETM+ satellite. Soil modeling involved using principal components of all sourced covariates, while principal components of the digital elevation model derivatives (only) were used for the climate variables (Webb et al., 2015).

3. Modeling of continuous variables was based on a regression kriging framework that entailed Cubist regression tree modeling (Quinlan, 1992) followed by model residual modeling (with variograms) and kriging. Spherical or exponential models were considered only. Visual criteria of the global variogram of residuals were used to determine whether regression kriging should be pursued or not. Otherwise regression modeling was used only. Categorical variable modeling entailed either the fitting of binomial or ordinal logistic models, dependent on the nature of the target variable information.

4. Prediction uncertainties for continuous variables were quantified using an empirical approach as described in Malone et al. (2014) where the model errors within each partition of a Cubist model were used to form geographically specific error distributions (via leave-one-out cross validation) in order to estimate 90% prediction intervals. For categorical variables, the prediction probabilities were used as measures of uncertainty.

5. Validation statistics for continuous variables included the root mean square error (RMSE) and Lin’s concordance correlation coefficient (CCC; Lin, 1989). The prediction interval coverage probability (PICP, Shrestha & Solomatine, 2006) was used to evaluate the efficacy of the 90% prediction intervals. The PICP is the proportion of actual observations that are encapsulated by their prediction interval, and ideally will be equivalent the level of confidence associated with the prediction interval. For categorical variables, overall accuracy and kappa statistic (Congalton, 1991) were used. Validation statistics are reported for both calibration and validation datasets. The PICP is reported for the validation data set only.

All maps of the biophysical properties (together with uncertainties) used for the hazelnut LSA are supplied as Supplemental Information 1 to this research. Some details specific to each modeling variable are as follows.

Soil depth

Soil depth was modeled in a two-step procedure. The first step was binomial modeling of whether soil depth greater than 1.5 m or not. The rationale behind this was that soil coring was done to a maximum of 1.5 m depth, or depth-to-lithic contact, whichever was first. The second step involved regression kriging modeling of soil depth where soil depth was less than 1.5 m. The outputs from both steps were used for the LSA, which is discussed further on.

The two-step modeling procedure described in this study has been used previously by Gastaldi, Minasny & McBratney (2012) for mapping the occurrence and thickness of soil horizons. In that study, the first step involved modeling the occurrence of horizon classes, while the second step involved modeling the depth of the soil horizons.

Soil pH, EC and clay

A mass-preserving soil depth spline (Bishop, McBratney & Laslett, 1999) was used to harmonise observed soil profile data in order to impute data for the 0–15 cm depth intervals for all locations. Spatial modeling was performed using regression kriging framework.

Soil drainage

Digital soil mapping of soil drainage was carried out as in the method described by Kidd et al. (2014a) and Kidd et al. (2014b). Essentially this method codifies the descriptive soil drainage classification as detailed in The National Committee on Soil and Terrain (2009) into an ordinal value. The numerical classes were then spatially modeled as a continuous variable using regression kriging.

Stoniness

Soil stoniness or percentage of coarse fragments greater than 2 mm was modeled in a two-step procedure, much for the same reasons as for soil depth. Here the first step was a binomial logistic regression model which considered the presence vs. absence of coarse fragments. The second step entailed an ordinal logistic regression of the coarse fragment incremental percentage classes (6) as described in The National Committee on Soil and Terrain (2009).

Incidence of frost

The method for estimating the probability of frost incidence at each location is described in Webb et al. (2014) and Webb et al. (2015). Spatial modeling occurred as a two-step procedure because many sites had no incidence of frost occurrence. Subsequently the first step entailed a binomial logistic regression of presence vs. absence of frost occurrence followed by the second step of regression kriging of frost occurrence probability using the sites where information relating to this was recorded.

Temperature, rainfall and chill hour requirements

These were all spatially modeled as continuous variables. A multiple linear regression model was used for rainfall modeling, due to the small number of rainfall observations. The validation statistics reported are from a leave-one-out-cross validation. The PICP for the validation of prediction uncertainties is based on the calibration data. Cubist regression models were used for temperature and chill hours, with residual kriging being appropriate for chill hours only.

Approach 1. Land suitability assessment without considering prediction uncertainties

Raster maps of each of the hazelnut LSA variables were interrogated pixel-by-pixel from which an assessment of suitability was derived using the most limiting factor approach. For variables that were modeled as a two-step procedure (soil depth, stoniness, frost incidence), a positive condition invoked the interrogation of the map from second part of the two-step procedure. Using the soil depth variable as an example, if a pixel on the prediction map recorded a positive score (likely presence that depth to lithic contact is less than 1.5 m) contact, the second map was interrogated to estimate the depth to the lithic contact.

Approach 2. Land suitability assessment in consideration of prediction uncertainties

The basis of the approach is a simulation of possible realisations (pixel-by-pixel) of the input variables before assessing the suitability using the most limiting factor approach. For the continuous variables the prediction intervals were sampled with an assumption of normal distribution upon each realisation. For the categorical variables, sampling from either the binomial or multinomial distributions of the prediction probabilities was performed. In consideration of variables subject to the two-step model procedure, the second condition was invoked where a presence or occurrence of the phenomenon in question was found for a given sample. If the condition was invoked, the probability distribution of the second variable was sampled.

There were two issues to consider in this approach. The first was an issue of computation where compared to approach 1, the computation time can increase dramatically with increasing number of realisations. Besides using multi-core compute facilities, Latin Hypercube Sampling (LHS) was used to sample from the multivariate distributions of the continuous variables. LHS sampling is a more efficient sampling approach to random sampling if the objective is to ensure the multivariate distributions are sampled entirely (Pebesma & Heuvelink, 1999).

The second issue was one of rationality in terms of the multivariate information generated from each realisation. The rationality here was the maintenance of correlations between the LSA variables. Subsequently a modified LHS was used where the correlation matrix of the sample is corrected to match that of the data that is being sampled. The correction was made through the use of the Huntington–Lyrintzis algorithm (Huntington & Lyrintzis, 1998). Implementation of correlation matrix corrected LHS was performed using of the LHScorcorr function in the pse R package (Chalom & Prado, 2014).

In this study 100 realisations were made in order to estimate probabilities of each suitability class at each map pixel. The probabilities were estimated as the number of occurrences of each suitability class divided by the total number of realisations. Further analysis of these outputs is presented in the results.

Software

The entire workflow of this study was implemented in the R scripting language (R Core Team, 2015) except for the creation of map products which were created using ArcGIS (ESRI Inc.). Specifically, custom functions were built for the hazelnut LSA. The raster (Hijmans, 2015) package with associated rgdal (Bivand, Keitt & Rowlingson, 2015) and sp (Bivand, Pebesma & Gomez-Rubio, 2013) packages were used for handling and manipulating all GIS processes. Cubist models were implemented from the Cubist (Kuhn et al., 2014) package and the multinomial and ordinal logistic modeling were carried out using the nnet and MASS packages respectively (Venables & Ripley, 2002). Variogram modeling entailed the usage of geostatistical functions from gstat (Pebesma, 2004). Many of the processes, including applying models spatially and applying the LSA in the presence of prediction uncertainties were distributed across multiple compute nodes (8) to improve the computation efficiency. This was done using specialist R packages parallel (R Core Team, 2015) and doParallel (Revolution Analytics & Weston, 2014).