Using Random Forest for Delineation of Dimension Stone Qarries. Mapping of the Santa Eulália Plutonic Complex (Portugal)

1. Introduction

Dimension stones, i.e., natural rocks susceptible of being quarried with the purpose of producing a variety of construction and decorative materials with appropriate extent and shape, need to correspond to rigorous criteria that include size, durability, strength, resistance to polishing and final aesthetic features (e.g., [1]).

Granites are commonly found worldwide with a large variety of textures that makes them a preferred raw material for being quarried. In particular, coarse and medium grained granites are prone to be polished by cost-efficient techniques, and the variety of colours (that necessarily depends on the mineralogical attributes) give them added visual value. The quality of granites used for architectural purposes plays an important role in its valorisation as a dimensional stone, however, the availability of exploitable types according to the global market demands may be challenging, either by the absence of the specified dimension stones in the region or by the lack of knowledge of the resources.

Thematic geological mapping is thus an indispensable step in the characterization of resources, as it should integrate the most important information, such as the identification (e.g., textural variety or different rock types), assess the quality and delimit the structural geometry of the target stone. This is especially important when considering coloured granites that are rather uncommon, such as the reddish and dark coloured (more valorised as ornamental stones), when compared to the typical grey varieties.

In Portugal, granite quarrying is a common practice but there are few exploitation centres that host reddish granites and other darker granitoid stones [2]. Therefore, achieving proper delineation of the extent of such uncommon rocks in areas where outcrops are scarce, near the limits, or where knowledge about the internal structure of the massif is low, could benefit from novel approaches (like Machine Learning and Artificial Intelligence methods) to ensure the stone continuity, which would then be useful for resource valorisation and support industry exploiting these types of stones.

In this work, we propose the use of Random Forest, a supervised Machine Learning (ML) method for classification that has proved to be reliable for thematic mapping, e.g., [3,4,5], as a first approach to attempt to determine distinct types of magmatic rocks (diverse types of “pinkish/reddish” and “greyish” granites, as much as gabbro-diorite rocks) that embody the Santa Eulália Plutonic Complex (SEPC). All the different magmatic facies of the SEPC have been exploited to some extent; nowadays there are fewer exploitations, most of them for ornamental stones, but also for aggregates. This activity left several exposed quarry sites, visible from satellite. Therefore, the SEPC is an excellent work-case scenario to test and apply the ML methods as a support for geological mapping in this context.

2. Regional Framework

2.1. Geodynamic Setting

The Iberian Massif is the southwestern segment of the European Variscan Belt, and comprises the structures emplaced during the Palaeozoic opening and closure of the Rheic Ocean, followed by the continental collision that formed Pangea - also referred to as the Variscan Wilson Cyle (e.g., [6]). The major stages of the Iberian Massif are the Cambrian-Ordovician rifting processes, with marine sedimentation overlaying previous Ediacaran formations and rift-related volcanism (e.g., [7]), followed by the Ordovician-Lower Devonian drifting phase (e.g., [8]), and finally the mid-Devonian subduction and continental collision processes, majorly represented by large magmatism that persisted until Permo-Carboniferous times ([6,9]). The processes related with the compressive phases imposed a tectono-stratigraphic zonation throughout the Iberian Massif, with different tectono-magmatic and metamorphic characteristics according to distance from the Rheic suture zone (e.g., [10]).

The Ossa-Morena Zone (OMZ), located north of the suture zone (in the southern region of the Iberian Massif), hosts several collision-related magmatic bodies. Large intrusions, mostly Middle Mississippian-Lower Pennsylvanian in age, outcrop in the OMZ western and southern borders: the Veiros, Vale de Maceiras and Campo Maior gabbroic plutons [11,12,13], the Beja Igneous Complex gabbroic rocks and the Évora Massif granitoids [9,14,15,16]. Post-kinematic granitoid bodies are less common and mostly occur on the northern and central regions of the OMZ during Upper Pennsylvanian-Cisuralian, such as the Ervedal, Fronteira and Elvas granites, and the large Santa Eulália Plutonic Complex (SEPC; [17,18,19,20,21]). All these lithological units of the Iberian Massif reflect the continental collision effects, that are translated by NNE-SSW and E-W brittle to brittle-ductile strike-slip faults, well expressed in local and regional fracturing patterns (referred as the late-Variscan deformation event [22]).

The SEPC spawns in the northern OMZ region (Figure 1) along with other collision-related stocks, cutting the Ediacaran-Cambrian metasedimentary rocks and rift-related magmatic bodies that occur in the Alter do Chão-Elvas Sector (at southwest) and the Tomar-Badajoz-Córdoba Shear Zone (in the northeast).

2.2. The Geology of the Santa Eulália Plutonic Complex

The SEPC is a granitoid ring-complex emplaced in the northern Ossa-Morena Zone region, as mentioned. It is a late- to post-Variscan magmatic batholith (area ≈ 400 km²) with an ellipsoidal W-E major axis shape that crosscuts the previous pre-collisional NW-SE trending structures from Ediacaran-Cambrian times (c.f. Figure 1). This sub-circular complex consists of an assemblage of intrusive units (ranging from several granite facies to diorite and gabbro) with metasedimentary enclaves that can be correlated with the country rocks. A mappable lithological discernment has been assumed, separating an inner unit of “greyish” granites (informally designated as G1 group), from an external ring of “pink-reddish” granites (G0 group) that envelop numerous enclaves of a diorite and gabbro rocks (M group) and country rocks xenoliths [18,21,24].

Granite facies from both the G1 and G0 groups have been further subdivided in three major zones each based on textural varieties (grain size and colour) that implicate the aesthetic interest for the industry (Figure 2). Likewise, the G1 group is divided in three concentrical belts of grey granite facies with coarse (G1 facies), medium (G2 facies) and fine grain size (G3 facies) changing respectively from the outer to the central, and the G0 group in three areas of distinct fabric around the massif that go from “coarse-grained pink to red” (in the west), “coarse-grained grey to pink” (in the northeast) and “coarse-grained pink to red” (in the southeast) granites. Despite all the divisions, there are authors that further subdivide each of the G0 group facies in smaller patches of the same unit but with a slight variety on the grain size and shade [2,24].

Works regarding the fracturing pattern of the SEPC showed two major directions of sub-vertical to vertical joints (N30° - N45° and N90° - N115°) [26,27] congruent with the patterns of the late Variscan deformation event [22]. The major directions cut all the facies, though structures ranging N320°- N330° derived from late-magmatic and magmatic-hydrothermal processes, especially evidenced in the G1 group area (quartz lodes, veins and greisen structures [28,29].

The inner region of the SEPC is characterized by a typical granite landscape (Figure 3) of whitish-grey granites from the G1 group, commonly outcropping as boulders or rounded blocks, and often forming tors. The outcrops show some particular features such as the development of micaceous schlieren (e.g., Figure 3b) and delta structures indicative of ascending melt flow (compatible with magnetic lineaments published in literature [19]). Gravimetric data show a thick continuity beneath the surface and the probable deep crustal G1 group melt origin, towards east [19,21]. Zircon U-Pb radiometric data provided a consistent crystallization age of 302 ± 2 Ma [25] for the rocks of this group. Despite the limits shown in the available regional geological maps found in literature [2,18,21,24], the spatial relationship between the G1 group and rest of the SEPC is still not well defined.

Contrastingly, the external ring is extremely heterogenous regarding its lithological units. Not only granites from the G0 group outcrop with the most notorious textural variation throughout the region (Figure 4), but also the hosted enclave type varies from gabbro-diorite masses to country rocks xenoliths according to the segment they occur, with an apparent constraint regarding the host G0 facies (Figure 2).

In the western segment, the medium-grained “pink to red” granite is the prevalent facies (Figure 4b), covering the largest area. The presence of xenoliths is very common in this part, either as small irregular enclaves within the granite or as large blocks of country rocks sequences affected by contact metamorphism, mostly hornfels and calc-silicate rocks [30,31,32] geometrically disposed as metasedimentary roof pendants above the western granite facies (Figure 4d). The roof pendants comprise rocks from at least four lithostratigraphic units, namely the Ediacaran and Cambrian schistose formations, the Carbonate Formation, a quartzite bar, and meta-mafic rocks. Enclaves of M group rocks are rarer in this part of the SEPC.

In the northeastern region granites are coarse-grained and have a “pink-greyish” colour (Figure 4c). Despite the common occurrence of xenoliths as small mesoscale enclaves, no roof pendants have been mapped in this area. However, large sites of M group rocks have been delimited as the only macroscale enclaves for this segment.

The southeastern most part of the SEPC is dominated by coarse-grained “pink to red” granite and is the only facies of the G0 group spatially associated with both enclaves of the gabbro-diorite M group and with roof pendant xenoliths.

The increase of the reddish intensity on all granite facies from the G0 group is assumed to be due to the impregnation of iron oxides in feldspars [2] which provide textural aspect that range from Figure 4c to the granite aspect from Figure 4b.

The “pink-reddish” granites of the SEPC have been dated with zircon U-Pb radiometric methods which yielded crystallization ages of 297 ± 4 Ma [20], and 301 ± 0.9 Ma [25], as well as Rb-Sr isochron of mica concentrates that indicate an age of 290 ± 5 Ma [17].

The gabbro-diorite rocks from the M group outcrop as minor masses and lenses within the northern and eastern segments of the external aureole. The limit between rocks from the M group and the G0 granites rarely outcrops, and is characterized by a sharp contact, with apparent emplacement of the granites over the gabbro (Figure 5a). The gabbro-diorite association commonly shows complex magmatic structures of mingling and mixing/unmixing, suggesting different coexisting melts that formed either by different magma injection events, replenishment, or by immiscibility (Figure 5b,c). Available geochronological data determined 303 ± 3 Ma [20] and 307 ± 2.5 Ma [25] zircon U-Pb ages for these rocks.

The close overlapping ages determined for the rocks of the SEPC and the contrasting geochemical and petrographic features found especially between the G1 and G0 groups induced a debate regarding their nature and genetic association that still remains open, of either the G0 and G1 derived from a primeval melt, or there are distinct sources of liquid. In any case the M group rocks are considered as autoliths present in the external facies. But regardless the interpretation of petrogenesis, and independently of the origin and contemporaneity of the melts that gave origin to these rocks, the SEPC consists of an important target for dimension stone exploitation, and all granite and gabbro facies have been exploited to some extent in the past. A list of exploitable facies of the SEPC can be found in Table 1.

3. Results

This chapter presents the EDA and the RF modelling results, divided by the datasets. The latter is further sectioned by the two model predictions, firstly to the Test dataset, and then to the SEPC. In both instances, the models were subjected to the performance metrics, namely confusion matrices (see supplementary materials), accuracy (overall, producer’s, and user’s), and Kappa coefficient, proceeded by data rasterization as final output.

3.1. Exploratory Data Analysis (EDA)

3.1.1. SEPC Lithological Atlas

EDA is essential to assess data balance and band correlation, providing a strong support for data validation claims. A better understanding of the Subset areas can be achieved through characterization of each component, including the spectral signatures, for each group. Figure 8 contains areas representative of each lithology (for complete list, see supplementary materials). Here, the masking process of removing vegetation, water and roads/ cities are filled with white pixels (Masked). Corine Land Cover (CLC18) reveals no relation between land cover uses and the lithotypes, and as such it was only used to identify and remove the urban fabric (not represented here). Likewise, the spectral signatures (SS) do not reveal any clear difference between the groups.

3.1.2. Descriptive Statistics

A comparison of the results between the two datasets (Complete and Subset) to the dataset before masking the pixels (Original data) was made, assessing for pixel proportion (and respective frequency) (Table 2), and correlation matrices (Figure 9). Lastly, Table 3 indicates the minimum, median, and maximum values of each band, grouped by dataset.

Table 2 shows that G0 is the largest group across all groups, representing ~57% of the total number of pixels, followed by G1 group, averaging ~37%; M group has a small representation, with less than 6% average. Across datasets, the difference between the groups is smaller than 3%, indicating that group frequency and proportion is maintained.

Based on the correlation matrices, four band groups can be distinguished, common in all datsets: i) B1 is always a separated group, ii) Bands 2 to 5, representing visible and very near infrared (VNIR), iii) Bands 6 to 9, representing near (NIR) to shortwave infrared (SWIR), and iv) Bands 11 and 12, representing the higher SWIR wavelenghts.

3.2. Complete Dataset Models

3.2.1. Complete Testing Sets Predictions

The testing dataset predictions are fundamental to prevent overfitting and provide a general assessment to model performance. Confusion matrices were produced, which compare the classed labels to the predicted outcome, and allow for the extraction of OA, UA and PA, that validate the model along with the Kappa coefficient (Table 4).

Table 4 reveals a positive correlation between the percentage of training data used and the validation metrics, with Kappa coefficient classifying these predictions as moderate (10%, 30%, and 50% training sample sizes) to substantial (70% and 90% training sample sizes) according to [41] proposed scale. At 10% the UA is approximately the same as the OA, but has a higher growth rate, attributed to the M group, with more correct predictions than G0 and G1 groups (Figure 10a, see supplementary materials for individual accuracy classes for each model), albeit only reaching “good prediction” results. PA shows the lowest scores off all accuracy metrics, attributed to the poor (10%) to moderately poor (30%, 50%, and 70%) M group percentages (Figure 10b). However, the PA’s increased growth is also attributed to the same group, as G0 and G1 have a less steep growth incline.

3.1.2. Complete SEPC Group Predictions

The models were used to predict the SEPC, i.e., the Complete dataset before splitting into train and testing datasets (Table 5). OA shows good predictions, increasing from 75,94% to almost perfect at 98,33% accuracy. UA’s lowest accuracy is ~75%, by G0 group (with 10% training sample size), closely followed by G1 group, and then M group (Figure 11a). As the training sample size increases, the difference between the three groups is reduced, with all three having very good to perfect UA predictions. PA is the lowest accuracy metric, with a more pronounced difference in the 10% and 30% training sample sizes, attributed by the M group (Figure 11b). The Kappa coefficient values show a great increase in model confidence, from moderate (10%), substantial (30% and 50%) to almost perfect (70% and 90%).

The UA and PA for both testing and SEPC group set predictions reveal that the M group has greater UA scores and lower PA scores than the other groups. This is indicative of a large number of M group predictions matching the real M group labels, but for a small number of M group pixels accurately predicted as such.

3.3. Subset Dataset Models

3.3.1. Subset Testing Sets Predictions

As with the Complete dataset, the Subset models were used to predict a test dataset, essential for data validation and to provide an assessment on model performance, including the calculation of accuracy and Kappa coefficient metrics (Table 6).

Table 6 shows an increase in all accuracy and Kappa coefficient values, as the percentage of data used to train the models is increased. UA is higher than OA, due to the very good classification of M group rocks, ranging from 81 to 93% (Figure 12a), while G0 and G1 rocks share approximately the same percentage as the OA. PA increase in the different models can be attributed to the larger correct number of pixels in the M group, from 44.32% in the 10% Model, to 73.53% in the 90% Model (Figure 12b). The G0 shows the smallest growth, with an average increase in accuracy of 1.68% in each step. However, this group’s predictions are the closest to a near perfect prediction record.

The Kappa coefficient is between 0.52 and 0.74, classified as moderate (10 and 30% models) to substantial (50, 70, and 90% models).

3.2.2. Subset SEPC Group Predictions

Following model validation by the test dataset, the predictions on the study area was performed, classifying the soil and rock pixels of the SEPC as either G0, G1, or M group rocks.

Table 7 indicates a moderate OA ranging between 67.95 and 70.84%, with a poorer performance than test dataset. Kappa coefficient is also lower, classified as fair (10 and 30% models) to moderate (50, 70, and 90% models). Both UA and PA steadily increase, albeit with at a slower rate compared to the test dataset. Contrarily to the test dataset, the UA value for M group is consistently low, between 38% and 45%, balanced by the G0 and G1 group classifications, that are ~70% across all models (Figure 13a). The PA shows very good classifications on G0 group (~81%), moderately good classifications for G1 group (~58%), and poor to moderately poor for M group (~27%), showing the variable difficulty in classifying and identifying the different groups (Figure 13b).

Comparing both testing and SEPC group set predictions, the observed behaviour for G0 and G1 is the same, with approximately equal UA percentages, and PA values that grow approximately at the same rate. However, it is noteworthy that while M group’s behaviour for the Testing set is like the one encountered in the Complete model predictions, the same is not found when predicting the SEPC, and the UA is far below the remaining groups.

For the SEPC group sets, the predictions were rasterized into thematic maps following the new classifications (Figure 14a,b, see supplementary materials for all figures). Based on the total number of pixels used in training and predicting the SEPC groups, the Subset model with 90% training size sample (which equates to approximately 11.93% of the total number of pixels in the SEPC) was compared to the Complete model with 10% training sample size, allowing to match pixels with the same classification predicted (coloured by respective group colour) and those that do not match (coloured in blue) (Figure 14c). This last figure is also overlayed by the geological map, informing which regions were misclassified.

4. Discussion

This work explores an established workflow for remote sensing data classification [38,39], including the removal of vegetation, water bodies, urban centres and roads. The RF algorithm trained two models (Complete and Subset datasets), with 10%, 30%, 50%, 70%, and 90% of data, used to predict the test and the SEPC groups. Results were validated using confusion matrices, supported with accuracy and Kappa coefficient metrics.

4.1. Test Dataset

The test dataset was used to verify the validity of the trained models, using the remainder of the data that was not used for training. The results indicate a good overall accuracy, without reaching perfect predictions, inferring a well-trained model without overfitting. Furthermore, the high UA percentage indicates good reliability on the predicted classifications, as they are mostly correctly classifying the real features. This is especially visible for the M group, where the UA is at least 5% higher than the other classes and reaches 95.70% on the 90% model. It is noteworthy that most of the incorrectly identified M group classifications fall within the G0 group, that envelop these rocks. The PA reveals that G0 group, the pink-reddish granite, are the best classified rocks in the test dataset, with an accuracy above 86%. Conversely, the M group shows low PA, detailing that most of the M group is not classified as such. However, this is also the group with the highest increase in accuracy as the train set is increased, showing how important it is to use a minimum quantity of pixels for each class to be correctly classified. As such, the differences in UA and PA for the train dataset can be used to infer that the G0 and G1 groups are either relatively homogeneous on their respective spectral signature, or that their intra-variational differences are sufficiently well represented within the confines of the chosen training areas. This does not seem to be true for the M group rocks, requiring a larger sample size to be correctly classified, as observed for the 50%, 70%, and 90% models. This group is composed by a mix of gabbro-diorite rocks, and as the diorites are usually found as masses at outcrop scale, is unlikely to interfere with the spectral signature.

4.2. SEPC Dataset

The main interest in training the models was to assess how well the RF classification algorithm could classify the SEPC groups using regions of training. The performance metrics indicate an overall moderate performance in accuracy, slightly increasing as a larger sample size is used. The average UA is lower than the average OA, but in individual accuracies, both G0 and G1 groups share approximately the same percentage than the OA, and the M group unit is responsible for the lower UA, contrarily to the observed dataset predictions.

The individual assessment of the classes for PA clearly shows that G0 predictions mostly correctly classify G0 areas mapped in the geological maps, moderately predicted the G1 group, and have a poor performance classifying the actual M group. Furthermore, both UA and PA have only a small increment of accuracy when increasing model percentage, indicating that the models, using this configuration of selected areas as training data, have reached a plateau of the feasible accuracy. This limitation is imposed by the total number of pixels considered in the dataset, since at most (90% training sample size) it only represents about 11% of the total SEPC dataset.

The combination of UA and PA leads the authors to believe that the G0 granites can be correctly classified using this methodology, however the M group gabbro-diorites need further refinement, either by increasing the number of pixels used as training, or by using in field knowledge to capture more detailed intra-variational differences.

When comparing the total number of pixels used, the OA accuracy obtained from the 90% training sample size is inferior to the Complete 10% training size model (70.84% to 75.30%). The thematic map that identifies the difference between the models (as blue pixels in Figure 14c), overlayed with the geological map, infers that G0 group has the most reliable prediction using the training sets. The G1 group has a salt and pepper effect, which may denote a possible lack of representation in all intra-variational differences for each lithology. The M group has the least number of concordant pixels between the models, as verified by the UA (Figure 12a), contrary to the trend observed for all other predictions. Nonetheless, the authors believe that the moderate OA could bode well for the lithological classification, since the limits between G0 and the rest of the units are not well defined, and the limits between M group and G1 rarely outcrops, indicating that some of the defined limits may not be entirely correct. Given the difficulty in classifying the last group, a higher resolution approach should be considered, such as drone-based spectral ground-truthing techniques, that better separate each class, refining the spectral signatures.

5. Final Remarks and Future Insights

The SEPC continues to be actively exploited for dimension stone quarries, sought out for physical characteristics. As such, the need for accurately mapping the lithological distribution, in an efficient and cost-effective way, becomes paramount for decision-makers who want to open or further progress their quarries. Despite the differences found across the models, the results indicate that using RF in known training areas is a useful “first approach” tool for delineating the limits of the major groups defined in the SEPC. We observed large scale predicted maps (based on Sentinel-2 bands) that provided overall good predictions for the G0 group, moderately good predictions for the G1 group and poor predictions for the M group, and a general outline compatible with the established geological maps of the region. The raster interpretations are supported by the validation metrics, and both UA and PA indicate good classification percentages for G0 group, but poor for the M group. Kappa coefficient confirms that the models have a fair to moderate strength of agreement, confirming its usefulness as a free-cost tool for first limit delineations.

However, we only used satellite images limited by the spatial resolution (10m pixel size); the next steps embrace a high-resolution approach, in which the results obtained herein will be confronted more robust datasets from several key-areas within the SEPC. This will complement the results and conclusions presented here and includes:

- Multispectral drone imaging to quickly capture the spectral signatures of each selected area, representing each group and boundary zones.

- Multispectral analysis of representative rock and soil samples. These will be used to train the predictive model (i.e., calibration).

- Application of RF and other machine learning algorithms to the selected areas to assess its ability for lithological classification and final ground-truth validation.