Assessing the effectiveness of ground truth data to capture landscape variability from an agricultural region using Gaussian simulation and geostatistical techniques.

Predictive modeling with remotely sensed data requires an accurate representation of spatial variability by ground truth data. In this study, we assessed the reliability of the size and location of ground truth data in capturing the landscape spatial variability embedded in the Airborne Visible Infrared Imaging Spectrometer-Next Generation (AVIRIS-NG) hyperspectral image in an agricultural region in Anand, India. We derived simulated spectral vegetation and soil indices using Gaussian simulation from AVIRIS-NG image for two point-location datasets, (1) ground truth points from adaptive sampling and (2) points from conditional Latin Hypercube Sampling (cLHS). We compared values of the simulated image indices against the actual image indices (measured) through the analysis of mean absolute errors. Modeling the variogram of the measured indices with the hyperspectral image in high spatial resolution (4m), is an effective way to characterize the spatial heterogeneity at the landscape level. We used geostatistical techniques to analyze the shapes of experimental variograms in order to assess whether or not the ground truth points, when compared against the cLHS-derived points, captured the spatial structures and variability of the studied agricultural area using measured indices. In addition, we explored the capability of the variogram by running tests in different point sample sizes. The ground truth and cLHS datasets were able to derive equivalent values for field spatial variability from image indices, according to our findings. Furthermore, this research presents a methodology for selecting spectral indices and determining the best sample size for efficiently replicating spatial patterns in hyperspectral images.

Introduction

Remote sensing (RS) technology offers advantages in high-resolution soil property mapping and analysis of spatial soil data across multiple farms (Robert 2002; Seelan et al., 2003; Brisco et al., 1998; Duffera et al., 2007). For predictive modeling of target soil variables and validation of modeled data, utilizing the capability of RS in soil mapping frequently begins with extensive ground truth data and information collection (Plourde and Congalton, 2003; Ge et al., 2006). However, in practice, intensive data collection takes time and is often costly, necessitating the development of an optimal sampling method to create representative ground truth (GT) soil data (Mulla, 2013).

Field sampling strategies are essential in generating reliable remotely sensed products (Stehman and Czaplewski, 1998; Mandal and Ghosh, 2000). The ideal sample design is one that is statistically sound, simple to apply, economical, produces trustworthy estimates, and minimizes sampling errors (Stehman, 2001; de Gruijter et al., 2006). There are a number of soil sampling methods suggested in literature that minimize the errors, optimize the sample size, and identify the representative sampling locations. The vast literature on spatial sampling frameworks can be broadly categorized into 1) design-based (probability) (DB) strategy and 2) model-based (prediction) (MB) strategy (Brus and de Gruijter, 2011).

The DB sampling can be a simple random, a systematic, a stratified random or an adaptive sampling (Cochran, 1977; Stehman and Czaplewski, 1998). For a simple landscape, a probability based stratified random sampling approach may be sufficient (Miklos et al., 2010; Arrouays et al., 2017). However, for complex landscapes and terrains, adaptive DB sampling may produce better representatives of target variables (Graniero and Robinson, 2003). In the MB sampling, the relationship between the sample and population properties are quantified by the model, and the model prescribes ideal sample locations (Stevens, 2006). For example, the conditioned Latin hypercube Sampling (cLHS) is a frequently used MB sampling design in digital soil mapping (Gao et al., 2016). The choice between the two sampling strategies depends on the goals of the study, the availability of legacy and auxiliary data, availability of local resources, and local accessibility (Biswas and Zhang, 2018). A comparison study by (Wang et al. (2012)) elaborated that the DB strategy is more suited for confronting the “how much” questions (estimating global quantities, frequency distribution of the target variable), while the MB is more on the “where” questions (predict values at unknown locations). Lesch (2005) and Spöck (2012) reviewed the details of both frameworks.

Whether the DB or MB sampling is used, it is generally recognized that collecting more samples gives a better reconstruction and more reliable mapping of the target soil property (Dicks and Lo, 1990; Vasat et al., 2012). For example, to generate high resolution soil organic carbon (SOC) maps, (Guo et al. (2006)) had to rely on finer-scale inventories and increased sample size to estimate carbon values. The authors acknowledged that the approach was both labor and time-intensive, and sometimes redundant. Because available resources and, in some cases, site accessibility influence soil sample size, determining the optimal representative sample size and locations could be difficult (Kumar et al., 2016). A sample size identification method was proposed by An et al. (2018) where sizes of representative samples were selected differently based on soil property variation and environmental covariates. However, the spatial relationships between a soil property and the environmental covariates are complex and various. This complexity could result in uncertainties in soil property estimates.

Often the source of uncertainty stems from the suboptimal field sampling design (Dicks and Lo 1990; Ramsey, 1998). The optimum number, spacing, and frequency of samples play crucial roles and, if not carefully evaluated, could contribute to large errors (due to distribution heterogeneities and spatial trends) in estimates of target variable (Van der Perk et al., 2008). Thus, an assessment and quantification of uncertainties from GT data measurements following a GT campaign for predictive modeling is critical (Mulla 2016). In sampling error analysis, these uncertainties are often expressed as variances, , estimated by SD2 (the square of the standard deviation) (Jarvis et al., 2004). For decades, the accuracy of field sampling for soil properties has been of concern because representative samples may not provide sufficient information (Cameron et al., 1971; Buscaglia and Varco, 2003; Kurfürst et al., 2011). Buscaglia and Varco (2003) highlighted that spatial variability within agricultural fields, for example, influences trueness and precision of predictions based from a bulk of GT samples collected.

Generally, there has been little attention to in-depth analysis of the collected GT soil samples and its ability to effectively capture the variability in the landscape of interest. Most RS literature deal with improving sampling technique rather than on quantifying the errors and uncertainty introduced as a result of suboptimal sampling before developing predictive models (Zorzi et al., 2002; Minasny et al., 2013; Eurachem, 2019). To describe the uncertainty of measurements is to quantify the closeness or similarities of the simulated RS image indices from the GT data and the measured values from the imagery itself (Bell, 2001; Kurfürst et al., 2011). There are existing studies on spatial similarities and differences in soil properties from GT samples, for example, soil pH (Laslett et al., 1987), soil zinc (Eze and Kumahor, 2019), soil phosphorous ((Boluwade and Madramootoo, 2013)), and soil carbon (Zhang and McGrath, 2004). These studies resulted in spatially mapped distribution of soil properties and quantified uncertainties and correlations in spatial models (Eze and Kumahor, 2019).

The need for hyperspectral and high-spatial resolution satellite imagery continues to increase because of its potential for producing reliable maps of soil properties (Denis et al., 2014; Franceschini et al., 2015). In addition, information derived from remote sensing including crop and soil-related spectral indices, could provide sufficient information that is highly correlated with the target farm soil properties and characterize the spatial heterogeneity of the landscapes. The spatial heterogeneity of the terrain is particularly crucial when selecting a sample technique that captures the spatial scale of the target property while simultaneously maximizing field collection resources (Stein and Ettema, 2003). In this paper we assessed the reliability of the GT data locations collected during a field campaign in Anand, India, in effectively capturing the variability observed in the airborne Visible Infrared Imaging Spectrometer-Next Generation (AVIRIS-NG) data over the study region. We quantified the uncertainty in the landscape by (1) analyzing experimental variograms and evaluating the underlying spatial structures of image indices generated from the cLHS and GT datasets and (2) comparing absolute errors between simulated and measured values. By utilizing Gaussian simulation and geostatistical technique analysis, we determined the optimal number of samples to effectively capture the variability in the study area.

Materials and methods

Study area

This study is located in the Anand District in the state of Gujarat, India (Figure 1). The site covers approximately 216 square kilometers (21,600 ha). The investigated site represents farmlands with diverse agricultural management and land-use systems.

Figure 1

(a) Map showing the location of the study site in the Anand District in Gujarat, India and, (b) sample distributions of the 1,000 points derived from cLHS sampling method (red) and those collected from the farm fields (yellow). The gray patches are urban areas. Crops in the sampled fields were either (c) sparse or (d) non-existent.

General flow of the study

Figure 2 provides the overall flow of the adopted methodology, which is divided into three major steps:

Figure 2

A three-step methodology that includes the data collection, analysis, and interpretation of results using two sets of data points from the cLHS and GT.

Data collection: This included the preprocessing of AVIRIS-NG and GT datasets by subsetting the image to the study area and formatting the GT data.

Data analysis: This included deriving image indices, running an R algorithm and generating cLHS points, obtaining measured values for GT and cLHS points by overlaying them over the image indices, producing simulated image indices using the GT and cLHS points, computing the mean absolute error through comparison between the simulated and measured values, and applying variogram using the GT and cLHS points in order to assess whether or not the datasets behaved similarly in capturing the spatial variability in the study area through the image indices;

Interpretation of results: This included the assessment of the similarities or differences between the GT and cLHS points, and determining the optimal number of samples that effectively capture the spatial variability in the study area.

We executed the steps by using not only one sample size, but four different sample sizes: 100, 400, 700, and 1000 points.

AVIRIS-NG hyperspectral image

For this paper, we used the high-resolution airborne hyperspectral imagery AVIRIS-NG that was collected over India in March 2018 as part of a cooperative mission between National Aeronautics and Space Administration (NASA) and Indian Space Research Organization (ISRO). The AVIRIS-NG has a spatial resolution of 4 m and samples 425 spectral bands over the 380 nm–2510 nm range at approximately 5 nm spectral resolution. Geometrically corrected level-2 reflectance data for the study area was downloaded from the NASA server (https://aviris-ng.jpl.nasa.gov).

Ground truth dataset

Field campaign took place during the summer of 2019 when soil conditions were mostly dry at the study area in the Anand District in Gujarat, India. Together with our collaborators from Anand Agricultural University (AAU), we collected a total of 1000 soil samples (Figure 1) in a span of two weeks through intensive adaptive field-or farm-level sampling in order to detect differences between sampling strategies. The planning phase was the most important part of the fieldwork. Due to accessibility and complexity of the local road networks in the study area, we prepared detailed daily driving routes before heading to the field and to avoid unnecessary incidents. The sampling design was modified in real time as data collection continued, in order to improve the representativeness of the dataset. We ensured that the collection of samples followed the required sampling protocols or techniques (Jahn et al., 2006; Massawe et al., 2016). In order to sample points from non-urban areas only, we ran an algorithm to remove settlements and urban areas from the image. We also collected data about crop types and agricultural practices by interacting with local farmers.

Image indices

Image indices are an effective way to distinguish geographic characteristics. We chose our spectral indices (Table 1) based on two prior studies that mapped agricultural crops, soil, and tillage using AVIRIS (Salas and Subburayalu, 2019) and AVIRIS-NG (Salas and Subburayalu, 2019). The Photochemical Reflectance Index (PRI) (Gamon and Serrano, 1997) was the most accurate predictor of crop primary productivity for wheat, legumes, and eggplant. Over the last two decades, PRI has gotten a lot of interest for remote sensing mapping of various vegetation functional types (Zhang et al., 2016). The Cellulose Absorption Index (CAI) (Daughtry 2001) was sensitive to decreasing crop cover and residue and performed well in spectrally separating tobacco from other crops (Salas and Subburayalu, 2019). In mapping eggplant, tobacco, and legumes, the Moment Distance Index Normalized (MDIN) was discovered to be a useful predictor (Salas and Subburayalu, 2019). The Normalized Difference Vegetation Index (NDVI) (Tucker, 1979) has already been demonstrated to provide useful information on plant growth patterns, physiological state, and biochemistry (Roberts et al., 2012; Burai et al., 2015). The modified soil-adjusted vegetation index (MSAVI) reduces the impact of soil background on vegetation sensitivity (Qi et al., 1994). Chen et al. (2004) proposed using remote sensing data to map seasonal bare soil patches using the bare soil index (BSI). Finally, using experimental variograms, we evaluated the underlying spatial features of the agricultural study site using all the derived image indices.

Table 1

Spatial heterogeneity is described using image indices. Spectral bands were averaged to represent NIR (750–850 nm), red (600–700 nm), green (500–600 nm), and blue (400–500 nm) for broadband indices. The represents the AVIRIS-NG's target wavelength.

Index	Equation	Source
NDVI	Red−NIRRed+NIR	Tucker (1979)
CAI	100[(λ2030+λ2210)−λ21002]	Daughtry (2001)
PRI	λ529−λ580λ529+λ580	Gamon and Serrano (1997)
MDIN	MDRPMDLPMD_RP=∑i=λ_RPλ_LP(ρ_i^2+(λ_RP−i)^2)^0.5MD_LP=∑i=λ_LPλ_RP(ρ_i^2+(i−λ_LP)^2)^0.5	Salas and Subburayalu (2019)
MSAVI	2NIR+1−(2NIR+1)^2−8(NIR−Red)2	Qi et al. (1994)
BSI	(Red+Blue)−Green(Red+Blue)+Green	Chen et al. (2004)

Image spatial heterogeneity

With variogram (Curran and Atkinson, 1998), we quantified and assessed the spatial heterogeneity of the study area from the image indices using the locations and values of the cLHS and GT points. Variogram provides quantitative tools to measure spatial autocorrelation simply by characterizing its parameters to quantify the distinctive spatial properties of the image indices. Current literature used variogram in characterizing spatial patterns of vegetation cover using medium spatial resolution (Garrigues et al., 2006; Ding et al., 2014) and high spatial resolution (Nijland et al., 2009; Gómez et al., 2012) image indices.

The examination of the autocorrelation typically observed in geographic data, where data values from sites close to each other are more comparable than data values from sites far apart, is included in variogram analysis (Isaaks and Srivastava, 1989). A measure of a variable's correlation with its spatial location is called spatial autocorrelation structure and is depicted in Figure 3. The experimental variogram is created by calculating the average squared difference in data values for every pair of data sites along a specified direction and graphing these values against the distances between the data pairs. This could be expressed using Eq. (1):where x is a data location, h is a lag vector, Z(x) is the data value at location x, N is the number of data pairs spaced a distance and direction h units apart.

Figure 3

The typical variogram and the important parameters (range, sill, nugget) for spatial autocorrelation analysis.

Afterwards, the experimental variogram is utilized to fit a theoretical variogram model with known mathematical features (Garrigues et al., 2006). Variogram values can be used in a variety of calculations, including estimating at unsampled places. In this work, an exponential model is used.

The key parameters of the modeled variogram include the following: range, the distance beyond which data are no longer spatially correlated; (b) nugget, the level of random variation within the data; and (c) sill: the magnitude of variation present at which the variogram levels off at the range. The value of the variogram at a distance of 0 is 0 by definition (i.e., data points for the same location are identical). The nugget effect is the jump discontinuity that denotes small-scale spatial disparities within the fields at scales lower than the pixel size. Chilès and Delfiner (1999) and Isaaks and Srivastava (1989) used a variogram to describe the nature of the data's geographic variability.

Conditional Latin Hypercube Sampling (cLHS)

The conditional Latin Hypercube Sampling (cLHS) is a prominent model-based sampling method that employs a stratified random technique to optimize sample selection using continuous and/or categorical covariates as input (Minasny and McBratney 2006). The cLHS is frequently utilized in digital soil mapping and soil property prediction applications (Ng et al., 2018). When compared to Monte Carlo sampling, the cLHS has proven to improve the sampling scheme and reduce the computational overhead (Yin et al., 2011). In a benchmark study on comparison of sampling techniques, Santos and Beck (2015) found that, while Importance Sampling and Subset Simulation were efficient sampling techniques in Monte Carlo Simulation, the use of Latin Hypercube Sampling had a significant and positive influence for all sampling techniques.

We used the cLHS approach to select sample sizes – 100, 400, 700, and 1000 samples – from each input image index. We computed the sample sizes by doubling the confidence interval, starting at 2, then 4, and 8 in order to capture the potential variabilities in a gradient field. The cLHS procedure follows these steps (Minasny and McBratney 2006): Given N sites with ancillary data (W), select n sample sites (n ⪡ N) so that the sampled sites w form a Latin hypercube in the feature space, or the multivariate distribution of W is maximally stratified. For k continuous variables (e.g., image index), each component of W is divided into n (e.g., sample size = 1000) equally probable strata based on their distributions, and w is a sub-sample of W.

Sequential Gaussian simulation

Delbari et al. (2009) emphasized that the main idea behind sequential Gaussian simulation (SGS) is that the conditional distribution of the observed variable could be used to simulate successive grid points. SGS has been used to explore the geographical variability and scarcity of soil moisture (Zhang et al., 2017), geochemical anomaly (Chen et al., 2013), soil carbon release (Teixeira et al., 2011), geographical patterns of soil organic matter (Chai et al., 2007), and soil salinity (Zhao et al., 2018).

SGS works in the following order: (a) normalize the observed values, (b) compute the experimental variogram using the normalized values, (c) improve the sampling grid by defining a random path through all nodes, (d) pick a random node in the sampled grid, (e) use kriging to calculate the estimated value and its variance, (f) a normal distribution obtained from the data in step 5, (g) determine the cumulative distribution function, (h) remove a random observed value of this function; this value represents the investigated attribute in the same location as it was removed, (i) reiterate steps 4–8 until you have visited all of the points, (j) reiterate steps 4–9 until all of the realizations have been completed, and then (k) transform the obtained values for the original data's mean and variance. SGS creates a mean map of the n realizations of the examined index after all of the processes are completed. Goovaerts (1998) has a comprehensive account of the SGS process.

Using the cLHS and GT datasets, we created simulated image indices and compared descriptive statistics. The values of the simulated image indices were then compared to the measured image indices using mean absolute errors (MAE) analysis.

Comparison between simulated vs measured

The MAE is a useful measure widely reported for model evaluations in the field of geosciences (Gallaun et al., 2010; Heiden et al., 2012). MAE has been used to represent average difference, instead of the average error, when no set of estimates is known to be the most reliable. Calculation of MAE is fairly straightforward. It involves summing the magnitudes of the errors to obtain the ‘total error’ and then dividing the total error by the number of observations (Eq. 2).where y is the actual value of the index and is the predicted value, and n is the number of reference sample sizes.

We compute the MAE as the absolute value of the difference between the simulated value and the actual value. Smaller MAE values could only mean that the actual and simulated images capture the same spatial variability in the study area.

Results

The summary statistics of the measured indices (NDVI, PRI, CAI, MDIN, MSAVI, BSI) derived from the model-based cLHS points and design-based GT points are shown in Table 2. We found that the GT dataset captured a wider range of mean values compared with the cLHS points. For example, NDVI ranged from 0.383 to 0.401 vs. 0.439 to 0.414 for GT and cLHS, respectively. Overall, the sample sizes of test points (1000, 700, 400, and 100) revealed that the descriptive statistics (mean, maximum, minimum, standard deviation) obtained by using design-based GT samples could capture the same statistics of NDVI, PRI, CAI, MDIN, MSAVI, BSI images vis-à-vis to cLHS set (Table 2).

Table 2

Measured 5-index (NDVI, PRI, CAI, MDIN, MSAVI, BSI) statistics from sample points derived by the cLHS method and those collected in the field (GT).

Image Index/Sample Size	Mean		Max		Min		SD
	cLHS	GT	cLHS	GT	cLHS	GT	cLHS	GT
NDVI1000	0.439	0.383	0.827	0.818	0.201	0.170	0.174	0.182
NDVI700	0.428	0.385	0.841	0.818	0.201	0.170	0.172	0.181
NDVI400	0.427	0.383	0.830	0.818	0.206	0.170	0.178	0.185
NDVI100	0.414	0.401	0.803	0.818	0.202	0.171	0.161	0.187
PRI1000	-0.050	-0.026	0.005	0.008	-0.109	-0.108	0.023	0.033
PRI700	-0.050	-0.059	0.091	0.008	-0.108	-0.108	0.024	0.032
PRI400	-0.052	-0.027	0.004	0.008	-0.113	-0.106	0.024	0.033
PRI100	-0.053	-0.029	0.001	0.002	-0.111	-0.101	0.024	0.033
CAI1000	-0.013	-0.009	0.029	0.045	-0.037	-0.034	0.011	0.012
CAI700	-0.012	-0.010	0.025	0.046	-0.037	-0.035	0.012	0.012
CAI400	-0.015	-0.010	0.012	0.033	-0.032	-0.035	0.008	0.012
CAI100	-0.011	-0.011	0.012	0.020	-0.028	-0.029	0.010	0.011
MDIN1000	0.375	0.445	0.968	0.992	0.004	0.008	0.205	0.235
MDIN700	0.377	0.446	0.976	0.992	0.004	0.008	0.214	0.235
MDIN400	0.374	0.437	0.925	0.992	0.019	0.112	0.198	0.221
MDIN100	0.388	0.395	0.894	0.988	0.031	0.024	0.202	0.239
MSAVI1000	0.649	0.617	0.996	0.986	0.004	0.014	0.278	0.244
MSAVI700	0.649	0.617	0.996	0.986	0.004	0.014	0.270	0.244
MSAVI400	0.638	0.610	0.996	0.986	0.016	0.014	0.272	0.240
MSAVI100	0.662	0.628	0.996	0.996	0.040	0.142	0.266	0.230
BSI1000	0.263	0.259	0.365	0.364	0.048	0.030	0.071	0.102
BSI700	0.266	0.260	0.368	0.364	0.015	0.030	0.069	0.101
BSI400	0.267	0.255	0.362	0.364	0.046	0.030	0.073	0.105
BSI100	0.271	0.274	0.363	0.363	0.087	0.030	0.063	0.084

Mean = mean index values; Max = maximum index value; Min = minimum index value; SD = standard deviation of index values.

We compared the spatial structures of the image indices by constructing experimental variograms of cLHS and GT point datasets. Figure 4 showed the model fits for the experimental variograms of 1000 sample size. All model fits using NDVI, PRI, CAI, MDIN, and MSAVI indicated that there was always some nugget variance present. Except for PRI, we observed comparable and relatively similar variogram patterns for cLHS and GT datasets, substantiating that both datasets quantified similar spatial characteristics of the study area. While we detected less spatial autocorrelation with MSAVI, BSI, and MDIN, nonetheless, these comparisons of variography results for cLHS and GT samples indicated that both datasets captured similar spatial structures from the image indices.

Figure 4

Sample experimental variograms of (a) NDVI (b) PRI (c) CAI (d) MDIN (e) MSAVI (f) BSI images for 1000 sample size from cLHS and GT datasets. Variograms use an exponential model.

Among indices, MDIN image displayed the most similarities between the two datasets (Figure 4d). For BSI, the major difference between the two datasets was on the amount of variance and range detected (Figure 4f). All variograms used an exponential model (Table 3), where the semivariogram increases gradually and forms a smooth curve when it reaches the sill.

Table 3

Fitted variogram models of GT and cLHS samples, characterizing the spatial variations inherent in the image indices NDVI, PRI, CAI, MDIN, MSAVI, and BSI.

Image Index/Sample Size	Model	GT				cLHS
		Parameters
		C0	C0+C	C0/C0+C	a(m)	C0	C0+C	C0/C0+C	a(m)
NDVI1000	Exp	0.41	0.82	0.50	425	0.03	0.95	0.03	168
NDVI700	Exp	0.41	0.80	0.51	385	0.38	1.04	0.36	553
NDVI400	Exp	0.42	0.77	0.55	304	0.46	0.99	0.46	574
NDVI100	Exp	0.45	0.70	0.64	511	0.50	0.96	0.52	686
PRI1000	Exp	0.05	0.51	0.09	52	0.18	1.23	0.14	156
PRI700	Exp	0.07	0.50	0.14	63	0.34	1.00	0.34	9139
PRI400	Exp	0.06	0.45	0.13	42	0.57	0.96	0.59	9542
PRI100	Exp	0.07	0.37	0.19	79	0.71	0.86	0.82	9761
CAI1000	Exp	0.62	0.83	0.74	610	0.85	1.19	0.71	639
CAI700	Exp	0.63	0.80	0.78	629	0.87	1.18	0.73	675
CAI400	Exp	0.61	0.78	0.78	581	0.87	1.15	0.75	791
CAI100	Exp	0.61	0.77	0.79	647	0.89	1.10	0.80	685
MDIN1000	Exp	0.81	0.94	0.86	844	0.60	0.98	0.61	185
MDIN700	Exp	0.81	0.93	0.87	816	0.69	1.10	0.62	271
MDIN400	Exp	0.82	0.90	0.91	880	0.67	0.96	0.69	768
MDIN100	Exp	0.83	0.89	0.93	934	0.68	0.94	0.72	871
MSAVI1000	Exp	0.62	0.85	0.72	409	0.79	0.94	0.84	400
MSAVI700	Exp	0.60	0.84	0.72	375	0.77	0.87	0.88	522
MSAVI400	Exp	0.66	0.83	0.80	562	0.80	0.89	0.89	456
MSAVI100	Exp	0.76	0.94	0.81	2977	0.81	1.00	0.81	1023
BSI1000	Exp	0.70	0.99	0.70	387	0.95	1.15	0.82	22
BSI700	Exp	0.63	0.94	0.67	955	0.93	1.12	0.83	153
BSI400	Exp	0.70	0.92	0.76	1111	0.99	1.19	0.83	918
BSI100	Exp	0.71	1.10	0.64	1689	0.93	1.01	0.92	1186

Exp = Exponential Model (γ = C0, h = 0; γ = C0 + C [1- e–], h < 0).

Except for PRI, the GT samples showed comparable nugget effect with the cLHS in the image indices, with NDVI having the least C0 difference (average of 0.12 for all sample sizes) (Table 3). PRI had the largest C0 difference of 0.38 average for all sample sizes. The average sill values (C0+C, which is the plateau value that a variogram stops changing) for cLHS (NDVI = 0.98; CAI = 1.15; MDIN = 0.77; MSAVI = 0.92; BSI = 1.12) and GT samples (NDVI = 0.77; CAI = 0.79; MDIN = 0.91; MSAVI = 0.73; BSI = 0.99) were relatively similar, values ranged between 0.73 and 1.15. The observed similarities between cLHS and GT samples was a validation that both datasets captured similar variabilities of the image indices.

Between indices, NDVI and MDIN had comparable variability in grayscale values as indicated by their sill values (Table 3). NDVI, however, had greater overall variability compared with MDIN and CAI as shown in the nugget to sill ratio. In terms of range values, CAI and NDVI, indicated similarities in terms of variabilities at distances <800 m and of distances at which spatial correlation were present.

We also ran the normal QQ plots for the cLHS and GT samples and most image indices indicated univariate normality of the datasets. We observed a main departure from the 45-degree line for few samples with high values of NDVI (0.6–0.8). Additionally, they correctly assessed the variability in prediction, with values of the standard errors being close to the values of the root mean squared prediction errors (Tables 4 and 5).

Table 4

Prediction error statistics for various sampling sizes applied to the image indices.

Samples	Root Mean Sq.		Ave. Std. Error		Mean Standardize
	CLHS	GT	cLHS	GT	cLHS	GT
NDVI1000	0.173	0.167	0.167	0.165	-0.022	-0.0038
NDVI700	0.175	0.174	0.165	0.173	-0.0009	-0.0006
NDVI400	0.179	0.178	0.174	0.175	0.0030	-0.0065
NDVI100	0.162	0.182	0.166	0.181	0.0405	-0.0038
PRI1000	0.023	1.014	0.023	0.494	0.0001	0.00001
PRI700	0.022	1.003	0.022	0.490	0.0004	0.0021
PRI400	0.024	0.994	0.024	0.496	0.0001	-0.0122
PRI100	0.024	1.060	0.025	0.472	0.0001	0.0411
CAI1000	0.011	0.013	0.011	0.011	-0.0011	-0.0261
CAI700	0.012	0.013	0.011	0.011	-0.0003	-0.0276
CAI400	0.008	0.010	0.008	0.011	-0.0004	-0.0300
CAI100	0.010	0.019	0.010	0.011	0.015	0.0106
MDIN1000	0.205	0.200	0.205	0.255	-0.002	-0.004
MDIN700	0.215	0.209	0.231	0.237	0.006	-0.004
MDIN400	0.196	0.202	0.196	0.241	-0.005	-0.007
MDIN100	0.209	0.204	0.205	0.224	0.003	0.018
MSAVI1000	0.279	0.234	0.282	0.227	0.003	0.015
MSAVI700	0.266	0.235	0.264	0.226	0.017	0.015
MSAVI400	0.267	0.233	0.251	0.226	0.049	0.014
MSAVI100	0.264	0.210	0.267	0.218	0.012	0.037
BSI1000	0.070	0.100	0.075	0.098	-0.025	0.020
BSI700	0.069	0.099	0.070	0.098	0.010	0.026
BSI400	0.073	0.103	0.072	0.102	0.0001	0.040
BSI100	0.063	0.102	0.064	0.102	-0.005	0.039

Table 5

Statistics from simulated image indices using sample points from cLHS and GT datasets.

Simulated	Mean		Max		Min		SD
	cLHS	GT	cLHS	GT	cLHS	GT	cLHS	GT
NDVI1000	0.43	0.38	0.45	0.40	0.40	0.37	0.005	0.006
NDVI700	0.40	0.39	0.41	0.62	0.39	0.26	0.003	0.063
NDVI400	0.39	0.35	0.41	0.73	0.38	0.33	0.003	0.079
NDVI100	0.39	0.38	0.42	0.40	0.39	0.37	0.004	0.005
PRI1000	-0.05	-0.01	-0.04	-0.01	-0.05	-0.02	0.0002	0.0004
PRI700	-0.05	-0.01	-0.05	-0.01	-0.05	-0.01	0.0002	0.0003
PRI400	-0.05	-0.01	-0.05	-0.01	-0.05	-0.01	0.0003	0.0003
PRI100	-0.05	-0.02	-0.05	-0.02	-0.05	-0.03	0.0001	0.0009
CAI1000	-0.01	-0.01	-0.01	-0.01	-0.01	-0.01	0.0002	0.0003
CAI700	-0.01	-0.01	-0.01	-0.01	-0.01	-0.01	0.0003	0.0002
CAI400	-0.01	-0.01	-0.01	-0.01	-0.02	-0.01	0.0001	0.0002
CAI100	-0.01	-0.01	-0.01	-0.009	-0.01	-0.01	0.0002	0.0004
MDIN1000	0.36	0.41	0.40	0.41	0.34	0.40	0.006	0.003
MDIN700	0.37	0.34	0.39	0.35	0.34	0.32	0.007	0.005
MDIN400	0.36	0.39	0.38	0.40	0.34	0.38	0.003	0.005
MDIN100	0.36	0.39	0.37	0.40	0.35	0.37	0.000	0.007
MSAVI1000	0.64	0.64	0.68	0.65	0.62	0.63	0.008	0.004
MSAVI700	0.67	0.64	0.70	0.65	0.66	0.63	0.005	0.005
MSAVI400	0.65	0.63	0.68	0.65	0.63	0.63	0.008	0.004
MSAVI100	0.71	0.65	0.73	0.66	0.70	0.64	0.004	0.005
BSI1000	0.26	0.26	0.27	0.27	0.25	0.25	0.002	0.003
BSI700	0.28	0.30	0.29	0.29	0.28	0.29	0.001	0.001
BSI400	0.29	0.29	0.29	0.29	0.29	0.28	0.000	0.001
BSI100	0.28	0.28	0.29	0.29	0.28	0.27	0.001	0.002

Mean = mean values; Max = maximum value; Min = minimum value; SD = standard deviation of values.

We used Gaussian simulation (e.g., 1000 simulations) to derive simulated NDVI, PRI, CAI, MDIN, MSAVI, and BSI images using the cLHS and GT datasets. We compared the values of the simulated image indices against the measured image indices through the analysis of mean absolute errors (MAE). In Figure 5, the spatial heterogeneity between the image indices followed the same pattern for larger number of samples using the semivariogram parameters obtained from the AVIRIS-NG. The MAEs of cLHS dataset were observed to have a linear correlation with GT dataset (r2 > 0.85) in four out of six indices (NDVI, r2 = 0.88; PRI, r2 = 0.94; CAI, r2 = 0.67; MDIN, r2 = 0.99; MSAVI, r2 = 0.17; and BSI, r2 = 0.94). Figure 5 showed that larger number of samples have higher accuracy rate for both cLHS and GT datasets. The MAE also indicated that the GT dataset has lesser simulation errors than the cLHS when using the NDVI image.

Figure 5

The mean absolute error (MAE) comparing between simulated and measured (a) NDVI (b) PRI (c) CAI (d) MDIN (e) MSAVI and (f) BSI values for the different sampling sizes using the cLHS and GT datasets.

Discussion

Below, we addressed several issues related to the quantification of spatial heterogeneity in order to assess the reliability of the GT dataset that was collected using an adaptive field-or farm-level sampling.

Characterization of spatial heterogeneity from images

This work characterized and quantified the spatial heterogeneity of farmlands with diverse agricultural management and land-use systems from the modeling of the variogram of high spatial resolution image indices (NDVI, CAI, MDIN, MSAVI, BSI, and PRI). From our exploratory analysis, we observed spatial dependence between image indices and point locations for cLHS and GT datasets. We expected this observation as previous studies have used spectral indices to facilitate in the quantification of spatial structures from imagery, such as of soil properties and landscape changes (Oliver and Webster, 2014; Lin et al., 2009). Image indices extracted from the AVIRIS-NG showed how spatial variability, as quantified by the variogram parameters, behaved similarly for the cLHS and GT datasets, with minimal differences in values. The difference in variability, such as the observed lower sill values in the GT dataset and the slightly higher nugget-sill ratio for NDVI, MDIN and CAI, is explained by the influence of mixed farm features during the selection of fields to sample. Such features, when sampled much closer, could influence the variability and the degree of spatial correlation within objects in a farm field. For instance, some samples were collected from locations that had higher moisture content and others were from fields with newly-planted crops. This mixture of field characteristics could have caused the contrasting index values for NDVI, MDIN, and CAI. In the case of BSI, the apparent difference in sill values detected by the cLHS and GT datasets may have been caused by our effort to diversify the GT dataset by collecting them in various farm management practices and crop stages, as per our observation driving around the study area.

In higher sample sizes, indices such as NDVI and MDIN showed similar behaviors in characterizing spatial heterogeneity based from fitted variogram models of the GT and cLHS samples, which provided an understanding that both indices may have captured the same pattern of spatial structures. Maps for the simulated images for the maximum 1000 and minimum 100 sample sizes are shown in Figures 6 and 7, respectively. In addition, NDVI and MDIN exhibited comparable effective ranges, which means that when these indices are used for sampling schemes, they would require the same sampling interval. The larger range for NDVI, CAI, and MDIN, especially for the GT measurements, may have been caused by variability between fields, likely related to differences in tillage practices or field features. Results obtained by Dvorakova et al. (2020) utilizing CAI derived from a hyperspectral image, attributed the larger range to the differences in residue cover during the period when the farms were being prepared for seeding.

Figure 6

Sample simulated image indices based on the maximum 1000 sampling size for GT and cLHS datasets. Along with the images, histogram plots are displayed.

Figure 7

Sample simulated image indices based on the minimum 100 sampling size for GT and cLHS datasets. Along with the images, histogram plots are displayed.

Between CAI and MDIN, the former showed a much stronger nugget-to-sill ratio than the latter for the GT dataset. This could be because CAI demonstrated a strong correlation with plant litter and crop residue that the index distinguished spatially (Nagler et al., 2000; Daughtry, 2001). MDIN, in contrast, is more effective in differentiating crops (Salas et al., 2020) but less in tilled crops (Salas and Subburayalu, 2019). Interestingly, while fitted variogram models of the GT and cLHS samples for CAI, NDVI, and MDIN showed shared patterns, PRI behaved otherwise. Our findings are in agreement with Rahman et al. (2003) which showed the local variance of PRI to be much lower than NDVI.

Number of samples needed to quantify variability

We computed the MAE as the absolute value of the difference between the simulated value and the actual value. For example, when we took the difference between the actual NDVI and the simulated NDVI, this resulted to the MAE of the NDVI images. Smaller MAE values could only mean that the actual and simulated images captured the same spatial variability of NDVI in the study area.

The larger the number of samples, the higher is the accuracy rate for both cLHS and GT data sets. The comparison between the values of MAE from the cLHS and GT data could be used to derive optimum sampling schemes for the study area. For example, in Figure 5, we repeated the process for various sampling sizes and plotted the results. The number of samples that reduced the MAE and reflected the variability of samples for both datasets was around 700. This number had manifested in the six spectral indices that we used to characterize the landscape spatial heterogeneity. Additionally, it also around this number of samples that the MAE curve tended to flatten out or lessen for both datasets. Because the spatial variability could vary based on the surface characteristics and imaging resolution utilized (Silveira et al., 2018), this number should only be used as an example for our study area and not as a guideline for other areas.

Reliability of ground truth data

Variograms have been widely used to characterize the spatial variations inherent in remote sensing images (Woodcock et al., 1988; Wallace et al., 2000; Hamada et al., 2019). Our variogram results implied that the set of GT samples was sufficient and reliable, and was capable of capturing the spatial variance components of field-level sampling using image indices. Results from cLHS and GT datasets rendered fairly similar variogram trends of spatial heterogeneity of the studied agricultural fields. As an example, the higher sill values for GT data correspond to similar higher sill values for cLHS data. In addition, both datasets exhibited similar fashions of decreasing variance as sample size decreases from 1000 to 100 in 88% of the image indices. In Table 2, the cLHS and GT datasets also showed comparability as the measured 6-index statistics analyzed from both datasets were comparable for all tests of sample sizes (1000, 700, 400, and 100).

Key implication

Field spatial variability is an important factor in determining small but significant changes in soil property content (e.g. SOC) as a result of farm management practices. The geostatistical analysis presented above showed that the adaptive field sampling design could be reliably used to estimate the spatial heterogeneity as explained by the hyperspectral image indices. Despite the fact that the GT and cLHS datasets showed similar patterns based on model fits for the experimental variograms, the cLHS yielded larger estimates of spatial heterogeneity. Because the sample size provided appropriate coverage of the study area, the cLHS technique allowed for better discrimination of structural elements in the analyzed agricultural fields. However, based on our familiarity in the field, most of the points derived from cLHS were unrealistic. They were either inaccessible or hard to reach by motor vehicles for sampling due to the complexity of the local road networks in the study area. In practice, this is a significant challenge because a sampling design necessitates both an effective model approach and the precise determination of possible georeferenced sample locations.

While cLHS has demonstrated that it can define the underlying spatial landscape features, we have shown that the adaptive field- or farm-level GT points we collected could produce similar findings. The implication of this comparison between approaches is that, in rural areas where accessibility is unknown and sample locations from model-based procedures such as cLHS may be rendered ineffective, an alternative soil sampling approach such as adaptive field sampling could be practical and, in theory, could provide a higher level of assurance of landscape feature spatial characterization. Also, modifying sample locations depending on ancillary information provided by farmers during fieldwork could also provide an adequate balance and effective sampling site selection (Trobia, 2011). When dealing with complex landscapes and terrains, however, design-based GT sampling could be problematic and may require special attention to generate better approximations of target elements (Graniero and Robinson, 2003).

Conclusion

In the mapping of a target soil property, capturing the spatial variation in the landscape of interest from GT samples using an optimal field sampling design is critical. Our results from image indices indicated that the GT and cLHS datasets could arrive at comparable values for field spatial variability. By using the design-based adaptive sampling in the collection of the GT dataset, specific sampling points were delineated at target areas during fieldwork that contributed to the optimization of sample size and locations and increase certainty in results. Further, this study provides a framework that would assist in the selection of spectral indices and the optimal sample size that efficiently replicate spatial patterns in hyperspectral images.

More research into the relationship between spatial resolution of image indices and spatial heterogeneity will likely be required to better comprehend the discrepancies between the GT and cLHS datasets. More importantly, these investigations may provide useful information and methodologies for future soil sampling and modeling of agricultural landscapes with substantial heterogeneities and complexities.

Declarations

Author contribution statement

Eric Ariel L. Salas: Conceived and designed the experiments; Performed the experiments; Analyzed and interpreted the data; Wrote the paper.

Sakthi Kumaran Subburayalu: Conceived and designed the experiments; Contributed reagents, materials, analysis tools or data.

Brian Slater, Rucha Dave, Parshva Parekh: Performed the experiments; Contributed reagents, materials, analysis tools or data.

Kaiguang Zhao: Contributed reagents, materials, analysis tools or data.

Bimal Bhattacharya: Contributed reagents, materials, analysis tools or data.

Funding statement

This work was supported by the , under award No(s) 80NSSC17K0653 P00001 for the joint NASA and Indian Space Research Organization AVIRIS-NG Campaign in India. It was also supported by / through Central State University Evans-Allen Research Program Fund Number NI201445XXXXG018-0001.

Data availability statement

The authors do not have permission to share data.

Declaration of interests statement

The authors declare no conflict of interest.

Additional information

No additional information is available for this paper.