Final project.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tree cover classification on Azuero Peninsula, Panama <br>\n",
    "## GEOS 505: Research Computing in Earth Sciences <br>\n",
    "### Cristina Barber <br>\n",
    "#### Introduction and Problem Statement<br>\n",
    "\n",
    "Second-growth forests can provide ecosystem services including carbon storage, biodiversity conservation, and human livelihood. Consequently, restoring ecosystem function to millions of hectares of deforested land is a global conservation priority and forest restoration at landscape-scales have increasing prominence in international and national environmental policies. These large-scale goals are complicated by uncertainty in forest recovery rates across heterogeneous landscapes. For example, some sites may require active land management to avoid stagnation of succession (forest growth process), while others may not. This unpredictability of successional trajectories means that restoration is often an expensive investment with high failure rates. Predictive models could enable ecologists to forecast successional trajectories over large areas and long time periods, thus improving restoration success. \n",
    "Seed sources is a good example of an important common landscape feature which abundance and composition are heterogeneous across the landscape and with a high predictive capacity on recruitment. A widely used variable representing seed sources and with important implications for microclimatic and soil conditions is forest cover.  Still, studies looking at distance to forest cover have obtained a contradictory effect of this variable on recruitment. Whereas some studies have found forest cover to be a good predictor for forest recovery that can improve the diversity and structure characteristics of secondary forest others have not found significant effects. Chapter 1 of my thesis objective is to predict recruitment based on remote sensing imagery that can capture landscape characteristics. Since the surrounding forest could have great effects on recruitment, it is a variable that I would like to include in my model. Therefore, the objective of this project is to obtain a layer of forest cover on 23.000ha that can be included in the model. To achieve this goal, I will use June 05, 2018 data from cube satellite and Normalized Difference Vegetation Index (NDVI) data from Landsat 8 for this project. Using the information from both data layers I parameterized a model capable of detecting the forest cover in each pixel using NDVI. I expect that pixels with higher NDVI values will have higher forest cover percentage.<br>\n",
    "#### Data and Methods <br>\n",
    "*Data*<br>\n",
    "Cube satellite data belongs to the planet project and has a 3 m resolution in the visual spectra that will allow classifying forest and non-forest areas. NDVI data will be a Landsat image obtained from google earth engine at 30 m pixel spatial resolution. Since cube spatial resolution is very accurate, classifying the forest cover percentage for the 23.000 ha of my study area will be very computationally intensive, so I will create a model from a smaller sample of the area to upscale forest percentage cover to the whole landscape.\n",
    "<br>\n",
    "*Methods*<br>\n",
    "First we imported the needed libraries into python and set the working directory\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import os\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "import cv2\n",
    "import gdal\n",
    "from sklearn.cluster import KMeans\n",
    "from sklearn import linear_model\n",
    "from sklearn.metrics import mean_squared_error, r2_score\n",
    "from sklearn.metrics import explained_variance_score\n",
    "os.getcwd()\n",
    "os.chdir('C:\\\\Users\\\\Cristina\\\\Documents\\\\planet_order_197131\\\\NDVi_clipaligned')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After importing the needed libraries we will load the NDVI image and explore it"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "NDVI = cv2.imread('NDVItrue_instideplanet.tif',-1)\n",
    "cv2.imshow('image',NDVI)\n",
    "cv2.waitKey(0)\n",
    "cv2.destroyAllWindows()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and then load the CubeSAt image and rearange the layers to match a better structure to work with them and do the clustering."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gfid = gdal.Open('planet_insideplanet2.tif')\n",
    "planet = gfid.ReadAsArray()\n",
    "planet = np.float32(planet)\n",
    "\n",
    "planet = np.swapaxes(planet,0,2)\n",
    "planet = np.swapaxes(planet,0,1)\n",
    "\n",
    "small = cv2.resize(planet, (0,0), fx=0.2, fy=0.2)\n",
    "cv2.imshow('indice 0',small[:,:,:3]*2)\n",
    "cv2.waitKey(0)\n",
    "cv2.destroyAllWindows()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next step is to  randomly selected 100 pixels from the NDVI image and its overlapping 10x10 pixel arrays from CubeSat. In order to avoid the random point to fall into water, explored the image to obtain which is the rwo of pixels closer to the water but still on land to explude the next rwos from the random selection. The random points were selected on the NDVI image and then the 10 x 10 pixels array from the CubeSat image locations matched."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Visualization of the NDVI image to look for the threshold\n",
    "NDVI_copy = cv2.cvtColor(NDVI.copy(), cv2.COLOR_GRAY2BGR)\n",
    "\n",
    "NDVI_copy[342,:] = [122,0,0] \n",
    "\n",
    "NDVI_copy = cv2.resize(NDVI_copy, (0,0), fx=2, fy=2)\n",
    "\n",
    "cv2.imshow('image',NDVI_copy)\n",
    "cv2.waitKey(0)\n",
    "cv2.destroyAllWindows()\n",
    "# Select the random point in the NDVI image and extract irs values\n",
    "size = 100\n",
    "x=np.random.randint(low=0, high=342, size=size)\n",
    "y=np.random.randint(low=0, high=NDVI.shape[1], size=size)\n",
    "\n",
    "for index in range(x.shape[0]):\n",
    "    print('index: {}, x: {}, y: {}'.format(index,x[index],y[index]))\n",
    "\n",
    "NDVI_training=list()\n",
    "for i in range(x.shape[0]):\n",
    "    NDVI_training.append(NDVI[x[i],y[i]])\n",
    "\n",
    "    print(NDVI_training[-1])\n",
    "    \n",
    "# Locate the CubeSat pixels arrays and extract the values of the bands\n",
    "x_planet_i=x*10\n",
    "x_planet_e=x_planet_i+10\n",
    "y_planet_i=y*10\n",
    "y_planet_e=y_planet_i+10\n",
    "\n",
    "planet_training=list()\n",
    "for i in range(x.shape[0]):\n",
    "    planet_training.append(planet[x_planet_i[i]:x_planet_e[i],y_planet_i[i]:y_planet_e[i]])\n",
    "\n",
    "NDVI_copy = cv2.cvtColor(NDVI.copy(), cv2.COLOR_GRAY2BGR)\n",
    "\n",
    "for indox in range(x.shape[0]):\n",
    "    NDVI_copy[x[index],y[index]] = [0,0,255] \n",
    "\n",
    "# Show the selected pixels on the NDVI image\n",
    "NDVI_copy[340,:] = [122,0,0] \n",
    "\n",
    "NDVI_copy = cv2.resize(NDVI_copy, (0,0), fx=2, fy=2)\n",
    "\n",
    "cv2.imshow('image',NDVI_copy)\n",
    "cv2.waitKey(0)\n",
    "cv2.destroyAllWindows()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once the random pixels were seleccted I calculated the percentage of forest on each 10x10 pixels array from CubeSat. To do this I used the clustering K-mean tool in Python from sci-kit-learn and applied it to all the selected CubeSat pixels, in total 10,000, using the three visible light bands."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Visualiza the 3 visible light bands values os the CubeSat image \n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111, projection='3d')\n",
    "\n",
    "for pl in planet_training:\n",
    "    value_list = pl.reshape((-1,4))\n",
    "    for index in range(value_list.shape[0]):\n",
    "        ax.scatter(value_list[index,0], value_list[index,1], value_list[index,2])\n",
    "\n",
    "ax.set_xlabel('R')\n",
    "ax.set_ylabel('G')\n",
    "ax.set_zlabel('B')\n",
    "\n",
    "plt.show()\n",
    "\n",
    "# Transform the list into a big matrix to allow the clustering to work in all areas at the same time \n",
    "# so forest and non forest are classified using the same criteria in all pixels\n",
    "planet_training_v=planet_training.copy()\n",
    "planet_training_vs=planet_training_v[0]\n",
    "for selected_image_index in range(1,len(planet_training_v)):\n",
    "    planet_training_vs=np.concatenate((planet_training_vs,planet_training_v[selected_image_index]))\n",
    "    \n",
    "# Clustering into two classes\n",
    "X = planet_training_vs.reshape((-1,4))[:,:3]\n",
    "kmeans = KMeans(n_clusters=2, random_state=0).fit(X)   \n",
    "\n",
    "#Show how the pixels where classified \n",
    "fig = plt.figure()\n",
    "    ax = fig.add_subplot(111, projection='3d')\n",
    "    \n",
    "    for index in range(X.shape[0]):\n",
    "        \n",
    "        colr = 'b'\n",
    "        if kmeans.labels_[index] == 1:\n",
    "            colr = 'y'\n",
    "        \n",
    "        ax.scatter(X[index,0], X[index,1], X[index,2], c = colr)\n",
    "    \n",
    "    for index in range(kmeans.cluster_centers_.shape[0]):\n",
    "        \n",
    "        ax.scatter(kmeans.cluster_centers_[index,0], kmeans.cluster_centers_[index,1], kmeans.cluster_centers_[index,2], c = 'r')\n",
    "        \n",
    "    ax.set_xlabel('R')\n",
    "    ax.set_ylabel('G')\n",
    "    ax.set_zlabel('B')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then I calculated the number of pixels classified as forest on each 10x10 CubeSat pixels array, being this number the forest percentage cover."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "percent=[]\n",
    "for n in range(len(c_pixels)):    \n",
    "    percent.append(1-c_pixels[n].sum())    \n",
    "for index in range(x.shape[0]):   \n",
    "    print('percentaje: {}, NDVI: {}'.format(percent[index],NDVI_training[index]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From the total sample size, we selected half of it as training data set and the other half as a test dataset. Using a glm I fitted a model using the training data set in which the predictor was the NDVI and the response forest percentage."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot both variables to determine the glm to use\n",
    "plt.scatter(percent,NDVI_training)    \n",
    "plt.show()\n",
    "\n",
    "# Set parameters to divide the data set into test and traning\n",
    "mid=int(len(percent)/2)\n",
    "end=int(len(percent))\n",
    "\n",
    "#Fit the regression\n",
    "reg = linear_model.BayesianRidge()\n",
    "reg.fit(np.asarray(percent[:mid]).reshape(-1,1),np.asarray(NDVI_training[:mid]).reshape(-1,1))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, I used the model to predict the forest cover percentage in each 30 m pixel."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "NDVI_complete=np.asarray(NDVI.copy())\n",
    "NDVI_com_vec=NDVI_complete.reshape(-1)\n",
    "percent_all_pred=[]#np.asarray(reg.predict([[percent_test[0]]]))\n",
    "for i in range(len(NDVI_com_vec)):\n",
    "    percent_all_pred.append(reg.predict([[NDVI_com_vec[i]]]))\n",
    "    \n",
    "percent_all_pred1=np.asarray(percent_all_pred)\n",
    "percent_all_pred1=percent_all_pred1.reshape(447,451)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I validated this model predicting tree cover in the test data set and comparing the results with tree cover calculated using the clustering segmentation. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "NDVI_test=np.asarray(NDVI_training.copy()[mid:end])\n",
    "percent_y_pred=[]#np.asarray(reg.predict([[percent_test[0]]]))\n",
    "for i in range(len(NDVI_test)):\n",
    "    percent_y_pred.append(reg.predict([[NDVI_test[i]]]))\n",
    "np.asarray(percent_y_pred)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Results\n",
    "*Clustering*<br>\n",
    "The data used for the clustering showed a high correlation of reflectance from the three visible spectrum bands. The clustering algorithm clustered areas with a high reflectance of the three visible spectrums as non-forest. Most of the non-forest areas are bare ground soil or grass, which have a higher reflectance on the visible spectrum, whereas vegetation absorbs most of the blue and red light and reflect part of the green one. The results of the segmentation were visually inspected to evaluate the capacity of the clustering to properly classify the land cover."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Shape the kmeans result from the clustering in the pixels again\n",
    "kmean_pixels=kmeans.labels_.copy()\n",
    "c_pixels=list()\n",
    "for i in range(len(planet_training)):\n",
    "    c_pixels.append(kmean_pixels[(i*100):(i*100)+100].reshape((10,10)))\n",
    "    \n",
    "# Visualize in red the pixel that is being display in the NDVI image####   \n",
    "    NDVI_copy = cv2.cvtColor(NDVI.copy(), cv2.COLOR_GRAY2BGR)\n",
    "    \n",
    "    NDVI_copy[x[selected_image_index],y[selected_image_index]] = [0,0,255] \n",
    "    \n",
    "    NDVI_copy = cv2.resize(NDVI_copy, (0,0), fx=2, fy=2)\n",
    "    \n",
    "    cv2.imshow('image',NDVI_copy)\n",
    "    #cv2.waitKey(0)\n",
    "    #cv2.destroyAllWindows()\n",
    "    \n",
    "# Subset of planet image centered ont the segment area to compare with the matrix \n",
    "# representing the clustering we formated above and find out if forest is clustered\n",
    "# as 1 or 0 and asses the quality of the clustering\n",
    "\n",
    "for  selected_image_index in range(len(planet_training)):\n",
    "selecter_image_index=0\n",
    "    x_planet_i=(x[selected_image_index]*10)-50\n",
    "    x_planet_e=((x[selected_image_index]*10)+10)+50\n",
    "    y_planet_i=(y[selected_image_index]*10)-50\n",
    "    y_planet_e=((y[selected_image_index]*10)+10)+50\n",
    "    small_copy_crop = planet.copy()[x_planet_i:x_planet_e,y_planet_i:y_planet_e,:3]\n",
    "    small_copy_crop[50,:] = [122,0,0]\n",
    "    small_copy_crop[-50,:] = [122,0,0]\n",
    "    small_copy_crop[:,50] = [122,0,0]\n",
    "    small_copy_crop[:,-50] = [122,0,0]\n",
    "    small_copy_crop = cv2.resize(small_copy_crop, (0,0), fx=5, fy=5)\n",
    "    \n",
    "    print('image number: {}'.format(selected_image_index))\n",
    "    print(c_pixels[selected_image_index])\n",
    "    cv2.destroyAllWindows()\n",
    "\n",
    "    cv2.imshow('indice 0',small_copy_crop*5)\n",
    "    cv2.waitKey(0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*General linear model*<br>\n",
    "The glm model intercept was -1.16 with an increase of forest percentage of 2,29 for each NDVI unit increase (Fig 1.). When comparing the forest percentage predicted using the test data set on the model with the real data obtained with the clustering segmentation, I obtained a variance score of 0,45."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "reg.coef_\n",
    "reg.intercept_\n",
    "\n",
    "plt.scatter(NDVI_training,percent)\n",
    "X_plot = np.linspace(0.4,0.8,100)\n",
    "plt.plot(X_plot,reg.intercept_+(reg.coef_*X_plot))\n",
    "plt.show()\n",
    "\n",
    "# The coefficients\n",
    "print('Coefficient: \\n', reg.coef_)\n",
    "# The mean squared error\n",
    "print(\"Mean squared error: %.2f\"\n",
    "      % mean_squared_error(np.asarray(percent[mid:end]).reshape(-1,1), np.asarray(percent_y_pred)))\n",
    "# Explained variance score: 1 is perfect prediction\n",
    "print('Variance score: %.2f' % r2_score(np.asarray(percent[mid:end]).reshape(-1,1), np.asarray(percent_y_pred)))\n",
    "print('Explained variance score: %.2f' % explained_variance_score(np.asarray(percent[mid:end]).reshape(-1,1), np.asarray(percent_y_pred),multioutput='uniform_average')) \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Map of forest cover on 23,000 ha*<br>\n",
    "The obtained image from predicting forest cover in the whole area was visually examined to confirm that the classification of forest percentage corresponded with the expected results. The new image seemed to appropriately classify areas with no forest in a more categorical way and more informative than the NDVI image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cv2.imshow('image',percent_all_pred1)\n",
    "cv2.waitKey(0)\n",
    "cv2.destroyAllWindows()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Next steps<br>\n",
    "The chosen methods for clustering and creation of the model were relatively simple and provided a good ratio of prediction capacity. Still, using a beta regression instead of a linear model will probably improve the predictions.  The following step other than improve the model using the a beta regression is to use the obtained layer inside a model to predict forest recovery"
   ]
  }
 ],
 "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.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}