{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\GIDA2\\Anaconda3\\lib\\site-packages\\geopandas\\_compat.py:88: UserWarning: The Shapely GEOS version (3.4.3-CAPI-1.8.3 r4285) is incompatible with the GEOS version PyGEOS was compiled with (3.8.1-CAPI-1.13.3). Conversions between both will be slow.\n",
      "  shapely_geos_version, geos_capi_version_string\n"
     ]
    }
   ],
   "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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epsg:32611\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x28cb7d97f08>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 9016 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#open point shapefile\n",
    "pointData = gpd.read_file('Shp/pointData.shp')\n",
    "print(pointData.crs)\n",
    "pointData.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "EPSG:32611\n",
      "1\n"
     ]
    }
   ],
   "source": [
    "#open raster file\n",
    "ndviRaster = rasterio.open('Rst/ndviImage.tiff')\n",
    "print(ndviRaster.crs)\n",
    "print(ndviRaster.count)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x28cb90cf248>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 628656 characters]\n",
      "text/plain": [
       "<Figure size 864x864 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#show point and raster on a matplotlib plot\n",
    "fig, ax = plt.subplots(figsize=(12,12))\n",
    "pointData.plot(ax=ax, color='orangered')\n",
    "show(ndviRaster, ax=ax)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "239728.82443269365 3990037.7658296037\n",
      "249272.44019424845 3995708.972054402\n",
      "247098.89419619803 3987464.4872342106\n",
      "254893.6798443789 3991661.679506308\n",
      "255893.01133773543 3997207.9692944367\n",
      "257491.94172710588 3984091.7434441326\n",
      "243001.63507343628 3978745.319954675\n",
      "238054.94418132148 3978370.5706446664\n"
     ]
    }
   ],
   "source": [
    "#extract xy from point geometry\n",
    "for point in pointData['geometry']:\n",
    "    print(point.xy[0][0],point.xy[1][0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Point correspond to row, col: 708, 332\n",
      "Raster value on point 0.11 \n",
      "\n",
      "Point correspond to row, col: 519, 650\n",
      "Raster value on point 0.34 \n",
      "\n",
      "Point correspond to row, col: 794, 578\n",
      "Raster value on point 0.47 \n",
      "\n",
      "Point correspond to row, col: 654, 837\n",
      "Raster value on point 0.26 \n",
      "\n",
      "Point correspond to row, col: 469, 871\n",
      "Raster value on point 0.33 \n",
      "\n",
      "Point correspond to row, col: 906, 924\n",
      "Raster value on point 0.44 \n",
      "\n",
      "Point correspond to row, col: 1084, 441\n",
      "Raster value on point 0.26 \n",
      "\n",
      "Point correspond to row, col: 1097, 276\n",
      "Raster value on point 0.36 \n",
      "\n"
     ]
    }
   ],
   "source": [
    "#extract point value from raster\n",
    "for point in pointData['geometry']:\n",
    "    x = point.xy[0][0]\n",
    "    y = point.xy[1][0]\n",
    "    row, col = ndviRaster.index(x,y)\n",
    "    print(\"Point correspond to row, col: %d, %d\"%(row,col))\n",
    "    print(\"Raster value on point %.2f \\n\"%ndviRaster.read(1)[row,col])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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
}
