Quick Soil Sampling Data Summary With Python

Geologists often use soil sampling methods in mineral exploration to gather information about element anomalies in the soil of a certain area. The output of this activity is a database containing coordinates and various element concentrations from the sample taken at that point. Before diving deep into the data and doing any fancy analysis, it is always a good idea to familiarize ourselves with the data by looking into the summary statistics and the general and spatial relationship of the data as a whole. By using Python we can write a simple code to get the summary statistics and some useful graphs to summarize the data in just a few minutes.

We will use pandas, numpy, and matplotib libraries for this code, so let’s begin with importing our library

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

We can use pandas to read our dataset and put it into a variable called data. We can then use head() function from pandas to look at the first 5 rows of the data.

data = pd.read_csv('soil_data_example.csv')
data.head()

Note: for this example, we use data with dummy coordinates.

Summary statistics

By using pandas, getting summary statistics is as simple as this:

pd.set_option('display.max_columns', None)
data.iloc[:,2:].describe()

For practicality, pandas by default will not show the full range of rows and columns of our dataset. The first line of the code tells pandas to show the full range of the columns. In the second line, we use iloc[:, 2:] attribute to ignore the first two columns of our dataset because we don’t need summary statistics of the x and y coordinate (unless for some reason you do).

Histograms

Next, we will draw the histogram of each element’s concentration in the area. Because of the logarithmic nature of element concentration in the soil, we will plot the log concentration value into the histogram. Since we are dealing with a lot of 0’s, we will add a very small value to all the concentration data to make it easier with the logarithmic transformation. There are various ways to deal with zero values in the dataset, but for this example, let’s just add 0.0001 to all concentration data.

data.iloc[:,2:] = np.log(data.iloc[:,2:]+0.0001)

np.log is a natural logarithm with base e if we want logarithm transformation with base 10 we can use np.log10. Now if we call data.head() we will see that that the concentration of each element already transformed into its logarithmic value.

Note: code above will return an error if there are still zero values in the datasets. There are various ways to deals with 0 and NaN value in our dataset and it was up to you which method you will use for your case. You can also ignore this line if you don’t want to do any logarithmic transformation for your dataset

data.head()

Next, we will use matplotlib to define our canvas and axes and then loop into each element in our dataset to plot each histogram.

fig, ax = plt.subplots(nrows=7, ncols=5, figsize=(20,8), sharey=True)
plt.subplots_adjust(top=3)


row = 0
col = 0
for element in data.iloc[:,2:].columns:
    x=list(data[element])
    ax[row,col].hist(x, bins=35)
    ax[row,col].set_title("log " + element)
    if col<4:
        col+=1
    else:
        col=0
        row+=1
        
plt.show()

The code above can be separated into two sections. The first section of the code will define the figure and axis canvas for our plot. Since we are dealing with 35 variables we will divide our subplots into 7 rows and 5 columns.

We can check the length of our variable with:

'''
it is getting annoying to keep typing data.iloc[:,2:] 
to get our elements so let's just put it into a
variable called grade
'''
grade = data.iloc[:,2:]
    
len(grade.columns)

35

The second section of the code defines the loop behavior that will plot each histogram into predefined subplots.

Correlation Matrix

Our next figure will be a correlation matrix. This is a very useful and neat graph that will tell us the correlation coefficient of one variable with another variable within the dataset. The correlation matrix value ranges from -1 to 1, where:

• -1 = perfect negative correlation
• 0 = no correlation
• 1 = perfect positive correlation.

Pandas come in handy with its corr() function to calculate a correlation matrix of our dataset. Combine it with matplotlib functions such as: matshow() to draw the matrix into a figure; cmap() to denote it with colormap; and colorbar() to draw the legend, and we can get this informative graph.

cor_matrix = grade.corr()
fig = plt.figure(figsize=(20, 18))
ax = fig.add_subplot(111)
cax = ax.matshow(cor_matrix,cmap='coolwarm', vmin=-1, vmax=1)
ticks = np.arange(0,len(grade.columns),1)

fig.colorbar(cax)
ax.set_xticks(ticks)
plt.xticks(rotation=90)
plt.title('Correlation Matrix')
ax.set_yticks(ticks)
ax.set_xticklabels(grade.columns)
ax.set_yticklabels(grade.columns)

plt.show()

Map plots

With the same concept as the histogram plot, we can make quick and dirty map plots of our points to get a glimpse of the spatial correlation of each element. For this purpose, we use scatter() function with our UTM coordinates as its x and y value and denote the concentration value of each element to a color map using cmap(). We do this using for loop so it can loop to each element in our grade data.

fig, ax = plt.subplots(figsize=(20,50), nrows=12, ncols=3)

rows=0
cols=0 
for element in grade:
    cax = ax[rows,cols].scatter(data['EASTING'], data['NORTHING'], c=grade[element], cmap='coolwarm', s=10, alpha=0.9)
    ax[rows,cols].set_xticklabels('')
    ax[rows,cols].set_yticklabels('')
    ax[rows,cols].set_title(element)
    ax[rows,cols].axis('equal')
    fig.colorbar(cax, ax=ax[rows,cols])
    if cols<2:
        cols+=1
    else:
        cols = 0
        rows += 1
ax[11,2].set_axis_off()

plt.show()

And that’s it. Your soil sampling data summary in just about a minute. You can download the full code on my GitHub page and reuse it for your next project. Hope it can save you the trouble and you can be more focused on your real geologist work. Have fun exploring!