Soil_Classification.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Soil Classification using Deep Learning\n",
    "\n",
    "Classification of soils can be an important capability in many fields including geotechnical engineering, agricultural science, and construction. While many classification schemes exist, a popular and intuitive way of classifying soils is based on soil grain size, such as gravel, silt, sand, etc. This is a rather simple way of classifying soil because there are firm cutoffs between each classification. For example, gravel is any grain over 2 mm in diameter, sand is between $\\frac{1}{16}$ mm to 2 mm, and silt/clay is smaller than $\\frac{1}{16}$ mm. See [additional information on soil classification](https://apmonitor.com/pds/index.php/Main/SoilClassification).\n",
    "\n",
    "<img align=left width=400px src='https://apmonitor.com/pds/uploads/Main/soil_classification.png'>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Objective\n",
    "In this exercise, computer vision makes a classification model of soils based on photos of gravel, sand, and silt. A convolutional neural network framework is used classify photos of soil.\n",
    "\n",
    "### Key Concepts:\n",
    "#### Data Augmentation\n",
    "Data Augmentation is a practice that is commonly used in image classification that enlarges the existing image dataset by creating new data sets from the existing data. This augmentation is most commonly seen in operations like shifting the image, mirroring the image horizontally or vertically, zooming in on the image, and so on. These augmentations then create new images and thus enhance the number of images to train on. Read more about this in an [article from Jason Brownlee](https://machinelearningmastery.com/how-to-configure-image-data-augmentation-when-training-deep-learning-neural-networks/).\n",
    "\n",
    "#### Convolutional Neural Network (CNN)\n",
    "A convolutional neural network (CNN) is a popular model for image classification. CNNs distinguish characteristics of an image by looping through the pixel values in an image and calculating the dot product of a those pixels with a filter `kernel` matrix. These kernel matrices accentuate different aspects of the image such as vertical and horizontal lines, curvature, etc. For more information on CNNs, see [APMonitor - Computer Vision with Deep Learning](https://apmonitor.com/pds/index.php/Main/VisionDeepLearning)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "training_data_directory = 'train'\n",
    "test_data_directory = 'test'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Setup"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Import the following Python modules. Remember to use ```pip``` to install any packages you are missing. For example, if you get an error: ```ModuleNotFoundError: No module named 'cv2'```, add a new cell by clicking the \"+\" button above and run the following command in another cell: ```!pip install opencv-python```."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "yF35Ku0Ms-O6"
   },
   "outputs": [],
   "source": [
    "import os\n",
    "import re\n",
    "import cv2\n",
    "import time\n",
    "import shutil\n",
    "import zipfile\n",
    "import urllib.request\n",
    "import numpy as np\n",
    "from PIL import Image\n",
    "from os import listdir\n",
    "from os.path import isfile, join\n",
    "from random import randrange\n",
    "import matplotlib.pyplot as plt\n",
    "from tensorflow.keras.models import load_model\n",
    "from tensorflow.keras.preprocessing import image\n",
    "from tensorflow.keras.preprocessing.image import ImageDataGenerator\n",
    "from tensorflow.keras.models import Sequential\n",
    "from tensorflow.keras.layers import Dense, Activation, Flatten\n",
    "from tensorflow.keras.layers import Conv2D, MaxPooling2D"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Data\n",
    "\n",
    "The data is separated into directories test and train. The test folder and train folder contain subdirectories corresponding to the possible soils. The labeled tree structure of the folders assists with marking the photos for training and testing."
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "├───test\n",
    "│   ├───Gravel\n",
    "│   ├───Sand\n",
    "│   └───Silt\n",
    "└───train\n",
    "    ├───Gravel\n",
    "    ├───Sand\n",
    "    └───Silt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# download soil_photos.zip\n",
    "file = 'soil_photos.zip'\n",
    "url = 'http://apmonitor.com/pds/uploads/Main/'+file\n",
    "urllib.request.urlretrieve(url, file)\n",
    "\n",
    "# extract archive and remove soil_photos.zip\n",
    "with zipfile.ZipFile(file, 'r') as zip_ref:\n",
    "    zip_ref.extractall('./')\n",
    "os.remove(file)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Import the data into the Python session. The first step is to process the images into a format that 1) makes the data readable to the model, and 2) provides more training material for the model to learn from. For example, the ```training_data_processor``` variable scales the data so that it can be a model input, but also takes each image and augments it so that the model can learn from multiple variations of the same picture. It flips it horizontally, rotates it, shifts it, and more so that the model learns from the soil photo rather than the orientation or size."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "executionInfo": {
     "elapsed": 23726,
     "status": "ok",
     "timestamp": 1618065873113,
     "user": {
      "displayName": "Peter Van Katwyk",
      "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14GjhknOaA2vSO9txrYBB1U01ff4hmNe3wj5xrjP9Nzg=s64",
      "userId": "16032387266776692701"
     },
     "user_tz": 360
    },
    "id": "topbCwEWsr3q",
    "outputId": "a4116caa-a435-41f2-d2d8-71ca1e307d8c"
   },
   "outputs": [],
   "source": [
    "# Initiate data processing tools\n",
    "training_data_processor = ImageDataGenerator(\n",
    "    rescale = 1./255,\n",
    "    horizontal_flip = True,\n",
    "    zoom_range = 0.2,\n",
    "    rotation_range = 10,\n",
    "    shear_range = 0.2,\n",
    "    height_shift_range = 0.1,\n",
    "    width_shift_range = 0.1\n",
    ")\n",
    "\n",
    "test_data_processor = ImageDataGenerator(rescale = 1./255)\n",
    "\n",
    "# Load data into Python\n",
    "training_data = training_data_processor.flow_from_directory(\n",
    "    training_data_directory,\n",
    "    target_size = (256, 256),\n",
    "    batch_size = 32,\n",
    "    class_mode = 'categorical',\n",
    ")\n",
    "\n",
    "testing_data = test_data_processor.flow_from_directory(\n",
    "    test_data_directory,\n",
    "    target_size = (256 ,256),\n",
    "    batch_size = 32,\n",
    "    class_mode = 'categorical',\n",
    "    shuffle = False\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Model Building"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Build the CNN model. The following cell sets parameters for model building. This includes the number of convolutional layers, fully connected dense layers, the number of nodes in each layer, and the number of training epochs. For more information on these parameters and Convolutional Neural Networks in general, see [APMonitor - Computer Vision with Deep Learning](https://apmonitor.com/pds/index.php/Main/VisionDeepLearning). The `layer_size` parameter has a large influence on the training speed, accuracy, and size of the model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "executionInfo": {
     "elapsed": 3288116,
     "status": "ok",
     "timestamp": 1618026147600,
     "user": {
      "displayName": "Peter Van Katwyk",
      "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14GjhknOaA2vSO9txrYBB1U01ff4hmNe3wj5xrjP9Nzg=s64",
      "userId": "16032387266776692701"
     },
     "user_tz": 360
    },
    "id": "9-jCm9Y3Ehbs",
    "outputId": "c99cfeb8-dd44-4da7-e0b7-e97ba1455875"
   },
   "outputs": [],
   "source": [
    "# choose model parameters\n",
    "num_conv_layers = 2\n",
    "num_dense_layers = 1\n",
    "layer_size = 32\n",
    "num_training_epochs = 20\n",
    "MODEL_NAME = 'soil'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initiate model variable\n",
    "model = Sequential()\n",
    "\n",
    "# begin adding properties to model variable\n",
    "# e.g. add a convolutional layer\n",
    "model.add(Conv2D(layer_size, (3, 3), input_shape=(256,256, 3)))\n",
    "model.add(Activation('relu'))\n",
    "model.add(MaxPooling2D(pool_size=(2, 2)))\n",
    "\n",
    "# add additional convolutional layers based on num_conv_layers\n",
    "for _ in range(num_conv_layers-1):\n",
    "    model.add(Conv2D(layer_size, (3, 3)))\n",
    "    model.add(Activation('relu'))\n",
    "    model.add(MaxPooling2D(pool_size=(2, 2)))\n",
    "\n",
    "# reduce dimensionality\n",
    "model.add(Flatten())\n",
    "\n",
    "# add fully connected \"dense\" layers if specified\n",
    "for _ in range(num_dense_layers):\n",
    "    model.add(Dense(layer_size))\n",
    "    model.add(Activation('relu'))\n",
    "\n",
    "# add output layer\n",
    "model.add(Dense(3))\n",
    "model.add(Activation('softmax'))\n",
    "\n",
    "# compile the sequential model with all added properties\n",
    "model.compile(loss='categorical_crossentropy',\n",
    "                optimizer='adam',\n",
    "                metrics=['accuracy'],\n",
    "                )\n",
    "\n",
    "# use the data already loaded previously to train/tune the model\n",
    "model.fit(training_data,\n",
    "            epochs=num_training_epochs,\n",
    "            validation_data = testing_data)\n",
    "\n",
    "# save the trained model\n",
    "model.save(f'{MODEL_NAME}.h5')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Model Testing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The model has now been trained and has been saved to the computer (in the same folder as this notebook). The last line of the printed output above contains the accuracy for both the training data and the validation, or testing data. Pay attention to the ```val_accuracy``` value on the last line, which corresponds to the accuracy on images that the model was not trained on. This is the best measure of the estimator performance on other images.\n",
    "\n",
    "The following function, ```make_prediction```, takes the file path to a soil photo as an input and outputs the classification that the model predicted. From the test folder, copy the filepath into the test_image_filepath variable in the second cell. Try out a few photos and see how the model preforms. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 286
    },
    "executionInfo": {
     "elapsed": 3080,
     "status": "ok",
     "timestamp": 1617992171243,
     "user": {
      "displayName": "Peter Van Katwyk",
      "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14GjhknOaA2vSO9txrYBB1U01ff4hmNe3wj5xrjP9Nzg=s64",
      "userId": "16032387266776692701"
     },
     "user_tz": 360
    },
    "id": "Chpl4ZrSPtHc",
    "outputId": "cb9bd334-7945-4a5d-d765-6abd8f0f0eec",
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "def make_prediction(image_fp):\n",
    "    im = cv2.imread(image_fp) # load image\n",
    "    plt.imshow(im[:,:,[2,1,0]])\n",
    "    img = image.load_img(image_fp, target_size = (256,256))\n",
    "    img = image.img_to_array(img)\n",
    "\n",
    "    image_array = img / 255. # scale the image\n",
    "    img_batch = np.expand_dims(image_array, axis = 0)\n",
    "    \n",
    "    class_ = [\"Gravel\", \"Sand\", \"Silt\"] # possible output values\n",
    "    predicted_value = class_[model.predict(img_batch).argmax()]\n",
    "    true_value = re.search(r'(Gravel)|(Sand)|(Silt)', image_fp)[0]\n",
    "    \n",
    "    out = f\"\"\"Predicted Soil Type: {predicted_value}\n",
    "    True Soil Type: {true_value}\n",
    "    Correct?: {predicted_value == true_value}\"\"\"\n",
    "    \n",
    "    return out"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "test_image_filepath = test_data_directory + r'/Sand/0.jpg'\n",
    "print(make_prediction(test_image_filepath))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Soil Classification Percentages"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Soils are not perfectly homogenous (all one soil type). Soils are often a blend of types and may be better represented using a percentage. For example, the cell below shows a test photo that is labeled \"Silt\":"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "percentage_photo = test_data_directory + r\"/Silt/5.jpg\"\n",
    "im = cv2.imread(percentage_photo) # load image\n",
    "plt.imshow(im[:,:,[2,1,0]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This photo has some gravel, however, there appears to be sand, gravel, and silt all in one photo. In order to better classify this image, create a proportion of each label. For example, this photo might be 30% Gravel, 20% Sand, and 50% Silt. \n",
    "\n",
    "Subdivide the photo into many smaller segments or squares and train a classifier on the smaller squares. Loop through the photo  to classify each small square and take a proportion of the squares pertaining to Gravel, Sand, and Silt. This proportion is then converted into a percentage of the respective soil type."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run these cells below to subdivide the existing training photos into smaller segments. The new directory will be called train_divided and test_divided.  \n",
    "\n",
    "Note: This will take time (~30 seconds to 2 minutes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def split_images(image_dir, save_dir):\n",
    "    classification_list = ['Gravel', 'Sand', 'Silt']\n",
    "    for classification in classification_list:\n",
    "        folder = image_dir + '/' + classification + '/'\n",
    "        save_folder = save_dir + '/' + classification + '/'\n",
    "        files = [f for f in listdir(folder) if isfile(join(folder, f))]\n",
    "\n",
    "        for file in files:\n",
    "            if '.ini' in file:\n",
    "                continue\n",
    "            fp = folder + file\n",
    "            img = cv2.imread(fp)\n",
    "            h,w,c = img.shape\n",
    "            im_dim = 64\n",
    "            # For cropping images\n",
    "            for r in range(0,img.shape[0],im_dim):\n",
    "                for c in range(0,img.shape[1],im_dim):\n",
    "                    cropped_img = img[r:r+im_dim, c:c+im_dim,:]\n",
    "                    ch, cw, cc = cropped_img.shape\n",
    "                    if ch == im_dim and cw == im_dim:\n",
    "                        write_path = f\"{save_folder + str(randrange(100000))}img{r}_{c}.jpg\"\n",
    "                        cv2.imwrite(write_path,cropped_img)\n",
    "                    else:\n",
    "                        pass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "try:\n",
    "    parent = training_data_directory.replace('train', '')\n",
    "    dirs = ['train_divided', 'test_divided']\n",
    "    class_ = [\"Gravel\", \"Sand\", \"Silt\"]\n",
    "    for dir in dirs:\n",
    "        os.mkdir(os.path.join(parent, dir))\n",
    "        for classification in class_:\n",
    "            os.mkdir(os.path.join(parent, dir, classification))\n",
    "\n",
    "    # split training images\n",
    "    split_images(image_dir=training_data_directory,\n",
    "                save_dir=training_data_directory.replace('train', 'train_divided'))\n",
    "    # split test images\n",
    "    split_images(image_dir=test_data_directory,\n",
    "                save_dir=test_data_directory.replace('test', 'test_divided'))\n",
    "except FileExistsError:\n",
    "    pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load Model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If the number of input pixels changes, a new model should be trained or the image rescaled to the original training dimensions. For the purposes of this demo, retraining is not needed because the photos are upscaled to 1024x1024 before subsampling 256x256 blocks (16) in each image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model_fp = os.getcwd()+'/'+'soil.h5'\n",
    "print(model_fp)\n",
    "model = load_model(model_fp)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Classify Image"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use the loaded model to classify a test image. The function `classify_images` takes an image and a model and loops through each 256x256 square. It classifies each square and adds to the counter to create the fractional soil prediction. The function outputs the proportion of each soil type that was classified."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def classify_images(image_fp, model):\n",
    "    classes = ['Gravel', 'Sand', 'Silt']\n",
    "    gravel_count = 0\n",
    "    sand_count = 0\n",
    "    silt_count = 0\n",
    "\n",
    "    img = cv2.imread(image_fp)\n",
    "    img = cv2.resize(img,(1024,1024))\n",
    "    im_dim = 256\n",
    "\n",
    "    for r in range(0, img.shape[0], im_dim):\n",
    "        for c in range(0, img.shape[1], im_dim):\n",
    "            cropped_img = img[r:r + im_dim, c:c + im_dim, :]\n",
    "            h, w, c = cropped_img.shape\n",
    "            if h == im_dim and w == im_dim:\n",
    "                classification = model_classify(cropped_img, model)\n",
    "                if classification == classes[0]:\n",
    "                    gravel_count += 1\n",
    "                elif classification == classes[1]:\n",
    "                    sand_count += 1\n",
    "                elif classification == classes[2]:\n",
    "                    silt_count += 1\n",
    "            else:\n",
    "                continue\n",
    "    total_count = gravel_count + sand_count + silt_count\n",
    "    proportion_array = [gravel_count / total_count, sand_count / total_count, silt_count / total_count]\n",
    "    return proportion_array\n",
    "\n",
    "\n",
    "def model_classify(cropped_img, model):\n",
    "    classes = ['Gravel', 'Sand', 'Silt']\n",
    "    image_array = cropped_img / 255.\n",
    "    img_batch = np.expand_dims(image_array, axis=0)\n",
    "    prediction_array = model.predict(img_batch)[0]\n",
    "    first_idx = np.argmax(prediction_array)\n",
    "    first_class = classes[first_idx]\n",
    "    return first_class\n",
    "\n",
    "def classify_percentage(image_fp):\n",
    "    start = time.time()\n",
    "    out = classify_images(image_fp=image_fp, model=model)\n",
    "    finish = str(round(time.time() - start, 5))\n",
    "    \n",
    "    im = cv2.imread(image_fp) # load image\n",
    "    plt.imshow(im[:,:,[2, 1, 0]])\n",
    "\n",
    "    print(f'''---\n",
    "Percent Gravel: {round(out[0] * 100, 2)}%)\n",
    "Percent Sand: {round(out[1] * 100, 2)}%)\n",
    "Percent Silt: {round(out[2] * 100, 2)}%)\n",
    "Time to Classify: {finish} seconds\n",
    "---''')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "classify_percentage(image_fp=percentage_photo)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Does the percentage classified match the photo?  \n",
    "\n",
    "Both approaches to soil classification have advantages. For example, classifying the entire image is faster. However, percentage classification is more computationally expensive but provides more information and may be a better representation of the true classification."
   ]
  }
 ],
 "metadata": {
  "accelerator": "TPU",
  "colab": {
   "authorship_tag": "ABX9TyNhoK7FQHGoNgSPJqWd8qt+",
   "collapsed_sections": [],
   "name": "SoilClassification.ipynb",
   "provenance": []
  },
  "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}