A Fusion of Feature-Oriented Principal Components of Multispectral Data to Map Granite Exposures of Pakistan

Abstract

Despite low spatial resolutions, thermal infrared bands (TIRs) are generally more suitable for mineral mapping due to fundamental tones and high penetration in vegetated areas compared to shortwave infrared (SWIR) bands. However, the weak overtone combinations of SWIR bands for minerals can be compensated by fusing SWIR-bearing data (Sentinel-2 and Landsat-8) with other multispectral data containing fundamental tones from TIR bands. In this paper, marble in a granitic complex in Mardan District (Khyber Pakhtunkhwa) in Pakistan is discriminated by fusing feature-oriented principal component selection (FPCS) obtained from the ASTER, Landsat-8 Operational Land Imager (OLI), Thermal Infrared Sensor (TIRS) and Sentinel-2 MSI data. Cloud computing from Google Earth Engine (GEE) was used to apply FPCS before and after the decorrelation stretching of Landsat-8, ASTER, and Sentinel-2 MSI data containing five (5) bands in the Landsat-8 OLI and TIRS and six (6) bands each in the ASTER and Sentinel-2 MSI datasets, resulting in 34 components (i.e., 2 × 17 components). A weighted linear combination of selected three components was used to map granite and marble. The samples collected during field visits and petrographic analysis confirmed the remote sensing results by revealing the region’s precise contact and extent of marble and granite rock types. The experimental results reflected the theoretical advantages of the proposed approach compared with the conventional stacking of band data for PCA-based fusion. The proposed methodology was also applied to delineate granite deposits in Karoonjhar Mountains, Nagarparker (Sindh province) and the Kotah Dome, Malakand (Khyber Pakhtunkhwa Province) in Pakistan. The paper presents a cost-effective methodology by the fusion of FPCS components for granite/marble mapping during mineral resource estimation. The importance of SWIR-bearing components in fusion represents minor minerals present in granite that could be used to model the engineering properties of the rock mass.

1. Introduction

Conventional lithological mapping is time-consuming, costly, and limited due to the practicalities of hilly terrains. Remote sensing (RS) can be used as a cost and time-effective means for lithological mapping of mineral deposits [1,2,3]. The suitability of data for lithological mapping depends on the rock type in question, data availability, and spatial, spectral, and radiometric resolutions.

The Landsat-8 Operational Land Imager (OLI) consists of 6 visible and near-infrared (VNIR), two shortwave infrared (SWIR), and one panchromatic sensors with 30 m, 30 m, and 15 m resolutions, respectively. The Landsat-8 Thermal Infrared Sensor (TIRS) has two TIR bands of a 100 m spatial resolution available at the 30 m resolution in Google Earth Engine (GEE). In comparison, the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) has 14 bands, i.e., 3 VNIR bands with a 15 m resolution, 6 SWIR bands with a 30 m resolution, and 5 TIR bands with a 90 m resolution [4]. However, the radiometric response of SWIR bands from the ASTER is unreliable after 2008 due to a technical fault [5]. The Sentinel-2 Multispectral Instrument (MSI) has a relatively higher (20 m) spatial resolution for SWIR bands, but no TIR bands [6,7]; therefore, further studies are needed to explore its capabilities for mineral exploration [8]. The Sentinel-2 MSI has a resolution of 10–20 m, but it lacks TIR bands, whereas the Landsat-8 TIRS has TIR bands, but its resolution is 100 m. However, GEE can convert bands to the exact spatial resolution by resampling, e.g., 100 m TIR bands are available as 30 m by default in its database. The TIR emissivity for igneous rocks decreases with the high silica content, while a corresponding shift in absorption to a higher wavelength band is observed in the mineral groups due to types of dominant silicates [4]. The minimum TIR emissivity of granite at 9 µm is due to the vibration absorption of the Si–O silicate structure, which shifts to 11 µm in mafic rocks having a tectosilicate structure [4,9,10]. However, TIR satellite sensors have low temporal, spatial, and spectral resolutions, and the spectral responses in these bands may contain a mixture of features [11]. Vegetation or soil cover may cause high variation in granite spectral response, making it challenging to map the deposits [4,12]. High-spatial-resolution SWIR bands show spectral absorption for different granite compositions (e.g., muscovite and epidote) and calcite in 2.2 μm regions; a complete absorption is only observed for calcite near the 2.3 μm regions [4,10]. The fusion of TIR with other high-resolution SWIR datasets can provide insight into understanding different features [4,11,13]. Integrating results from multiple sources, e.g., the Landsat-8, the Sentinel-2, and the ASTER, by a suitable criterion would generally improve the outcome compared to single-source data [14].

Analytical techniques for RS are grouped into data-driven and knowledge-based techniques [15]. Several techniques, such as band ratios [16,17], band ratio matrix transformation (BRMT) [18], spectral angle mappers (SAM) [19,20], decorrelation stretching (DS) [17,21], principal components analysis (PCA) [17,22,23], directed principal component analysis (DPCA) [24], feature-oriented principal components selection (FPCS) also known as Crosta technique [25,26,27], spatially weighted PCA [28,29], directed PCA [30,31], and other variants of PCA [32], have been reported for mapping different lithological [33] and environmental [34] features. Recently, Support Vector Machine (SVM), PCA, and independent component analysis (ICA) have been used to map gold-bearing granite-greenstone rocks using AVIRIS-NG hyperspectral and ASTER data [35].

This research presents a fusion strategy of principal components (PCs), representing effective but low-resolution TIR bands loading with relatively less responsive but high-spatial-resolution Sentinel data. This research used Landsat-8, ASTER, and Sentinel-2 data to devise a strategy for developing a surface outcrop map at regional and local scales in Marble, granite-bearing vegetated regions in Shewa Shahbazghari Mardan District of Khyber Pakhtunkhwa, Pakistan. Such low-cost, fast and reliable mapping strategies can be used in any granite field and maybe combined further with geophysical data for three-dimensional (3D) resource quantification and modelling. The main objectives of this study are as follows: (i) to apply the feature-oriented principal components selection (FPCS) technique to decorrelation-stretched and raw Sentinel-2, ASTER, and Landsat-8 datasets for generating granite and marble outcrop maps; (ii) to select the three best PCs by the Crosta technique, supported by the following: (1) field validation; (2) known granite and marble mining leases; and (3) petrographic investigations of collected samples; (iii) to devise a weighting criterion for the fusion of the selected components to generate a false-color composite (FCC) image, differentiating between marble and granite in the region; (iv) to further extend the technique for mapping extents of granite in the southern (Thar District in the Province of Sindh) and northern (Malakand District of Khyber Pakhtunkhwa Province) Pakistan to establish further applicability.

2. Geology of Shewa Shahbazghari and Ambela Granite Complexes

The Shewa-Shahbaz Garhi complex is an isolated outcrop (longitude: 72.10′ E, latitude: 34.10′ to 34.30′ N) dominantly comprised of alkali granites with a subordinate amount of basic injections (Figure 1) [36,37]. This region also contains an outcrop of Nowshera formation that mainly contains sandy dolomite, limestone, and marble with calcareous quartzite and calcareous argillite. This formation has varying grain sizes and is rich in magnesium carbonates (MgCO₃). The other lithology reported from the complex includes riebeckite gneiss and porphyritic micro-granite [36]. These rocks constitute some of the significant components of the intra-plate magmatism known as the Peshawar Plain Alkaline Igneous Province (PPAIP) complex. The exposure of the PPAIP complex extends from the Tarbela region in the east towards the Pakistan–Afghanistan border (Loe Shilman area) in the west.

In addition to the Shewa Shahbazgarhi complex, the other exposures of the PPAIP complex include the Ambela granitic complex (AGC), the Warsak complex, and the Tarbela alkaline complex. The Shewa Shahbazgarhi complex’s granites are white to grey, ranging in grain size from medium to coarse [36]. Primarily, they consist of 50–55% alkali feldspar, 2–3% plagioclase, and 28–33% quartz with accessory minerals not exceeding 15% in total. Most abundant among the accessory is biotite (<8%), along with muscovite, sphene, zircon, garnet, and epidote [38]. Based on mineralogical composition, the basic dikes that intrude in these granites can be divided into two lithologies [37]: (a) dolerite, which is mostly fresh, consisting of plagioclase, clinopyroxene, orthopyroxene, and minor olivine; and (b) amphibolite, consisting of amphibole and plagioclase with a minor concentration of chlorite, epidote, and clays.

The Ambela granitic complex, exposed in the northeast of the Shewa complex, is divided into three groups based on their modal mineralogy, age, and field relations [39]. Group 1 includes granites, alkali granites, and micro-porphyries intruded by group 2 rocks containing quartz-syenites, syenites, carbonatites, and fenites. The basic injection constitutes group 3 rocks and cross-cuts of both group 1 and 2 rocks, representing the latest magmatic episode in the Ambela complex [39]. The granite of the Ambela complex has the same texture, color, and compositional characteristics as the granite of the Shewa-Shahbazgarhi complex.

The second study area of the Kotah Dome in Malakand is north of Shewa Shahbazgarhi, predominantly consisting of the auger and flaser granite and granodiorite gneiss intruded by tourmaline granite, meta gabbro, amphibolite, and pegmatite [40,41]. Quartz, feldspar, biotite, and muscovite are the main constituents, while accessories include garnet, chlorite, and amphibole in the igneous, metamorphic complex [41]. Thirdly, the granite of the Nagarparker area is in the 19-km-long and 305-m-high Karoonjhar Mountains, extending across from the Aravalli range of India. This granite is dominated by quartz, two feldspars (plagioclase and alkali) that occur as larger grains, and other minerals include biotite, iron oxide, chlorite, sphene, and epidote [42]. The ilmenite and Ca-plagioclase give pink to red and grey colors to the quartz diorite + diorite (pink) and monzogranite (grey granite), respectively [43].

3. Materials and Methods

3.1. PCA

PCA is a well-known technique that enhances artefact discriminations by accumulating all image data bands [44]. Eigenvector and eigenvalue pairs are derived from the eigendecomposition of the data’s covariance matrix to project it along the principal axes. Image data are multiplied with the eigenvector matrix to obtain PCs. The eigenvector indicates the amount of information retained from each band. A target surface is highlighted by dark or bright pixels in the relevant PC image. The PCs obtained from PCA contain the information from all bands in varying proportions. Thus, components retaining most spectral information of a specific lithological unit can be identified from eigen analysis by the Crosta technique, also known as FPCS [45]. The Crosta technique identifies spectral-information components by visualizing individual component images, indicating a particular mineral as bright/dark at a specified location [46]. This technique presents dark and bright pixels in a PC according to the magnitude and sign of the eigenvector loadings.

3.2. DS

The DS technique stretches a PC image to improve the contrast for different lithologies in the PC image and suppress the noise [47]. Similarly, a Gaussian stretch can be performed on data to analyze eigenpairs for interpreting PC images. In DS, the transformation matrix given in Equation (1) is derived using eigenvalues and eigenvector matrices obtained by the covariance matrix’s eigendecomposition. The transformation function can be written as:

$$T={\left(v\right)}^{t}s v,$$

(1)

where $t$ represents the transpose, $s$ is the rotation matrix, and $v$ is the eigenvector matrix obtained by the eigendecomposition of the bands’ covariance matrix. This transformation $T$ is needed to stretch the pixel data along the principal axes. Such that,

$$s=\frac{1}{\sqrt{\lambda }} .$$

(2)

where the parameter $\lambda$ is the eigenvalue matrix. The bands are stretched with the following mathematical operation:

$$D{S}^i=T · P^i.$$

(3)

where $P^i$ is the $i^{th}$ principal axis and $D{S}^i$ is the $i^{th}$ stretched band. This transformation matrix is used to stretch the data, i.e., applying “rotation” on each band along its principal axis, thus decreasing band-to-band correlation. This stretching of the bands’ pixel data improves the contrast between different target minerals. Although the technique significantly improves the contrast in the images, the information is still scattered across all bands, and at most, three bands can be visualized in an FCC. Stretching the results from FPCS results by the DS algorithm can improve its lithological discernibility, making visual interpretation easier.

3.3. Mapping in Dense Vegetation by the Fusion of FPCS Components from Different Data Sources Using GEE

Frequently used multispectral data sources, such as the ASTER, the Landsat-8 OLI, and the TIRS, and Sentinel-2 MSI collections, are publicly available as part of the Google Cloud public data program [48,49,50]. The satellite data are stored as pixels rather than images in GEE [50], which means that one can acquire data from the GEE data catalogue at the pixel level and not at the image level. Mosaicking can be performed at a pixel level, e.g., retaining the best pixel with the minimum cloud cover among the same image collected at different times. Furthermore, cloud storage allows web-based access for processing, applying various algorithms and analysis by utilizing Google distributed computing architecture with computing efficiency exponentially higher than that of standard desktop computers [49].

A pixel in RS images holds the mixed spectral response from the earthly objects present in the corresponding location. The mixture gets more complicated, when the spatial resolution is low. Due to inconsistent findings and contradictory observations in TIR imageries, it is difficult to interpret and understand TIR imageries in highly vegetated areas. Therefore, high-resolution SWIR bands perform better than lower-resolution TIR bands in a medium-to-low vegetation region [11]. Integrating multiple datasets can help compensate for the lower resolution of TIR and the shallow penetration of SWIR for granitic rock mapping [7,51]. Such an operation is generally known as a fusion [52], broadly categorized as pixel-level (pan-sharpening before processing), object-level (integrating features from pixels), or product-level fusion (after processing of individual data sources) [53]. The application of PCA and the product-level fusion of selective PCs in the lower-dimension space have been reported for computational efficiency [54]. Similar to the weighted FCCs of different band ratios [55], PCs can be assigned a unique weight and linearly combined to get weighted FCCs.

3.3.1. Pre-Processing and Selection of Suitable Multispectral Bands

First, filtered images with the least vegetation, noise, clouds, and shadows were acquired from the earth engine database within the ROI shown in Figure 2. This figure is the base map for evaluating the results that indicated the leased areas for marble/granite quarries and low, medium, and dense vegetation areas. Landsat-8 data were collected from the Level-1 collection 1 Tier 1 (precision and terrain corrected) suite of Landsat-8 OLI and TIRS images scaled and calibrated at the sensor level. The SWIR bands of the Sentinel-2 MSI L2A data with the bottom of atmosphere reflectance in cartographic geometry were atmospherically corrected using Sen2Cor processor and the PlanetDEM digital elevation model (DEM). True-color composites (TCCs) and FCCs were used for filtering data of any discontinuities, cloud, noise, vegetation, and shadows in the image. The median value for each pixel was obtained from all the images of these datasets over the past five years for selecting pixels with the minimum vegetation, clouds, and other variations in a pixel at a different time from Landsat-8 and Sentinel-2 MSI data.

GEE’s “Clip” function was used to crop the area of interest to avoid mixing the spectral response of granitic rock bodies with the surrounding residential areas and agricultural fields. This helped process only the desired location in the study region, which allowed faster computational time and better results.

Granite and marble bodies showed significant spectral variations in 0.9 µm, 1.4 µm, and 2.3 µm high-spatial-resolution NIR and SWIR regions [4,10]. Granite-associated minerals such as muscovite showed discriminative spectral signatures responses in the SWIR range between 1.400 µm and 2.200 µm [13]. On the other hand, alkaline rocks emit larger wavelength radiations, specifically the TIR range of 8–12 µm [56]. Therefore, the corresponding bands, i.e., 0.88 µm, 1.600 µm, and 2.430 µm, of the three datasets (the Landsat-8 OLI, the ASTER, and the Sentinel-2 MSI) were used to carry out further investigations (see

Table 1. Relevant near infrared (NIR) and shortwave infrared (SWIR) bands of each dataset.

Dataset		Bands	Spectral Range	Resolution	Acquisition Time
Landsat-8 OLI [59]	LANDSAT/LC08/C01/T1_RT	B5 (NIR)	0.82 to 0.88 µm	30 m	Mosaicked median-pixel data acquired from 2015 to 2020
		B6 (SWIR 1)	1.57 to 1.65 µm
		B7 (SWIR 2)	2.11 to 2.29 µm
Landsat-8 TIRS [59]		B10 (TIR 1)	10.6–11.19	100 m resampled to 30 m
		B11 (TIR 2)	11.50–12.51
ASTER [5]	ASTER/AST_L1T_003	B1 (Green)	0.520–0.600 µm	15 m	Mosaicked median-pixel data acquired from 2015 to2020
		B2 (Red)	0.630–0.690 µm
		B3N (NIR)	0.78 to 0.86 µm
		B10 (TIR 1)	8.125–8.475 µm	90 m
		B11 (TIR 2)	8.475–8.825 µm
		B12 (TIR 3)	8.925–9.275 µm
Sentinel-2 MSI [60]	COPERNICUS/S2_SR	B2 (Blue)	0.490 µm	10 m	Mosaicked median-pixel data acquired from 2015 to 2020
		B3 (Green)	0.560 µm
		B4 (Red)	0.665 µm
		B8 (NIR)	0.842 µm	20 m
		B11 (SWIR 1)	1.57 to 1.65 µm
		B12 (SWIR 2)	2.11 to 2.29 µm

).

Calcite (a main mineral constituent of marble) shows diagnostic spectral absorptions in the SWIR region of electromagnetic radiation. Since the SWIR bands of ASTER are no longer functional and have previously produced poor results [57]; hence, the SWIR bands from the Landsat-8 OLI and Sentinel-2 data were utilized. The 15 m VNIR (B01, B02, and B3N) and the 90 m TIR (B10, B11, and B12) bands of the ASTER, the NIR (B5), SWIR (B6 and B7) of the Landsat-8 OLI, and TIR (B10 and B11) bands of the Landsat-8 TIRS, and the visible (B2, B3, and B4), NIR (B8), and SWIR (B11 and B12) bands of the Sentinel-2 were selected. The red-green-blue (RGB) bands of the Landsat-8 OLI were not used due to their lower resolution. Based on the spectral response of granite (in SWIR and TIR bands) and marble (1.6 μm and 2.3 μm) [58], suitable bands for FPCS were identified, as shown in Table 1. The schematic representation of the detailed methodology is also shown in Figure 3.

The bands’ spatial resolution variation was automatically handled by nearest neighbor resampling while processing the scale provided (i.e., 15 m in this study). The TCCs from Landsat-8 OLI data had the least shadows, but a relatively greater vegetation cover than Sentinel-2 TCCs. The Sentinel-2 images were level 2A-corrected; therefore, they were rich in contrast, free of clouds and had less vegetation than the Landsat-8 OLI images. Since satellites have different revisit times, due to varying weather conditions, the images of Sentinel-2 reported higher vegetation than the Landsat-8 ones.

3.3.2. Application of PCA

All the bands were mean-centered, i.e., subtracting the corresponding band mean from the band values before deriving the covariance matrix and eigen analysis. Finally, FPCS was applied separately to these three datasets (ASTER, Sentinel-2, and Landsat-8) in the raw and stretched (using DS) form to obtain PC images. The results were obtained by applying FPCS to both the stretched and raw data; for all the three datasets, 34 PCs were generated, i.e., 10 from the Landsat-8 OLI and the Landsat-8 TIRS, 12 from the ASTER, and 12 from the Sentinel-2 data sources with the same spatial resolution, by using nearest neighbor resampling.

3.3.3. Fusion of the Selected Crosta Components

Both the raw and stretched data were used to incorporate more information, as a high degree of stretching means the loss of information for minor pixels [42]. The proposed criteria assigned an equal weight to all the datasets involved in the final FCC. However, the final FCC was supposed to contain two components from the same dataset, and then, the weight assignment was based on the information of a particular mineral in a component (i.e., granite and marble in this case). After the assignment of weights, components were linearly combined to derive a final FCC. The three best PCs, retaining significant granitic spectral features, were selected and linearly combined after assigning weights to the components. Each component image was used to validate the known locations of four granite leases and marble quarries of Shewa-Shahbaz Garhi and Ambela complexes. The known mining leases of granite in these regions were marked as ML-1, ML-2, ML-3, and ML-4, while the marble mining leases were marked marble mining leases (MMLs) in Figure 2. The criterion for the assignment of weight for each component is presented as:

$$FCC=w_{1}P{C}_d^i+w_{2}P{C}_d^j+w_{3}P{C}_d^k,$$

(4)

where $d$ can be any combination of the chosen components among the FPCS applied to raw or stretched Sentinel-2, ASTER, Landsat-8 OLI, and Landsat-8 TIRS data; $w_1, w_2, and w_3$ are the weights assigned to generate an FCC for the chosen components $P{C}_d^i$ $P{C}_d^j$, and $P{C}_d^k;$ $i, j$, and $k$ indicate the $PC$ of dataset $d$ subject to the following constraints:

$$w_1+w_2+w_3=3.$$

(5)

If the chosen components are from the different datasets, equal weights are assigned to each component, i.e., $w_1=w_2=w_3=1$. If any two components are from the same dataset, then the sum of weights of these two components is half, i.e., 1.5, and the third component is assigned as the remaining 1.5 as per Equation (5). This 1.5 weight is split among the two components from same dataset, such that the Crosta component with more explicit indications about the targeted mineral is given a weight between 2/3 and 3/4 while the other is given the remaining weight. This is conducted to incorporate more information from the component about the mineral/rock type. In this way, the weights can be distributed among components to produce the best results for different rock types and data.

The proposed methodology was also extended to map the well-known granitic rock bodies of Karoonjhar Mountains in Nagarparker (Sindh Province of Pakistan) and the Kotah Dome, Malakand (Khyber Pukhtunkhwa Province of Pakistan).

4. Results

The application of DS on the selected bands and the PCs of the three datasets reduced the correlation among the SWIR bands and the NIR bands. Resultantly, the PCs represented more contrasting features in the data, as seen from the scatter plots (shown in Figure 4) of PC1, PC2, and PC3 of the Landsat-8 OLI data. FPCS performed after DS showed relatively better results than the FPCS on raw data; therefore, the final FCC had two components from stretched data, i.e., Sentinel-2 and Landsat-8, and one component from raw data, i.e., Landsat-8. The eigen analysis of all data sources revealed a correlation between the eigenvector loadings and the spectral ranges being absorbed/reflected in the PCs. The spectral characteristics of granite were used to identify suitable PCs by analyzing the amount of information from the SWIR and TIR bands in each component and by the visual comparison of the PCs and their associated eigenvector loadings. The spectral information about the vegetation, granite, and marble was obtained through GEE’s ui.chart feature. Since the number of bands in the multispectral data was limited, the spectral profile showed smooth variations (as it interpolated the spectral profile between two known points) in Figure 5. The plots of different mineral components of granite compositions of case studies from the USGS spectral library are shown in Figure 6. The eigenvector loadings from the application of FPCS to the stretched and unstretched (raw) Sentinel-2, Landsat-8 OLI, and TIRS in ASTER (shown in Figure 7, Figure 8 and Figure 9) are reported in

Table 2. Eigenvector loadings after applying feature-oriented principal component selection (FPCS) to raw (Raw) and stretched (Str) Landsat-8 Operational Land Imager (OLI) and Landsat-8 Thermal Infrared Band Sensor (TIRS) data. The highlighted components were used in generating the final image shown in Figure 10a.

Landsat-8		NIR Band 5 (850–880 nm)	SWIR1 Band 6 (1570–1660 nm)	SWIR2 Band 7 (2110–2290 nm)	TIR1 Band 10 (10,600–11,190 nm)	TIR2 Band 11 (11,500–12,510 nm)
PC1	Raw	0.095	−0.125	−0.243	0.675	−0.679
	Str	0.561	0.621	0.470	0.230	0.163
PC2	Raw	0.858	0.178	0.120	0.308	0.351
	Str	0.805	−0.262	−0.484	−0.176	−0.133
PC3	Raw	0.127	−0.675	−0.707	−0.053	−0.163
	Str	−0.020	0.247	0.258	−0.792	−0.496
PC4	Raw	0.329	−0.591	−0.643	−0.358	0.030
	Str	0.191	−0.696	0.692	0.007	−0.005
PC5	Raw	−0.362	−0.385	−0.116	0.565	0.623
	Str	−0.009	0.021	0.011	0.538	−0.843

,

Table 3. Eigenvector loadings after applying FPCS on the raw (Raw) and stretched (Str) ASTER LIT data.

	ASTER		Visible and Near-Infrared (VNIR) Band 1 (520–600 nm)	Red Band 2 (630–690 nm)	NIR Band 3N (780–860 nm)	TIR1 Band 10 (8120–8470 nm)	TIR2 Band 11 (8470–8820 nm)	TIR 3 Band 12 (8920–9270 nm)
ASTER	PC1	Raw	−0.117	−0.127	−0.064	−0.505	−0.575	−0.616
		Str	0.337	−0.257	−0.075	−0.611	0.661	−0.061
	PC2	Raw	−0.624	−0.736	−0.117	−0.033	0.122	0.196
		Str	−0.877	0.148	−0.097	−0.422	0.094	−0.110
	PC3	Raw	0.014	0.128	−0.981	0.122	0.042	−0.067
		Str	0.246	0.683	0.367	−0.473	−0.233	0.255
	PC4	Raw	−0.102	−0.126	0.120	0.756	0.036	−0.621
		Str	0.236	0.278	−0.647	−0.152	−0.280	−0.590
	PC5	Raw	0.030	0.007	0.017	−0.395	0.807	−0.437
		Str	−0.027	0.607	−0.119	0.445	0.647	−0.016
	PC6	Raw	−0.765	0.640	0.072	−0.009	0.016	−0.001
		Str	0.000	0.000	0.650	0.060	0.059	−0.756

and

Table 4. Eigenvector loadings after applying FPCS to raw (Raw) and stretched (Str) Sentinel-2 MSI data. The highlighted components were used in generating the final image shown in Figure 10a.

	Sentinel-2		Blue Band 2 (490 nm)	Green Band 3 (560 nm)	Red Band 4 (665 nm)	NIR Band 8 (842 nm)	SWIR1 Band 11 (1570–1650 nm)	SWIR2 Band 12 (2115–2290 nm)
Sentinel-2	PC1	Raw	0.558	−0.387	−0.187	0.280	0.481	0.441
		Str	−0.211	−0.272	−0.363	−0.253	−0.573	−0.598
	PC2	Raw	0.494	0.674	0.114	0.303	0.199	−0.395
		Str	0.056	−0.001	0.160	−0.929	−0.050	0.325
	PC3	Raw	0.278	−0.434	0.183	−0.516	0.257	−0.607
		Str	0.501	0.525	0.468	−0.005	−0.333	−0.380
	PC4	Raw	0.257	−0.392	−0.122	0.492	−0.644	−0.329
		Str	0.098	0.122	−0.184	−0.266	0.737	−0.572
	PC5	Raw	−0.106	0.095	−0.941	−0.104	0.130	−0.257
		Str	0.317	0.485	−0.764	−0.017	−0.122	0.255
	PC6	Raw	0.538	0.212	−0.133	−0.557	−0.481	0.326
		Str	−0.769	0.632	0.077	−0.042	−0.019	−0.027

.

Several PCs that mostly retained topographical or vegetation responses and no discriminative spectral response of granite and marble were ignored. Key FPCS components were selected from Table 2, Table 3 and Table 4 based on the absolute loadings of critical bands responding to different associated minerals.

4.1. FPCS of the Landsat-8 Data

Generally, granite (feldspar + quartz) has major absorption features in the TIR region [61]. In PC3 of the stretched Landsat-8 data, the TIR bands showed significantly low emissivity, i.e., negative loadings for the TIR1 band (−0.792) and minutely higher emissivity for the TIR2 band (−0.496). At the same time, it has positive loadings in the SWIR1 (0.247) and SWIR2 (0.258) bands, showing a mild reflection in this region due to the presence of associated quartz and plagioclase. Granite is shown in dark pixels due to the indicative absorption characteristics in TIR1 and TIR2 bands (Figure 7).

PC3 of the raw Landsat-8 dataset had major contributions from the SWIR1 (−0.675) and SWIR2 (−0.707) bands and showed marble in bright pixels due to the diagnostic absorption features of calcite associated with this band. Vegetation cover and granite were not distinctive, as the NIR (0.127) and TIR1 (0.053) bands had low contributions.

PC4 of the raw Landsat-8 dataset showed granite regions in mild, bright pixels due to the low emissivity indicated by negative loadings in the TIR1 band (−0.358). In contrast, marble showed higher absorption in the SWIR2 band (−0.643) and was dark due to the absorption of calcite in the SWIR2 band. Vegetation was reported as bright pixels due to the reflection of the NIR band (0.329).

In PC3 of the stretched Landsat-8 data, the TIR bands showed significantly low emissivity for the TIR1 band (−0.792) and relatively higher in the TIR2 band (−0.496), while positive loadings in the SWIR1 (0.247) and SWIR2 (0.258) bands were observed due to the reflection of quartz and orthoclase.

PC3 of the raw Landsat-8 dataset had significant contributions from the SWIR2 band (−0.707) and showed marble in bright pixels due to the diagnostic absorption features of calcite associated with marble. Granite and vegetation cover were not distinctive, as the NIR (0.127), and TIR1 (0.053) bands had low contributions.

PC4 of the raw Landsat-8 dataset showed granite regions in mild, bright pixels due to the low emissivity indicated by negative loadings in the TIR1 band (−0.358), while marble showed higher absorption in the SWIR2 band (−0.643) and was dark due to the diagnostic absorption features of calcite in the SWIR2 band. Vegetation was reported as bright pixels due to the reflection of the NIR band (0.329).

4.2. FPCS of the ASTER Data

PC5 of the stretched ASTER data had positive loadings (0.445) and (0.647) for the TIR1 and TIR2 bands, respectively, showing granite regions in the bright color corresponding to the significantly low emissivity of granite in the TIR2 band. Furthermore, the opposite signs of the red and NIR bands indicated that the negative loading (−0.119) in the NIR band represented reflection and the positive loading (0.607) indicated the absorption of vegetation in the red band corresponding to the spectral characteristics of vegetation. PC6 of the raw data had major contributions from the red band with a positive loading (0.640) showing vegetation in dark pixels due to high absorption. It also showed granite-specific low emissivity characteristics in the TIR2 band (0.016) and mild emissions in the TIR1 and TIR3 bands.

4.3. FPCS of the Sentinel-2 MSI Data

PC4 of the raw Sentinel-2 data had the maximum negative loadings of (−0.644) and (−0.329) from the SWIR1 and SWIR2 bands showing granite in medium-toned pixels due to the absorption of the associated muscovite in both bands. It had the maximum positive loadings (0.492) from the NIR band and a negative (−0.122) red band corresponding to vegetation’s reflection and absorption features in these bands. PC3 of the stretched data has negative loadings of −0.333 and −0.380 from the SWIR1 and SWIR2 bands, respectively, due to the absorption of calcite and reflection in RGB bands. PC3 of the raw data had positive loadings of 0.278 and 0.257 for the blue and SWIR1 bands, respectively, indicating the reflection features of quartz and orthoclase; however, it is indicated as medium-toned darker pixels due to epidote and muscovite absorption in the SWIR2 band (−0.607).

4.4. The Fusion of the Selected Sentinel-2 MSI and Landsat-8 Components

The validation of the results with granite (locations ML-1–ML-4), MMLs, and field visit locations (P1–P7) of the Shewa Shahbazgarhi and Ambela complexes are shown in Figure 10. The final map (shown in Figure 10a) was obtained by the FCCs of a weighted combination of PC3 obtained from the FPCS of the stretched Sentinel-2 MSI data, PC4 of the raw Landsat-8, and PC3 of the stretched Landsat-8 datasets. PC3 of the stretched Sentinel-2 MSI component was assigned a weight of 1.5, and PC3 of Landsat-8 stretched datasets had more matching with the target (granite/marble) than the raw Landsat-8 data. Therefore, it was assigned ¾ of the remaining 0.5 weight, i.e., 1.125, and PC4 of the raw Landsat-8 datasets was weighted as 0.375. Figure 10b,c shows the results obtained by applying PCA on stacking all bands of the raw and stretched ASTER, Sentinel-2, and Landsat-8 OLI and TIR sensors mentioned in Table 1. The results from the proposed method in Figure 10a showed a good correlation with ground truth locations compared to Figure 10b,c.

A binarized normalized difference vegetation index (NDVI) was calculated and overlaid on the final map for comparison, which showed pixels containing of >0.5 NDVIs in black and pixels containing less than 0.5 NDVIs in white color (as seen in Figure 11). Granite having no vegetation cover was visible in sky-blue and sea green. The dark reddish black color at the MML region showed the marble mining areas, while the dark sky-blue color near P4, P5, and ML-3 indicated the granite with scree, clay, and disturbed land area, respectively, due to weathering. The results from our proposed methodology (Figure 10a) correctly identified all granite and marble known locations compared to the FPCS of the stacked raw (Figure 10b) and stretched (Figure 10c) bands, identifying 5 and 10 out of the 11 locations, respectively.

Furthermore, information on granite mining leases was acquired from the Directorate General of Mines and Minerals (DGMM) of Khyber Pakhtunkhwa Province in Pakistan (

Table 5. Validation of results after they were compared with known locations of granite Leases (locations ML-1–ML-4, MMLs, and field survey points P1–P7 in Figure 2 of Shewa Shahbazghari and Ambela complexes.

S. No	Techniques	Field Points (Ground Truth)							Mining Leases (Ground Truth)					Accuracy
		P1 (G)	P2 (G)	P3 (G)	P4 (G)	P5 (G)	P6 (M)	P7 (M)	ML1 (G)	ML2 (G)	ML3 (G)	ML4 (G)	MML (M)
1	Proposed weighted linear combination of FPCS components	G	G	G	G	G	M	M	G	G	G	G	M	100%
2	FCC of PCs obtained after PCA applied to 17 raw bands from 3 data sources	VL	VL	S	G	G	G	G	G	G	G	VL	G	44%
3	FCC of PCs obtained after PCA applied to 17 stretched bands from 3 data sources	VL	VL	G	G	G	M	M	G	G	G	G	M	82%

G, granite; M, marble; V_L, light vegetation; V_D, dense vegetation.

) and the field validation of the area under study to traverse the target regions using the existing geological map [38]. The field features such as the color, texture, mineralogy, and weathering grade of rocks were correlated with the maps generated via the RS techniques discussed above. Dense to mild vegetation was indicated in a pink to the red gradient (as shown in Figure 10a), which correlated with the vegetation intensity on the rocks exposures (as shown in Figure 11). The magnified views of the final map presented in Figure 10a, with the details of the study area is shown in Figure 12. Detailed petrographic and textural examination characterized the rocks such as granites and marble (Figure 13), supporting the RS analytical results.

The petrographic study in Shewa Shahbazgarhi also highlighted the presence of potassium feldspar and silica-rich granite rock bodies and calcite-rich marble rock bodies in the Shewa Shahbazgarhi region. It indicated that the stacked FPCS could highlight small differences for geological mapping in the study area. The application of the methodology in Malakand and Karoonjhar Mountains showed the effectiveness of this technique to provide economic geological maps of an area.

The results from the proposed methodology were compared with the PCA of stacked bands from all three datasets. PC2, PC3, PC4, and PC9 of the raw stacked data had significant information from the SWIR and TIR spectral ranges to discriminate granite and marble. Likewise, PC1–PC9 of the stretched stacked data provided useful information about the targeted rock bodies. However, each PC being a mixture of 17 bands made it difficult to select a specific PC retaining distinctive spectral response of granite and marble. Therefore, the FPCS application on bands from the same data sources separately is preferred due to a 100% accuracy based on the field survey and the known mining leases of granite and marble. On the other hand, stacking all the bands from all three datasets resulted in 50% and 75% accuracies in raw and stretched bands because of spectral mixing.

Table 6. Eigenvalues of the selected band ranges/center (in µm) of different data sources and components in Malakand region.

Data Sources Selected Component	Description of Associated Bands of the Data Sources Associated Eigen Loadings of Selected Component
Sentinel-2	Blue B2 (0.49)	Green B3 (0.56)	Red B4 (0.665)	NIR B8 (0.842)	SWIR1 B11 (1.570–1.65)	SWIR2 B12 2.115–2.29)
PC5 of the stretched data	0.1631	−0.0168	−0.1254	−0.0206	−0.2979	−0.9317
ASTER	VNIR B1 (0.52–0.6)	VNIR B2 (0.63–0.69)	NIR B3N (0.78–0.86)	TIR 1 B10 (8.12–8.47)	TIR 2 B11 (8.47–8.82)	TIR 3 B12 (8.92–9.27)
PC2 of the raw data	−0.2282	−0.2717	−0.2819	−0.1059	0.4852	0.2368
Landsat-8	NIR B5 (0.85–0.88)	SWIR1 B6 (1.57–1.66)	SWIR2 B7 (2.11–2.29)	TIR1 B10 (10.6–11.19)	TIR2 B11 (11.5–12.51)
PC5 of the stretched data	−0.4988	−0.4413	−0.0721	0.5150	0.5347

and

Table 7. Eigenvalues for the selected band ranges/center (in µm) and the components of Nagarparkar granite.

Data Sources Selected Component	Description of Associated Bands of the Data Sources Associated Eigen Loadings of Selected Component
Sentinel-2	Blue B2 (0.49)	Green B3 (0.56)	Red B4 (0.665)	NIR B8 (0.842)	SWIR1 B11 (1.570–1.65)	SWIR2 B12 (2.115–2.29)
PC2 of the raw data	−0.1879	−0.1615	−0.3477	0.8675	0.1512	−0.2055
PC5 of the stretched data	0.2935	0.2873	−0.2607	0.6071	0.5116	0.3644
ASTER	VNIR B1 (0.52–0.6)	VNIR B2 (0.63–0.69)	NIR B3N (0.78–0.86)	TIR 1 B10 (8.12–8.47)	TIR 2 B11 (8.47–8.82)	TIR 3 B12 (8.92–9.27)
PC5 of the stretched data	−0.3434	−0.1360	−0.4276	−0.3374	−0.6024	−0.4515

present the eigen loadings of the selected FPCS components for mapping granite outcrops of the Malakand and Nagarparkar in Thar region. PC2 of ASTER raw data reported granite in bright due to positive loadings of 0.4825 and 0.2368 in the TIR1 and TIR2 bands, respectively, as given in Table 6, showing the low emissivity of granite against opposite (negative) loadings on B1, B2, and B3N. PC5 of the stretched Landsat-8 data showed the maximum positive loadings of 0.515 and 0.534 in B10 (TIR1) and B11 (TIR2) bands, respectively, corresponding to the significantly low emissivity (shown in bright which can be seen in dark after the negation of the component) due to granite rock in the thermal wavelengths. On the other hand, this component had opposite-sign (negative) loadings for the SWIR1 band (−0.44) corresponding to the reflection of granite-associated quartz in the SWIR1 band and no information about the SWIR2 band. PC5 of the stretched Sentinel-2 data showed negative loading (−0.9317) in the SWIR2 band, mainly corresponding to the high absorption of muscovite, epidote, and amphibole in the SWIR2 band.

Figure 14 and Figure 15 present the results from the application of the proposed methodology to the granite outcrops of the Kotah Dome in Malakand (Khyber Pakhtunkhwa, Pakistan) and Karoonjhar (Nagarparkar, Sindh) graqz8nitic rock body surrounded by sandy hills of the Thar Desert.

The Nagarparker region had no considerable vegetation cover over the granite rock bodies; therefore, the ASTER performed relatively better in highly vegetated areas. In the case of PC5 of the stretched ASTER data, the TIR2 band (−0.6024) showed low emissivity as a distinctive feature of granite minerals indicated by negative loadings, as shown in Table 7. Similarly, relatively higher emissivity in the TIR1 (−0.3374) and TIR3 (−0.4515) bands was observed due to the high reflection of quartz countering absorption of all the other minerals. PC2 of the raw Sentinel data had a positive loading (0.151) in the SWIR1 band and a negative loading in the SWIR2 band (−0.2055) due to mild absorption tones of muscovite. The negation of PC5 of the Sentinel-2 stretched data would represent negative loadings (absorption) due to the positive loadings of the NIR (0.60) and SWIR1 (0.51) bands and the reflection of the red band (−0.26) due to iron oxide.

5. Discussion

In Shewa Shahbazgarhi area of Khyber Pakhtunkhwa Province of Pakistan components, SWIR bands (Sentinel-2 and Landsat-8) were also helpful in mapping marble outcrop in the area. Absorption was indicated in the SWIR2 band due to calcite_, and both components of Landsat-8 indicated low emissivity in the TIR1 band due to the presence of orthoclase and quartz. PC3 of the stretched Landsat-8 data is granite, indicating components representing the reflection in the SWIR1 and SWIR2 bands and the absorption in TIR bands due to quartz, orthoclase, and biotite in granite. PC4 of the raw Landsat-8 dataset was marble, indicating components representing the high absorption of calcite in the SWIR2 band. In contrast, PC3 of the stretched Sentinel-2 data represented the weak absorption in the SWIR1 band due to biotite and muscovite (Figure 15) and a relatively high absorption of calcite in the SWIR2 band.

The Kotah Dome, Malakand, rich in quartz, orthoclase (K-feldspar), biotite, and muscovite with garnet, chlorite, and amphibole accessories [60,61] was mapped by equally weighted PC5 components of the stretched Landsat-8 and Sentinel-2, and PC2 of the raw ASTER data. PC5 of the stretched Sentinel-2 data with negative loadings (−0.9317) in the SWIR2 band represented the high absorption of muscovite, epidote, and amphibole in the SWIR2 band. While PC5 of the stretched Landsat-8 represented quartz and orthoclase due to its reflection in the NIR bands and absorption in the TIR bands. ASTER components represented the mixed response of all granite-associated minerals quartz, orthoclase, biotite, and muscovite due to their significantly low emissivity in the TIR2 band and the high reflection in the SWIR bands.

Nagarparker pink-colored granites is rich in orthoclase + plagioclase and quartz, with small amounts of iron oxide, biotite, titanite, chlorite, and epidote [42]. PC2 of the raw (weighted as 0.5) Sentinel-2 data showed the highest loading in the NIR bands, representing the reflection of granite-associated minerals (epidote, muscovite, and amphibole) and their absorption in the SWIR2 band. PC5 of the stretched (weighted as 1.0) Sentinel-2 data contained the highest loading in the NIR bands and relatively lower loadings in SWIR bands due to the mixed reflection of granite-bearing minerals. PC5 of the stretched ASTER data (weighted as 1.5) represented the corresponding absorption features of quartz, plagioclase, limonite, orthoclase, and amphibole in the TIR bands.

The Karoonjhar and Kotah Dome has lesser vegetation cover as compared to Shewa Shahbazgarhi; therefore, the ASTER performed relatively better to identify granitic rock bodies in those areas due to the low emissivity in TIR bands. However, in the presence of vegetation in Shewa Shahbazgarhi, the TIR bands were less effective due to their lower spatial, temporal, and spectral resolutions [11]. Therefore, the fusion of the TIR-bearing (Landsat-8) and SWIR-bearing data (Sentinel-2) were productive. Additionally, the SWIR bands in selective components with the TIR bands (Landsat-8) may represent reflection of orthoclase (potassium feldspar) quartz and the absorption of accessory amounts of amphibole, sphene, zircon, garnet, and mica in the absence of TIR bands as in Sentinel-2 data [62,63].

In all the three case studies, granite low emissivity was indicated in TIR bands (8.8 µm), while the carbonate absorption of the SWIR2 band was prominent in most of the components during the mapping of Shewa Shahbazghari granite. The low emissivity in TIR bands was observed due to granite-bearing quartz, orthoclase, and plagioclase in PC4 of the raw Landsat-8 band and PC3 of the stretched Landsat-8 (TIR1) band. The SWIR1 band generally indicates the reflection features of main granite components such as quartz and plagioclase, while the absorption is shown in the SWIR2 band mainly due to epidote and muscovite. The granite absorption in the TIR bands was captured by PC2 of the raw ASTER (TIR2) band and PC5 of the stretched Landsat-8 (TIR1) band during the mapping of the Kotah Dome in Malakand region. Absorption due to amphiboles, chlorites, epidote, tourmaline, and muscovite in the SWIR2 band was reported by PC5 of the stretched Landsat-8 and the stretched Sentinel-2 data due to OH/metal–OH [64] generally and Mg–OH combinations, particularly in case of muscovite [4]. Nagarparkar granite was delineated by PC2 of the raw ASTER component-bearing TIR band and PC5 of the stretched Sentinel-2 band representing high reflection for amphibole, quartz, and plagioclase in NIR and SWIR1 bands and relatively lesser reflection in the SWIR2 band due to limonite and muscovite.

The importance of the fusion of FPCS components-bearing datasets containing TIR, SWIR, and/or SWIR + TIR bands has been highlighted. One Sentinel-2 and two Landsat-8 components in the case of high-vegetation Shewa Shahbazghari granite; one ASTER with two Sentinel-2 components for mild-vegetation Malakand granite and in the case of non vegetated Nagarparkar one component from the ASTER, Landsat-8, and Sentinel-2 data each were used to map granite. In all cases, fundamental low emissivity was observed in the TIR bands (ASTER TIR2 or Landsat-8 TIR1), and the absorption of the SWIR bands (Sentinel-2 or Landsat-8) by granite accessory minerals bands played a role. Due to the high spatial resolution and absorption of the SWIR2 band by accessory minerals, Sentinel-2-bearing components were useful to fuse selective components from ASTER and/or Landsat-8 datasets containing TIR bands.

FPCS can be further used by selecting a limited number of bands from the high-resolution hyperspectral data for mapping dimension stones. The proposed methodology can also be used for data annotation [65] in supervised techniques for lithological mapping. Mineralogy is also used to predict geotechnical properties [38,66,67]; therefore, SWIR2 bands/components identifying accessory minerals could play an important role in modelling geotechnical parameters via mineralogy through RS. The proposed mapping strategy can be used as a cheap source for mineral resource modelling combined with shallow seismic geophysics to model the deposit features in three-dimensions.

6. Conclusions

• A time-efficient methodology for fusing selective FPCS components from stretched and raw Landsat-8 OLI + Landsat-8 TIRS, ASTER, and Sentinel-2 MSI is presented to map granite and marble of Shewa Shahbazghari area. The process took 147.14 s on the GEE Cloud computing platform to carry out PCA in three multispectral data sources’ raw and stretched forms.
• A weighted linear combination of PCs was carried out based on each component’s granite and marble indicating features to enhance granitic spectral features while suppressing others. The results were validated by matching all the exposed granite and marble quarries successfully.
• The TIR bands of the ASTER had a 90 m resolution that did not perform well in densely vegetated areas; therefore, a combination of Landsat-8 and Sentinel-2 components discriminated the granite/marble rock types better in Shewa Shahbazgarhi due to much better resolutions. Furthermore, the non-availability of granite absorption features from the TIR bands in Sentinel-2 MSI datasets was overcome by Landsat-8 TIR bands. As a result, the high-resolution SWIR bands from Sentinel-2 and the medium-to-low-resolution TIR bands from Landsat-8 TIRS helped obtain the final FCC map for discriminating granite from marble.
• The weighted PCs extrapolation to other granite outcrops across Pakistan showed the technique’s usefulness to obtain results with ease and accuracy. In less vegetated regions, i.e., the Kotah Dome, Malakand, and non-vegetated Karoonjhar Nagarparkar, the ASTER component was helpful due to the granite-specific TIR2 band. The advantage of Sentinel-2 components was highlighted in all case studies, which shows that high-spatial-resolution SWIR bands support TIR bands from the ASTER and/or Landsat in mapping granite and accessory minerals. The information about accessory granite minerals due to their spectral absorption in the SWIR2 band is incorporated in granite through the high-spatial-resolution Sentinel-2 components bearing high-resolution SWIR2 bands.
• Based on the analysis of eigenvector loadings, the low emissivity in Landsat-8 TIR1 and the reflection of the SWIR1 (1570 nm–1660 nm) and SWIR2 (2110 nm–2290 nm) bands were observed due to orthoclase and quartz in granite. The SWIR2 components from Landsat-8 and Sentinel-2 indicated absorption characteristics of granite accessory minerals, i.e., muscovite, epidote, and amphibole. Similarly, marble showed absorption features in the SWIR2 (2110–2290 nm) band of Landsat-8 OLI and Sentinel-2.
• The methodology can be further used for high-end RS applications:

  a.

  Predicting geotechnical properties through mineralogy using freely available multispectral or commercial high-resolution datasets (Worldview3, Quick bird, Spot Eye, etc.).

  b.

  Unsupervised classification for better data annotation before using machine learning/deep learning algorithms.