Soil parent material context

Nearly the entire extent of Saskatchewan's soils is formed on unconsolidated transported parent materials, predominantly the result of glacial processes (Anderson and Cerkowniak 2010). The most common soil parent materials in Saskatchewan are glacial tills (often enriched with limestone-rich bedrock), followed by glaciolacustrine and fluvial deposits (Table 1). Due to the overall slope of the land and drainage to the northeast, large glacial lakes formed that caused extensive glaciolacustrine deposits associated with proglacial lakes (Anderson and Cerkowniak 2010). Additionally, high-energy streams from the melting glaciers were responsible for the creation of fluvial deposits along deltas that formed where streams entered glacial lakes. Minor areas of eolian landscapes are also present, which formed due to post-glacial wind erosion. Minor areas of bedrock, eolian, peat, and colluvium are also present. The existing parent material map, from Saskatchewan's detailed soil survey polygons, is provided in Fig. 1.

Fig. 1.

Overview maps of imagery and parent material in Saskatchewan. (a) True colour red-green-blue median surface reflectance image (Landsat 7 Bands 3–2-1) for Saskatchewan from July and August from 1999 to 2021. The black points indicate the location of the National Pedon Database sample locations used for validation. The squares with yellow diagonal lines indicate the locations of the finer spatial resolution example maps. (b) Existing Saskatchewan soil parent material map. Coordinates are in UTM Zone 13 N NAD83.

Table 1.

Distribution of parent material types in the detailed soil survey polygons of Saskatchewan (SKSIS Working Group 2018) and the National Pedon Database parent material class counts (Agriculture and Agri-Food Canada 2016).

Remote sensing data

Remote sensing data for the agricultural regions of Saskatchewan (Fig. 1) were acquired using Google Earth Engine (Gorelick et al. 2017) to be used as environmental covariates in the parent material disaggregation models. Landsat 7 Tier 1 surface reflectance data were acquired from 1999 to 2021, with several band indices calculated (Table 2). The median anthocyanin reflectance index (ARI), which is associated with anthocyanin pigments in plant foliage (Gitelson et al. 2001), was calculated for data from the months of July and August. This index was included as it has been useful for discriminating vegetation differences in Alberta for peatland mapping (DeLancey et al. 2019). The median canopy response salinity index (CRSI) was also calculated using data from the months of July and August, as it has been used to distinguish soil differences associated with salinity (Scudiero et al. 2015). The median normalized difference vegetation index (NDVI) was calculated for the months of July and October, separately. These months were selected as July corresponds to peak photosynthetic activity, and only non-arable cropland would have any photosynthetic activity in October. Non-arable cropland would be expected to be preferentially located on sandier parent materials or those with steeper slopes. The standard deviation of NDVI using data from May to October was also determined. Level 1 C data from Sentinel 2 were acquired for the months of May to October from 2015 to 2020 to generate the median red-edge inflection index (REIP). The REIP was included in the study, as it has been an important parameter for discriminating vegetation differences in Alberta (DeLancey et al. 2019). Each parameter was median focal filtered using a 3 × 3 window size and exported from Google Earth Engine at 30 m spatial resolution. The Google Earth Engine scripts for generating these covariates are available on GitHub (Sorenson 2021a).

Table 2.

Complete list of features included in the analysis prior to feature selection.

Data from the Sentinel 1 platform were also retrieved from Google Earth Engine. Processing of the Sentinel 1 data was done using the methodology in Hird et al. (2017). The median vertical–vertical (VV) and vertical–horizontal (VH) backscatter was determined using data from May to October from 2015 to 2020. To focus more on broader landscape trends, rather than fine-scale variation, the data were median focal filtered with a 9 × 9 window size and exported at a 30 m spatial resolution, rather than the default 10 m, to match the resolution of the Landsat 7 data. The Google Earth Engine script for generating these covariates is available on GitHub (Sorenson 2021a).

Bare soil composite imagery

BSCI was generated using Google Earth Engine and the methodology from Sorenson et al. (2021). Atmospherically corrected Tier 1 surface reflectance data for Landsat 5 were retrieved for the months of July and August from 1984 to 1995. Landsat 5 data from this time period, rather than more recent Sentinel-2 or Landsat 7 and 8 data, were used because the total number of acres in Saskatchewan under conventional tillage peaked in 1970 and declined to only 18% of farmland by 2006 (Statistics Canada 2015). Therefore, more recent imagery is heavily affected by crop residues. The months were limited to July and August to ensure that crop residues from previous years had time to decompose, thereby reducing the influence of crop residues on BSCI. Clouds, shadows, and low-quality pixels were filtered using the Landsat quality assessment band (pixel_qa). The Landsat 5 data had seven bands: Band 1 corresponding to blue light (450–520 nm), Band 2 corresponding to green light (520–600 nm), Band 3 corresponding to red light (620–690 nm), Band 5 (1550–1750 nm), Band 6 (10400–12500 nm), and Band 7 (2080–2350 nm). Band 6 was not used in the analysis.

Image collections were filtered using NDVI, normalized burn ratio (NBR₂), normalized difference water index (NDWI) (Demattê et al. 2018), and normalized difference index 7 (NDI7) (Sorenson et al. 2021). Pixels were kept as likely representing the bare soil surface if the NDVI value was less than 0.3 (Demattê et al. 2020), NBR₂ was less than 0.1 (Demattê et al. 2020), NDWI was less than 0.5 (Du et al. 2016), and NDI7 was less than 0 (Kokaly et al. 2017). Additionally, any region mapped as pasture or grassland in Agriculture and Agri-Food Canada's 2019 Annual Crop Inventory was masked (Agriculture and Agri-Food Canada 2019). While 2019 does not match the years of the bare soil imagery collections, both land uses tend to be managed for long term. Following filtering and masking, the median reflectance values for each pixel were calculated, and then spatially filtered using a circular 10 × 10 median focal filter. The focal filtering was done to account for the 300 m location uncertainty reported in the National Pedon Database and to smooth out fine-scale variation more likely attributable to factors other than parent material. The final filtered composite image was then aggregated to a spatial resolution of 30 m. Additional details and band ratio equations are available in Sorenson et al. (2021) and the Google Earth Engine script is available on GitHub (Sorenson 2021b). Example results for the BSCI are provided in Fig. 2, which includes true colour red-green-blue images of the bare soil surface.

Fig. 2.

Imagery for the four finer spatial resolution example areas. The first row corresponds to the northeast example area, the second row is the northwest example area, the third row is the southwest example area, and the fourth row is the southeast example area. The figures in column (a) are true colour red-green-blue median surface reflectance image (Landsat 7 Bands 3-2-1) for Saskatchewan from July and August from 1999 to 2021. The figures in column (b) are true colour red-green-blue median reflectance bare soil pixel composite images (Landsat 5 Bands 3-2-1) for Saskatchewan from 1984 to 1995. This represents a true colour image of the exposed soil surface. The white areas in column (b) correspond to areas either without bare soil pixels present or permanent pasture areas that have been masked. Coordinates are in UTM Zone 13 N NAD83.

Terrain attributes

Terrain attributes were determined using the 30 m digital surface model from the Advanced Land Observation Satellite (ALOS) (Japan Aerospace Exploration Agency 2015) retrieved using Google Earth Engine. Two terrain attributes, terrain ruggedness index (TRI) and the standard deviation of elevation, were calculated using a range of median focal filtering window sizes and calculation window sizes (Table 2). Median focal filtering was used, rather than mean, as it is less sensitive to the effects of outliers. These attributes were selected, as they represent coarser scale patterns in landscape variability rather than specific fine-scale patterns related to factors such as local landscape morphometry, hydrological characteristics, and landscape exposure. While these finer-scale attributes were documented to improve predictive parent material mapping at 100 m in British Columbia (Heung et al. 2014; Bulmer et al. 2016), the characteristics of prairie landscapes make these attributes challenging to accurately calculate with currently available DEMs in Saskatchewan. Satellite-derived DEMs with spatial resolutions of 30 m, like the ALOS model used here, lack the horizontal and vertical accuracy to capture fine-scale variation in the relatively low-relief Canadian Prairies. Specifically, the ALOS 30 m DEM has a reported vertical root mean square error of 1.78 m (Caglar et al. 2018), which is greater than the relief differences driving soil property differences in some hummocky landscapes in Saskatchewan (Landi et al. 2004). Additionally, 30 m spatial resolution is not fine enough to consistently characterize slope and topographic variation in many Saskatchewan landscapes (Landi et al. 2004). However, currently these are the best freely available DEM data for Saskatchewan. Therefore, this study focused on using terrain attributes that accounted for general degrees of variability more broadly.

TRI was calculated using the system for automated geoscientific analyses (SAGA) (Conrad et al. 2015). The TRI was calculated with either a window size of 10 or 20 (Table 2), and with the following median focal filters applied to the DEM prior to processing: no filter, 3 × 3, 5 × 5, and 9 × 9 (Table 2). Applying a 3 × 3 median focal filter to the TRI results with no filter on the DEM was also tested. The standard deviation of elevation was calculated using Google Earth Engine with window sizes of 3 × 3, 5 × 5, 9 × 9, and 21 × 21 following 3 × 3 median focal filtering of the DEM (Table 2). The standard deviation of elevation was also calculated using window sizes of 21 × 21 and 101 × 101 for a 9 × 9 median focal-filtered DEM. Examples of the terrain features determined as important environmental covariates during model development are illustrated in Fig. 3. TRI was not calculated using a 101 × 101 window because of the computational time required compared to standard deviation. Multiple scales were considered as multi-scale digital terrain analysis has been shown to improve predictive soil mapping results (Behrens et al. 2010).

Fig. 3.

Terrain attribute imagery for the three terrain parameters selected as important predictors by the feature selection process. The first row provides imagery for the portion of Saskatchewan with conventional soil survey maps at 1:100000. The bottom row is the imagery for the southwest example area. The figures in column (a) are the ALOS DEM with a 3 × 3 median filter applied and the standard deviation of elevation calculated with a 21 × 21 window. Figs. in column (b) are the ALOS DEM with a 9 × 9 median filter applied and the standard deviation of elevation calculated with a 21 × 21 window. The figures in column (c) are the ALOS DEM with a 9 × 9 median filter applied and the standard deviation of elevation calculated with a 101 × 101 window. Coordinates are in UTM Zone 13 N NAD83.

Model development

Predictive models were developed using the detailed soil survey polygons as training data, which were mapped at a 1:100000 scale. The disaggregation approach was similar to Møller et al. (2019) where a single set of area-proportionally sampled synthetic training points were generated and a random forest model was trained. This approach was selected for the computational efficiency (Møller et al. 2019). Synthetic training data were developed by randomly assigning points within each polygon on an area proportional basis. The number of points generated for each polygon was equal to the natural logarithm of the polygon area divided by the area of the smallest polygon in the data set. This number was then multiplied by 5 for the final number of points per polygon. The rationale for this approach was that it ensured multiple training points are generated for each polygon and that larger polygons have more points generated to reduce the likelihood that they are underrepresented in the training data.

Points were assigned to the parent material class based on the first deposition layer described for the polygon that they were generated in, with the exception of the lacustrine over till class, which reflects both the first and second deposition layers. Polygons with parent material values of undifferentiated indicated that there is too much variability in parent materials at the mapping scale of the detailed soil survey polygons to assign a parent material class, so data from these polygons were not used for model development. As each polygon has a single component, the single component was used to assign the training data label. Multi-component polygons in the parent material maps were labelled as undifferentiated and they were excluded during training data creation. In total, 5786772 synthetic training points were generated. The scripts for generating the training data from the parent material polygons are available on GitHub (Sorenson 2021c).

The parent material classes for the synthetic training points generated from the detailed soil survey polygons were simplified. Fluvial eolian and fluvial lacustrine parent materials were grouped as fluvial because of the similar characteristics of these materials. Additionally, minor parent materials with very similar characteristics were removed, specifically glaciolacustrine and glaciofluvial. Fluvial eolian and fluvial lacustrine were grouped as fluvial. Fen and sphagnum peat were grouped as a single peat class. Lacustro-till were relabelled as lacustrine due to its infrequent occurrence and high similarity to lacustrine deposits. Residual and undifferentiated bedrock were combined into a bedrock class. The final parent material classes were: bedrock, colluvium, eolian, fluvial, lacustrine, lacustrine over till, peat, and till.

To improve the computational efficiency of the model building, the total points were randomly subset to three sets with 1000000 training points from the original 5786772. Separate training subsamples were generated for the three following models: (1) the BSCI-only model, which includes only training points from locations where BSCI is available and includes BSCI as environmental covariates, (2) the exhaustive model, which includes training points from all locations, regardless if BSCI is available and includes BSCI as environmental covariates, and (3) the BSCI-excluded model, which also includes training points from all locations but does not include BSCI as environmental covariates. Where BSCI imagery was not available, a null value of 0 was set for model development purposes.

Following creation of the training data, random forest models were built using the ranger package in R (Wright and Ziegler 2017). As class imbalances exist in the training data, the case weights term in the ranger model was used to balance across classes during model building. The case weights term in the ranger package allows for probabilities to be assigned to each training point that determines the likelihood a point will be sampled during individual tree building in the random forest. As a result, each tree of the random forest can be built with balanced classes, and the overall random forest model will therefore be built on class balanced data even if the original data set has imbalanced classes. Case weights were set so that each parent material class had an equal chance of being subsampled for each individual tree built in the random forest. An initial model was built using all environmental covariates listed in Table 2 to determine feature importance based on Gini index. Features were then selected for the final models using a backward feature selection process, where random forest models were built using sequentially less features. For each iteration, the feature importances were recalculated, and the least important feature was removed. The final features selected were determined based on where the out-of-bag classification error was minimized. In total, the BSCI-only model used 16 features, the exhaustive model used 20 features, and the BSCI-excluded model used 11 features (Table 3). Final random forest models were built for each scenario with the following model hyperparameters: probability set to true, extra trees as the split rule, and case weights set so that subsampling probability was equalized across parent material classes to reduce the influence of class imbalances on the model. The scripts for generating the random forest models and resulting parent material maps are available on GitHub (Sorenson 2021c).

Table 3.

Final features used for each of the three predictive models along with feature importance values for the final model.

Validation data

Data from the National Pedon Database (Table 1) were used as model validation data (Agriculture and Agri-Food Canada 2016). Historically, the samples were collected to represent modal soils from dominant soil types from a region and were collected by horizon with analyses performed on the homogenized horizon samples. The location accuracy for these samples is low, with specific locations estimated relative to land marks by the soil surveyors. The estimated location uncertainty is approximately ±300 m. The 300 m uncertainty is the location uncertainty reported for the data in the National Pedon Database. These soil pedon observations were used to inform the soil survey parent material mapping and so they are not completely independent from the maps used as training data, but the point observations themselves were not used in the model training. The soil parent material classes were simplified for this study based on the original parent material classes in the survey database using the same approach as the training data.

The distribution of parent material classes amongst pedon validation data did not match the distribution of parent material classes reflected by the detailed soil survey polygons (Table 1). This was likely because the pedon data were collected to provide example profiles for the range of soil associations present in Saskatchewan, not to sample parent material classes on a proportional basis. The most frequent parent material in the pedon validation data was fluvial, followed by lacustrine, till, and then lacustrine over till (Table 1). Minor amounts of eolian, bedrock, and colluvium parent materials were present. No peat data are present in the validation data. For the detailed soil survey polygon data, the most common parent material class was till, followed by fluvial and then lacustrine (Table 1).

Model validation

Model performance was evaluated based on overall model accuracy, Cohen's Kappa, which measures inter-rater reliability and considers agreement occurring by chance, and specific class producer accuracy. Accuracy and Kappa values in this study are reported on a scale of 0–1. Producer accuracy refers to the number of correctly predicted observations of a soil class per the total number of observations of that class (Malone et al. 2017), and is referred in the results simply as class accuracy.

Uncertainty estimates

Prediction confidence was determined using the confusion index (Burrough et al. 1997; Odgers et al. 2011; Brungard et al. 2015), which is calculated with the following equation:

Model performance

Overall, the BSCI-only model was the best performing model. The accuracy for the most probable parent material was 0.52 and Kappa was 0.35 (Table 4). The model performance for the exhaustive model was very similar, but slightly lower, with an overall accuracy of 0.48 and Kappa of 0.30 (Table 5). The model performance for the BSCI-excluded model was lowest, with an accuracy of 0.38 and a Kappa value of 0.19 (Table 6). Predictive results were comparable between the two models that included BSCI environmental covariates (Fig. 4). However, it is important to note that the BSCI-only model could only generate predicted maps for those regions where BSCI imagery exists, which is a significant limitation (Fig. 5).

Fig. 4.

Parent material prediction results for each of the four example areas illustrated in Fig. 1. The first row corresponds to the northeast example area, the second row is the northwest example area, the third row is the southwest example area, and the fourth row is the southeast example area. The figures in column (a) are the prediction results from the BSCI-only model. The figures in column (b) are the prediction results from the exhaustive model. Figs. in column (c) are the prediction results from the BSCI-excluded model. The figures in column (d) are the detailed soil survey polygons displayed by parent material. Coordinates are in UTM Zone 13 N NAD83.

Fig. 5.

Parent material prediction results for the agricultural region of Saskatchewan. The figures in column (a) are the prediction results from the BSCI-only model, column (b) are the prediction results from the exhaustive model, and column (c) are the prediction results from the BSCI-excluded model. Coordinates are in UTM Zone 13 N NAD83.

Table 4.

Confusion matrix for the independent validation results of the BSCI-only model.

Table 5.

Confusion matrix for the independent validation results of the exhaustive model.

Table 6.

Confusion matrix for the independent validation results of the BSCI-excluded model.

In comparison, classification accuracy for regional-scale parent material mapping in British Columbia had an overall accuracy value of 0.779 and Kappa values of 0.697 (Heung et al. 2014). That study focused on topographic derivatives that appear to work well in British Columbia. The greater relief in British Columbia improves results with coarser-scale DEMs with higher vertical error, compared to that in the Canadian Prairies where relief changes related to changes in soil properties can be shallower than the DEM vertical error (Landi et al. 2004). The better results in British Columbia compared to this Saskatchewan study might be explained by the less spatially extensive study area mapped in the British Columbia study (5472 km² in British Columbia versus 298000 km² in Saskatchewan). A second British Columbia study targeting a more extensive study area (945000 km²) had comparable accuracies to this study with values of 0.411 or 0.415, depending on whether a balanced or constrained approach was used for predicting parent material (Bulmer et al. 2016)

Cross-class confusion for the BSCI-only model showed variance in performance across the parent material classes. Due to the lack of validation points, the mapping performance for colluvium and peat classes cannot be evaluated as part of this study (Table 4). Eolian and bedrock classes had few validation points resulting in less certainty in the overall predictions for those classes. The bedrock class observations were predicted across the parent material classes. Predictive accuracy for the eolian parent material class was the highest with a class accuracy of 0.75 (Table 4). Fluvial and till parent materials had the majority of the validation points classified correctly, with the primary misclassification for fluvial being lacustrine (Table 4). Interestingly, an area mapped as till in the southwest corner of the southeast example area was mapped as fluvial deposits rather than till even though it was mapped as till in the original detailed soil polygons (Fig. 4). Without validation data in this area, the accuracy of this mapping cannot be confirmed. However, some of this area could include tills reworked by fluvial processes as they are adjacent to a large glacial lake that was present in this region of Saskatchewan (SKSIS Working Group 2018).

Till was primarily misclassified as lacustrine over till, which is the parent material with the most similar characteristics, and surface expression influenced by the underlying Till. Predictions for lacustrine parent material observations were less than 50% correct, with a class accuracy of 0.47. There was a fairly even distribution of misclassifications across fluvial, lacustrine over till, and till observations. Lacustrine over till observations were classified correctly 25% of the time, with the class being primarily misclassified as lacustrine. This is likely because of the importance of bands 5 and 7, which are related to clay content (Sorenson et al. 2021), to the model, and the similar surface textures between these two parent materials.

Cross-class confusion for the exhaustive model (Table 5) was similar to the BSCI-only model. Overall, this model had slightly lower class accuracies for the eolian, fluvial, and till parent material classes. There was a slight increase in class accuracies for the lacustrine and lacustrine over till classes. There were major changes in the performances per class for the BSCI-excluded model. The accuracies of all classes dropped (Table 6). This was particularly the case for the eolian class. This is likely because of the lack of bare soil imagery Bands 5 and 7 to identify the lower clay content areas associated with eolian parent materials. Based on these shifts in overall accuracy, Kappa, and individual class accuracies, the utility of the BSCI-excluded model is low.

Feature importance

Overall, the most important feature for the BSCI-only model was the Sentinel-1 VH backscatter (Table 3). A possible explanation for the importance of this feature is that VH backscatter has been shown to be related to aboveground biomass (Laurin et al. 2018). By comparison, the VV backscatter was less important (11th most important feature). VV is sensitive to soil moisture, but does vary with vegetation characteristics as well (Vreugdenhil et al. 2018). Particularly under water-limiting conditions, finer textured soils have increased yields and associated biomass (He et al. 2013), and parent materials having higher average clay contents, and therefore biomass, are likely the reason Sentinel-1 VH backscatter was an important feature. Additionally, increased biomass would be associated with areas managed as pasture, which is also associated with coarser parent materials such as fluvial deposits.

The second most important feature was the standard deviation of 9 × 9 median focal filtered elevation with a 101 × 101 focal window. This terrain feature was the most smoothed and coarsest-scale feature of the terrain attributes (Fig. 3), suggesting that the coarser-scale terrain attributes were more useful for parent material mapping in Saskatchewan, at least given the ALOS DEM used in this study. The coarse smoothing resulted in a terrain feature associated with larger-scale landscape heterogeneity. This likely was helpful in differentiating between glacial till with heterogeneous hummocky landforms compared to level or undulating lacustrine landscapes that would have low amounts of terrain variability.

Optical remote sensing features were then the next most important features, with the median ARI being one of the most important of these features (Table 3). The ARI is a function of anthocyanin pigments in foliage (Gitelson et al. 2001) and is likely corresponding to variation in plant composition, which would be closely related to management practices. Coarser-textured hummocky parent materials, such as fluvial parent materials, are more likely to be managed as pasture compared to annually cropped lacustrine deposits. The standard deviation of NDVI was the most important NDVI feature (third-most important feature overall). Areas with higher NDVI standard deviations indicated more changes in NDVI over the growing season. This is expected to be influenced by factors such as management practices, which, in turn, are influenced by soil and parent material type. The standard deviation of NDVI has been shown to be a useful remote sensing feature for land use monitoring and land classification (Becker et al. 2021). The median October NDVI and Median July NDVI were less important as the eleventh and fourteenth most important features, respectively (Table 3).

While not the most important features, BSCI features were reasonably important: the fifth and eighth most important features were the BSCI bands 5 and 7 (Table 3), which correspond to shortwave infrared bands. The importance of these bands could be due to increased absorption in these wavelengths with higher clay contents considering they were found to be the most important bands for mapping clay content in Saskatchewan with BSCI (Sorenson et al. 2021). Clay has absorption features within Landsat 5’s shortwave bands range (Rossel et al. 2010). However, these bands cover a wide portion of the shortwave infrared portion of the electromagnetic spectrum that may reflect other soil properties such as soil organic carbon, so their relationship to clay content is not certain. The sixth, ninth, and tenth most important features were the visible light bands for the BSCI. These bands were important predictors for SOC prediction in Saskatchewan (Sorenson et al. 2021), along with the near-infrared band which was the 12th most important feature. Silt and clay content, which are affected by parent material type, has been documented to be positively related to soil organic carbon concentration in Saskatchewan (Plante et al. 2006).

Standard deviation of a 3 × 3 median-filtered elevation with a 21 × 21 focal window was the 13th most important feature, which would be reflective of finer-scale terrain patterns compared to the 9 × 9 filtered 101 × 101 standard deviation (Fig. 3). Including a feature that characterizes the terrain patterns at different scales likely helped distinguish different parent materials as the contrasts amongst some parent materials were greater at coarser scales (i.e., lacustrine versus lacustrine over till) and others were greater at finer scales (i.e., lacustrine versus fluvial). While inclusion of finer-scale terrain attributes did improve the results, the coarsest-scale terrain attribute was more important (Table 3), at least for mapping parent material in Saskatchewan with a 30 m DEM, and in conjunction with other remote sensing features. While TRI has value quantifying topographic heterogeneity in other prediction soil mapping studies (Brungard et al. 2015), the simple standard deviation of elevation was determined to be more important by the random forest models in this study for the BCSI-only mapping.

The most important features for the exhaustive model were very similar to the BSCI-only model (Table 3). There were, however, differences in terms of the relative importance of each feature. The importance of Sentinel-1 VV backscatter was much greater for the exhaustive model than for the BSCI-only model. The likely reason for this is that because this model included areas where bare soil composite data were not present, the model needed to rely on other feature to distinguish differences amongst parent material classes. There was also an additional terrain feature included in the final selected feature set, particularly different focal window sizes for the standard deviation of elevation and two TRI features. In contrast with the BSCI-only model, BSCI band 3 was not included in the final model. Overall, BSCI was less important and bands 5 and 7 were found to be the most important BSCI bands. The reduction in importance of the BSCI bands could be partially explained by the fact that a significant portion of the training data would not have variance in bare soil values.

The median REIP and CRSI were also important features for the exhaustive model. The REIP is an approximation of the near-infrared red inflection point in vegetation spectra. The location of REIP is highly correlated with foliar chlorophyll content, photosynthetic activity, and a good predictor of leaf area index in wheat (Herrmann et al. 2010). CRSI has been documented to be related to salinity in other contexts (Scudiero et al. 2015). Salinity is affected by many factors with drainage and hydraulic conductivity being two important factors (Eilers et al. 1997). Both factors vary within the different parent material classes in this study, which explains why they were less important for parent material mapping compared to other features.

For the BSCI-excluded model, Sentinel 1 backscatter, standard deviation of elevation, and NDVI attributes were most important, with the standard deviation of NDVI as most important NDVI feature followed by the median October NDVI (Table 3). As bare soil composite data were not present, the model found features related to land use such as the standard deviation of NDVI to be the best as separating parent material classes. Sentinel 1 VH backscatter was the most important feature followed by the standard deviation of 9 × 9 median filtered elevation with a 101 × 101 focal window, and then Sentinel 1 VV backscatter. Other important features were ARI, REIP, the standard deviation of 3 × 3 median filtered elevation with a 21 × 21 focal window, median July NDVI, and then the standard deviation of 9 × 9 median filtered elevation with a 21 × 21 focal window, followed by the CRSI. For all three models, the coarsest scale terrain attribute was the most important terrain attribute, as the greatest differences amongst parent material classes were apparent in the larger-scale terrain variability-associated features. Other studies have found local and landscape morphometric features along with hydrological characteristics useful (Heung et al. 2017). With a higher quality DEM, these features would also likely be useful in the Canadian Prairies as well; however, with the 30 m DEM available, the coarser scale properties were more useful in the context of this study. This could be because the scale of terrain variation in the Canadian prairies is such that a 30 m DEM cannot consistently detect all slope positions in the landscape, and the coarser scale terrain analysis has enough of the upper and lower slope positions in the DEM to better characterize overall variability if not specific slope positions at a given point.