{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# DSI Summer Workshops Series \n",
"## July 19, 2018\n",
"\n",
"\n",
" Peggy Lindner
\n",
" Center for Advanced Computing & Data Science (CACDS)
\n",
" Data Science Institute (DSI)
\n",
" University of Houston
\n",
" plindner@uh.edu\n",
" \n",
" Please make sure you have Jupyterhub running and all the required python modules installed. Data for this and other tutorials can be found in the github repsoitory for the Summer 2018 DSI Workshops https://github.com/peggylind/Materials_Summer2018\n",
" \n",
" presenting:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Facies classification using Machine Learning\n",
"\n",
"#### Brendon Hall, [Enthought](https://www.enthought.com/)\n",
"\n",
"This notebook demonstrates how to train a machine learning algorithm to predict facies from well log data. The dataset we will use comes from a class excercise from The University of Kansas on [Neural Networks and Fuzzy Systems](http://www.people.ku.edu/~gbohling/EECS833/). This exercise is based on a consortium project to use machine learning techniques to create a reservoir model of the largest gas fields in North America, the Hugoton and Panoma Fields. For more info on the origin of the data, see [Bohling and Dubois (2003)](http://www.kgs.ku.edu/PRS/publication/2003/ofr2003-50.pdf) and [Dubois et al. (2007)](http://dx.doi.org/10.1016/j.cageo.2006.08.011). \n",
"\n",
"The dataset we will use is log data from nine wells that have been labeled with a facies type based on oberservation of core. We will use this log data to train a support vector machine to classify facies types. Support vector machines (or SVMs) are a type of supervised learning model that can be trained on data to perform classification and regression tasks. The SVM algorithm uses the training data to fit an optimal hyperplane between the different classes (or facies, in our case). We will use the SVM implementation in [scikit-learn](http://scikit-learn.org/stable/modules/svm.html).\n",
"\n",
"First we will [explore the dataset](#Exploring-the-dataset). We will load the training data from 9 wells, and take a look at what we have to work with. We will plot the data from a couple wells, and create cross plots to look at the variation within the data. \n",
"\n",
"Next we will [condition the data set](#Conditioning-the-data-set). We will remove the entries that have incomplete data. The data will be scaled to have zero mean and unit variance. We will also split the data into training and test sets.\n",
"\n",
"We will then be ready to [train the SVM classifier](#Training-the-SVM-classifier). We will demonstrate how to use the cross validation set to do [model parameter selection](#Model-parameter-selection).\n",
"\n",
"Finally, once we have a built and tuned the classifier, we can [apply the trained model](#Applying-the-classification-model-to-new-data) to classify facies in wells which do not already have labels. We will apply the classifier to two wells, but in principle you could apply the classifier to any number of wells that had the same log data."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Exploring the dataset\n",
"\n",
"First, we will examine the data set we will use to train the classifier. The training data is contained in the file `facies_vectors.csv`. The dataset consists of 5 wireline log measurements, two indicator variables and a facies label at half foot intervals. In machine learning terminology, each log measurement is a feature vector that maps a set of 'features' (the log measurements) to a class (the facies type). We will use the pandas library to load the data into a dataframe, which provides a convenient data structure to work with well log data."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib inline\n",
"import pandas as pd\n",
"import numpy as np\n",
"import matplotlib as mpl\n",
"import matplotlib.pyplot as plt\n",
"import matplotlib.colors as colors\n",
"from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
"\n",
"from pandas import set_option\n",
"set_option(\"display.max_rows\", 10)\n",
"pd.options.mode.chained_assignment = None\n",
"\n",
"filename = 'dataJuly19th/facies_vectors.csv'\n",
"training_data = pd.read_csv(filename)\n",
"training_data"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Remove a single well to use as a blind test later."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"blind = training_data[training_data['Well Name'] == 'SHANKLE']\n",
"training_data = training_data[training_data['Well Name'] != 'SHANKLE']"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This data is from the Council Grove gas reservoir in Southwest Kansas. The Panoma Council Grove Field is predominantly a carbonate gas reservoir encompassing 2700 square miles in Southwestern Kansas. This dataset is from nine wells (with 4149 examples), consisting of a set of seven predictor variables and a rock facies (class) for each example vector and validation (test) data (830 examples from two wells) having the same seven predictor variables in the feature vector. Facies are based on examination of cores from nine wells taken vertically at half-foot intervals. Predictor variables include five from wireline log measurements and two geologic constraining variables that are derived from geologic knowledge. These are essentially continuous variables sampled at a half-foot sample rate. \n",
"\n",
"The seven predictor variables are:\n",
"* Five wire line log curves include [gamma ray](http://petrowiki.org/Gamma_ray_logs) (GR), [resistivity logging](http://petrowiki.org/Resistivity_and_spontaneous_%28SP%29_logging) (ILD_log10),\n",
"[photoelectric effect](http://www.glossary.oilfield.slb.com/en/Terms/p/photoelectric_effect.aspx) (PE), [neutron-density porosity difference and average neutron-density porosity](http://petrowiki.org/Neutron_porosity_logs) (DeltaPHI and PHIND). Note, some wells do not have PE.\n",
"* Two geologic constraining variables: nonmarine-marine indicator (NM_M) and relative position (RELPOS)\n",
"\n",
"The nine discrete facies (classes of rocks) are: \n",
"1. Nonmarine sandstone\n",
"2. Nonmarine coarse siltstone \n",
"3. Nonmarine fine siltstone \n",
"4. Marine siltstone and shale \n",
"5. Mudstone (limestone)\n",
"6. Wackestone (limestone)\n",
"7. Dolomite\n",
"8. Packstone-grainstone (limestone)\n",
"9. Phylloid-algal bafflestone (limestone)\n",
"\n",
"These facies aren't discrete, and gradually blend into one another. Some have neighboring facies that are rather close. Mislabeling within these neighboring facies can be expected to occur. The following table lists the facies, their abbreviated labels and their approximate neighbors.\n",
"\n",
"Facies |Label| Adjacent Facies\n",
":---: | :---: |:--:\n",
"1 |SS| 2\n",
"2 |CSiS| 1,3\n",
"3 |FSiS| 2\n",
"4 |SiSh| 5\n",
"5 |MS| 4,6\n",
"6 |WS| 5,7\n",
"7 |D| 6,8\n",
"8 |PS| 6,7,9\n",
"9 |BS| 7,8\n",
"\n",
"Let's clean up this dataset. The 'Well Name' and 'Formation' columns can be turned into a categorical data type. "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"training_data['Well Name'] = training_data['Well Name'].astype('category')\n",
"training_data['Formation'] = training_data['Formation'].astype('category')\n",
"training_data['Well Name'].unique()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"These are the names of the 10 training wells in the Council Grove reservoir. Data has been recruited into pseudo-well 'Recruit F9' to better represent facies 9, the Phylloid-algal bafflestone. \n",
"\n",
"Before we plot the well data, let's define a color map so the facies are represented by consistent color in all the plots in this tutorial. We also create the abbreviated facies labels, and add those to the `facies_vectors` dataframe."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# 1=sandstone 2=c_siltstone 3=f_siltstone \n",
"# 4=marine_silt_shale 5=mudstone 6=wackestone 7=dolomite\n",
"# 8=packstone 9=bafflestone\n",
"facies_colors = ['#F4D03F', '#F5B041','#DC7633','#6E2C00',\n",
" '#1B4F72','#2E86C1', '#AED6F1', '#A569BD', '#196F3D']\n",
"\n",
"facies_labels = ['SS', 'CSiS', 'FSiS', 'SiSh', 'MS',\n",
" 'WS', 'D','PS', 'BS']\n",
"#facies_color_map is a dictionary that maps facies labels\n",
"#to their respective colors\n",
"facies_color_map = {}\n",
"for ind, label in enumerate(facies_labels):\n",
" facies_color_map[label] = facies_colors[ind]\n",
"\n",
"def label_facies(row, labels):\n",
" return labels[ row['Facies'] -1]\n",
" \n",
"training_data.loc[:,'FaciesLabels'] = training_data.apply(lambda row: label_facies(row, facies_labels), axis=1)\n",
"training_data.describe()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This is a quick view of the statistical distribution of the input variables. Looking at the `count` values, most values have 4149 valid values except for `PE`, which has 3232. In this tutorial we will drop the feature vectors that don't have a valid `PE` entry."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"PE_mask = training_data['PE'].notnull().values\n",
"training_data = training_data[PE_mask]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's take a look at the data from individual wells in a more familiar log plot form. We will create plots for the five well log variables, as well as a log for facies labels. The plots are based on the those described in Alessandro Amato del Monte's [excellent tutorial](https://github.com/seg/tutorials/tree/master/1504_Seismic_petrophysics_1)."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def make_facies_log_plot(logs, facies_colors):\n",
" #make sure logs are sorted by depth\n",
" logs = logs.sort_values(by='Depth')\n",
" cmap_facies = colors.ListedColormap(\n",
" facies_colors[0:len(facies_colors)], 'indexed')\n",
" \n",
" ztop=logs.Depth.min(); zbot=logs.Depth.max()\n",
" \n",
" cluster=np.repeat(np.expand_dims(logs['Facies'].values,1), 100, 1)\n",
" \n",
" f, ax = plt.subplots(nrows=1, ncols=6, figsize=(8, 12))\n",
" ax[0].plot(logs.GR, logs.Depth, '-g')\n",
" ax[1].plot(logs.ILD_log10, logs.Depth, '-')\n",
" ax[2].plot(logs.DeltaPHI, logs.Depth, '-', color='0.5')\n",
" ax[3].plot(logs.PHIND, logs.Depth, '-', color='r')\n",
" ax[4].plot(logs.PE, logs.Depth, '-', color='black')\n",
" im=ax[5].imshow(cluster, interpolation='none', aspect='auto',\n",
" cmap=cmap_facies,vmin=1,vmax=9)\n",
" \n",
" divider = make_axes_locatable(ax[5])\n",
" cax = divider.append_axes(\"right\", size=\"20%\", pad=0.05)\n",
" cbar=plt.colorbar(im, cax=cax)\n",
" cbar.set_label((17*' ').join([' SS ', 'CSiS', 'FSiS', \n",
" 'SiSh', ' MS ', ' WS ', ' D ', \n",
" ' PS ', ' BS ']))\n",
" cbar.set_ticks(range(0,1)); cbar.set_ticklabels('')\n",
" \n",
" for i in range(len(ax)-1):\n",
" ax[i].set_ylim(ztop,zbot)\n",
" ax[i].invert_yaxis()\n",
" ax[i].grid()\n",
" ax[i].locator_params(axis='x', nbins=3)\n",
" \n",
" ax[0].set_xlabel(\"GR\")\n",
" ax[0].set_xlim(logs.GR.min(),logs.GR.max())\n",
" ax[1].set_xlabel(\"ILD_log10\")\n",
" ax[1].set_xlim(logs.ILD_log10.min(),logs.ILD_log10.max())\n",
" ax[2].set_xlabel(\"DeltaPHI\")\n",
" ax[2].set_xlim(logs.DeltaPHI.min(),logs.DeltaPHI.max())\n",
" ax[3].set_xlabel(\"PHIND\")\n",
" ax[3].set_xlim(logs.PHIND.min(),logs.PHIND.max())\n",
" ax[4].set_xlabel(\"PE\")\n",
" ax[4].set_xlim(logs.PE.min(),logs.PE.max())\n",
" ax[5].set_xlabel('Facies')\n",
" \n",
" ax[1].set_yticklabels([]); ax[2].set_yticklabels([]); ax[3].set_yticklabels([])\n",
" ax[4].set_yticklabels([]); ax[5].set_yticklabels([])\n",
" ax[5].set_xticklabels([])\n",
" f.suptitle('Well: %s'%logs.iloc[0]['Well Name'], fontsize=14,y=0.94)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Placing the log plotting code in a function will make it easy to plot the logs from multiples wells, and can be reused later to view the results when we apply the facies classification model to other wells. The function was written to take a list of colors and facies labels as parameters. \n",
"\n",
"We then show log plots for wells `SHRIMPLIN`. "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"make_facies_log_plot(\n",
" training_data[training_data['Well Name'] == 'SHRIMPLIN'],\n",
" facies_colors)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In addition to individual wells, we can look at how the various facies are represented by the entire training set. Let's plot a histgram of the number of training examples for each facies class."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"#count the number of unique entries for each facies, sort them by\n",
"#facies number (instead of by number of entries)\n",
"facies_counts = training_data['Facies'].value_counts().sort_index()\n",
"#use facies labels to index each count\n",
"facies_counts.index = facies_labels\n",
"\n",
"facies_counts.plot(kind='bar',color=facies_colors, \n",
" title='Distribution of Training Data by Facies')\n",
"facies_counts"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This shows the distribution of examples by facies for the 3232 training examples in the training set. Dolomite (facies 7) has the fewest with 141 examples. There are also only 185 bafflestone examples. Depending on the performance of the classifier we are going to train, we may consider getting more examples of these facies.\n",
"\n",
"Crossplots are a familiar tool in the geosciences to visualize how two properties vary with rock type. This dataset contains 5 log variables, and scatter matrix can help to quickly visualize the variation between the all the variables in the dataset. We can employ the very useful [Seaborn library](https://stanford.edu/~mwaskom/software/seaborn/) to quickly create a nice looking scatter matrix. Each pane in the plot shows the relationship between two of the variables on the x and y axis, with each point is colored according to its facies. The same colormap is used to represent the 9 facies. "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"scrolled": false
},
"outputs": [],
"source": [
"#save plot display settings to change back to when done plotting with seaborn\n",
"inline_rc = dict(mpl.rcParams)\n",
"\n",
"import seaborn as sns\n",
"sns.set()\n",
"sns.pairplot(training_data.drop(['Well Name','Facies','Formation','Depth','NM_M','RELPOS'],axis=1),\n",
" hue='FaciesLabels', palette=facies_color_map,\n",
" hue_order=list(reversed(facies_labels)))\n",
"\n",
"#switch back to default matplotlib plot style\n",
"mpl.rcParams.update(inline_rc)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Conditioning the data set\n",
"\n",
"Now we extract just the feature variables we need to perform the classification. The predictor variables are the five wireline values and two geologic constraining variables. We also get a vector of the facies labels that correspond to each feature vector."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"scrolled": true
},
"outputs": [],
"source": [
"correct_facies_labels = training_data['Facies'].values\n",
"\n",
"feature_vectors = training_data.drop(['Formation', 'Well Name', 'Depth','Facies','FaciesLabels'], axis=1)\n",
"feature_vectors.describe()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Scikit includes a [preprocessing](http://scikit-learn.org/stable/modules/preprocessing.html) module that can 'standardize' the data (giving each variable zero mean and unit variance, also called *whitening*). Many machine learning algorithms assume features will be standard normally distributed data (ie: Gaussian with zero mean and unit variance). The factors used to standardize the training set must be applied to any subsequent feature set that will be input to the classifier. The `StandardScalar` class can be fit to the training set, and later used to standardize any training data."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from sklearn import preprocessing\n",
"\n",
"scaler = preprocessing.StandardScaler().fit(feature_vectors)\n",
"scaled_features = scaler.transform(feature_vectors)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"feature_vectors"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Scikit also includes a handy function to randomly split the training data into training and test sets. The test set contains a small subset of feature vectors that are not used to train the network. Because we know the true facies labels for these examples, we can compare the results of the classifier to the actual facies and determine the accuracy of the model. Let's use 20% of the data for the test set."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from sklearn.cross_validation import train_test_split\n",
"\n",
"X_train, X_test, y_train, y_test = train_test_split(\n",
" scaled_features, correct_facies_labels, test_size=0.2, random_state=42)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Training the SVM classifier\n",
"\n",
"Now we use the cleaned and conditioned training set to create a facies classifier. As mentioned above, we will use a type of machine learning model known as a [support vector machine](https://en.wikipedia.org/wiki/Support_vector_machine). The SVM is a map of the feature vectors as points in a multi dimensional space, mapped so that examples from different facies are divided by a clear gap that is as wide as possible. \n",
"\n",
"The SVM implementation in [scikit-learn](http://scikit-learn.org/stable/modules/generated/sklearn.svm.SVC.html#sklearn.svm.SVC) takes a number of important parameters. First we create a classifier using the default settings. "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from sklearn import svm\n",
"\n",
"clf = svm.SVC()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we can train the classifier using the training set we created above."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"clf.fit(X_train,y_train)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now that the model has been trained on our data, we can use it to predict the facies of the feature vectors in the test set. Because we know the true facies labels of the vectors in the test set, we can use the results to evaluate the accuracy of the classifier."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"predicted_labels = clf.predict(X_test)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We need some metrics to evaluate how good our classifier is doing. A [confusion matrix](http://www.dataschool.io/simple-guide-to-confusion-matrix-terminology/) is a table that can be used to describe the performance of a classification model. [Scikit-learn](http://scikit-learn.org/stable/modules/generated/sklearn.metrics.confusion_matrix.html) allows us to easily create a confusion matrix by supplying the actual and predicted facies labels.\n",
"\n",
"The confusion matrix is simply a 2D array. The entries of confusion matrix `C[i][j]` are equal to the number of observations predicted to have facies `j`, but are known to have facies `i`. \n",
"\n",
"To simplify reading the confusion matrix, a function has been written to display the matrix along with facies labels and various error metrics. See the file `classification_utilities.py` in this repo for the `display_cm()` function."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from sklearn.metrics import confusion_matrix\n",
"from classification_utilities import display_cm, display_adj_cm\n",
"\n",
"conf = confusion_matrix(y_test, predicted_labels)\n",
"display_cm(conf, facies_labels, hide_zeros=True)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The rows of the confusion matrix correspond to the actual facies labels. The columns correspond to the labels assigned by the classifier. For example, consider the first row. For the feature vectors in the test set that actually have label `SS`, 23 were correctly indentified as `SS`, 21 were classified as `CSiS` and 2 were classified as `FSiS`.\n",
"\n",
"The entries along the diagonal are the facies that have been correctly classified. Below we define two functions that will give an overall value for how the algorithm is performing. The accuracy is defined as the number of correct classifications divided by the total number of classifications."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def accuracy(conf):\n",
" total_correct = 0.\n",
" nb_classes = conf.shape[0]\n",
" for i in np.arange(0,nb_classes):\n",
" total_correct += conf[i][i]\n",
" acc = total_correct/sum(sum(conf))\n",
" return acc"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As noted above, the boundaries between the facies classes are not all sharp, and some of them blend into one another. The error within these 'adjacent facies' can also be calculated. We define an array to represent the facies adjacent to each other. For facies label `i`, `adjacent_facies[i]` is an array of the adjacent facies labels."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"adjacent_facies = np.array([[1], [0,2], [1], [4], [3,5], [4,6,7], [5,7], [5,6,8], [6,7]])\n",
"\n",
"def accuracy_adjacent(conf, adjacent_facies):\n",
" nb_classes = conf.shape[0]\n",
" total_correct = 0.\n",
" for i in np.arange(0,nb_classes):\n",
" total_correct += conf[i][i]\n",
" for j in adjacent_facies[i]:\n",
" total_correct += conf[i][j]\n",
" return total_correct / sum(sum(conf))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"print('Facies classification accuracy = %f' % accuracy(conf))\n",
"print('Adjacent facies classification accuracy = %f' % accuracy_adjacent(conf, adjacent_facies))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Model parameter selection\n",
"\n",
"The classifier so far has been built with the default parameters. However, we may be able to get improved classification results with optimal parameter choices.\n",
"\n",
"We will consider two parameters. The parameter `C` is a regularization factor, and tells the classifier how much we want to avoid misclassifying training examples. A large value of C will try to correctly classify more examples from the training set, but if `C` is too large it may 'overfit' the data and fail to generalize when classifying new data. If `C` is too small then the model will not be good at fitting outliers and will have a large error on the training set.\n",
"\n",
"The SVM learning algorithm uses a kernel function to compute the distance between feature vectors. Many kernel functions exist, but in this case we are using the radial basis function `rbf` kernel (the default). The `gamma` parameter describes the size of the radial basis functions, which is how far away two vectors in the feature space need to be to be considered close.\n",
"\n",
"We will train a series of classifiers with different values for `C` and `gamma`. Two nested loops are used to train a classifier for every possible combination of values in the ranges specified. The classification accuracy is recorded for each combination of parameter values. The results are shown in a series of plots, so the parameter values that give the best classification accuracy on the test set can be selected.\n",
"\n",
"This process is also known as 'cross validation'. Often a separate 'cross validation' dataset will be created in addition to the training and test sets to do model selection. For this tutorial we will just use the test set to choose model parameters."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"#model selection takes a few minutes, change this variable\n",
"#to true to run the parameter loop\n",
"do_model_selection = True\n",
"\n",
"if do_model_selection:\n",
" C_range = np.array([.01, 1, 5, 10, 20, 50, 100, 1000, 5000, 10000])\n",
" gamma_range = np.array([0.0001, 0.001, 0.01, 0.1, 1, 10])\n",
" \n",
" fig, axes = plt.subplots(3, 2, \n",
" sharex='col', sharey='row',figsize=(10,10))\n",
" plot_number = 0\n",
" for outer_ind, gamma_value in enumerate(gamma_range):\n",
" row = int(plot_number / 2)\n",
" column = int(plot_number % 2)\n",
" cv_errors = np.zeros(C_range.shape)\n",
" train_errors = np.zeros(C_range.shape)\n",
" for index, c_value in enumerate(C_range):\n",
" \n",
" clf = svm.SVC(C=c_value, gamma=gamma_value)\n",
" clf.fit(X_train,y_train)\n",
" \n",
" train_conf = confusion_matrix(y_train, clf.predict(X_train))\n",
" cv_conf = confusion_matrix(y_test, clf.predict(X_test))\n",
" \n",
" cv_errors[index] = accuracy(cv_conf)\n",
" train_errors[index] = accuracy(train_conf)\n",
"\n",
" ax = axes[row, column]\n",
" ax.set_title('Gamma = %g'%gamma_value)\n",
" ax.semilogx(C_range, cv_errors, label='CV error')\n",
" ax.semilogx(C_range, train_errors, label='Train error')\n",
" plot_number += 1\n",
" ax.set_ylim([0.2,1])\n",
" \n",
" ax.legend(bbox_to_anchor=(1.05, 0), loc='lower left', borderaxespad=0.)\n",
" fig.text(0.5, 0.03, 'C value', ha='center',\n",
" fontsize=14)\n",
" \n",
" fig.text(0.04, 0.5, 'Classification Accuracy', va='center', \n",
" rotation='vertical', fontsize=14)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The best accuracy on the cross validation error curve was achieved for `gamma = 1`, and `C = 10`. We can now create and train an optimized classifier based on these parameters:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"clf = svm.SVC(C=10, gamma=1) \n",
"clf.fit(X_train, y_train)\n",
"\n",
"cv_conf = confusion_matrix(y_test, clf.predict(X_test))\n",
"\n",
"print('Optimized facies classification accuracy = %.2f' % accuracy(cv_conf))\n",
"print('Optimized adjacent facies classification accuracy = %.2f' % accuracy_adjacent(cv_conf, adjacent_facies))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"[Precision and recall](https://en.wikipedia.org/wiki/Precision_and_recall) are metrics that give more insight into how the classifier performs for individual facies. Precision is the probability that given a classification result for a sample, the sample actually belongs to that class. Recall is the probability that a sample will be correctly classified for a given class.\n",
"\n",
"Precision and recall can be computed easily using the confusion matrix. The code to do so has been added to the `display_confusion_matrix()` function:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"display_cm(cv_conf, facies_labels, \n",
" display_metrics=True, hide_zeros=True)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To interpret these results, consider facies `SS`. In our test set, if a sample was labeled `SS` the probability the sample was correct is 0.8 (precision). If we know a sample has facies `SS`, then the probability it will be correctly labeled by the classifier is 0.78 (recall). It is desirable to have high values for both precision and recall, but often when an algorithm is tuned to increase one, the other decreases. The [F1 score](https://en.wikipedia.org/wiki/Precision_and_recall#F-measure) combines both to give a single measure of relevancy of the classifier results.\n",
"\n",
"These results can help guide intuition for how to improve the classifier results. For example, for a sample with facies `MS` or mudstone, it is only classified correctly 57% of the time (recall). Perhaps this could be improved by introducing more training samples. Sample quality could also play a role. Facies `BS` or bafflestone has the best `F1` score and relatively few training examples. But this data was handpicked from other wells to provide training examples to identify this facies.\n",
"\n",
"We can also consider the classification metrics when we consider misclassifying an adjacent facies as correct:\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"display_adj_cm(cv_conf, facies_labels, adjacent_facies, \n",
" display_metrics=True, hide_zeros=True)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Considering adjacent facies, the `F1` scores for all facies types are above 0.9, except when classifying `SiSh` or marine siltstone and shale. The classifier often misclassifies this facies (recall of 0.66), most often as wackestone. \n",
"\n",
"These results are comparable to those reported in Dubois et al. (2007)."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Applying the classification model to the blind data\n",
"\n",
"We held a well back from the training, and stored it in a dataframe called `blind`:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"blind"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The label vector is just the `Facies` column:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"y_blind = blind['Facies'].values"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can form the feature matrix by dropping some of the columns and making a new dataframe:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"well_features = blind.drop(['Facies', 'Formation', 'Well Name', 'Depth'], axis=1)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we can transform this with the scaler we made before:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"X_blind = scaler.transform(well_features)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now it's a simple matter of making a prediction and storing it back in the dataframe:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"y_pred = clf.predict(X_blind)\n",
"blind['Prediction'] = y_pred"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's see how we did with the confusion matrix:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"cv_conf = confusion_matrix(y_blind, y_pred)\n",
"\n",
"print('Optimized facies classification accuracy = %.2f' % accuracy(cv_conf))\n",
"print('Optimized adjacent facies classification accuracy = %.2f' % accuracy_adjacent(cv_conf, adjacent_facies))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We managed 0.75 using the test data, but it was from the same wells as the training data. This more reasonable test does not perform as well..."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"display_cm(cv_conf, facies_labels,\n",
" display_metrics=True, hide_zeros=True)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"...but does remarkably well on the adjacent facies predictions. "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"display_adj_cm(cv_conf, facies_labels, adjacent_facies,\n",
" display_metrics=True, hide_zeros=True)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def compare_facies_plot(logs, compadre, facies_colors):\n",
" #make sure logs are sorted by depth\n",
" logs = logs.sort_values(by='Depth')\n",
" cmap_facies = colors.ListedColormap(\n",
" facies_colors[0:len(facies_colors)], 'indexed')\n",
" \n",
" ztop=logs.Depth.min(); zbot=logs.Depth.max()\n",
" \n",
" cluster1 = np.repeat(np.expand_dims(logs['Facies'].values,1), 100, 1)\n",
" cluster2 = np.repeat(np.expand_dims(logs[compadre].values,1), 100, 1)\n",
" \n",
" f, ax = plt.subplots(nrows=1, ncols=7, figsize=(9, 12))\n",
" ax[0].plot(logs.GR, logs.Depth, '-g')\n",
" ax[1].plot(logs.ILD_log10, logs.Depth, '-')\n",
" ax[2].plot(logs.DeltaPHI, logs.Depth, '-', color='0.5')\n",
" ax[3].plot(logs.PHIND, logs.Depth, '-', color='r')\n",
" ax[4].plot(logs.PE, logs.Depth, '-', color='black')\n",
" im1 = ax[5].imshow(cluster1, interpolation='none', aspect='auto',\n",
" cmap=cmap_facies,vmin=1,vmax=9)\n",
" im2 = ax[6].imshow(cluster2, interpolation='none', aspect='auto',\n",
" cmap=cmap_facies,vmin=1,vmax=9)\n",
" \n",
" divider = make_axes_locatable(ax[6])\n",
" cax = divider.append_axes(\"right\", size=\"20%\", pad=0.05)\n",
" cbar=plt.colorbar(im2, cax=cax)\n",
" cbar.set_label((17*' ').join([' SS ', 'CSiS', 'FSiS', \n",
" 'SiSh', ' MS ', ' WS ', ' D ', \n",
" ' PS ', ' BS ']))\n",
" cbar.set_ticks(range(0,1)); cbar.set_ticklabels('')\n",
" \n",
" for i in range(len(ax)-2):\n",
" ax[i].set_ylim(ztop,zbot)\n",
" ax[i].invert_yaxis()\n",
" ax[i].grid()\n",
" ax[i].locator_params(axis='x', nbins=3)\n",
" \n",
" ax[0].set_xlabel(\"GR\")\n",
" ax[0].set_xlim(logs.GR.min(),logs.GR.max())\n",
" ax[1].set_xlabel(\"ILD_log10\")\n",
" ax[1].set_xlim(logs.ILD_log10.min(),logs.ILD_log10.max())\n",
" ax[2].set_xlabel(\"DeltaPHI\")\n",
" ax[2].set_xlim(logs.DeltaPHI.min(),logs.DeltaPHI.max())\n",
" ax[3].set_xlabel(\"PHIND\")\n",
" ax[3].set_xlim(logs.PHIND.min(),logs.PHIND.max())\n",
" ax[4].set_xlabel(\"PE\")\n",
" ax[4].set_xlim(logs.PE.min(),logs.PE.max())\n",
" ax[5].set_xlabel('Facies')\n",
" ax[6].set_xlabel(compadre)\n",
" \n",
" ax[1].set_yticklabels([]); ax[2].set_yticklabels([]); ax[3].set_yticklabels([])\n",
" ax[4].set_yticklabels([]); ax[5].set_yticklabels([])\n",
" ax[5].set_xticklabels([])\n",
" ax[6].set_xticklabels([])\n",
" f.suptitle('Well: %s'%logs.iloc[0]['Well Name'], fontsize=14,y=0.94)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"compare_facies_plot(blind, 'Prediction', facies_colors)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Applying the classification model to new data\n",
"\n",
"Now that we have a trained facies classification model we can use it to identify facies in wells that do not have core data. In this case, we will apply the classifier to two wells, but we could use it on any number of wells for which we have the same set of well logs for input.\n",
"\n",
"This dataset is similar to the training data except it does not have facies labels. It is loaded into a dataframe called `test_data`."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"well_data = pd.read_csv('dataJuly19th/validation_data_nofacies.csv')\n",
"well_data['Well Name'] = well_data['Well Name'].astype('category')\n",
"well_features = well_data.drop(['Formation', 'Well Name', 'Depth'], axis=1)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The data needs to be scaled using the same constants we used for the training data."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"X_unknown = scaler.transform(well_features)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Finally we predict facies labels for the unknown data, and store the results in a `Facies` column of the `test_data` dataframe."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"#predict facies of unclassified data\n",
"y_unknown = clf.predict(X_unknown)\n",
"well_data['Facies'] = y_unknown\n",
"well_data"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"well_data['Well Name'].unique()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can use the well log plot to view the classification results along with the well logs."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"make_facies_log_plot(\n",
" well_data[well_data['Well Name'] == 'STUART'],\n",
" facies_colors=facies_colors)\n",
"\n",
"make_facies_log_plot(\n",
" well_data[well_data['Well Name'] == 'CRAWFORD'],\n",
" facies_colors=facies_colors)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Finally we can write out a csv file with the well data along with the facies classification results."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"well_data.to_csv('well_data_with_facies.csv')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## References\n",
"\n",
"Amato del Monte, A., 2015. Seismic Petrophysics: Part 1, *The Leading Edge*, 34 (4). [doi:10.1190/tle34040440.1](http://dx.doi.org/10.1190/tle34040440.1)\n",
"\n",
"Bohling, G. C., and M. K. Dubois, 2003. An Integrated Application of Neural Network and Markov Chain Techniques to Prediction of Lithofacies from Well Logs, *KGS Open-File Report* 2003-50, 6 pp. [pdf](http://www.kgs.ku.edu/PRS/publication/2003/ofr2003-50.pdf)\n",
"\n",
"Dubois, M. K., G. C. Bohling, and S. Chakrabarti, 2007, Comparison of four approaches to a rock facies classification problem, *Computers & Geosciences*, 33 (5), 599-617 pp. [doi:10.1016/j.cageo.2006.08.011](http://dx.doi.org/10.1016/j.cageo.2006.08.011)"
]
}
],
"metadata": {
"anaconda-cloud": {},
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.5"
}
},
"nbformat": 4,
"nbformat_minor": 1
}