{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Import requires libraries\n",
    "This tutorial uses just Python core libraries, Scipy (Numpy, Matplotlib and others) and Rasterio"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "#from osgeo import ogr, gdal, osr\n",
    "import rasterio\n",
    "from rasterio.transform import Affine\n",
    "import numpy as np\n",
    "import os\n",
    "import matplotlib.pyplot as plt\n",
    "from rasterio.plot import show"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Setting up the input and output files\n",
    "We declare the path for the raster bands and the output NDVI and land cover change rasters. The land cover change contour shape path is also defined."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Input Raster and Vector Paths\n",
    "\n",
    "#Image-2019\n",
    "path_B5_2019=\"../Image20190203clip/LC08_L1TP_029047_20190203_20190206_01_T1_B5_clip.TIF\"\n",
    "path_B4_2019=\"../Image20190203clip/LC08_L1TP_029047_20190203_20190206_01_T1_B4_clip.TIF\"\n",
    "\n",
    "#Image-2014\n",
    "path_B5_2014=\"../Image20140205clip/LC08_L1TP_029047_20140205_20170307_01_T1_B5_clip.TIF\"\n",
    "path_B4_2014=\"../Image20140205clip/LC08_L1TP_029047_20140205_20170307_01_T1_B4_clip.TIF\"\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Output Files\n",
    "\n",
    "#Output NDVI Rasters \n",
    "path_NDVI_2019 = '../Output/NDVI2019.tif'\n",
    "path_NDVI_2014 = '../Output/NDVI2014.tif'\n",
    "path_NDVIChange_19_14 = '../Output/NDVIChange_19_14.tif'\n",
    "\n",
    "#NDVI Contours\n",
    "contours_NDVIChange_19_14 = '../Output/NDVIChange_19_14.shp'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Open the Landsat image bands with GDAL \n",
    "In this part we open the red (Band 4) and near infrared NIR (Band 5) with commands of the GDAL library and then we read the images as matrix arrays with float numbers of 32 bits."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[    0, 15510, 14779, ..., 18204, 18722, 18515],\n",
       "       [    0, 11765, 13832, ..., 15842, 16966, 17672],\n",
       "       [    0, 10027, 11470, ..., 14321, 14948, 15572],\n",
       "       ...,\n",
       "       [    0, 14584, 13352, ..., 15260, 14624, 13359],\n",
       "       [    0, 14569, 13616, ..., 15065, 13242, 11065],\n",
       "       [    0,     0,     0, ...,     0,     0,     0]], dtype=uint16)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Check bands\n",
    "B5_2019 = rasterio.open(path_B5_2019)\n",
    "B5_2019.read(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Open raster bands\n",
    "B5_2019 = rasterio.open(path_B5_2019)\n",
    "B4_2019 = rasterio.open(path_B4_2019)\n",
    "B5_2014 = rasterio.open(path_B5_2014)\n",
    "B4_2014 = rasterio.open(path_B4_2014)\n",
    "\n",
    "#Read bands as matrix arrays\n",
    "B52019_Data = B5_2019.read(1)\n",
    "B42019_Data = B4_2019.read(1)\n",
    "B52014_Data = B5_2014.read(1)\n",
    "B42014_Data = B4_2014.read(1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compare matrix sizes and geotransformation parameters\n",
    "The operation in between Landsat 8 bands depends that the resolution, size, src, and geotransformation parameters are the same for both images. In case these caracteristics do not coincide a warp, reproyection, scale or any other geospatial process would be necessary."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+init=epsg:32613\n",
      "+init=epsg:32613\n",
      "SRC OK\n"
     ]
    }
   ],
   "source": [
    "print(B5_2014.crs.to_proj4())\n",
    "print(B5_2019.crs.to_proj4())\n",
    "if B5_2014.crs.to_proj4()==B5_2019.crs.to_proj4(): print('SRC OK')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(610, 597)\n",
      "(610, 597)\n",
      "Array Size OK\n"
     ]
    }
   ],
   "source": [
    "print(B5_2014.shape)\n",
    "print(B5_2019.shape)\n",
    "if B5_2014.shape==B5_2019.shape: print('Array Size OK')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[652500.0, 29.98324958123953, 0.0, 2166000.0, 0.0, -30.0]\n",
      "[652500.0, 29.98324958123953, 0.0, 2166000.0, 0.0, -30.0]\n",
      "Geotransformation OK\n"
     ]
    }
   ],
   "source": [
    "print(B5_2014.get_transform())\n",
    "print(B5_2019.get_transform())\n",
    "if B5_2014.get_transform()==B5_2019.get_transform(): print('Geotransformation OK')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Get geotransformation parameters\n",
    "Since we have compared that rasters bands have same array size and are in the same position we can get the some spatial parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "597 610\n"
     ]
    }
   ],
   "source": [
    "#geotransform = B5_2014.get_transform()\n",
    "#originX,pixelWidth,empty,finalY,empty2,pixelHeight=geotransform\n",
    "cols =  B5_2014.width\n",
    "rows =  B5_2014.height\n",
    "print(cols,rows)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Affine(29.98324958123953, 0.0, 652500.0,\n",
       "       0.0, -30.0, 2166000.0)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#transform = Affine.from_gdal(originX,pixelWidth,empty,finalY,empty2,pixelHeight)\n",
    "transform = B5_2014.transform\n",
    "transform"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "CRS.from_epsg(32613)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "projection = B5_2014.crs\n",
    "projection"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compute the NDVI and store them as a different file\n",
    "We can compute the NDVI as a matrix algebra with some Numpy functions. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "ndvi2014 = np.divide(B52014_Data - B42014_Data, B52014_Data + B42014_Data,where=(B52014_Data - B42014_Data)!=0)\n",
    "ndvi2014[ndvi2014 > 1] = 1 #for water bodies\n",
    "\n",
    "ndvi2019 = np.divide(B52019_Data - B42019_Data, B52019_Data+ B42019_Data,where=(B52019_Data - B42019_Data)!=0)\n",
    "ndvi2019[ndvi2019 > 1] = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x19d3a36cfd0>"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 298548 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#plot the array\n",
    "plt.imshow(ndvi2014*100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def saveRaster(dataset,datasetPath,cols,rows,projection):\n",
    "    rasterSet = rasterio.open(datasetPath,\n",
    "                             'w',\n",
    "                             driver='GTiff',\n",
    "                             height=rows,\n",
    "                             width=cols,\n",
    "                             count=1,\n",
    "                             dtype=np.float32,\n",
    "                             crs=projection,\n",
    "                             transform=transform,)\n",
    "    rasterSet.write(dataset,1)\n",
    "    rasterSet.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "saveRaster(ndvi2014.astype('float32'),path_NDVI_2014,cols,rows,projection)\n",
    "\n",
    "saveRaster(ndvi2019.astype('float32'),path_NDVI_2019,cols,rows,projection)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot NDVI Images\n",
    "We create a function to plot the resulting NDVI images with a colobar"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "def plotNDVI(ndviImage,cmap):\n",
    "    src = rasterio.open(ndviImage,'r')\n",
    "    fig, ax = plt.subplots(1, figsize=(12, 12))\n",
    "    show(src, cmap=cmap, ax=ax)\n",
    "    ax.set_xlabel('Este')\n",
    "    ax.set_ylabel('Norte')\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plotNDVIContour(ndviImage,cmap):\n",
    "    src = rasterio.open(ndviImage,'r')\n",
    "    fig, ax = plt.subplots(1, figsize=(12, 12))\n",
    "    show(src, cmap=cmap, contour=True, ax=ax)\n",
    "    ax.set_xlabel('Este')\n",
    "    ax.set_ylabel('Norte')\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 1523416 characters]\n",
      "text/plain": [
       "<Figure size 1200x1200 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plotNDVI(path_NDVI_2014,'YlGn')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 1484132 characters]\n",
      "text/plain": [
       "<Figure size 1200x1200 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plotNDVIContour(path_NDVI_2014,'YlGn')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 1569600 characters]\n",
      "text/plain": [
       "<Figure size 1200x1200 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plotNDVI(path_NDVI_2019,'YlGn')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create a land cover change image\n",
    "\n",
    "We create a land cover change image based on the differences in NDVI from 2019 and 2014 imagery. The image will be stored as a separate file. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.        ,  0.03547079,  0.03795126, ...,  0.03159038,\n",
       "         0.03800274,  0.02669263],\n",
       "       [ 0.        ,  0.03765402,  0.05892348, ..., -0.00586917,\n",
       "         0.01892243,  0.0195063 ],\n",
       "       [ 0.        ,  0.03218401,  0.03739506, ..., -0.06377307,\n",
       "        -0.0306757 ,  0.03894268],\n",
       "       ...,\n",
       "       [ 0.        , -0.01659706,  0.01140219, ...,  0.00272184,\n",
       "         0.02552676,  0.06663901],\n",
       "       [ 0.        ,  0.00599447, -0.00407705, ...,  0.05750174,\n",
       "         0.00467577,  0.08448431],\n",
       "       [ 0.        ,  0.        ,  0.        , ...,  0.        ,\n",
       "         0.        ,  0.        ]])"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ndviChange = ndvi2019-ndvi2014\n",
    "ndviChange = np.where((ndvi2014>-999) & (ndvi2019>-999),ndviChange,-999)\n",
    "ndviChange"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "saveRaster(ndviChange.astype('float32'),path_NDVIChange_19_14,cols,rows,projection)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 1217816 characters]\n",
      "text/plain": [
       "<Figure size 1200x1200 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plotNDVI(path_NDVIChange_19_14,'Spectral')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
