Depositional Facies Belt Mapping with SVM in Python + some improvization
Say that you have well data with facies label in it, and you meant to make the facies map out of it. Well, we can do in a conventional way, but this article meant to show you an alternative solution in making Facies Belt Map using Python + SVM.
Facies Belt Map is one of the vital input in Reservoir Modelling in subsurface that made by geologist. The data required to make this map is point data consist of coordinate (x,y), and Facies Code. I will try to give you the basic example on how to make this map.
What is Facies
There are ambiguity in defining what Facies actually is in geological terminology. Facies has been described as “.. as a way to distinguish bodies of rock units in terms of physical characteristic, composition, formation, or various other attributes..” (Source: SEG Wiki). In sedimentological terminology Facies can be described as “…a unit of rock that deposited in same depositional environment…”. In this article, we’re gonna use Facies terminology that being used by Sedimentologist.
Facies bands are zones that characterized by same depositional process. Each bands possess same depositional process, thus the sediment product supposed to be same in term of sedimentary texture and structure.
What is SVM?
SVM (Support Vector Machine) is a set of supervised learning methods used for classification, regression and outlier detection. This algorithm has developed in 1992 and proven as one of the best classification solution at that time.
SVM algorithm will draw a vector line that best separates two (or more) objects in a n-dimensional dataset, and the plane dictates which category the data falls into.
Into Python
There are four simple steps in order to make facies belt map using SVM.
- Load the data form source. The data supposed to consist x,y coordinate as well as facies label in Integer data format.
- Create classification Model, and fit the data into the model. This one is very simple
- Create meshgrid that covers spatial extend selected area. And reshape the Meshgrid from 2d array into 1d array.
- predict the result using Predict method, and visualize the result
- Make the facies map “More Geology Looks”
Data Loading
import the required package.
import pandas as pd
import numpy as npimport matplotlib.pyplot as pltfrom sklearn.svm import SVC
the data we’re about to use is Toy data that consist 30 points of X, Y coordinate and lablels:
x = [0.38436421, 0.82736184, 0.61237846, 0.39413618, 1. ,
0.90879493, 0.49511322, 0.30618844, 0.27687253, 0.75244339,
0.85993507, 0.61237846, 0.4788266 , 0.68078226, 0.08794774,
0.78827395, 0.6449517 , 0.4234521 , 0.31270309, 0.18241014,
0.68403959, 0.37459224, 0.0814331 , 0.57980522, 0.21172605,
0. , 0.46253998, 0.15635154, 0.56677435, 0.46579731,
0.2182407 ]
y = [1. , 0.89160839, 0.85664336, 0.82867133, 0.81118881,
0.74125873, 0.73076924, 0.72027972, 0.67482517, 0.64685314,
0.64335664, 0.61888113, 0.58391608, 0.55244755, 0.51398602,
0.48601398, 0.45104895, 0.41608392, 0.40209791, 0.36013986,
0.36013986, 0.35664336, 0.33916084, 0.2937063 , 0.2902098 ,
0.24825175, 0.16433567, 0.16083917, 0.13636364, 0.00699302,
0. ]
label = [2, 0, 1, 2, 0, 0, 2, 2, 2, 1, 0, 1, 1, 1, 2, 0, 0, 1, 2, 2, 0, 1,2, 0, 2, 2, 0, 1, 0, 0, 1]plt.scatter(x,y,c=label) #lets visualize the data
The data consist of well aranged points that showing orientation preference in distribution. We can easily draw the line that discriminate the labels for each facies with pen and pencil, but that’s not our goal at this moment.
Classification Model and Train the Model
We will use SVC (Support Vector Classification) function from Sklearn to perform SVM classification. The kerenel that being used for this model is ‘poly’ with degree of polinomial of 1 (it is similar with linear, but resulting very different output, i dont know why). Lets set up the model:
clf = SVC(probability=True, degree=1, kernel='poly')to train the data into the model, we need to merge the data X Y column into an array with dimension of 2 with np.c_ function. Then we will be able to train the data using fit method. we can perform scoring for this model by using score function.
xy = np.c_[x,y]
clf.fit(xy,label)
clf.score(xy,label)Create Meshgrid that Covers Spatial Extent selected Area
Meshgrid is a method to create a 2d array that consist of xy location in the extend that we has been defined. See this documentation to get a better understanding what meshgrid is.
min_x = np.min(x)
min_y = np.min(y)max_x = np.max(x)
max_y = np.min(y)xx = np.linspace(min_x,max_x,100)
yy = np.linspace(min_y,max_y,100).T
xx, yy = np.meshgrid(xx, yy)
Now we have xx and yy variable that store the coordinate for our meshgrid. What actually is xx, yy store? We can visualize xx and yy by these following code:
plt.subplot(121)
plt.imshow(xx, origin='lower')
plt.subplot(122)
plt.imshow(yy, origin='lower')
Predict the and Visualize the Result
To predict the result, we need to merge and conver xx and yy array into an array with dimension of 2 with np.c_. Then we will be able to put them for prediction:
xx = np.linspace(0,1,100)
yy = np.linspace(0,1,100).T
xx, yy = np.meshgrid(xx, yy)
xxyy = np.c_[xx.flatten(),yy.flatten()]grid_shape = xx.shapepred = clf.predict(xxyy)
prob = clf.predict_proba(xxyy)
prediction value are stored in pred variable, and the predicton probability are stored in prob variable. To visualize the prediction result we can use:
plt.imshow(pred.reshape(grid_shape), extent=[0,1,0,1], origin='lower') #origin = ‘lower’ args added in order to inverse y cartesian space in imshow
plt.scatter(x,y,c=label, edgecolors='black')
prob variable store the probability value for each pixel and for each label. We will be able to visualize the porbability for each label by using imshow in prob data.
plt.figure(figsize=(10,10))
for i in np.unique(label):
plt.subplot(1,len(np.unique(label)),i+1)
plt.imshow(prob.T[i].reshape(grid_shape), extent=[0,1,0,1], origin=’lower’)
plt.scatter(x,y,c=label, edgecolors=’black’)
plt.title(‘Prob for Fac Code: ‘+ str(i))
Make the facies map “More Geology Looks”
OK.. now we have our map ready. But it doesnt looks geologicaly acceptable!! the change is too sharp, can we make it more “transitionnal”? sure. Here’s the code:
fac_transsitional = np.array([np.random.choice(a=(0,1,2),p=(i)) for i in prob])
plt.imshow(fac_transsitional.reshape(grid_shape), origin='lower')
Lets compare the histogram between Non Transsitional Prediction, Transsitional Pred, and original dataset:
the code for histogram generation are:
plt.hist((pred, fac_transsitional,label), bins=range(0,4), density=True, label = (“non Transsitional Prediction”,”Transsitional Prediction”, “Orignial Data”))
plt.xticks((0.5, 1.5, 2.5),(0,1,2))
plt.xlabel(‘Facies’)
plt.ylabel(‘Normalized Value’)
plt.legend()The Use Case
1 This method is suitable for the field where fulfil the following condition:
- Consist of a lot of well data, and each well already has facies label in it.
- In the area where facies band are exist. It is possible to use this method to map a spotted facies model such as carbonate system, deep water lobe system, but we need to fine tune the SVM kernel type, and the degree parameter.
2 This map can be used as conditioning map in making facies object modelling, using python of other software that available in market.
3 I recomend you not to use this map as final product. This not work of Magic. Always have gologist QC’d the result, ask their judgement wheteher the map has ‘geological make sense’ or not. In case it is not geologically appropriate,we need to Redraw the Map with paper and pencil (or other software that available on the market).
Happy coding. Im very happ to answer you Question..!!
