{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#import required libraries\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import geopandas as gpd\n",
    "import rasterio\n",
    "from rasterio.plot import show\n",
    "from skimage.feature import match_template\n",
    "import numpy as np\n",
    "from PIL import Image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#open point shapefile\n",
    "pointData = gpd.read_file('Shp/pointData.shp')\n",
    "print('CRS of Point Data: ' + str(pointData.crs))\n",
    "\n",
    "#open raster file\n",
    "palmRaster = rasterio.open('Rst/palmaOrthoTotal_14cm.tif')\n",
    "print('CRS of Raster Data: ' + str(palmRaster.crs))\n",
    "print('Number of Raster Bands: ' + str(palmRaster.count))\n",
    "print('Interpretation of Raster Bands: ' + str(palmRaster.colorinterp))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#show point and raster on a matplotlib plot\n",
    "fig, ax = plt.subplots(figsize=(18,18))\n",
    "pointData.plot(ax=ax, color='orangered', markersize=100)\n",
    "show(palmRaster, ax=ax)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#selected band: green\n",
    "greenBand = palmRaster.read(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#extract point value from raster\n",
    "surveyRowCol = []\n",
    "for index, values in pointData.iterrows():\n",
    "    x = values['geometry'].xy[0][0]\n",
    "    y = values['geometry'].xy[1][0]\n",
    "    row, col = palmRaster.index(x,y)\n",
    "    print(\"Point N°:%d corresponds to row, col: %d, %d\"%(index,row,col))\n",
    "    surveyRowCol.append([row,col])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# number of template images\n",
    "print('Number of template images: %d'%len(surveyRowCol))\n",
    "# define ratio of analysis\n",
    "radio = 25"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#show all the points of interest, please be careful to have a complete image, otherwise the model wont run\n",
    "fig, ax = plt.subplots(1, len(surveyRowCol),figsize=(20,5))\n",
    "\n",
    "for index, item in enumerate(surveyRowCol):\n",
    "    row = item[0]\n",
    "    col = item[1]\n",
    "    ax[index].imshow(greenBand)\n",
    "    ax[index].plot(col,row,color='red', linestyle='dashed', marker='+',\n",
    "     markerfacecolor='blue', markersize=8)\n",
    "    ax[index].set_xlim(col-radio,col+radio)\n",
    "    ax[index].set_ylim(row-radio,row+radio)\n",
    "    ax[index].axis('off')\n",
    "    ax[index].set_title(index)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Match the image to the template\n",
    "\n",
    "listaresultados = []\n",
    "\n",
    "templateBandList = []\n",
    "for rowCol in surveyRowCol:\n",
    "    imageList = []\n",
    "    row = rowCol[0]\n",
    "    col = rowCol[1]\n",
    "    #append original band\n",
    "    imageList.append(greenBand[row-radio:row+radio, col-radio:col+radio])\n",
    "    #append rotated images\n",
    "    templateBandToRotate = greenBand[row-2*radio:row+2*radio, col-2*radio:col+2*radio]\n",
    "    rotationList = [i*30 for i in range(1,4)]\n",
    "    for rotation in rotationList:\n",
    "        rotatedRaw = Image.fromarray(templateBandToRotate)\n",
    "        rotatedImage = rotatedRaw.rotate(rotation)\n",
    "        imageList.append(np.asarray(rotatedImage)[radio:-radio,radio:-radio])\n",
    "    #plot original and rotated images\n",
    "    fig, ax = plt.subplots(1, len(imageList),figsize=(12,12))\n",
    "    for index, item in enumerate(imageList):\n",
    "        ax[index].imshow(imageList[index])\n",
    "    #add images to total list\n",
    "    templateBandList+=imageList"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# match the template image to the orthophoto\n",
    "matchXYList = []\n",
    "for index, templateband in enumerate(templateBandList):\n",
    "    if index%10 == 0:\n",
    "        print('Match template ongoing for figure Nº %d'%index)\n",
    "    matchTemplate = match_template(greenBand, templateband, pad_input=True)\n",
    "    matchTemplateFiltered = np.where(matchTemplate>np.quantile(matchTemplate,0.9996))\n",
    "    for item in zip(matchTemplateFiltered[0],matchTemplateFiltered[1]):\n",
    "        x, y = palmRaster.xy(item[0], item[1])\n",
    "        matchXYList.append([x,y])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot interpreted points over the image\n",
    "fig, ax = plt.subplots(figsize=(20, 20))\n",
    "matchXYArray = np.array(matchXYList)\n",
    "ax.scatter(matchXYArray[:,0],matchXYArray[:,1], marker='o',c='orangered', s=100, alpha=0.25)\n",
    "show(palmRaster, ax=ax)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# cluster analysis\n",
    "from sklearn.cluster import Birch\n",
    "brc = Birch(branching_factor=10000, n_clusters=None, threshold=2e-5, compute_labels=True)\n",
    "brc.fit(matchXYArray)\n",
    "birchPoint = brc.subcluster_centers_\n",
    "birchPoint"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot clustered points\n",
    "fig = plt.figure(figsize=(20, 20))\n",
    "ax = fig.add_subplot(111)\n",
    "ax.scatter(birchPoint[:,[0]],birchPoint[:,[1]], marker='o',color='orangered',s=100)\n",
    "show(palmRaster, ax=ax)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# save xy to a csv file\n",
    "np.savetxt(\"Txt/birchPoints.csv\", birchPoint, delimiter=\",\")"
   ]
  }
 ],
 "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.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
