Next Article in Journal
A Data Fusion Method for Generating Hourly Seamless Land Surface Temperature from Himawari-8 AHI Data
Previous Article in Journal
Nonlinear Unmixing via Deep Autoencoder Networks for Generalized Bilinear Model
Previous Article in Special Issue
Towards Streamlined Single-Image Super-Resolution: Demonstration with 10 m Sentinel-2 Colour and 10–60 m Multi-Spectral VNIR and SWIR Bands
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Discrimination of Rock Units in Karst Terrains Using Sentinel-2A Imagery

1
Croatian Geological Survey, Sachsova 2, 10000 Zagreb, Croatia
2
Faculty of Geodesy, University of Zagreb, Kačićeva 26, 10000 Zagreb, Croatia
3
Faculty of Science, University of Zagreb, Horvatovac 102a, 10000 Zagreb, Croatia
*
Author to whom correspondence should be addressed.
Remote Sens. 2022, 14(20), 5169; https://0-doi-org.brum.beds.ac.uk/10.3390/rs14205169
Submission received: 26 August 2022 / Revised: 26 September 2022 / Accepted: 12 October 2022 / Published: 15 October 2022
(This article belongs to the Special Issue Open Access Satellite Imagery Processing and Applications)

Abstract

:
We explored the potential incorporation of Sentinel-2A imagery for rock unit determination in the Croatian karst region dominated by carbonate rocks. The various lithological units are potential sources of both stone aggregates and dimension stone, and their spatial distribution is of high importance for mineral resource management. The presented approach included the preprocessing and processing of existing analog data (geological maps), Sentinel-2A satellite images and the United States Geological Survey spectral indices, all in combination with ground truth data. Geological mapping and digital processing of legacy maps using the K-means and random forest algorithm reduced the spatial error of the geometry of geological boundaries from 100 m and 300 m to below 100 m. The possibility of discriminating individual lithological units based on spectral analysis and discriminant function analysis was also examined, providing a tool for evaluating the geological potential for mineral resources. Despite the challenges posed by the lithological homogeneity of karst terrain, the results of this study show that the use of spectral signature data derived from Sentinel-2A satellite images can be successfully implemented in such terrains for the enhancement of existing geological maps and mineral resources exploration.

1. Introduction

Remote sensing methods are very limited regarding subsurface depth. Therefore, they are often used as support for or an addition to surface geological research [1]. The use of multispectral satellite data in geological mapping dates back to the 1970s with the establishment of the ERTS (Earth Resources Technology Satellite) system by NASA (National Aeronautics and Space Administration), later renamed the Landsat Program [1,2]. Regarding the study area, the application of multispectral satellite image processing dates back to the middle of the 1970s [3,4,5]. Most of the published studies in the multidisciplinary scientific discipline of remote sensing–geology–mining deal with terrains consisting of lithologies with substantial differences in spectral signature, especially when it comes to the field combination of various sedimentary, metamorphic and igneous rocks [6,7]. In addition, many of these studies refer to the exploitation of ores which commonly represent distinct mineralogical outliers within their respective terrains [7,8,9,10]. The mineral resources discussed in this paper are dimension stone, stone aggregates, raw materials in the cement industry, gravel and sand.
Karst regions are specific terrains characterized by uniform mineralogical composition, which poses a challenge for carbonate unit discrimination by remote sensing techniques. Most efforts to study karst terrains by remote sensing techniques have concentrated on the mapping of distinct karst landforms, hydrogeology studies and monitoring environmental change [11,12]. Here, we examined the possibility of carbonate lithological unit discrimination (on a medium and large scale) for geological mapping and mineral resource exploration purposes using European Space Agency (ESA) Sentinel mission multispectral imagery in one of the most distinct karst terrains in the world, the Dinarides. Sentinel-2A satellite images were chosen because of the positive relationship availability and spatial and radiometric resolution possibilities that exist according to available papers. The potential of Sentinel-2A sensors in geological research was discussed in [13]. Machine learning (ML) methods, which are based on the theory of statistics to build mathematical models [14], were chosen as appropriate to be tested in these regions. The methods chosen were unsupervised classification using the K-means algorithm and supervised classification using the random forest (RF) algorithm [15,16,17].
The Croatian karst region is characterized by large areas of bedrock surface exposure. The legacy geological maps (Basic Geological Map of the Yugoslavia to the scale 1:100,000) are based on direct interpretations of field and areal photogrammetry data. These maps are used in spatial planning for evaluating the regional mineral potential and the management practices related to mineral resources. Since the maps are used for spatial plans at a scale of 1:25,000, there is a need to improve their accuracy.
The chosen study area in the Dalmatia region of Croatia is distinct as it hosts not only carbonate units of the Mesozoic Adriatic Carbonate Platform [18] but also calciclastic units of a foreland basin fill [19] and a Quaternary cover largely derived from their erosion.
In this study, we used detailed geological map data (up to 1:10,000), which we produced by field work using precise global navigation satellite systems (GNSSs). Here, available data from the spectral library at the United States Geological Survey (USGS) were applied in addition to digital satellite images. The normalized difference vegetation index (NDVI) was also used to distinguish rocks from vegetation, as well as digital orthophoto images (DOF) of the State Geodetic Administration of the Republic of Croatia (SGA). All obtained data were analyzed with geographic information system (GIS) tools and processed by discriminant function analysis (DFA, [20]).
In this study, we tested whether basic ML algorithms can be used during geological/lithological mapping of karst terrains. In addition, the possibility of improving the geometry of existing geological boundaries using the ML method and the reduction of the time required to conduct mineral resource surveys in karst terrains were tested. The deviation in the geometry of the geological boundaries, i.e., their accuracy, due to geological mapping methods used in the past do not meet the standard of modern spatial planning or of locating potential mineral exploration and exploitation targets (deposits). Therefore, the aim was to reduce such deviations as much as possible, preferably at the beginning of the project.
The assumption was that one of the possible ways to reduce deviations in the geometry of geological boundaries is using a combination of available satellite imagery, machine learning and classical fieldwork methods. Remote sensing data processing using ML algorithms at the medium-scale map level in various studies has demonstrated that such an approach can contribute to the accuracy enhancement of existing geological maps [11,21,22,23,24,25]. Geological mapping in combination with satellite imagery processing presents an approach for spectral signature location reading at the large-scale map level that has not been applied to the karst area of the Dinarides. Based on this method, spectral signatures were read and statistically processed, and lithological units in karst terrain were discriminated.

2. Materials and Methods

2.1. Study Area

The study area encompasses a part of the eastern Adriatic coastal region with an area of 8982 ha, situated in the central part of Northern Dalmatia, Croatia (Figure 1a). The relief includes the southern slopes of the Promina Mountains (Figure 1b). The study area is covered by the Basic Geological Map of the Yugoslavia (OGK SFRJ), sheet Šibenik, at a scale of 1:100,000 [26,27] and by the field geological map of the OGK SFRJ, sheet Bribirske Mostine, at a scale of 1:25,000 (Figure 1c). The terrain proved to be suitable due to the great exposure of outcrops, lack of vegetation and the geological potential of strategically important mineral resources for the Republic of Croatia. This area is part of the Croatian karst region of the Dinarides (Figure 1a), which is represented by a thick, largely carbonate succession deposited on carbonate platforms of different ages, types and paleogeographic settings from the Upper Carboniferous to the Eocene [18]. The oldest deposits on the Bribirske Mostine field geological map (Figure 1c) are rudist limestones (Figure 1c, Senonian, K2,3), in which the percentage of CaCO3 ranges from 95% to 99%. They range from microcrystalline to recrystallized limestones and, rarely, are dolomitized bioaccumulated limestones. Senonian rudist limestones (SRL) are used as dimension stone and/or as stone aggregates and are important mineral resources [28].
Paleogene deposits include Eocene foraminiferal limestones and calciclastic deposits of the Eocene and Lower Oligocene “Promina Beds” [27]. Eocene foraminiferal limestones (Figure 1c, E1,2, FL) are well-layered, light-gray microcrystalline limestones containing up to 95–98% CaCO3 and have significant potential for exploitation as dimension stone and/or stone aggregates [28]. The next three lithological units are part of the “Promina Beds”. The first is made up of fine-grained detrital to marly limestones with traces of marl and carbonate conglomerates (Figure 1c, E2,3, LMC). Farther upward in the succession, two units consist of dominant carbonate conglomerates and alternations of these conglomerates with limestone and marls, as well as platy limestones (Figure 1c, aE3 (CO), bE3 (PL)). The lithological members of the “Promina Beds” can be exploited as raw materials in the cement industry as dimension stone and/or stone aggregates [28]. Quaternary sediments (Figure 1c,d (CL) and j (SLS)) consist mostly of gravel, sand and terra rosa. The colluvial (deluvial) deposits (Figure 1c,d, CL) are essentially deposits of weathered flysch, mostly gravels and sands of different grain sizes. These Quaternary sediments have geological potential for the exploitation of sand and gravel [28].

2.2. Datasets

The field geological map of the OGK SFRJ, sheet Bribirske Mostine, 1:25,000 (Figure 1c), a satellite image obtained with a multispectral sensor on Sentinel-2 satellites (MSI S-2A), the USGS spectral library data on limestones [30] and geological ground truth field data were used. Level-2A data using ortho-corrected bottom-of-atmosphere (BOA) data with subpixel multispectral notation were used [31]. During December 2021, two sites (Subsets 1 and 2) were geologically mapped. Subset 1 (Figure 2a) covers an area of 100 ha, where 36 observation points were recorded. Four (4) different lithological units were extracted: SRL, FL, CL and SLS (Figure 2a). The working principle was the same for Subset 2 (Figure 2b), covering an area of 134 ha. In this case, 35 geological field observation points were recorded, and two lithological units were singled out. These were PL and CO (Figure 2b).

2.3. Imagery Processing

At a medium scale, a digital form of the Bribirske Mostine geological map was created and transformed into the HTRS/96TM system. A zip file with title ID T33TWJ of 896,68 MB size was downloaded using the Copernicus Open Access Hub. The ingestion date was 2 April 2020. Resampling was used, the goal of which was to reduce the product to a single spatial resolution [15]. The resampling algorithm was bilinear, the downsampling method was mean and the flag downsampling method was first (Sentinel Applications Platform (SNAP) [32]). According to the SNAP transformations, the product was transformed from the WGS84 to the HTRS96/TM geographic system. First, the K-means algorithm was used [33,34]. Within the K-means classification, all bands of MSI S-2A were used, and a composite of bands 2, 3, 4, 8 and 11, all for 5, 10, 15 and 20 classes, was used separately.
The K-means algorithm is one of the simplest and most popular unsupervised classification algorithms, the aim of which is to divide M points in N dimensions into K clusters so that the within-cluster sum of squares is minimized. A matrix of M points in N dimensions and a matrix of K initial clusters are required as algorithm input data. Sets of objects are grouped in such a way that the objects in the same cluster (group) are as similar as possible. Similarity is determined by the distance between two data points, and the Euclidean distance is used to measure this distance [35,36]. The K-means algorithm objective is reducing the intra-cluster distance and increasing the distance between clusters [33,34,37,38]. This algorithm is currently used in numerous remote sensing projects [39,40,41]. For instance, a novel method based on a combination of spectral indices and K-means was developed for automatic land cover classification (ALCG, [42]).
The RF algorithm was subsequently used [43] to carry out the following important functions. The first was to serve as a means to conduct an impartial revision of the initial geological datasets, while the second was to serve as a basis for creating the final geological map [44,45]. In supervised classification, the user determines the input data, in this case, pixel values. After that, the algorithm uses the spectral signatures of the training sites and applies them to the entire image. It is a regression algorithm that works by creating a forest with decision trees, the object of which is to create a prediction model, the results of which are derived from a combination of different input variables. At the same time, the larger the number of trees in the forest, the more accurate the results [15,16,43,46,47]. The result of such a process, in this case, was a classified lithological map with unique and separable spectral signatures of each class. A review of machine learning in processing remote sensing data for mineral exploration was given by [48].
In this study, the input data for the RF classification were a K-means and DOF-derived geological map of Bribirske Mostine. Based on this map, 5200 circle polygons were selected to serve as ground truth data in RF classification. These polygons were placed on all lithological members, as shown in Table 1. It can be seen that a variant analysis was made for colluvium and platy limestones. The exact position of the input circle polygons is shown in Figure 1b.
The total number of circle polygons was determined in accordance with the surface distribution of each individual lithological member. The shape of the polygon was defined as a circle with a 10 m radius (20 m diameter). Such a radius was chosen due to the spatial resolution of MSI S-2A of 10 m. Due to possible deviations in the geometry of boundaries, the polygons were arranged in such a way that they were not positioned directly on the geological boundary. In the RF classification, the MSI S-2A band combination was used as in the K-means classification. Two variants of RF classification were considered. An account of the two variants of the RF classification is illustrated in the Results chapter. It is clearly marked with a red hachure.
Regarding large-scale maps (Figure 3), data collected from geological field mapping were used as ground truth data, i.e., validation data for Subset 1 and Subset 2 sites. The distribution of limestones in the study area was obtained by processing the USGS spectral library dataset. The NDVI was used to discriminate vegetation from the other forms of land cover [49]) and as a backup tool for evaluating the position of points or pixels (pins) on which spectral records were determined and statistically analyzed.
File T33TWJ was imported into SNAP. Selection points (Figure 4 and Figure 5) based on NDVI (Figure 4a and Figure 5a), the RGB composite (Figure 4b and Figure 5b), USGS limestone index (Figure 4c and Figure 5c), DOF (Figure 4a and Figure 5a), and ground truth data generated in ArcMap followed (Figure 2 and Figure 3). The main aim of using a combination of these five types of backgrounds was the maximum reduction of vegetation influence on spectral signature readings. Based on the location of these points (pins), 40 spectral signatures of each lithological unit were read and statistically processed. Six lithological unit spectral signatures were singled out (Figure 3 and Figure 6).

2.4. Discriminant Function Analysis

After spectral analysis (Figure 4, Figure 5, Figure 6 and Figure 7), a DFA was conducted on large-scale maps to obtain more precise data on the differences in spectral signatures within lithological units. In this study, DFA was used to explore the differences in the spectral index data-bearing potential for improving classifications in a few distinct lithologies mapped in the karst area of the Dinarides. In methodological terms, DFA is all about finding differences between a number of originally observed (natural) groups with the least possible margin of error made during classification [20,50,51]. In geology, it is regularly used in geochemistry, environmental geology or bauxite studies [52,53]. In this study, the purpose was to construct a predictive discriminant function model (DFM) with the highest possible classification efficacy, taking into account six a priori outlined groups of mapped lithological units (see legend in Figure 2)—CO, PL, SRL, FL, CL and SLS—and 12 variables representing the spectral indices derived from satellite band imagery—B1, B2, B3, B4, B5, B6, B7, B8, B8A, B9, B11 and B12.

3. Results

Based on the K-means algorithm maps (Figure 8) and DOF (Figure 1b), a satellite geological map of the Bribirske Mostine area was derived (Figure 9). The GBs derived in this way were compared with the original GBs of the Bribirske Mostine area (Figure 1b). The first result of the comparison was that the clearly observed difference in the geometry of GBs was not linear (Figure 1b). The second result concerned differences in the interpretation of lithological units, which were sorted into four (4) levels of differences (LODs, Figure 9 and Figure 10). Additionally, in Figure 9a, the lithological unit marked “CL, colluvium, gravel and sand” is highlighted with hachures in red. It was clearly extracted using the K-means algorithm (Figure 8) but was not in the original geological map (Figure 1c). Subset 2 site was geologically mapped during December 2021 (Figure 2b and Figure 9b,c). Then, it was concluded that this is an area represented by a lithological unit marked “bE3, platy limestones, PL”, which follows the original geological map (Figure 1c). However, this area is covered with shallow soil with platy limestone (PL) fragments, as clearly shown in Figure 9d. By processing analog and K-means data results that correspond to ground truth data, apropos PL and CL exchanges were obtained, as shown in Figure 9d.
The results of digital satellite image processing using the RF algorithm (Figure 11) were presented as the relative difference between the original map and, thus, created geological maps. In all four cases (Figure 11a–d), the anomalies were clearly separable, in a very similar way to the K-means algorithm (Figure 8). In all four cases (Figure 11a–d), the lithological unit labeled “K2,3, rudist limestones (ex. Senonian, Maastrichtian, Campanian, Santonian, and Coniacian) SRL”, was not extracted, nor was the GB between it and the lithological unit marked “E1,2, foraminiferal limestones, FL”. In the first three cases (Figure 11a–c), in the southern part of the map, green-colored pixels were dissociated, which may indicate rudist limestones, while, in the fourth case, they were not. Based on ground truth data derived by geological field mapping, it was concluded that these were also anomalies caused by dense vegetation reflectance.
In the DFA of the investigated lithological units, six groups (K = 6) and twelve variables (descriptors, p = 12) were processed, resulting in a DFM of five discriminant functions (K − 1). Relevant multivariate tests were used before the analysis (Table 2), revealing extremely low associated probability (p = 0.000), a necessary precondition for the computation of discriminant functions (DFs). This breakdown process also included the univariate statistics of input data displayed as grouped boxplots for each band and lithological group (Table 3 and Figure 12). Results showed (p-level and eigenvalue, Table 2) that the first two DFs, explaining over 93% of the total variance, effectively configure the discriminant space (Figure 13a,b).
Essential in DFA, is derivation of the functional model (describing natural, geological processes) out of the structural model (mathematical structure) of the analyzed data. To this effect, discriminant loadings are especially helpful in validating how much of the overall discrimination is explained by each separate variable (band). Reading DFA in functional terms highlights the natural (geological) processes which are assumed accountable for the pattern of the structural model derived from spectral index data. Clearly, in building the best discrimination strategy, all variables with small discriminant loadings can be safely discarded for interpretation.
In this regard, comparison between the variable and group diagrams (Figure 13a,b) was the most helpful technique for structure-to-function relation and elucidation of the key reasons for group separation. In the process, affiliation between variables and associated groups was always interpreted considering their joint position along the corresponding discriminant axis in both diagrams (Figure 13a,b).
Regularly, the mutual relationship between the variables and group means (centroids) along the respective DF defines the geological (via remote sensing technology) nature of a computed mathematical model.
By virtue of variables (multispectral imagery bands) assuming the role of process descriptors, each DF was labeled geologically, contributing to the overall separation of groups (lithological units). Thus, the first discriminant function (DF1) was indubitably the greatest discriminator, clarifying more than 87% of the total variation in the analyzed model and, together with DF2, just over 93%. Despite their high statistical significance (p = 0.000), other functions could be safely excluded from further analysis owing to their low contribution to discrimination between groups (low eigenvalues) and low capability of inferring the relevant geological process.
The scatterplots of variable loadings and group means (centroids, Figure 13b) were bipolar, separating the FL/SRL from both the CL/PL and CO/SLS groups. However, a closer inspection revealed that all descriptor variables (bands) loaded on the positive side of DF1, a peculiar situation suggesting that the discrimination between FL/SRL and other groups posted on the negative pole of DF1 on their respective diagram was based solely on the “inverse image”, that is, the “absence” of relevant descriptor variables. In practical terms, this means that the respective group pairs (CL/PL and CO/SLS) were characterized entirely by decreased values in all descriptor variables regarding FL/SRL. Additionally, B1 was the variable with the highest discriminant potential, while the rest was crammed in several clusters, possibly indicating a more parsimonious scheme of band selection in further investigations. These include B2/B11/B12, then B8/B8A and, finally, B6/B7/B9 as obvious variable clusters. The B3, B4 and B5 bands can be seen as “transitory” variables communicating between the aforementioned clusters (Figure 13). Eventually, owing to its highlighted position in the discriminant model, DF1 clearly underscored the central theme of the analysis—separation of distinct lithological units mapped in the karst area of the Dinarides based on selected spectral indices (bands) derived from satellite imagery.
This was established by the high discrimination of both limestone formations (SRL and FL) compared to other lithological units marked by increased values of all spectral indices (bands), particularly of B1 (followed by B2, B11 and B12). Further insight into spectral differentiation of the investigated lithological groups was not impeded by the inferior discriminant potential of DF2 (Figure 13). Allocating the abovementioned “inverse image” properties (the fourth quadrant variable clustering) to DF2, the latter was predominantly involved with isolating SRL from the FL group (Figure 13b). As with DF1, separation revolved around the “absence/presence” polarized rapport contrasting SRL with FL via the generally higher spectral indices of the former, particularly B5 but also B6/B7/B9, B4 and B1 and others most distant from the axis intersection, as seen from the scatterplot of variable loadings (Figure 13a) depicting the group distances.

4. Discussion

Existing geological maps are the basic input data for the creation of the mineral resources potentiality maps, which are, in turn, used in the construction of spatial plans as basic documents for the physical planning of a certain administrative areas (state, region, county). In addition to geological data, an integral part of spatial plans is the data of other territory users. In most cases, they are regularly updated and roughly correspond to the time frame of the spatial plan construction. However, updating the geological data of a wider area is demanding and requires a significant expenditure of time and financial resources.
Therefore, the possibility of finding an efficient way of updating existing geological maps was considered, with the aim of obtaining more precise data on the geology of the area. Further in the discussion, the facts about the existing geological maps that are used as input data are commented on, followed by a discussion in two parts. The first part refers to the applicability of ML algorithms for the correction of geological boundaries as the basic elements of a geological map, in this case, at the scale 1:25,000. The second part refers to the use of spectral signatures of rocks and statistics in extracting lithological members on detailed maps at the scale 1:10,000. It is discussed whether such methods can be used in more precise determination of mineral resources potential. The karst area of the Dinarides was chosen due to the availability of rock outcrops and lack of vegetation but also the presence of mineral resources that are strategically important for the Republic of Croatia.
The creation of the OGK SFRJ in the 1960s and 1970s provided a basic geological map which is the necessary basis for further geological research and planning, including research of mineral resources [54]. This map is based on a chronostratigraphic division, and, despite certain disadvantages in correlation with maps of neighboring countries or some neighboring sheets, it is still extensively used in research related to mineral resources. Its use in combination with digital GIS technology and numerous new, precisely obtained field data (GNSS, real-time kinematic (RTK), etc.) has revealed certain deviations from the recorded field data. These deviations are primarily related to the geometry of the GB. As mentioned before, there has been a lack of previous research on lithological unit discrimination within karst terrains using remote sensing methods. Therefore, there was an initiative by the authors to initiate these studies. Conversely, technology development, and, more recently, the development of machine learning algorithms, mathematics and statistics in general, has intensified on a global scale. However, due to a lack of previous research on the mapping of karst terrains, here, we initially explored basic K-means and RF tools. Following these results, efforts will certainly be made to further implement other remote sensing methods and their mutual combination.
By combining the results collected with the K-means algorithm and DOF, four (4) LODs were identified. It is important to emphasize that these were determined LODs in geological maps without further quantification. Ideally, because of the number of LODs, it would be best to conduct the remapping of the entire area. However, this is unlikely to occur in the foreseeable future, and, according to the data offered by modern methods of remote sensing, it is considered unnecessary. By qualitatively “reprocessing” the existing maps, it is possible to obtain very accurate data. In [45], it was concluded that K-means could be used in mapping lithology and affirmed that this method can be implemented in existing map enhancement. The results of the RF classification showed that, in all cases, the anomalies were dissociated in a similar way. It is very likely that this was applied to surface objects with more significant differences in composition compared to rocks or vegetation. The next result was that the lithological unit marked as “K2,3, rudist limestones, SRL” was also not dissociated, nor was its GB with the lithological unit marked “E1,2 foraminiferal limestones, FL”. Although they were formed under different temporal and spatial conditions, due to the very similar percentage of CaCO3, they are very difficult to dissociate. In the fourth case (Figure 11d), RF did not dissociate green pixels, as it did in the first three cases shown in Figure 11a–c. In this case, the map shown in Figure 11d displays data very similar to the actual field situation compared with the previous three listed. Regarding RF, Figure 11d is considered the most representative for the area in terms of further exploration of mineral resources. If alternatives with colluvium were excluded, other geological relations were very satisfactory in this research stage. They could be used as a basis for the potential investigation of mineral resources. The effectiveness of ML techniques for classifying geological units based on multispectral maps was written about in [55]. With ML processing methods, existing geological boundaries can be corrected as basic elements of geological maps, which means that, in this way, geological maps can also be completed and can be used in delineating areas favorable for mineral exploration and exploitation through the creation of adequate mineral resources maps with accuracy that can be used in regional spatial plans at a scale of 1:25,000.
The raster-derived data based on the USGS limestone index agreed well with the field data. The reason is very likely due to the high content of CaCO3 in rudist (SRL), foraminiferal (FL) and platy limestone (PL). In this case, the NDVI was used as an additional tool for the general separation of vegetation from other land cover types (mostly rocks). Being obtained by a reliable sensor, it was considered useful for this type of use [56]. In DFA, the classification efficiency has often demonstrated higher efficiency in discrimination than various statistical tests [20]. In this regard, classification rates of previously outlined lithostratigraphic units can be straightforwardly examined, comparing the mathematically “predicted” (calculated) and original (“natural”) group membership of single cases. Table 4 reveals that the DFM was structured with very high efficacy exceeding 80%. However, for some groups, such as CL and PL, overlapping was considerable, which was also evident from their low squared Mahalanobis distance (3.33, Table 5) and, visually, from the respective group centroid diagrams (Figure 13b). In the abovementioned case, barely 22 of the 40 CL members were allocated to that specific group, whereas the remainder was lost to PL (11), SLS (6) and CL groups, collapsing the classification rate to a modest 55%. Given the PL group, the classification accuracy was slightly better and amounted to 72.5%. Here, 11 cases among the 40 were lost to other groups, including CO (6) and CL (5). This information indicates significant inhomogeneity of the CL lithological unit, the 45% allocations of which were unequivocally misclassified, showing remarkable cross-exchange with the PL group in particular (Table 4). Other groups were characterized by high classification rates, especially CO and SRL (95%).
In this paper, the authors tried to single out not only the types of limestone (Senonian rudist limestones (SRL) and Eocene foraminiferal limestones (FL), first type of discrimination) but also the sediments that are result of their weathering. The second type of discrimination refers to the “Promina Beds” sediments, which are known in the literature on karst terrains [19,57,58,59,60]. In this research, these sediments were represented by platy limestones (PL) and mainly carbonate conglomerates and limestones (CO). They were transported by strong currents which flowed from the raised Dinarides during the Eocene and Oligocene, eroded the carbonate platform sediments and deposited the weathered material in the land or marine sediment basins. The literature on the “Promina Beds” states that these sediments formed due to rapid changes in the river morphology, sudden fluctuations in flow and lack of stable banks. They also formed in marine conditions. This makes the lithological units PL and CO very similar, and it is often difficult to distinguish them in the field, as such sedimentary environments are extremely complex. The third type of lithological unit that has been singled out is Quaternary sediments (colluvium, gravel and sand (CL) and lake and swamp sediments, fine-grained sands and clays, SLS). They are also the weathering products of all previous lithological units and are usually additionally covered with a thin layer of soil. Therefore, they are also difficult to single out.
Using USGS spectral signatures and NDVI, the areas from which spectral signatures could be read were first selected, and these signatures were statistically processed. The obtained data indicated the possibility of distinguishing lithological members even within the complex karst terrain. This method is applicable in two ways. The first one refers to the validation of maps to the scale 1:25,000, while the second one refers to detailed geological mapping to the scale 1:10,000. This makes it possible to extract the potential of mineral resources in even more detail. DFA as a statistical method was also used as an additional tool for the validation of results, but also for proving the possibility of detailed separation of lithological members based on spectral signatures.

5. Conclusions

This paper emphasized that the rapid development of digital technologies enables research of mineral resources in karst terrains. Using the K-means algorithm, a lithological unit that was not dissociated in the original geological map was displayed at an early stage of the research. Its dissociation is important because its subsequently calculated quantities can significantly affect the geological potential of the study area. The K-means algorithm additionally showed its importance in outlier extraction. By processing the Sentinel-2A data, different lithological units with a carbonate composition were also extracted, but not without field validation and/or in combination with one of the supervised classification methods. In this case, the existing geological maps were finely, digitally processed, and the K-means algorithm with the DOF and field validation enabled their geometric correction. The GBs derived by ML algorithms were more precisely determined than before and can be used with greater precision when calculating geological potential. The supervised classification using the RF algorithm showed great potential in exploring mineral resources because it successfully divided the research area into given classes, although it was found to be mainly carbonate terrain. The supervised classification on raster datasets using USGS indices met its expectations exclusively as a support instrument, not for the absolute separation of lithological units. The NDVI also proved its value as a support tool. Finally, DFA showed its immeasurable potential in separating lithological units scanned by remote sensing and treated in the form of multispectral satellite data.
The use of ML algorithms made it possible to correct the geological boundaries by one order of magnitude, that is, the determined error in the geometry of the geological boundaries was reduced from the established 300 m to below 100 m. This is very significant for maps to the scale 1:25,000, since these data are directly used during spatial planning. In karst terrains without a significant percentage of vegetation, it is also possible to use spectral signatures. By integrating them into geological maps, the data can be displayed on an even more detailed scale, even to the scale of 1:10,000. This has not been the case so far. More precisely obtained geological data later result in more precise spatial plans, which is very important for other territory users but also for economic activities in a certain administrative unit. In addition to the above, the authors highlight the topic’s importance for the wider region, such as the Alpine–Dinaric karst area, Southeast Europe, etc.

Author Contributions

Conceptualization, N.G.; Formal analysis, N.G.; Funding acquisition, N.I.; Investigation, N.G. and M.G.; Methodology, M.G.; Project administration, N.G.; Software, M.G.; Supervision, M.G. and S.M.; Validation, S.M.; Visualization, B.L.-O.; Writing—original draft, N.G.; Writing—review and editing, M.G., S.M., B.L.-O., N.I. and Z.P. All authors have read and agreed to the published version of the manuscript.

Funding

This work was supported by the Croatian Geological Survey, Ministry of Science and Education of the Republic of Croatia, EIT Raw Materials (European Institute of Innovation and Technology) Horizon 2020 project RESEERVE—Mineral potential of the ESEE region (project no. 17029) and CSA Horizon 2020 project GeoTwinn (project no. 809943).

Data Availability Statement

Not applicable.

Acknowledgments

Thanks to Marko Copić, Tea Fluksi and Erli Kovačević Galović for fieldwork and language assistance. Support for this study was provided by the Croatian Geological Survey and the Ministry of Science and Education of the Republic of Croatia.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Gupta, R.P. Remote Sensing Geology, 3rd ed.; Springer: Berlin/Heidelberg, Germany, 2018; pp. 1–428. [Google Scholar] [CrossRef]
  2. Bishop, C.; Rivard, B.; de Souza Filho, C.; van der Meer, F. Geological Remote Sensing. Int. J. Appl. Earth Obs. Geoinf. 2018, 64, 267–274. [Google Scholar] [CrossRef]
  3. Oluić, M. Tektonska analiza graničnog područja Hrvatske i Slovenije izrađena na snimcima napravljenim iz satelita ERTS-1. Geol. Vje. 1975, 28, 87–96. [Google Scholar]
  4. Oluić, M. Neke karakteristike tektonike u graničnom području Crne Gore, Dalmacije i Hercegovine na osnovi Landsat snimaka. Geol. Vje. 1979, 24–25, 43–49. [Google Scholar]
  5. Oluić, M. Snimanje i Istraživanje Zemlje iz Svemira; Hrvatska akademija znanosti i umjetnosti: Zagreb, Republika Hrvatska, 2001; pp. 10–516. [Google Scholar]
  6. Tangestani, M.H.; Shayeganpour, S. Mapping a Lithologically Complex Terrain Using Sentinel-2A Data: A Case Study of Suriyan Area, Southwestern Iran. Taylor Fr. 2020, 41, 3558–3574. [Google Scholar] [CrossRef]
  7. Karimzadeh, S.; Tangestani, M.H. Potential of Sentinel-2 MSI Data in Targeting Rare Earth Element (Nd3+) Bearing Minerals in Esfordi Phosphate Deposit, Iran. Egypt. J. Remote Sens. Space Sci. 2022, 25, 697–710. [Google Scholar] [CrossRef]
  8. Ibrahim, E.; Barnabé, P.; Ramanaidou, E.; Pirard, E. Mapping Mineral Chemistry of a Lateritic Outcrop in New Caledonia through Generalized Regression Using Sentinel-2 and Field Reflectance Spectra. Int. J. Appl. Earth Obs. Geoinf. 2018, 73, 653–665. [Google Scholar] [CrossRef]
  9. Purwadi, I.; van der Werff, H.M.; Lievens, C. Targeting rare earth element bearing mine tailings on Bangka Island, Indonesia, with Sentinel-2 MSI. Int. J. Appl. Earth Obs. Geoinf. 2020, 88, 102055. [Google Scholar] [CrossRef]
  10. Shebl, A.; Abdellatif, M.; Hissen, M.; Abdelaziz, I.A.; Csámer, Á. Lithological mapping enhancement by integrating Sentinel 2 and gamma-ray data utilizing support vector machine: A case study from Egypt. Int. J. Appl. Earth Obs. Geoinf. 2021, 105, 102619. [Google Scholar] [CrossRef]
  11. Qi, X.; Zhang, C.; Wang, K. Comparing remote sensing methods for monitoring karst rocky desertification at sub-pixel scales in a highly heterogeneous karst region. Sci. Rep. 2019, 9, 13368. [Google Scholar] [CrossRef] [Green Version]
  12. Horvat, B.; Rubinić, J. Evaluating the Applicability of Thermal Infrared Remote Sensing in Estimating Water Potential of the Karst Aquifer: A Case Study in North Adriatic, Croatia. Remote Sens. 2021, 13, 3737. [Google Scholar] [CrossRef]
  13. Van der Meer, F.D.; van der Werff, H.M.A.; van Ruitenbeek, F.J.A. Potential of ESA’s Sentinel-2 for Geological Applications. Remote Sens. Environ. 2014, 148, 124–133. [Google Scholar] [CrossRef]
  14. Alpaydin, E. Introduction to Machine Learning, 4th ed.; The MIT Press: Cambridge, MA, USA; London, UK, 2020; p. 712. [Google Scholar]
  15. Cracknell, M.J. Machine Learning for Geological Mapping: Algorithms and Applications. Ph.D. Thesis, University of Tasmania, Hobart, Australia, 2014; pp. 1–275. [Google Scholar]
  16. Cracknell, M.J.; Reading, A.M. Geological Mapping Using Remote Sensing Data: A Comparison of Five Machine Learning Algorithms, Their Response to Variations in the Spatial Distribution of Training Data and the Use of Explicit Spatial Information. Comput. Geosci. 2014, 63, 22–33. [Google Scholar] [CrossRef] [Green Version]
  17. Cracknell, M.J.; Reading, A.M.; McNeill, A.W. Mapping Geology and Volcanic-Hosted Massive Sulfide Alteration in the Hellyer–Mt Charter Region, Tasmania, Using Random ForestsTM and Self-Organising Maps. Aust. J. Earth Sci. 2014, 61, 287–304. [Google Scholar] [CrossRef]
  18. Vlahović, I.; Tišljar, J.; Velić, I.; Matičec, D. Evolution of the Adriatic Carbonate Platform: Palaeogeography, Main Events and Depositional Dynamics. Palaeogeogr. Palaeoclimatol. Palaeoecol. 2005, 220, 333–360. [Google Scholar] [CrossRef]
  19. Zupanič, J.; Babić, L. Sedimentary Evolution of an Inner Foreland Basin Margin: Palaeogene Promina Beds of the Type Area, Mt. Promina (Dinarides, Croatia). Geol. Croat. 2011, 64, 101–119. [Google Scholar] [CrossRef]
  20. Rock, N.M.S. Numerical Geology. A source guide, glossary and selective bibliography to geological uses of computers and statistics. In Lecture Notes in Earth Sciences; Bhattacharji, S., Friedman, G.M., Neugebauer, H.J., Seilacher, A., Eds.; Springer: Berlin/Heidelberg, Germany, 1988; Volume 18, pp. 1–379. [Google Scholar]
  21. Rodriguez-Galiano, V.; Sanchez-Castillo, M.; Chica-Olmo, M.; Chica-Rivas, M. Machine Learning Predictive Models for Mineral Prospectivity: An Evaluation of Neural Networks, Random Forest, Regression Trees and Support Vector Machines. Ore Geol. Rev. 2015, 71, 804–818. [Google Scholar] [CrossRef]
  22. Lary, D.J.; Alavi, A.H.; Gandomi, A.H.; Walker, A.L. Machine Learning in Geosciences and Remote Sensing. Geosci. Front. 2016, 7, 3–10. [Google Scholar] [CrossRef] [Green Version]
  23. Othman, A.A.; Gloaguen, R. Integration of Spectral, Spatial and Morphometric Data into Lithological Mapping: A Comparison of Different Machine Learning Algorithms in the Kurdistan Region, NE Iraq. J. Asian Earth Sci. 2017, 146, 90–102. [Google Scholar] [CrossRef]
  24. Talukdar, S.; Singha, P.; Mahato, S.; Pal, S.; Liou, Y.A.; Rahman, A. Land-Use Land-Cover Classification by Machine Learning Classifiers for Satellite Observations-A Review. Remote Sens. 2020, 12, 1135. [Google Scholar] [CrossRef] [Green Version]
  25. Shirmard, H.; Farahbakhsh, E.; Heidari, E.; Pour, A.B.; Pradhan, B.; Müller, D.; Chandra, R. A Comparative Study of Convolutional Neural Networks and Conventional Machine Learning Models for Lithological Mapping Using Remote Sensing Data. Remote Sens. 2022, 14, 819. [Google Scholar] [CrossRef]
  26. Mamužić, P.; Korolija, B.; Majcen, Ž.; Borović, I.; Magaš, N.; Bojanić, L.; Božićević, S.; Babić, L.; Šimunić, A. Osnovna Geološka Karta SFRJ 1:100000, List Šibenik K 33-8. Inst. Geol. Istrž. Zagreb (1962-1965); Savezni geološki zavod: Beograd, Serbia, SFRJ, 1971. [Google Scholar]
  27. Mamužić, P.; Vrsalović, I.; Muldini-Mamužić, S.; Korolija, B.; Majcen, Ž.; Borović, I. Osnovna Geološka Karta SFRJ 1:100000: Tumač za List Šibenik K 33-8. Inst. Geol. Istrž. Zagreb (1966); Savezni geološki zavod: Beograd, Serbia, SFRJ, 1975. [Google Scholar]
  28. Marković, S. Hrvatske Mineralne Sirovine; Institut za geološka istraživanja: Zagreb, Croatia, 2002; pp. 1–544. [Google Scholar]
  29. Korbar, T. Orogenic Evolution of the External Dinarides in the NE Adriatic Region: A Model Constrained by Tectonostratigraphy of Upper Cretaceous to Paleogene Carbonates. Earth-Sci. Rev. 2009, 96, 296–312. [Google Scholar] [CrossRef]
  30. Kokaly, R.F.; Clark, R.N.; Swayze, G.A.; Livo, K.E.; Hoefen, T.M.; Pearson, N.C.; Wise, R.A.; Benzel, W.M.; Lowers, H.A.; Driscoll, R.L.; et al. USGS Spectral Library Version 7. U.S. Geol. Surv. Data Ser. 2017, 1035, 61. [Google Scholar] [CrossRef] [Green Version]
  31. Gascon, F.; Bouzinac, C.; Thépaut, O.; Jung, M.; Francesconi, B.; Louis, J.; Lonjou, V.; Lafrance, B.; Massera, S.; Gaudel-Vacaresse, A.; et al. Copernicus Sentinel-2A Calibration and Products Validation Status. Remote Sens. 2017, 9, 584. [Google Scholar] [CrossRef] [Green Version]
  32. Press, W.H.; Teukolsky, S.A.; Vetterling, W.T.; Flannery, B.P. Numerical Recipes in C: The Art of Scientific Computing, 2nd ed.; Cambridge University Press: New York, NY, USA, 1992; pp. 1–965. [Google Scholar]
  33. Hartigan, J.A. Clustering Algorithms; John Wiley & Sons: New York, NY, USA, 1975; pp. 1–347. [Google Scholar]
  34. Hartigan, J.A.; Wong, M.A. Algorithm AS 136: A k-means clustering algorithm. J. R. Stat. Soc. Ser. C-Appl. Stat. 1979, 28, 100–108. [Google Scholar] [CrossRef]
  35. Liberti, L.; Lavor, C. Euclidean Distance Geometry; Springe: Cham, Switzerland, 2017; 133p. [Google Scholar] [CrossRef]
  36. Maor, E. The Pythagorean Theorem: A 4000-Year History; Princeton University Press: Princeton, NJ, USA, 2019; 296p. [Google Scholar] [CrossRef]
  37. Jain, A.K. Data Clustering: 50 Years beyond K-Means. Pattern Recognit. Lett. 2010, 31, 651–666. [Google Scholar] [CrossRef]
  38. Richards, J.A. Remote Sensing Digital Image Analysis: An Introduction, 5th ed.; Springer: Berlin/Heidelberg, Germany, 2013; pp. 1–494. [Google Scholar] [CrossRef]
  39. Celik, T. Unsupervised Change Detection in Satellite Images Using Principal Component Analysis and κ-Means Clustering. IEEE Geosci. Remote Sens. Lett. 2009, 6, 772–776. [Google Scholar] [CrossRef]
  40. Vattani, A. K-Means Requires Exponentially Many Iterations Even in the Plane. In Proceedings of the twenty-fifth annual symposium on Computational geometry, Aarhus, Denmark, 8–10 June 2009; pp. 324–332. [Google Scholar] [CrossRef] [Green Version]
  41. Di Giuseppe, M.G.; Troiano, A.; Troise, C.; De Natale, G. K-Means Clustering as Tool for Multivariate Geophysical Data Analysis. An Application to Shallow Fault Zone Imaging. J. Appl. Geophys. 2014, 101, 108–115. [Google Scholar] [CrossRef]
  42. Gašparović, M.; Zrinjski, M.; Gudelj, M. Automatic Cost-Effective Method for Land Cover Classification (ALCC). Comput. Environ. Urban Syst. 2019, 76, 1–10. [Google Scholar] [CrossRef]
  43. Breiman, L. Random Forests. Mach. Learn. 2001, 45, 5–32. [Google Scholar] [CrossRef] [Green Version]
  44. Cracknell, M.J.; Reading, A.M. The Upside of Uncertainty: Identification of Lithology Contact Zones from Airborne Geophysics and Satellite Data Using Random Forests and Support Vector Machines. Geophysics 2013, 78, WB113–WB126. [Google Scholar] [CrossRef]
  45. Kuhn, S.; Cracknell, M.J.; Reading, A.M. Lithological Mapping in the Central African Copper Belt Using Random Forests and Clustering: Strategies for Optimised Results. Ore Geol. Rev. 2019, 112, 103015. [Google Scholar] [CrossRef]
  46. Harris, J.R.; Grunsky, E.C. Predictive Lithological Mapping of Canada’s North Using Random Forest Classification Applied to Geophysical and Geochemical Data. Comput. Geosci. 2015, 80, 9–25. [Google Scholar] [CrossRef]
  47. Darijani, M.; Farquharson, C.G.; Perrouty, S. A Random Forest approach to predict geology from geophysics in the Pontiac subprovince, Canada. Can. J. Earth Sci. 2022, 59, 489–503. [Google Scholar] [CrossRef]
  48. Shirmard, H.; Farahbakhsh, E.; Müller, R.D.; Chandra, R. A Review of Machine Learning in Processing Remote Sensing Data for Mineral Exploration. Remote Sens. Environ. 2022, 268, 112750. [Google Scholar] [CrossRef]
  49. Weier, J.; Herring, D. Measuring Vegetation (NDVI and EVI); NASA Earth Observatory: Washington, DC, USA, 2000. [Google Scholar]
  50. Davis, J.C. Statistics and Data Analysis in Geology, 2nd ed.; John Wiley & Sons: New York, NY, USA, 1986. [Google Scholar]
  51. Dillon, W.; Goldstein, M. Multivariate Analysis: Methods and Applications; John Wiley & Sons: New York, NY, USA, 1984. [Google Scholar] [CrossRef]
  52. Peh, Z.; Kovačević Galović, E. Geochemistry of Istrian Lower Palaeogene Bauxites—Is It Relevant to the Extent of Subaerial Exposure during Cretaceous Times? Ore Geol. Rev. 2014, 63, 296–306. [Google Scholar] [CrossRef]
  53. Hasan, O.; Miko, S.; Ilijanić, N.; Brunović, D.; Dedić, Ž.; Šparica Miko, M.; Peh, Z. Discrimination of Topsoil Environments in a Karst Landscape: An Outcome of a Geochemical Mapping Campaign. Geochem. Trans. 2020, 21, 1. [Google Scholar] [CrossRef] [PubMed]
  54. Plummer, C.; Carlson, D.; Hammersley, L. Physical Geology, 14th ed.; McGraw-Hill Education: New York, NY, USA, 2012. [Google Scholar]
  55. Kovačevič, M.; Bajat, B.; Trivič, B.; Pavlovič, R. Geological Units Classification of Multispectral Images by Using Support Vector Machines. In Proceedings of the International Conference on Intelligent Networking and Collaborative Systems, Barcelona, Spain, 4–6 November 2009; pp. 267–272. [Google Scholar] [CrossRef]
  56. Huang, S.; Tang, L.; Hupy, J.P.; Wang, Y.; Shao, G. A Commentary Review on the Use of Normalized Difference Vegetation Index (NDVI) in the Era of Popular Remote Sensing. J. For. Res. 2021, 32, 1–6. [Google Scholar] [CrossRef]
  57. Mrinjek, E. Conglomerate fabric and paleocurrent measurement in the braided fluvial system of the Promina Beds in northern Dalmatia (Croatia). Geol. Croat. 1993, 46, 125–136. [Google Scholar] [CrossRef]
  58. Mrinjek, E. Sedimentology and depositional setting of alluvial Promina Beds in northern Dalmatia, Croatia. Geol. Croat. 1993, 46, 243–261. [Google Scholar] [CrossRef]
  59. Babić, L.; Zupanić, J.; Kurtanjek, D. Sharply-topped alluvial gravel sheets in the Palaeogene Promina Basin (Dinarides, Croatia). Geol. Croat. 1995, 48, 33–48. [Google Scholar] [CrossRef]
  60. Babić, L.; Zupanič, J. Major events and stages in the sedimentary evolution of the Paleogene Promina basin (Dinarides, Croatia). Nat. Croat. Period. Musei Hist. Nat. Croat. 2007, 16, 215–232. [Google Scholar]
Figure 1. (a) Topographic map of the wider region of the eastern Adriatic coast with basic tectonic elements, modified according to [29]; (b) digital orthophoto (DOF) of the study area. White lines are geological boundaries (GBs) according to Figure 1c. GBs represented by black lines are derived by K-means and the DOF. Black dots are input points for RF classification. Magenta rectangle represents the Subset 1 site. Yellow rectangle represents the Subset 2 site; (c) field geological map of the OGK SFRJ, sheet Bribirske Mostine. Red rectangles in Figure 1b,c are areas in which satellite images were processed.
Figure 1. (a) Topographic map of the wider region of the eastern Adriatic coast with basic tectonic elements, modified according to [29]; (b) digital orthophoto (DOF) of the study area. White lines are geological boundaries (GBs) according to Figure 1c. GBs represented by black lines are derived by K-means and the DOF. Black dots are input points for RF classification. Magenta rectangle represents the Subset 1 site. Yellow rectangle represents the Subset 2 site; (c) field geological map of the OGK SFRJ, sheet Bribirske Mostine. Red rectangles in Figure 1b,c are areas in which satellite images were processed.
Remotesensing 14 05169 g001
Figure 2. (a) Lithological map of Subset 1; (b) lithological map of Subset 2. Black lines on both images are GBs derived from field mapping during December 2021. Red lines are GBs derived by ML algorithms. White lines are GBs from the original map. Differences in GB geometry were measured, and, at these sites, it amounted to between 30 m and 300 m.
Figure 2. (a) Lithological map of Subset 1; (b) lithological map of Subset 2. Black lines on both images are GBs derived from field mapping during December 2021. Red lines are GBs derived by ML algorithms. White lines are GBs from the original map. Differences in GB geometry were measured, and, at these sites, it amounted to between 30 m and 300 m.
Remotesensing 14 05169 g002
Figure 3. Geological mapping with DOF analysis, USGS data processing, RGB and NDVI flowchart.
Figure 3. Geological mapping with DOF analysis, USGS data processing, RGB and NDVI flowchart.
Remotesensing 14 05169 g003
Figure 4. Selection points (pins) for spectral signature reading regarding Subset 1 site; (a) DOF is in the background. Green area represents vegetation, and orange area represents rocks and other landforms, all based on NDVI analysis; (b) RGB composite based on MSI S-2A imagery; (c) raster dataset derived from USGS spectral index on limestones.
Figure 4. Selection points (pins) for spectral signature reading regarding Subset 1 site; (a) DOF is in the background. Green area represents vegetation, and orange area represents rocks and other landforms, all based on NDVI analysis; (b) RGB composite based on MSI S-2A imagery; (c) raster dataset derived from USGS spectral index on limestones.
Remotesensing 14 05169 g004
Figure 5. Selection points (pins) for spectral signature reading regarding Subset 2 site; (a) DOF is in the background. Green area represents vegetation, and orange area represents rocks and other landforms, all based on NDVI analysis; (b) RGB composite based on Sentinel 2-A imagery; (c) raster dataset derived from USGS spectral index on limestones.
Figure 5. Selection points (pins) for spectral signature reading regarding Subset 2 site; (a) DOF is in the background. Green area represents vegetation, and orange area represents rocks and other landforms, all based on NDVI analysis; (b) RGB composite based on Sentinel 2-A imagery; (c) raster dataset derived from USGS spectral index on limestones.
Remotesensing 14 05169 g005
Figure 6. Six lithological unit spectral signatures which were read according to flowchart shown in Figure 3.
Figure 6. Six lithological unit spectral signatures which were read according to flowchart shown in Figure 3.
Remotesensing 14 05169 g006aRemotesensing 14 05169 g006b
Figure 7. (a) MAX; (b) MIN; (c) average spectral signature values of each lithological unit.
Figure 7. (a) MAX; (b) MIN; (c) average spectral signature values of each lithological unit.
Remotesensing 14 05169 g007
Figure 8. K-means classification results. All colors in each of the figures from a to h refer to a particular number of classes obtained by the K-means algorithm; (a) utilizes five classes and bands 2, 3, 4, 8 and 11; (b) utilizes five classes and all bands; (c) utilizes 10 classes and bands 2, 3, 4, 8 and 11; (d) utilizes 10 classes and all bands; (e) utilizes 15 classes and bands 2, 3, 4, 8 and 11; (f) utilizes 15 classes and all bands; (g) utilizes 20 classes and bands 2, 3, 4, 8 and 11; (h) utilizes 20 classes and all bands.
Figure 8. K-means classification results. All colors in each of the figures from a to h refer to a particular number of classes obtained by the K-means algorithm; (a) utilizes five classes and bands 2, 3, 4, 8 and 11; (b) utilizes five classes and all bands; (c) utilizes 10 classes and bands 2, 3, 4, 8 and 11; (d) utilizes 10 classes and all bands; (e) utilizes 15 classes and bands 2, 3, 4, 8 and 11; (f) utilizes 15 classes and all bands; (g) utilizes 20 classes and bands 2, 3, 4, 8 and 11; (h) utilizes 20 classes and all bands.
Remotesensing 14 05169 g008aRemotesensing 14 05169 g008b
Figure 9. (a) K-means and DOF-derived geological map of Bribirske Mostine with all LODs; (b) clip from field geological map of the OGK SFRJ, sheet Bribirske Mostine (Figure 1c); (c) clip from K-means and DOF-derived map (Figure 9a); (d) ground truth data—shallow soil with platy limestone fragments is shown.
Figure 9. (a) K-means and DOF-derived geological map of Bribirske Mostine with all LODs; (b) clip from field geological map of the OGK SFRJ, sheet Bribirske Mostine (Figure 1c); (c) clip from K-means and DOF-derived map (Figure 9a); (d) ground truth data—shallow soil with platy limestone fragments is shown.
Remotesensing 14 05169 g009
Figure 10. Four derived LODs.
Figure 10. Four derived LODs.
Remotesensing 14 05169 g010
Figure 11. RF classification results; (a) RF variant 1 and bands 2, 3, 4, 8 and 11; (b) RF variant 1 and all bands; (c) RF variant 2 and Bands 2, 3, 4, 8 and 11; (d) RF variant 2 and all bands.
Figure 11. RF classification results; (a) RF variant 1 and bands 2, 3, 4, 8 and 11; (b) RF variant 1 and all bands; (c) RF variant 2 and Bands 2, 3, 4, 8 and 11; (d) RF variant 2 and all bands.
Remotesensing 14 05169 g011
Figure 12. Grouped boxplots for each band (B1B12) and lithological group (SRL, FL, CO, PL, CL, SLS).
Figure 12. Grouped boxplots for each band (B1B12) and lithological group (SRL, FL, CO, PL, CL, SLS).
Remotesensing 14 05169 g012
Figure 13. (a) Scatterplot of variable loadings, bands B1 to B12; (b) scatterplot of group means of six lithological units (SRL, FL, CO, PL, CL, SLS).
Figure 13. (a) Scatterplot of variable loadings, bands B1 to B12; (b) scatterplot of group means of six lithological units (SRL, FL, CO, PL, CL, SLS).
Remotesensing 14 05169 g013
Table 1. RF input polygon data for each lithological member.
Table 1. RF input polygon data for each lithological member.
No.LabelNameArea in km2Surface Percentage %No. of PolygonsPolygon ShapePolygon Diameter (m)
1j
SLS
Lake and swamp sediments, fine-grained sands and clays14.2615.25793Circle20
2d
CL
Colluvium, gravel and sand27.2429.171516Circle20
26.6328.511406
3aE3
CO
Mainly carbonate conglomerates and limestones20.9422.421.166Circle20
4bE3
PL
Platy limestones12.9113.82719Circle20
13.5214.47819
5E2,3
LMC
Limestones, marls and conglomerates in alternation16.3717.53911Circle20
6E1,2
FL
Foraminiferal limestones1.241.3369Circle20
7K2,3
SRL
Rudist limestones0.250.2714Circle20
8aAnomalies0.200.2212Circle20
9 Total area93.39
10Total percentage100
11Total polygons 5200
Table 2. Multivariate test regarding general significance of discrimination including tests of residual roots (DFs).
Table 2. Multivariate test regarding general significance of discrimination including tests of residual roots (DFs).
No. of variables12
Wilks’ lambda0.0024
Approximate F-ratio45.979
Degrees of freedom[60; 1048]
p-levelp = 0.000
DFEigen
value
Eigen
(%)
Eigen
cum.
Canon. RWilks’
λ
Chi2dfp-Level
128.37487.3687.360.9830.0021389.2600.000
21.8545.7193.070.8060.070611.8440.000
31.1013.3996.460.7240.200370.6300.000
40.8852.7299.180.6850.419199.8180.000
50.2650.82100.000.4580.79154.080.000
Table 3. Breakdown table of descriptive statistics for each band and lithological group.
Table 3. Breakdown table of descriptive statistics for each band and lithological group.
GroupBandNMeanSt. Dev.MinMaxQ1MedianQ3
CO 400.0490.0050.0410.0640.0450.0480.052
CL 400.0500.0050.0390.0660.0460.0500.053
FL 400.0980.0060.0850.1070.0940.0990.105
PLB1400.0530.0050.0450.0610.0500.0530.057
SLS 400.0420.0020.0390.0480.0400.0410.043
SRL 400.1070.0090.0860.1220.1020.1080.113
All Groups 2400.0660.0270.0390.1220.0460.0530.097
CO 400.0630.0100.0440.0940.0560.0640.069
CL 400.0670.0150.0450.1080.0570.0630.072
FL 400.1220.0100.1000.1470.1160.1230.129
PLB2400.0670.0110.0470.0920.0580.0660.075
SLS 400.0470.0060.0340.0610.0430.0460.049
SRL 400.1260.0150.0980.1550.1130.1250.138
All Groups 2400.0820.0330.0340.1550.0560.0680.115
CO 400.0840.0120.0600.1150.0760.0820.090
CL 400.0960.0190.0670.1630.0840.0910.103
FL 400.1420.0110.1170.1720.1360.1420.147
PLB3400.0950.0150.0690.1240.0840.0910.108
SLS 400.0680.0060.0580.0820.0630.0660.070
SRL 400.1490.0150.1180.1790.1370.1480.158
All Groups 2400.1060.0330.0580.1790.0780.0950.137
CO 400.0980.0160.0720.1340.0840.1000.109
CL 400.1210.0290.0730.2150.1020.1140.135
FL 400.1590.0120.1290.2010.1530.1580.166
PLB4400.1200.0200.0780.1620.1040.1180.136
SLS 400.0770.0090.0580.1030.0720.0750.079
SRL 400.1780.0190.1350.2220.1630.1780.191
All Groups 2400.1250.0390.0580.2220.0930.1180.158
CO 400.1280.0130.1080.1570.1150.1270.139
CL 400.1550.0260.1210.2460.1380.1490.164
FL 400.1860.0130.1490.2080.1810.1870.194
PLB5400.1610.0220.1220.2060.1430.1630.173
SLS 400.1220.0080.1080.1420.1160.1210.127
SRL 400.2120.0190.1640.2510.1960.2080.223
All Groups 2400.1600.0360.1080.2510.1290.1570.190
CO 400.1780.0160.1520.2270.1680.1730.188
CL 400.2250.0320.1750.3000.2000.2250.246
FL 400.2230.0150.1870.2540.2160.2240.234
PLB6400.2300.0360.1720.3270.2010.2330.247
SLS 400.2060.0260.1660.2760.1850.2040.219
SRL 400.2450.0170.1960.2810.2320.2460.257
All Groups 2400.2180.0330.1510.3270.1900.2200.241
CO 400.1980.0170.1730.2520.1860.1920.208
DE 400.2500.0360.1900.3260.2220.2500.280
FL 400.2460.0160.2060.2780.2390.2470.256
PLB7400.2520.0370.1850.3490.2200.2550.270
SLS 400.2340.0290.1900.3170.2110.2310.251
SRL 400.2680.0170.2190.3010.2550.2680.279
All Groups 2400.2410.0340.1730.3490.2120.2450.263
CO 400.2290.0230.1950.3210.2140.2250.240
CL 400.2760.0430.2080.3850.2440.2700.313
FL 400.2900.0180.2400.3340.2820.2920.299
PLB8400.2790.0440.2050.3930.2490.2750.304
SLS 400.2620.0340.2100.3630.2390.2540.287
SRL 400.3050.0170.2580.3470.2940.3040.317
All Groups 2400.2730.0390.1950.3470.2940.3040.317
CO 400.2300.0180.2040.2760.2190.2230.243
CL 400.2770.0370.2140.3580.2490.2800.304
FL 400.2980.0180.2520.3320.2900.3010.308
PLB8A400.2790.0370.2130.3760.2470.2820.297
SLS 400.2670.0290.2190.3510.2430.2660.284
SRL 400.3110.0170.2660.3510.2980.3110.320
All Groups 2400.2770.0370.2040.3760.2440.2820.305
CO 400.2240.0130.2050.2670.2160.2220.233
CL 400.2790.0310.2290.3410.2570.2830.295
FL 400.2860.0120.2540.3050.2830.2880.294
PLB9400.2800.0340.2230.3690.2500.2780.303
SLS 400.2690.0210.2360.3150.2520.2710.282
SRL 400.3030.0130.2680.3220.2960.3020.313
All Groups 2400.2740.0330.2050.3690.2450.2820.296
CO 400.2780.0260.2390.3390.2540.2740.303
CL 400.2990.0320.2470.3810.2740.3020.314
FL 400.4310.0280.3620.4990.4140.4320.450
PLB11400.3150.0280.2470.3620.2970.3240.336
SLS 400.2580.0130.2350.2910.2520.2570.267
SRL 400.4210.0200.3820.4590.4060.4210.436
All Groups 2400.3340.0720.2350.4990.2700.3130.407
CO 400.1910.0200.1600.2350.1750.1890.207
CL 400.2050.0340.1520.3120.1840.2030.212
FLB12400.3120.0220.2530.3660.2970.3110.325
PL 400.2130.0230.1680.2530.1990.2110.231
SLS 400.1580.0130.1380.1890.1510.1570.164
SRL 400.2990.0170.2710.3300.2850.2990.313
All Groups 2400.2300.0610.1380.3660.1780.2110.293
Note: St. Dev. = standard deviation; Min = minimum; Max = maximum; Q1 = lower quartile; Q3 = upper quartile; and Med = median.
Table 4. Classification matrix of six lithological units (CO, CL, FL, PL, SLS, SRL).
Table 4. Classification matrix of six lithological units (CO, CL, FL, PL, SLS, SRL).
GroupsCOCLFLPLSLSSRL%
Correct
CO380020095.00
CL1220116055.00
FL003600490.00
PL650290072.50
SLS300037092.50
SRL002003895.00
Total48273842434283.33
Observed classifications (rows)—samples passed on other groups; Predicted classifications (columns)—samples obtained from other groups; a priori defined groups have an equal number of 40 samples each.
Table 5. Squared Mahalanobis distances (SMD) between the investigated lithological groups (CO–SRL).
Table 5. Squared Mahalanobis distances (SMD) between the investigated lithological groups (CO–SRL).
COCLFLPLSLSSRL
CO0.0012.07116.2712.5712.49129.37
CL12.070.00130.903.339.57136.98
FL116.27130.900.00122.04142.8219.25
PL12.573.33122.040.0011.08128.68
SLS12.499.57142.8211.080.00152.59
SRL129.36136.9819.25128.68152.600.00
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Gizdavec, N.; Gašparović, M.; Miko, S.; Lužar-Oberiter, B.; Ilijanić, N.; Peh, Z. Discrimination of Rock Units in Karst Terrains Using Sentinel-2A Imagery. Remote Sens. 2022, 14, 5169. https://0-doi-org.brum.beds.ac.uk/10.3390/rs14205169

AMA Style

Gizdavec N, Gašparović M, Miko S, Lužar-Oberiter B, Ilijanić N, Peh Z. Discrimination of Rock Units in Karst Terrains Using Sentinel-2A Imagery. Remote Sensing. 2022; 14(20):5169. https://0-doi-org.brum.beds.ac.uk/10.3390/rs14205169

Chicago/Turabian Style

Gizdavec, Nikola, Mateo Gašparović, Slobodan Miko, Borna Lužar-Oberiter, Nikolina Ilijanić, and Zoran Peh. 2022. "Discrimination of Rock Units in Karst Terrains Using Sentinel-2A Imagery" Remote Sensing 14, no. 20: 5169. https://0-doi-org.brum.beds.ac.uk/10.3390/rs14205169

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop