Technological Advances in Mineral Exploration. The future with TerraEye (Part 1)
Mineral prospecting, the endeavor to discover valuable geological deposits on Earth, has undergone significant technological advancements in recent years. The methods have evolved from the traditional ‘boots on the ground’ approach to aerial surveys and satellite imaging. Among these, one of the most promising and less intrusive methods emerging is spectral analysis. Utilizing spectral data from Earth-observing satellites, this technique enables us to ‘see’ beyond what is visible to the naked eye, facilitating more accurate assessments of a terrain’s mineral composition. By integrating this data with existing mineralogical frameworks, we can identify potential raw material deposits, a process we aim to fully automate through TerraEye.
In this article, we will delve into the technical aspects of spectral analysis, with a particular focus on a widely recognized method: Spectral Angle Mapper (SAM). SAM proves effective in mapping certain minerals such as jarosite and alunite, which could serve as indicators when prospecting for porphyry copper for instance.
We will employ actual satellite images from Sentinel-2, along with spectral curves of the mixture of jarosite, to illustrate this method. Anticipate a fusion of theory, Python code, and visualizations to navigate you through this intriguing field of research.
Case Study: Unearthing the Secrets of the Atacama Desert
The Atacama Desert in Chile, renowned as one of the driest and most inhospitable regions on Earth, has long captivated geologists and mineral prospectors. Its unique geological attributes render it a fertile ground for mineral exploration, including the search for rare minerals with significant commercial potential. In this case study, our focal point is jarosite, a mineral that can act as a signpost for other valuable materials lurking below the surface.
Jarosite, a hydrous potassium iron(III) sulfate mineral, is typically found in acidic, sulfate-rich environments, like those formed by the oxidation of sulfide minerals. It’s a pivotal indicator mineral for the presence of water and acid-sulfate alteration, making it significant not merely on Earth but also in planetary exploration scenarios, such as on Mars.
The harmonization of Sentinel-2 data with the jarosite spectral signature is far from coincidental in this context. Upon examining the jar5osite spectral signature plot, it becomes apparent that the spectral bands showcasing distinctive features align extensively with the bands of the individual spectral channels of the Sentinel-2 data (denoted by black squares). We will incorporate the jarosite spectral signature for Sentinel-2 Multi-Spectral Instrument (MSI) sourced from the U.S. Geological Survey (USGS) as a benchmark for our spectral analysis.
However, the correlation between Sentinel-2 data and what we labelled as a ‘jarosite spectral signature’ deserves clarification. In reality, the reference spectrum used is not solely indicative of jarosite; it’s a composite signature featuring a mixture of Quartz, Goethite, and Jarosite. The rationale for using a mixed reference spectrum is rooted in the fact that each pixel in Sentinel-2 data covers an area of 100 square meters; for SWIR bands, this extends to 400 square meters. Given this spatial resolution, it’s unrealistic to expect a ‘pure’ jarosite spectrum in the data. Instead, the spectral signature will likely be a complex interplay of multiple minerals existing within that footprint.
At this point, it’s worth highlighting the role of a skilled geologist on the analytics team. Her/His knowledge and expertise are crucial for accurately identifying and evaluating the kinds of mineral mixtures most likely to be encountered, as well as determining which spectral signatures may be most informative about geological settings, thereby enhancing the integrity and reliability of spectral analysis.
Data obtaining
Initially, we retrieve the spectral signature data for jarosite, followed by unpacking and reading the file contents. Given that the original files consist of a header succeeded by a series of float numbers, we will bypass the header and import the numerical values into an array. You have the option to download the data independently or utilize the text provided below to generate a file with the spectral data.
S07SNTL2 Record=1840: Jarosite NMNH95074-1 (Na) BECKb AREF
4.5296192e-002
1.1021909e-001
2.3987867e-001
4.1108599e-001
4.7953874e-001
4.7384530e-001
3.9147010e-001
3.0794874e-001
2.6794386e-001
2.5157878e-001
6.4951229e-001
7.0593065e-001
5.1015931e-001As observed, the jarosite spectral signature has been prepared for 13 bands (L1C product). However, for further analysis, it’s imperative to remove the value for band 10 as this band is not present in the L2A products.
import numpy as np
def read_curve_from_txt(file_path: str, skip_header: int = 1) -> np.array:
"""
Reads a spectral curve from a text file and returns it as a NumPy array.
...
"""
with open(file_path, "r") as f:
lines = f.readlines()[skip_header:]
ref_signature = [float(line.strip()) for line in lines]
return np.array(ref_signature)
def remove_band_from_curve(curve: np.array, band_index_to_rm: int = 9):
"""
Removes a specific band value from a spectral curve array.
"""
return np.delete(curve, band_index)
# Read jarosite curve from txt file
jarosite_curve = read_curve_from_txt("./spectral_curves/S2_jarosit.txt")
# Remove band 10 value
jarosite_curve = remove_band_from_curve(jarosite_curve)Having prepared the first component, we now transition to preparing data from satellite observations. The multispectral data was procured via SentinelHub.
Sinergise has delved into this topic in their texts, which are worth exploring for further insight if necessary.
Now, we’ll proceed to read the Sentinel-2 data, prepare its RGB representation, and examine the original form of our data. Below are two plots depicted in Figure 2. The left plot presents the RGB representation of a selected Sentinel-2 image with a specified pixel (80,120). Conversely, the right plot illustrates the spectral signature for the selected pixel, juxtaposed with the spectral signature of the jarosite from the USGS.
from typing import Tuple
def plot_rgb_and_signatures(rgb_image: np.array, pixel_coordinates: tuple, pixel_curve: np.array, jarosite_curve: np.array) -> None:
# Create a figure and a 1x2 subplot grid
fig, axes = plt.subplots(1, 2, figsize=(15, 6))
# Plot the RGB image with the selected pixel marked in red
axes[0].imshow(rgb_image)
axes[0].scatter(pixel_coordinates[1], pixel_coordinates[0], color='red', s=50) # Mark the selected pixel
axes[0].set_title("RGB Image with Selected Pixel")
axes[0].axis('off')
# Plot the spectral curves
axes[1].plot(pixel_curve, label='Spectral Curve of Selected Pixel', color='blue')
axes[1].plot(jarosite_curve, label='Reference Jarosite Curve', color='green', linestyle='--')
axes[1].set_title("Spectral Curves")
axes[1].set_xlabel("Band")
axes[1].set_ylabel("Reflectance")
axes[1].legend()
# Display the plots
plt.tight_layout()
plt.show()
# Read the Sentinel-2 image
with rasterio.open("./S2_images/chile.tif") as src:
img_array = src.read()/10000 # convert DN into reflectance
profile = src.profile
# Create an RGB image from Sentinel-2 bands (assuming bands 4, 3, 2 for RGB)
rgb_image = img_array[1:4, :, :] # Select only RGB bands
rgb_image = img_array[[2, 1, 0], :, :] # Slicing to get bands 3, 2, 1
rgb_image = np.transpose(rgb_image, (1, 2, 0)) # Transpose to get height, width, bands order
# Normalize the RGB image to [0, 1]
rgb_image = rgb_image.astype(float)
rgb_image = (rgb_image - np.min(rgb_image)) / (np.max(rgb_image) - np.min(rgb_image))
pixel_coordinates = (80, 120) # Replace with the coordinates of your selected pixel
selected_pixel_curve = img_array[:, pixel_coordinates[0], pixel_coordinates[1]] # Extracting the spectrum of the selected pixel
plot_rgb_and_signatures(rgb_image, pixel_coordinates, selected_pixel_curve, jarosite_curve)
Data normalization
Normalization is a key step in data analysis that helps in comparing different sets of data on a uniform scale. In the context of spectral analysis, normalization is particularly important for eliminating differences in intensity and enabling the comparison of the shapes of spectral curves for different materials.
A common practice in normalization is employing the L2 (Euclidean) norm. This approach proves beneficial when the aim is to ensure a fair comparison among different spectral signatures. The L2 norm transforms the original spectral curve such that its length, within the vector space, is normalized to 1. This transformation allows for a comparison of curve shapes, irrespective of their reflectance levels.
The choice of L2 normalization is notably justified when the objective centres around comparing the shapes of spectral curves rather than their absolute reflectance values. For instance, in scenarios such as mineral identification, including jarosite detection, the spectral curve’s shape often emerges as a more crucial indicator compared to its reflectance. Hence, the L2 norm serves as an apt choice for normalization in such contexts.
# Normalize the jarosite curve
jarosite_curve_norm = jarosite_curve / np.linalg.norm(jarosite_curve)
# Reshape the image array for easier calculations
bands, height, width = img_array.shape
img_array_reshaped = img_array.reshape(bands, -1)
# Normalize each pixel's spectral curve
img_array_norm = np.linalg.norm(img_array_reshaped, axis=0)
img_array_reshaped_norm = img_array_reshaped / img_array_normSpectral Angle Mapper (SAM) method
The Spectral Angle Mapper (SAM) is a physically-based spectral classification technique that uses the n-dimensional angle to match pixels to reference spectra. The smaller the angle between the pixel vector and the reference, the better the match. The angle θ in n-dimensional space is calculated as:
where:
- vec{a} is the n-dimensional vector of the pixel
- vec{b} is the n-dimensional vector of the reference spectrum
Here’s a Python code snippet for implementing SAM. This example assumes that you have two NumPy arrays: img_array_reshaped for the set of pixels spectrum and jarosite_curve_norm for the reference spectrum of jarosite.
# Calculate the spectral anglecosine_similarity = np.dot(jarosite_curve_norm, img_array_reshaped_norm)
spectral_angle = np.arccos(cosine_similarity) # in radians
# Reshape back to the original image shape
spectral_angle = spectral_angle.reshape(height, width)
# Save the SAM result with the original georeferencing and metadata
profile.update(dtype=rasterio.float32, count=1)
with rasterio.open("SAM_results.tif", "w", **profile) as dst:
dst.write(spectral_angle.astype(rasterio.float32), 1)
To illustrate how SAM works, consider a satellite image from the Atacama Desert focused on an area suspected to contain jarosite. Below is the SAM classification map produced by applying the algorithm on each pixel in this region, comparing it with the known spectral curve of jarosite.
Areas marked in red/yellow have the smallest spectral angles, indicating a high likelihood of jarosite presence. On the other hand, areas in blue/green have larger spectral angles, making it less likely to find the jarosite.
Let’s examine how well the curve from the image matches the reference curve at the point where the SAM value is the lowest (highest probability of occurrence). First, let’s identify the coordinates of this point, and then use the previously defined function to plot the spectral signatures for the selected pixel.
# Determining the pixel position with the minimum SAM value
min_index_2d = np.unravel_index(np.argmin(spectral_angle, axis=None), spectral_angle.shape)
pixel_coordinates = min_index_2d
selected_pixel_curve = img_array[:, pixel_coordinates[0], pixel_coordinates[1]]
plot_rgb_and_signatures(rgb_image, pixel_coordinates, selected_pixel_curve, jarosite_curve)
Working with values in radians can be unintuitive, as we are generally more accustomed to other angular units. This issue can be addressed by converting radian values to cosine similarity values. Cosine similarity values range from -1 to 1, where 1 indicates identical vectors, 0 signifies orthogonal vectors, and -1 denotes diametrically opposed vectors. This range can sometimes be more intuitive to interpret than the 0 to π range for the spectral angle in radians.
# Convert into cosine similarity values
spectral_angle = np.cos(spectral_angle)
Figure 5. The radian output from the SAM algorithm is converted to cosine values for easier interpretation.
Through this case study, the SAM algorithm successfully identified multiple potential jarosite occurrences within the Atacama Desert, demonstrating its efficacy and potential for mineral prospecting.
The final step is to compare the obtained results with the geological findings from the La Polvora exploration project (2019). In areas where the SAM values indicate a high probability of jarosite presence, this has also been confirmed by field studies. However, there are also areas that show high potential for the presence of jarosite but have not been identified by geologists.
— -
This chapter provides a comprehensive view of the SAM method, from the theoretical basis to practical Python code and real-world application. By the end, the reader should have a solid understanding of how SAM works and how to implement it for mineral prospecting in satellite images.
The next part will describe the SFF method. Stay tuned!
Sources
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., Klein, A.J., 2017, USGS Spectral Library Version 7: U.S. Geological Survey Data Series 1035, 61 p., https://doi.org/10.3133/ds1035.
Kruse, F. A., A. B. Lefkoff, J. B. Boardman, K. B. Heidebrecht, A. T. Shapiro, P. J. Barloon, and A. F. H. Goetz. “The Spectral Image Processing System (SIPS) — Interactive Visualization and Analysis of Imaging spectrometer Data.” Remote Sensing of Environment 44 (1993): 145–163.
