SVM_classification.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import all libraries\n",
    "import rasterio\n",
    "from rasterio.plot import show\n",
    "import numpy as np\n",
    "import numpy.ma as ma\n",
    "import matplotlib.pyplot as plt\n",
    "from glob import glob\n",
    "import os\n",
    "import sys"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Support Vector Machine Classification\n",
    "https://scikit-learn.org/stable/modules/svm.html#classification\n",
    "\n",
    "https://stackoverflow.com/questions/43331510/how-to-train-an-svm-classifier-on-a-satellite-image-using-python#43386604\n",
    "\n",
    "#### This code has been tested with Landsat 4 TM, Landsat 8 OLI, Landsat 7 ETM+, Sentinel 2 MSI, and Sentinel 3 OLCI optical data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### open and read input raster file\n",
    "ras = 'D:/your_input_directory/your_input_raster.tif'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### read opened raster as numpy array\n",
    "with rasterio.open(ras, 'r')  as src:\n",
    "    show(src, cmap='viridis') # to visualize the raster\n",
    "    ras_arr = src.read(1)      # to read raster as numpy array\n",
    "### these steps are not mendatory\n",
    "### but can run to see the info of input raster\n",
    "    print(ras_arr.shape)       # to see the number of rows and columns\n",
    "    print(ras_arr)             # to print the array\n",
    "    print(np.min(ras_arr))     # to print minimum value of the array\n",
    "    print(np.max(ras_arr))     # to print maximum value of the array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### plot array without geocoordinate for \n",
    "### selecting supervised classes' (Y1:Y2, X1:X2)\n",
    "## based on their array coordinate\n",
    "plt.imshow(ras_arr, cmap=\"nipy_spectral\")\n",
    "### plot colorbar\n",
    "plt.colorbar(shrink=0.5)\n",
    "plt.clim(np.min(ras_arr), np.max(ras_arr))\n",
    "plt.gcf().set_size_inches(9, 7)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### make the supervised classes\n",
    "### for test case five classes has been defined\n",
    "### you can change the name of each class in 'class*'\n",
    "### according to your needs\n",
    "### class dictionary with keys and values\n",
    "classes = {'class1': 0, 'class2': 1, 'class3': 2, 'class4': 3, 'class5': 4} \n",
    "n_classes = len(classes) # assign the lenght of class dictionary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### build training data with data from define classes\n",
    "### this train data will supervise input raster based on SVC algorithm\n",
    "ras_svm = n_classes*np.ones(shape=(ras_arr.shape), dtype=np.int)\n",
    "### fixed the area of each class on the raster\n",
    "### technically on array\n",
    "### y1,y2, and x1,x2 must be different for each class\n",
    "### precision of y1,y2, and x1,x2 control the outcome of SVC\n",
    "ras_svm[y1:y2, x1:x2] = classes['class1']\n",
    "ras_svm[y1:y2, x1:x2] = classes['class2']\n",
    "ras_svm[y1:y2, x1:x2] = classes['class3']\n",
    "ras_svm[y1:y2, x1:x2] = classes['class4']\n",
    "ras_svm[y1:y2, x1:x2] = classes['class5']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### fix the shape and meta info for training and test set\n",
    "X = ras_arr.reshape([-1, 1]) \n",
    "y = ras_svm.ravel()\n",
    "train = np.flatnonzero(ras_svm < n_classes)\n",
    "test = np.flatnonzero(ras_svm == n_classes)\n",
    "\n",
    "### call and run SVC algorithm\n",
    "from sklearn.svm import SVC\n",
    "\n",
    "clf = SVC(gamma='auto') # gamma can be changed, for details see the Scikit-learn website\n",
    "clf.fit(X[train], y[train])\n",
    "y[test] = clf.predict(X[test])\n",
    "ras_svm = y.reshape(ras_arr.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### plot SVC array for quick look\n",
    "plt.imshow(ras_svm, cmap=\"Paired\") # you can use any cmap from matplotlib cmap list\n",
    "### plot colorbar\n",
    "plt.colorbar(shrink=0.5)\n",
    "plt.clim(np.min(ras_svm), np.max(ras_svm))\n",
    "print(ras_svm.shape)\n",
    "print(ras_svm)\n",
    "print(np.min(ras_svm))\n",
    "print(np.max(ras_svm))\n",
    "\n",
    "### save array as jpg file\n",
    "### make the output file directory\n",
    "fig_direc = 'D:/your_output_directory/CI'\n",
    "if not os.path.exists(fig_direc):\n",
    "    os.makedirs(fig_direc)\n",
    "### Save the figure as jpg format\n",
    "plt.savefig('D:/your_output_directory/your_output_raster.jpg', dpi=600) # dpi can be changed according to \n",
    "# needs but 600 is high enough and such good resolution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### save the output as GeoTiff file\n",
    "### Get metadata from input raster\n",
    "with rasterio.open(ras) as src:\n",
    "    meta = src.meta\n",
    "meta.update(dtype=rasterio.float32)\n",
    "### Create output folder and write output file in it as geotiff\n",
    "file = \"D:/your_output_directory/your_output_raster.tif\"\n",
    "os.makedirs(os.path.dirname(file), exist_ok=True)\n",
    "with rasterio.open(file, 'w', **meta) as dst:\n",
    "    dst.write(ras_svm.astype(rasterio.float32), 1)"
   ]
  }
 ],
 "metadata": {
  "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": 4
}