Mobile Proximal Sensing with Visible and Near Infrared Spectroscopy for Digital Soil Mapping

: The objective of this study was to estimate multiple soil property local regression models, confirm the accuracy of the predicted values using visible near-infrared subsurface diffuse reflectance spectra collected by a mobile proximal soil sensor, and show that digital soil maps predicted by multiple soil property local regression models are able to visualize empirical knowledge of the grower. The parent materials in the experimental fields were light clay, clay loam, and sandy clay loam. The study was conducted in Saitama Prefecture, Japan. To develop local regression models for the 30 chemical and 4 physical properties, a total of 231 samples were collected; to evaluate accuracy of prediction, 65 samples were collected. The local regression models were developed using 2nd derivative pretreatment by the Savitzky–Golay algorithm and partial least squares regression. The local regression models were evaluated using the coefficient of determination ( R 2 ), residual prediction deviation (RPD), range error ratio (RER), and the ratio of prediction error to interquartile range (RPIQ). The R 2 accuracy of the 34 local regression models was 0.81 or higher. In the predicted values for 65 unknown samples, the local regression models could ‘distinguish between high and low’ for 3 of the 34 soil properties, but were ‘not useful’ as absolute quantitative values for the other 31 soil properties. However, it was confirmed that the predicted values followed the transition in measured values, and thus that the developed 34 regression models could be used for generating digital soil maps based on relative quantitative values. The grower changed the ridge direction in the field from east–west to north–south just looking at the digital soil maps. 100 − g


Introduction
Proximal soil sensing (PSS), coupled with GNSS (global navigation satellite system) and visible and near-infrared (VNIR) spectroscopy, is a promising approach for detailed characterization of spatial soil heterogeneity. Given that none of the existing on-the-go soil sensors can provide all the soil information required for precision agriculture (including sustainable farming, environmental load reduction, and risk management), the combination of spectroscopy and chemometrics offers a reasonable alternative for characterizing soil composition [1]. In particular, using PSS in agricultural fields can dramatically enhance this technique by enabling innovative approaches based on appropriate field understanding which characterize local soil and environmental conditions in space and time, improving the efficiency of production to minimize environmental side effects and maximize farm incomes by increasing crop quality and/or yield [2].
Conventional PSS can facilitate the measurement and monitoring of the soil's chemical [3,4], physical [5], and biological [6] attributes to better understand their dynamics and interactions with the environment, while at the same time revealing their broader spatial heterogeneity through digital soil mapping. In addition, when high accuracy prediction values for soil properties related to environmental risks are obtained by PSS, they can be used to effectively monitor methanogen in paddy fields, as well as soil organic carbon, and offer a sound foundation for the adoption of optimal agronomic practices that also enable carbon sequestration and reduction in greenhouse gas emissions. Thus, sensing by PSS can help us to better articulate the potential of agricultural soil to meet the world's needs with regard to soil regeneration, soil contamination monitoring, food security, global climate change, and environmental sustainability, enabling the design and implementation of innovative management practices and the efficient pursuit of sustainable development goals (SDGs). To effectively pursue environmental conservation and sustainable agriculture, a better understanding of soil, based on high spatial-resolution information, is needed. In precision agriculture, also required is simultaneous detection of soil fertility (for crop production) and soil contamination, and environmental monitoring. Thus, the need for timely and accurate multiple soil property information is greater now than ever before.
VNIR spectroscopy is one means of enabling rapid, non-destructive, straightforward, and comparatively inexpensive determination of the relative quantitative value of multiple soil properties. One advantage of VNIR spectroscopy is that a single spectral dataset enables the simultaneous characterization of various soil properties using regression models. Numerous agricultural spectroscopic methods measure the diffuse reflectance spectra of oven-dried soil in the laboratory. In such cases, however, it is necessary to collect soil samples from the target field, ovendry and sieve them (usually in a 2-mm sieve), and finally dispose of the remaining soil, and there has been little advancement in this process [7][8][9][10][11][12][13][14][15][16]. However, these laboratory measurements using VNIR spectroscopy enable highly accurate predictions of multiple soil properties.
In order to improve the conventional measurement methodology, and to determine the relative quantitative values within a given field and between fields, a number of mobile proximal soil sensors have been developed by researchers which are capable of simultaneously recording subsurface diffuse reflectance spectra and global positioning data for digital soil mapping [17][18][19][20][21]. Current mobile soil sensors may have lower accuracy in prediction values than laboratory measurement values, but still provide benefits due to their intense spatial visualization. Therefore, there is ongoing research and development based on in-situ field measurement using on-the-go PSS, in areas such as variable-rate fertilization, the relationship between yield and environmental load, soil fertility, soil pollution, etc. [22][23][24][25]. In most of the previous studies, researchers focused on few soil properties they would like to estimate and created soil maps for variable rate fertilizer using PSS predicted values. Therefore, researchers have found it difficult to provide soil property information required by growers. The grower's requests were that if the spatial variability between and within fields could be understandable by using several soil maps, and if that variability was confirmed by the grower's experience, it could be used for farm management. Furthermore, in actual farm management, crops are cultivated for various purposes (e.g., organic farming, carbon farming, high quality crops, biodiversity, and environmental conservation), making digital soil maps depicting multiple soil properties even more valuable. In our previous studies [26][27][28][29][30], the development of regression models for several soil properties used in soil management was incomplete. Consequently, we desired to provide growers with digital soil maps of all soil properties for soil diagnostics. The aims of this study, then, were to: (i) estimate local regression models for multiple soil properties, using subsurface field-moist soil diffuse reflectance spectra; (ii) evaluate the accuracy of predicted values for unknown samples (as both absolute and relative quantitative values); (iii) demonstrate the effectiveness of newly-developed digital soil map software for our proximal soil sensor, which can generate digital soil maps immediately after measurement in agricultural fields; (iv) describe an instance of decision making by a grower using such digital soil maps.
In this study, a total of 34 local regression models for 30 chemical and 4 physical properties were estimated by 2nd derivative pretreatment and partial least squares regression coupled with full crossvalidation.

Mobile Proximal Sensor
A mobile proximal soil sensor (model SAS3000, Shibuya Seiki Co., Ltd., Ehime Pref., Japan) was used to measure visible and near-infrared subsurface diffuse reflectance spectra (VNIR-SDRS) in agricultural fields. The SAS3000 records several types of data simultaneously, to provide highresolution soil maps of agricultural fields. It is composed of an operational touch monitor, soil penetrator drive switch, display monitor, GNSS antenna (model GA530, Nikon-Trimble Co., Ltd., Tokyo, Japan), DGPS (differential global positioning system) receiver (model SPS351, Nikon-Trimble Co., Ltd., Tokyo, Japan), and control box ( Figure 1). The electrical power for the SAS3000 was supplied by a gasoline inverter generator (EU-16i, Honda Motor Co., Ltd., Minato-ku, Tokyo, Japan). An oil control valve in the tractor is connected to the hydraulic power system of the SAS3000, and is used to drive the chisel unit into the subsurface. The core devices of the SAS3000 for measuring the VNIR-SDRS include a tungsten halogen lamp (JCR15V 150W AL, Fuji Electric Lamp Industrial Co., Ltd., Tokyo, Japan), visible range spectrophotometer (C10083CAH, 320 to 1000 nm, 1 nm interval; Hamamatsu Photonics K. K., Shizuoka Pref., Japan), and near infrared range spectrophotometer (C9406GC, 900 to 1700 nm, 7 nm interval, Hamamatsu Photonics K. K., Shizuoka Pref., Japan). The SAS3000 was designed to collect subsurface VNIR-SDRS at depths of 0.05 to 0.30 m (with 0.05 m intervals) from the soil flattener to the soil surface, created by the double-tire of the SAS3000, which is adjusted by a gauge on the double-tire. In this study, the soil flattener depth was set at 0.15 m. The VNIR-SDRS measurement was recorded at 3 s intervals with a tractor speed of 0.28 m s −1 .

VNIR Spectra, Field-Moist Soil Collection, and Soil Analysis
The experimental sites were in rotational upland fields of taro (Colocasia esculenta (L.) Schott; cultivar: Dodare and Hasuba) in Sayama City (SC), Tokorozawa City (TC), Kawagoe City (KC), and Miyoshi Town (MT), Saitama Pref., Japan. To estimate local regression models for multiple soil properties, the experiment was conducted in a total of 24 fields, of 6 growers [G1, 8 fields (total 0.94 ha); G2, 3 fields (0.45 ha); G3, 2 fields (0.41 ha); G4, 7 fields (0.74 ha); G5, 1 field (0.23 ha); G6, 3 fields (0.35 ha); Figure 2]. A total of 231 VNIR-SDRS data and soil samples were collected at the measured VNIR-SDRS data points after harvesting the taro yam plants in combination with tilling, from 2017 to 2018. The parent materials in the experimental sites were distributed on light clay, clay loam, and sandy clay loam as soil texture (Atterberg method) ( Figure 3).  The SAS3000 sounds an alarm at the time of each data acquisition, counts the number of data items, and displays the VNIR-SDRS data number on each monitor. Each time a predetermined VNIR-SDRS data number (for example, at 20-data-number intervals in a single sensing line) was displayed on each monitor, a placemark stick was inserted into the soil surface beside the chisel unit, as a soil sampling point. After the SAS3000 had fully passed the placemark stick, we dug up the soil surface smoothed over by the soil flattener and packed the field-moist soil samples in sealable plastic bags, taking two sets of field-moist soil samples from each soil sampling point.
To estimate these local regression models for the SAS3000, two sets of 231 field-moist soil samples and their VNIR-SDRS were collected from the total of 24 fields of the six growers, from 2017 to 2018. In addition, in 2019, two sets of 65 field-moist soil samples and their VNIR-SDRS were collected from a total of 10 fields of 4 of the growers (G2, G3, G4, and G5), to confirm the accuracy of predictions for unknown samples using the respective estimated local regression models. Thus, in total, two sets of 296 field-moist soil samples and their VNIR-SDRS were collected, at the VNIR-SDRS data measurement points, after harvesting the taro plants in combination with tilling, from 2017 to 2019. One set was analyzed at our laboratory in Tokyo University of Agriculture and Technology (TUAT; Tokyo, Japan), and the other at the Agricultural Product Chemical Research Laboratory of the Federation of Tokachi Agricultural Cooperative Association (APCRL; Hokkaido, Japan) and the Oita Laboratory of Sumika Chemical Analysis Service, Ltd. (SCAS; Oita Pref., Japan). The 296 fieldmoist soil samples sent to TUAT and APCRL, respectively, were transported in a refrigerator car at less than 10 °C. At TUAT, MC was determined using the oven-dry method (110 °C, 24 h) [31]. The samples were then stored in sealable plastic bags at 5 °C until the end of MC measurement. After MC measurement, each oven-dry soil sample was sieved through a 2 mm sieve, and these 2-mm-sieved oven-dry soil samples were used to measure SOM by the ignition loss method, using a muffle furnace (750 °C, 3 h) [31]. MC and SOM measurements were conducted three times, and the average value was adopted as a dataset for the multivariate statistical analysis. At APCRL, the field-moist soil samples were oven-dried (60 °C, 24 h) and 2-mm-sieved, to measure the following 30 soil physical and chemical properties: pH, y1, PAC, EC, CEC, HR, CN, C-t, N-t, N-a, N-h, N-n, P-a, K, Ca, Mg, Na, B-s, Cu, Mn, Zn, CaMg, MgK, BSP, CSP, ESP, DD, S, SL, and CL. These soil properties were assessed using the soil analysis method of the Agricultural Research Division [32]. The 296 2-mm-sieved ovendry soil samples left by APCRL were sent to SCAS for analysis of SiO and Fe. The SiO solution was extracted using soil nutrient analysis method 14 [33] and the Fe solution was extracted using soil nutrient analysis method 16.3.4 [34]. These sample solutions were assessed with an inductivity coupled plasma optical emission spectrometer SPS5520ICP-OES (Agilent Technologies International Japan, Ltd., Tokyo, Japan).

Absorption Spectral Wavelength of Soil Properties
It has been confirmed that MC, HR, SiO, Fe, and CL have absorption wavelengths at 410 (SiO), 420 to 800 (Fe 3+ [35][36][37][38][39][40][41][42][43][44][45][46][47][48]. As the spectra measurement range of the SAS3000 was 350 to 1700 nm, it is theoretically possible to estimate regression models for all five soil properties (MC, HR, SiO, Fe and CL) using the device. However, the partial least squares regression (PLSR) spectra analysis range was set at 500 to 1600 nm, meaning that absorbance wavelengths of 350 to 500 nm and 1600 to 1700 nm could not be used in this analysis for regression model estimation. In addition, the absorption wavelengths of the other 29 soil properties (29 SPs) did not fall within this 500 to 1600 nm range. Thus, to estimate the soil property regression models for these 29 SPs, it was necessary to confirm the correlation between their measured values and those of the 5 SPs [49][50][51]. Specifically, when there was sufficient (0.4 or greater) measured value correlation between the above 5 SPs and the remaining 29 SPs, we were able to estimate regression models.

Spectra Data Processing and Partial Least Squares Regression Setting
To estimate the regression models, two spectrometer (visible and NIR range spectrophotometer) data sets were converted to a VINR-SDRS data sets, and an absorbance spectra data sets, using the following procedure, based on Kodaira and Shibusawa [27] and using Data Monitor Software (DMS; Shibuya Seiki Co., Ltd.). To transform the two spectrometer data sets to VNIR-SDRS data sets, the former was connected at 1000 nm, coupled with equal (5-nm-interval) spacing of the wave number resolution, with a wavelength range from 350 to 1700 nm, using the interpolation method. Then, the VNIR-SDRS data sets were converted to absorbance spectra data sets, using white reference spectra data (with a standard reflector: SRS-40-020; Labsphere, Inc., North Sutton, NH, USA) and dark reference spectra data (with light shielding). The white and dark reference spectra data were recorded once in each experimental field before commencing the VNIR-SDRS measurement [27].
Since absorbance spectra data are affected by uncertain noise (caused by factors such as the soil particle size, moisture variability, and light scatter), pretreatment of the spectra data can often improve the regression model accuracy, compared to using untreated absorbance spectra data. In addition, to enhance weak signals and reduce the noise and light scatter that would affect the baseline, smoothing and derivative mathematical pretreatment can be useful, depending on the parameter conditions. For example, the first derivative can eliminate baseline fluctuation, and the second can eliminate multiplicative fluctuation. In addition, the second derivative offers the following beneficial effects: (i) the peak waveforms can be separated from the broad absorption wavelengths; (ii) the peak wavelength values are enhanced more than with the first derivative; (iii) the peak wavelength value corresponds to the original absorbance wavelength (the peak waveform is inverted).
Thus, in this study, the second derivative coupled with smoothing (2DS-G) was applied [52]. MLR (multiple regression analysis), PCA (principal component analysis), PCR (principal component regression), NN (neural network) are among the multivariate analysis methods. In this paper, PLSR analysis [53][54][55] coupled with full (leave-one-out) cross-validation was applied, because: (i) even when there are fewer response variables than explanatory variables, it is possible to estimate a regression model; (ii) it enables orthogonalization between wavelength variables, to avoid multicollinearity.
In the field measurement using the SAS3000, it is impossible to remove all the gravel and crop residues, roots, creatures, etc. from the soil. Therefore, there are some outliers in absorbance spectra data sets, and it is necessary to exclude outliers from the absorbance spectra data sets. In this study, the selection criteria for identifying a given sample as an outlier was the largest residual variance value, based on the residual sample variance plot in U98 after the PLSR [56,57]. This process of absorbance spectra data outlier exclusion was repeated until the coefficient of determination (R 2 ) of the regression model was 0.82 or higher. In this study, we decided to exclude outliers up to a maximum of 115 samples (roughly half of the 231 total samples). In this study, the following 2DS-G and PLSR parameters were set using U98: 2DS-G setting was, (i) the wavelength range used for calculation was 350 to 1700 nm; (ii) the 'polynomial order' was second order; (iii) the 'number of smoothing points' was 3 to 41 (only odd numbers could be selected), and the guideline for selection was the 'number of smoothing points' with the highest R 2 .
PLSR setting was, (i) the wavelength range used for calculation was 500 to 1600 nm, with reference to the report by Kodaira and Shibusawa [27]; (ii) the 'number of PCs' (principal components; i.e., the number of PLSR factors) was 20; limited to this number, for each regression model, to prevent overfitting; (iii) the remaining settings were default settings.

Evaluation of Regression Models and Predicted Values
The regression models and predicted values for the unknown samples were evaluated using the following performance indicators: the R 2 value, the root mean square error (RMSE), residual prediction deviation (RPD), range error ratio (RER), and the ratio of prediction error to interquartile range (RPIQ). Generally, a model that performs well would have large values of R 2 , RPD, RER, and RPIQ, and a small RMSE value. These performance indicators were computed for each soil property according to the following equations: where is the measured value of the ith sample measured by conventional soil analysis, is the predicted value of the ith sample, is the average of the predicted values, N is the number of samples, and S.D. is the standard deviation of the observed soil property in the full cross-validation or prediction dataset; where the RER is the ratio of the range of the measured values of soil properties in the full crossvalidation set or the prediction set to the RMSE; where IQ is the interquartile range of the observed soil property in the full cross-validation or prediction dataset (IQ = Q3 − Q1, quantifying the spread in the central 50% of the data (Q2), with Q1 as the median of the first half, and Q3 the median of the second half). RPD and RPIQ (ratio of performance to IQ) are used in the spectroscopic literature to make data more comparable, accounting for the differences in the spread of the data. Their explanatory power in terms of classifying the model performance is the same as R 2 , but they are also useful for comparing datasets with different data spreads; with the additional advantage that RPIQ can be used for skewed data [57,58]. Previous studies defined criteria to assess the relative fit of the regression model using performance indicators, as following: R 2 > 0.90, 0.82 < R 2 ≤ 0.90, 0.66 < R 2 ≤ 0.82, 0.50 < R 2 ≤ 0.66, and R 2 ≤ 0.5; RPD > 3.0, 2.5 < RPD ≤ 3.0, 2.0 < RPD ≤ 2.5, 1.5 < RPD2.0, and RPD ≤ 1.5; RER > 20, 15 < RER ≤ 20; 10 < RER ≤ 15, 8 < RER ≤ 10, and RER ≤ 8 [59,60]. RPIQ follows the same classification format as RPD. In this study, based on these indicators, the regression model performance was classified into five categories: Category A is excellent, B is good (successful for RER), C indicates approximate quantitative prediction (moderately successful for RER), D indicates capability to distinguish between high and low (moderately useful for RER), and E means not useful (screening for RER).

Digital Soil Mapping Software
Traditionally, ArcGIS Ver.10.0 software (ESRI Japan Co., Chiyoda-ku, Tokyo, Japan) has been used to display digital soil maps, but parameter setting of ArcGIS is complicated and is both time consuming and very expensive [27]. Thus, in this study, the digital soil maps were prepared using Soil Map Viewer (SMV, Figure 4) for the SAS3000. The generated digital soil maps are displayed with colored dots representing the soil sensing points, meaning that no interpolation method (e.g., Inverse Distance Weighting, Kriging) is required. The features of SMV are: (i) classification of predicted values can be specified up to 5 levels; (ii) the boundary line of each field is displayed (boundary position data is required); (iii) statistics data of the mean, max, min, and coefficient of variation can be confirmed; (iv) the absorbance spectra, the absorbance spectra after pretreatment and predicted value at each measurement point can be confirmed; (v) a histogram of the selected area is displayed.
As a result, it is possible to display the 34 soil properties soil map and talk with growers on site immediately after the measurement, and to provide a sound analytical foundation for soil management. In particular, (iv) can support outlier determination as unusual spectra. SMV was jointly developed for the SAS3000 by Shibuya Seiki Co., Ltd. and the authors.    Table 1 shows the correlation between measured values. As the absorption wavelengths of CL, Fe, HR, MC, and SiO (5 SPs, red character) have been confirmed to fall within the range of 500 to 1600 nm, it is theoretically possible to estimate local regression models for these properties using PLSR. The measured values for BSP, Ca, CaMg, CEC, CSP, C-t, Cu, DD, MgK, N-h, N-t, P-a, PAC, pH, S, SL, SOM and y1 (18 SPs, green character) were confirmed to be in correlation with those of the 5 SPs; while those of B-s, CN, EC, ESP, K, Mg, Na, N-a, N-n and Zn (10 SPs, blue character) were confirmed to be in correlation with those of these 18 SPs. Consequently, it was confirmed that the 10 SPs had an indirect correlation with the 5 SPs, through their correlation with the 18 SPs. There was no confirmed correlation with measured values in the case of Mn. Therefore, with the exception of Mn, it was confirmed that local regression model estimation was possible for both the 10 SPs and the 18 SPs. For Mn, it is statistical estimate.

Evaluation and Accuracy of Local Regression Models
The accuracy of the respective local regression models estimated using PLSR is shown in Table  2. Without outlier exclusion, the smoothing points (S.P.) of 2DS-G for the absorbance spectra with the highest R 2 were 5 to 33. The average absorbance spectra of the 231 samples, after pretreatment by 2DS-G, are shown in Figure 6. The S.P. of B-s, Fe, MC, N-t, and P-a were 30 or more, meaning that greater noise removal was required than in the case of the other soil properties. On the other hand, the S.P. of CSP and EC were 10 or less, suggesting that it was necessary to utilize small spectral fluctuation. When outliers were excluded until R 2 exceeded 0.82, the number of PLSR factors (F) for the local regression models was 2 to 14, except in the case of Cu ( Table 2). The accuracy of the Cu local regression model was R 2 = 0.81, and F was 15 when the maximum outlier was 115. Although F for Bs was 2, the number of correlations between measured values was 3, so the result was inferred to be 'not underfitting' ( Table 1). As F for y1 and DD was 14, there was concern about overfitting; however, the result was inferred to be 'not overfitting', because the F selected in the previous study was 19 [61]. The local regression model R 2 accuracy of full cross-validation was classified as Category C for Cu, Category D for N-n, and Category B for the other local regression models. The accuracy of R 2 , RPD, and RER for full cross-validation was 0.65 to 0.84 for R 2 , 1.70 to 2.50 for RPD, and 7.15 to 14.6 for RER. In the full cross-validation RPD classification, the local regression models for BSP, CL, CEC, Cu, CSP, DD, EC, MgK, N-a, N-h, N-n, pH, SiO, and y1 were classified as Category D. The full cross-validation RER of CaMg, DD, EC, pH, N-a, and y1 was classified as Category E. Therefore, in order to estimate these local regression models, it was confirmed that the range of measured values should be expanded. Scatter plots of the three soil properties with R 2 > 0.82 (B-s, C-t, HR) and the four soil properties with low R 2 (BSP, Mn, N-n, SiO) are shown in Figure 7 as examples of the scatter plots of measured values versus predicted full cross-validation values. The three soil properties with high R 2 accuracy were distributed on the 1:1 line, but the four soil properties with low R 2 accuracy were distributed away from the 1:1 line.
PAC showed response variable correlation with CL, Cu, DD, MC, P-a, S, SiO, and SL (Table 1). Furthermore, PAC showed correlation between local regression model coefficients for BSP, CL, CN, Fe, K, MC, Mg, pH, and SiO (Table 1). PAC's response variable and local regression coefficient correlation showed a positive correlation with MC and SiO, and a negative correlation with CL. In particular, the correlation with 5 SPs is very important. Given these results, even if the PAC has no absorption wavelength, it is assumed that the PAC local regression model is reliable. The same applies to the other local regression models in Table1. In contrast, Mn showed no correlation with the measured values of any soil property (Table 1). However, in the local regression model coefficients, a negative correlation was confirmed with Ca, CaMg, CN, HR, P-a, and S; a positive correlation was confirmed with Mg, K, and Zn (Table 1). Thus, only nine coefficient correlations were confirmed, and the Mn local regression model was assumed to be of low reliability.

Evaluation and Accuracy of Predicted Values for Unknown Samples
The accuracy of the predicted values for the 65 unknown samples, using the 34 local regression models, is shown in Table 3. In the case of R 2 , RPD, and RPIQ, the local regression models for 10 soil properties (CaMg, CL, C-t, HR, MC, N-t, PAC, S, SOM and SiO) were classified as Category C or D, with the other 24 local regression models classified as Category E. The RER of HR, N-t, and PAC was Category D, while the others were Category E. The RER of HR, N-t, PAC was Category D; others were Category E. Therefore, in the overall dataset for this experiment, the predicted values for HR, N-t, and PAC, in the 65 unknown samples, using the local regression models, were considered to 'distinguish between high and low', as absolute quantitative values for each sample. The predicted values for properties other than HR, N-t, and PAC were evaluated to be 'not useful' as absolute quantitative values for each sample. To investigate why the accuracy of the predictions decreased, the transition in the measured and predicted values for the 34 soil properties in the 65 unknown samples is shown in Figure 8, where the two blue lines are the upper and lower predictable limits (predictable range), the black line is the measured value, and the red line is the predicted value. Large error is confirmed in the case of 22 soil properties (B-s, BSP, Ca, CEC, CN, CSP, Cu, EC, ESP, Fe, MC, Mg, K, MgK, Na, N-n, P-a, PAC, pH, SL, y1, and Zn), where the unknown samples fell outside the predictable range of the respective local regression models. This result suggests that each of these local regression models must be re-calibrated by PLSR, using a new dataset that includes the unknown samples outside the predictable range. The mean of the measured values (M.M.) and predicted values (M.P.), and error percentage (E.P.), are shown in Table 3. A total of 15 properties (Bs, BSP, CEC, Cu, EC, ESP, MgK, Na, N-a, N-h, N-n, P-a, SiO, y1, and Zn) showed an error of 10% or more, suggesting that bias adjustment is required. On the other hand, it was visually confirmed that the predicted values followed the transition in the measured values. Therefore, it was determined that these predicted values estimated by the local regression models could be used for soil maps displaying relative quantitative values. In particular, in this study, 13 soil properties (CL, CN, C-t, DD, Fe, HR, MC, N-t, PAC, S, SL, SiO, and SOM) could be used as relative quantitative values.

In-Situ Digital Soil Mapping
In precision agriculture, it is necessary to understand site-specific variability in agricultural fields, and growers require in-situ digital soil mapping of relative quantitative values. After soil sensing by the SAS3000, several soil property maps were provided at the field, using SMV, to a taro grower, who then focused on the soil variability of CL, N-t, and P-a in his field ( Figure 9). The soil variation among CL, N-t, and P-a in the field was divided into eastern and western regions. At that time, the grower noted that the variability was consistent with top soil dressing because, in the past, fresh red soil was brought in only in the eastern region. Additionally, the growth, quality, and yield of taro also differed between the east and west regions. In the grower's experience, the western region contained a large amount of gravel and small taros. To improve the variability in the yield, and the quality and efficiency of the farm work, one scenario for precision agricultural soil management in this field might involve site-specific fertilization through topdressing and basal dressing, top-soil dressing of the western region, changing the ridge direction the field, etc. In fact, the grower changed the ridge direction in the field from east-west to north-south, resulting in greater farm work efficiency (for example, top dressing, pest control, and harvesting), with no soil improvement costs. Note that, the area without dots in the southeastern part of the field does not indicate a breakdown of the SAS3000, but rather a region in which taro plants were stored underground, and which was not measured.

Regression Model Accuracy and Actual Operation
With the range of prediction for each soil property, the local regression models is shown as 'Range' in Table 2. In the soil diagnosis of upland fields in Japan, the reference value range is 10 to 30 mg 100 g −1 for P-a (according to the Agricultural Policy Division, [62] standard). In this case, the local regression model enables effective soil diagnosis. However, since the reference value range is 2.0 to 40 ppm for Zn, the Zn local regression model is not usable, which requires a new local regression model for Zn, by adding a dataset with multiple Zn measured data of 13 ppm and more (e.g., 14, 15..., 40 ppm or more). The same logic applies in the case of the other local regression models.

Performance Indicators and Regression Model
Parameters R 2 , RPD, RER, and RPIQ are performance indicators for evaluating the absolute quantitative accuracy of the predicted values, and it is important to note that the criteria for defining 'excellent' and 'not useful' models using such performance indicators are rather arbitrary, and there is no statistical basis for how these classification thresholds are determined [63]. Shree et al. [63] suggest that it is more important for a user to evaluate a model's performance with respect to the specific objectives of the given study, and we concur.

Digital Soil Map for Growers
In this study, CL, N-t, and P-a soil maps were used for grower decision making. The properties of soil maps required vary depending on the crops and varieties, the agricultural equipment available, and the philosophy and cultivation experience of the grower. The six growers suggested that it was important to understand the relative variability of multiple soil properties within a given field and between fields. As a specification from growers, mobile proximal soil sensor has been required that can provide various soil properties maps. Therefore, it is highly important, for sustainable farm management, that a mobile soil sensing system has the ability to provide multiple soil property maps.

Potential of a Mobile Proximal Soil Sensor
The results of this study, using the SAS3000, are particularly promising in suggesting that the method can be used for in-situ digital soil mapping, and rapidly, cheaply, and with less labor compared to conventional laboratory measurement, and the grower's use of the relative quantitative soil maps to modify his crop management suggests the possibility of using the SAS3000 as information technology for precision agriculture. As potential applications of the proposed method, the local regression models, and SMV for digital soil mapping can be used not only for variable-rate fertilization, crop selection, land use, and site-specific soil management, but also for effective monitoring of soil organic carbon (e.g., SOM and C-t), and to provide a sound analytical foundation for agronomic practices that enable carbon sequestration and reduction in greenhouse gas emissions [10,[64][65][66].

Conclusions
The dataset used in this study comprised soil analysis values and VINR soil reflectance spectra, coupled with GNSS data, from the rotational upland fields of six growers in Saitama Prefecture, Japan. The local regression models for 30 chemical and 4 physical soil properties, which showed an R 2 accuracy greater than 0.81, were estimated using 2nd derivative pretreatment and PLSR. The evaluation of full cross-validation, based on the R 2 and RPD performance indicators, revealed no 'excellent' or 'not useful' local regression model estimations, while six soil properties were classified as 'screening' by the RER indicator. In the performance indicator-based evaluation of the predicted values of 34 soil properties for 65 unknown samples, estimated using the local regression models, 3 of the 34 models could 'distinguish between high and low', but the other 31 local regression models were 'not useful', as absolute quantitative values. The study demonstrated that local regression model-based predictions for 13 of the 34 soil properties in the 65 unknown samples could provide relative quantitative values for the properties, based on the comparative transition in measured values and predicted values. The prediction error percentage (10% or less) for the 65 unknown samples confirmed that 20 soil properties could be successfully estimated as relative quantitative values, suggesting that the proposed method could generate useful digital soil maps displaying relative quantitative soil property values. Therefore, we concluded that the SAS3000 digital soil mapping equipment provided relative quantitative accuracy.
Author Contributions: M.K. and S.S. designed the study and performed the experiments; M.K. analyzed the results and wrote the article; S.S. contributed to the manuscript editing process. All authors have read and agreed to the published version of the manuscript.
Funding: This study was supported by grants from the Project of the NARO Bio-oriented Technology Research Advancement Institution (project ID16822280 of the special scheme project on regional developing strategy), and partially funded by the Ministry of Agriculture, Forestry and Fisheries, Japan.