{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Assignments Fall 2022\n",
    "\n",
    "\n",
    "This assignment is compulsory and need to be delivered before 17th of November, 2022 12pm, Nigeria time.\n",
    "\n",
    "Please use this jupyter file by adding cell-codes and cell-text with your solution. Save the file as *assignment_fall2022_name_surname.ipynb* and send as e-mail attachment to *g.amatulli@spatial-ecology.net* and *Erin.Stearns@gatesfoundation.org*.\n",
    "\n",
    "The exercises are based on lectures and materials that we covered during the lectures.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Enter in your */media/sf_LVM_shared/my_SE_data/exercise/*, save there your *assignment_fall2022_name_surname.ipynb* and start to solve the exercises by adding new cell-codes and cell-text."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1st excercise\n",
    "20% grading value\n",
    "\n",
    "Uisng the *geodata/dem/SA_elevation_mn_GMTED2010_mn_msk.tif* ocean value mask out the \n",
    "*geodata/vegetation/SA_tree_mn_percentage_GFC2013.tif* in such a way that the output *geodata/vegetation/SA_tree_mn_percentage_GFC2013_msk.tif* has as ocean value -1. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Driver: GTiff/GeoTIFF\n",
      "Files: geodata/dem/SA_elevation_mn_GMTED2010_mn_msk.tif\n",
      "Size is 5880, 8400\n",
      "Coordinate System is:\n",
      "GEOGCRS[\"WGS 84\",\n",
      "    DATUM[\"World Geodetic System 1984\",\n",
      "        ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n",
      "            LENGTHUNIT[\"metre\",1]]],\n",
      "    PRIMEM[\"Greenwich\",0,\n",
      "        ANGLEUNIT[\"degree\",0.0174532925199433]],\n",
      "    CS[ellipsoidal,2],\n",
      "        AXIS[\"geodetic latitude (Lat)\",north,\n",
      "            ORDER[1],\n",
      "            ANGLEUNIT[\"degree\",0.0174532925199433]],\n",
      "        AXIS[\"geodetic longitude (Lon)\",east,\n",
      "            ORDER[2],\n",
      "            ANGLEUNIT[\"degree\",0.0174532925199433]],\n",
      "    ID[\"EPSG\",4326]]\n",
      "Data axis to CRS axis mapping: 2,1\n",
      "Origin = (-83.000000000000000,14.000000000000000)\n",
      "Pixel Size = (0.008333333333333,-0.008333333333333)\n",
      "Metadata:\n",
      "  AREA_OR_POINT=Area\n",
      "  TIFFTAG_DATETIME=2021:05:12 09:34:03\n",
      "  TIFFTAG_DOCUMENTNAME=geodata/dem/SA_elevation_mn_GMTED2010_mn_msk.tif\n",
      "  TIFFTAG_SOFTWARE=pktools 2.6.7.6 by Pieter Kempeneers\n",
      "Image Structure Metadata:\n",
      "  COMPRESSION=DEFLATE\n",
      "  INTERLEAVE=BAND\n",
      "Corner Coordinates:\n",
      "Upper Left  ( -83.0000000,  14.0000000) ( 83d 0' 0.00\"W, 14d 0' 0.00\"N)\n",
      "Lower Left  ( -83.0000000, -56.0000000) ( 83d 0' 0.00\"W, 56d 0' 0.00\"S)\n",
      "Upper Right ( -34.0000000,  14.0000000) ( 34d 0' 0.00\"W, 14d 0' 0.00\"N)\n",
      "Lower Right ( -34.0000000, -56.0000000) ( 34d 0' 0.00\"W, 56d 0' 0.00\"S)\n",
      "Center      ( -58.5000000, -21.0000000) ( 58d30' 0.00\"W, 21d 0' 0.00\"S)\n",
      "Band 1 Block=5880x1 Type=Int16, ColorInterp=Gray\n",
      "    Computed Min/Max=-400.000,6736.000\n",
      "  NoData Value=-9999\n"
     ]
    }
   ],
   "source": [
    "# check the value range of the mask image\n",
    "! gdalinfo -mm geodata/dem/SA_elevation_mn_GMTED2010_mn_msk.tif"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Default software rendering mode (use -h if accelerated video card installed).\n",
      "Loading tools from /usr/bin/openev/tools/Tool_Export.py\n",
      "Loading tools from /usr/bin/openev/tools/Tool_ShapesGrid.py\n",
      "Loading tools from /usr/bin/openev/tools/Tool_DriverList.py\n"
     ]
    }
   ],
   "source": [
    "# open the file to be sure that the ocean is labeled as -9999 no data \n",
    "! /usr/bin/openev/bin/openev geodata/dem/SA_elevation_mn_GMTED2010_mn_msk.tif "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Driver: GTiff/GeoTIFF\n",
      "Files: geodata/vegetation/SA_tree_mn_percentage_GFC2013.tif\n",
      "Size is 5880, 8400\n",
      "Coordinate System is:\n",
      "GEOGCRS[\"WGS 84\",\n",
      "    DATUM[\"World Geodetic System 1984\",\n",
      "        ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n",
      "            LENGTHUNIT[\"metre\",1]]],\n",
      "    PRIMEM[\"Greenwich\",0,\n",
      "        ANGLEUNIT[\"degree\",0.0174532925199433]],\n",
      "    CS[ellipsoidal,2],\n",
      "        AXIS[\"geodetic latitude (Lat)\",north,\n",
      "            ORDER[1],\n",
      "            ANGLEUNIT[\"degree\",0.0174532925199433]],\n",
      "        AXIS[\"geodetic longitude (Lon)\",east,\n",
      "            ORDER[2],\n",
      "            ANGLEUNIT[\"degree\",0.0174532925199433]],\n",
      "    ID[\"EPSG\",4326]]\n",
      "Data axis to CRS axis mapping: 2,1\n",
      "Origin = (-83.000000000000000,14.000000000000000)\n",
      "Pixel Size = (0.008333333333333,-0.008333333333333)\n",
      "Metadata:\n",
      "  AREA_OR_POINT=Area\n",
      "  Author=giuseppe.amatulli@gmail.com using pktools\n",
      "  Input dataset=Global Forest Change 2000-2012 (Hansen 2013)\n",
      "  Input layer=Tree canopy cover for year 2000 (%)\n",
      "  Offset=0\n",
      "  Output=Mean of Tree cover (%)\n",
      "  Scale=0.01\n",
      "Image Structure Metadata:\n",
      "  COMPRESSION=DEFLATE\n",
      "  INTERLEAVE=BAND\n",
      "Corner Coordinates:\n",
      "Upper Left  ( -83.0000000,  14.0000000) ( 83d 0' 0.00\"W, 14d 0' 0.00\"N)\n",
      "Lower Left  ( -83.0000000, -56.0000000) ( 83d 0' 0.00\"W, 56d 0' 0.00\"S)\n",
      "Upper Right ( -34.0000000,  14.0000000) ( 34d 0' 0.00\"W, 14d 0' 0.00\"N)\n",
      "Lower Right ( -34.0000000, -56.0000000) ( 34d 0' 0.00\"W, 56d 0' 0.00\"S)\n",
      "Center      ( -58.5000000, -21.0000000) ( 58d30' 0.00\"W, 21d 0' 0.00\"S)\n",
      "Band 1 Block=5880x1 Type=UInt16, ColorInterp=Gray\n",
      "    Computed Min/Max=0.000,10000.000\n"
     ]
    }
   ],
   "source": [
    "# check the value range of the input image \n",
    "! gdalinfo -mm geodata/vegetation/SA_tree_mn_percentage_GFC2013.tif"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Check the data type that support a range from -1 to 10000**  \n",
    "see http://spatial-ecology.net/docs/build/html/GDAL/gdal_osgeo.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0...10...20...30...40...50...60...70...80...90...100 - done.\n"
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "pksetmask -co COMPRESS=DEFLATE -co ZLEVEL=9 -ot Int16 \\\n",
    "-m  geodata/dem/SA_elevation_mn_GMTED2010_mn_msk.tif -msknodata -9999 -nodata -1 \\\n",
    "-i  geodata/vegetation/SA_tree_mn_percentage_GFC2013.tif -o geodata/vegetation/SA_tree_mn_percentage_GFC2013_msk.tif"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Default software rendering mode (use -h if accelerated video card installed).\n",
      "Loading tools from /usr/bin/openev/tools/Tool_Export.py\n",
      "Loading tools from /usr/bin/openev/tools/Tool_ShapesGrid.py\n",
      "Loading tools from /usr/bin/openev/tools/Tool_DriverList.py\n"
     ]
    }
   ],
   "source": [
    "# open the file to be sure that the ocean is labeled as -1 no data \n",
    "! /usr/bin/openev/bin/openev geodata/vegetation/SA_tree_mn_percentage_GFC2013_msk.tif"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2nd exercise\n",
    "20% grading value\n",
    "\n",
    "Using the .tif files in *geodata/landsat_ct/*, search for a gdal command that is able to create a vector-tile file having a tile-polygon for each tile tif. Run the commands and visualize the vector-tile file in qgis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "071W_41N_med_B1.tif  072W_41N_med_B1.tif  073W_41N_med_B1.tif\r\n",
      "071W_42N_med_B1.tif  072W_42N_med_B1.tif  073W_42N_med_B1.tif\r\n",
      "072W_40N_med_B1.tif  073W_40N_med_B1.tif\r\n"
     ]
    }
   ],
   "source": [
    "! ls geodata/landsat_ct/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Default software rendering mode (use -h if accelerated video card installed).\n",
      "Loading tools from /usr/bin/openev/tools/Tool_Export.py\n",
      "Loading tools from /usr/bin/openev/tools/Tool_ShapesGrid.py\n",
      "Loading tools from /usr/bin/openev/tools/Tool_DriverList.py\n"
     ]
    }
   ],
   "source": [
    "# open all the tiles  \n",
    "! /usr/bin/openev/bin/openev geodata/landsat_ct/*.tif"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**First solution**\n",
    "This is the most simple and elegant solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Creating new index file...\n",
      "total 364\n",
      "-rwxrwx--- 1 root vboxsf 26579 Oct 12 17:57 071W_41N_med_B1.tif\n",
      "-rwxrwx--- 1 root vboxsf   361 Nov 24 10:21 071W_41N_med_B1.tif.aux.xml\n",
      "-rwxrwx--- 1 root vboxsf 26877 Oct 12 17:57 071W_42N_med_B1.tif\n",
      "-rwxrwx--- 1 root vboxsf   363 Nov 24 10:21 071W_42N_med_B1.tif.aux.xml\n",
      "-rwxrwx--- 1 root vboxsf 27168 Oct 12 17:57 072W_40N_med_B1.tif\n",
      "-rwxrwx--- 1 root vboxsf   361 Nov 24 10:21 072W_40N_med_B1.tif.aux.xml\n",
      "-rwxrwx--- 1 root vboxsf 25833 Oct 12 17:57 072W_41N_med_B1.tif\n",
      "-rwxrwx--- 1 root vboxsf   361 Nov 24 10:21 072W_41N_med_B1.tif.aux.xml\n",
      "-rwxrwx--- 1 root vboxsf 24464 Oct 12 17:57 072W_42N_med_B1.tif\n",
      "-rwxrwx--- 1 root vboxsf   361 Nov 24 10:21 072W_42N_med_B1.tif.aux.xml\n",
      "-rwxrwx--- 1 root vboxsf 27736 Oct 12 17:57 073W_40N_med_B1.tif\n",
      "-rwxrwx--- 1 root vboxsf   361 Nov 24 10:21 073W_40N_med_B1.tif.aux.xml\n",
      "-rwxrwx--- 1 root vboxsf 26047 Oct 12 17:57 073W_41N_med_B1.tif\n",
      "-rwxrwx--- 1 root vboxsf   362 Nov 24 10:21 073W_41N_med_B1.tif.aux.xml\n",
      "-rwxrwx--- 1 root vboxsf 25851 Oct 12 17:57 073W_42N_med_B1.tif\n",
      "-rwxrwx--- 1 root vboxsf   351 Nov 24 10:21 073W_42N_med_B1.tif.aux.xml\n",
      "-rwxrwx--- 1 root vboxsf  2106 Nov 24  2022 alltif_tile.dbf\n",
      "-rwxrwx--- 1 root vboxsf 98304 Nov 24 10:21 alltif_tile.gpkg\n",
      "-rwxrwx--- 1 root vboxsf   145 Nov 24  2022 alltif_tile.prj\n",
      "-rwxrwx--- 1 root vboxsf  1188 Nov 24  2022 alltif_tile.shp\n",
      "-rwxrwx--- 1 root vboxsf   164 Nov 24  2022 alltif_tile.shx\n"
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "rm -f geodata/landsat_ct/alltif_tile.*                                  \n",
    "gdaltindex geodata/landsat_ct/alltif_tile.shp  geodata/landsat_ct/????_???_med_B1.tif   ### perfect 100%%\n",
    "ls -l geodata/landsat_ct/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\t- 'VirtualXPath'\t[XML Path Language - XPath]\n",
      "\t- 'VirtualXPath'\t[XML Path Language - XPath]\n",
      "\t- 'VirtualXPath'\t[XML Path Language - XPath]\n",
      "\t- 'VirtualXPath'\t[XML Path Language - XPath]\n",
      "\t- 'VirtualXPath'\t[XML Path Language - XPath]\n",
      "\t- 'VirtualXPath'\t[XML Path Language - XPath]\n",
      "\t- 'VirtualXPath'\t[XML Path Language - XPath]\n",
      "\t- 'VirtualXPath'\t[XML Path Language - XPath]\n",
      "\t- 'VirtualXPath'\t[XML Path Language - XPath]\n"
     ]
    }
   ],
   "source": [
    "# open all the tif tiles and the vector-tile\n",
    "! qgis geodata/landsat_ct/*.tif geodata/landsat_ct/alltif_tile.gpkg  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Second soulution**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This solution is a bit complex. Anyway can be usefull to see how to combine different gdal/pktools comands to be able to get the same results. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0...10...20...30...40...50...60...70...80...90...100 - done.\n",
      "    Computed Min/Max=1.000,1.000\n",
      "0...10...20...30...40...50...60...70...80...90...100 - done.\n",
      "    Computed Min/Max=2.000,2.000\n",
      "0...10...20...30...40...50...60...70...80...90...100 - done.\n",
      "    Computed Min/Max=3.000,3.000\n",
      "0...10...20...30...40...50...60...70...80...90...100 - done.\n",
      "    Computed Min/Max=4.000,4.000\n",
      "0...10...20...30...40...50...60...70...80...90...100 - done.\n",
      "    Computed Min/Max=5.000,5.000\n",
      "0...10...20...30...40...50...60...70...80...90...100 - done.\n",
      "    Computed Min/Max=6.000,6.000\n",
      "0...10...20...30...40...50...60...70...80...90...100 - done.\n",
      "    Computed Min/Max=7.000,7.000\n",
      "0...10...20...30...40...50...60...70...80...90...100 - done.\n",
      "    Computed Min/Max=8.000,8.000\n",
      "0...10...20...30...40...50...60...70...80...90...100 - done.\n",
      "0...10...20...30...40...50...60...70...80...90...Creating output geodata/landsat_ct/alltif_tile_poligonize.gpkg of format GPKG.\n",
      "100 - done.\n"
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "# for each tile change all the pixel value to 1 or 2 or 3 \n",
    "n=1\n",
    "for file in geodata/landsat_ct/????_???_med_B1.tif; do \n",
    "filename=$(basename $file .tif)\n",
    "pkgetmask -min -1 -max 999999999 -data $n -co COMPRESS=DEFLATE -co ZLEVEL=9 -ot Byte -i $file -o geodata/landsat_ct/${filename}_byte.tif\n",
    "gdal_edit.py -a_nodata 0 geodata/landsat_ct/${filename}_byte.tif\n",
    "n=$(expr $n + 1)\n",
    "gdalinfo -mm geodata/landsat_ct/${filename}_byte.tif | grep Comp\n",
    "done\n",
    "\n",
    "# build a vrt in order to combine the created tif \n",
    "gdalbuildvrt -overwrite   geodata/landsat_ct/all_byte.vrt   geodata/landsat_ct/*_byte.tif\n",
    "\n",
    "# run polygonize \n",
    "rm -f geodata/landsat_ct/alltif_tile_poligonize.gpkg\n",
    "gdal_polygonize.py  geodata/landsat_ct/all_byte.vrt geodata/landsat_ct/alltif_tile_poligonize.gpkg\n",
    "rm -f geodata/landsat_ct/all_byte.vrt geodata/landsat_ct/*_byte.tif"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\t- 'VirtualXPath'\t[XML Path Language - XPath]\n",
      "\t- 'VirtualXPath'\t[XML Path Language - XPath]\n",
      "\t- 'VirtualXPath'\t[XML Path Language - XPath]\n",
      "\t- 'VirtualXPath'\t[XML Path Language - XPath]\n",
      "\t- 'VirtualXPath'\t[XML Path Language - XPath]\n",
      "\t- 'VirtualXPath'\t[XML Path Language - XPath]\n",
      "\t- 'VirtualXPath'\t[XML Path Language - XPath]\n",
      "\t- 'VirtualXPath'\t[XML Path Language - XPath]\n",
      "\t- 'VirtualXPath'\t[XML Path Language - XPath]\n",
      "\t- 'VirtualXPath'\t[XML Path Language - XPath]\n"
     ]
    }
   ],
   "source": [
    "! qgis geodata/landsat_ct/????_???_med_B1.tif geodata/landsat_ct/alltif_tile_poligonize.gpkg "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This solution is also fine but you do not have the tile name in the vector attribute page."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3rd excercise \n",
    "30% grading value\n",
    "\n",
    "\n",
    " * Open one file geodata/LST/LST_MOYDmax_month1.tif in qgis and create a vector file with 3 polygons, save it in LST folder. Any shape of the poligons are fine just be sure that is \"on top\" of the raster file. \n",
    " * Run a pktools command that is able to extract zonal statistic for one .tif file and save the results to a \" GPKG -raster,vector- (rw+vs): GeoPackage\" file format. \n",
    " * Build a for loop that is able to extract zonal statistic for all .tif files and save the GPKG. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\t- 'VirtualXPath'\t[XML Path Language - XPath]\n",
      "\t- 'VirtualXPath'\t[XML Path Language - XPath]\n",
      "\t- 'VirtualXPath'\t[XML Path Language - XPath]\n",
      "\t- 'VirtualXPath'\t[XML Path Language - XPath]\n",
      "\t- 'VirtualXPath'\t[XML Path Language - XPath]\n",
      "\t- 'VirtualXPath'\t[XML Path Language - XPath]\n",
      "\t- 'VirtualXPath'\t[XML Path Language - XPath]\n",
      "\t- 'VirtualXPath'\t[XML Path Language - XPath]\n",
      "\t- 'VirtualXPath'\t[XML Path Language - XPath]\n",
      "\t- 'VirtualXPath'\t[XML Path Language - XPath]\n"
     ]
    }
   ],
   "source": [
    "# open the tif file and drow the polingons save it as LST/polygons_zonal.gpkg\n",
    "! qgis geodata/LST/LST_MOYDmax_month1.tif geodata/LST/polygons_zonal.gpkg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\t- 'VirtualXPath'\t[XML Path Language - XPath]\n",
      "\t- 'VirtualXPath'\t[XML Path Language - XPath]\n",
      "processing layer polygons_zonal\n",
      "0...10...20...30...40...50...60...70...80...90...100 - done.\n"
     ]
    }
   ],
   "source": [
    "! rm -f geodata/LST/polygons_zonal_month1.gpkg\n",
    "! pkextractogr  -f \"GPKG\" -srcnodata -9999 -r mean \\\n",
    "-i geodata/LST/LST_MOYDmax_month1.tif -s  geodata/LST/polygons_zonal.gpkg -o geodata/LST/polygons_zonal_month1.gpkg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\t- 'VirtualXPath'\t[XML Path Language - XPath]\n",
      "INFO: Open of `geodata/LST/polygons_zonal_month1.gpkg'\n",
      "      using driver `GPKG' successful.\n",
      "\n",
      "Layer name: polygons_zonal\n",
      "Geometry: Polygon\n",
      "Feature Count: 3\n",
      "Extent: (28.159091, 0.336869) - (33.182828, 3.912879)\n",
      "Layer SRS WKT:\n",
      "GEOGCRS[\"WGS 84\",\n",
      "    DATUM[\"World Geodetic System 1984\",\n",
      "        ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n",
      "            LENGTHUNIT[\"metre\",1]]],\n",
      "    PRIMEM[\"Greenwich\",0,\n",
      "        ANGLEUNIT[\"degree\",0.0174532925199433]],\n",
      "    CS[ellipsoidal,2],\n",
      "        AXIS[\"geodetic latitude (Lat)\",north,\n",
      "            ORDER[1],\n",
      "            ANGLEUNIT[\"degree\",0.0174532925199433]],\n",
      "        AXIS[\"geodetic longitude (Lon)\",east,\n",
      "            ORDER[2],\n",
      "            ANGLEUNIT[\"degree\",0.0174532925199433]],\n",
      "    USAGE[\n",
      "        SCOPE[\"Horizontal component of 3D system.\"],\n",
      "        AREA[\"World.\"],\n",
      "        BBOX[-90,-180,90,180]],\n",
      "    ID[\"EPSG\",4326]]\n",
      "Data axis to CRS axis mapping: 2,1\n",
      "FID Column = fid\n",
      "Geometry Column = geom\n",
      "ID: Integer (0.0)\n",
      "mean: Real (0.0)\n",
      "OGRFeature(polygons_zonal):1\n",
      "  ID (Integer) = 2\n",
      "  mean (Real) = 32.2190872990031\n",
      "\n",
      "OGRFeature(polygons_zonal):2\n",
      "  ID (Integer) = 1\n",
      "  mean (Real) = 37.5875141960014\n",
      "\n",
      "OGRFeature(polygons_zonal):3\n",
      "  ID (Integer) = 3\n",
      "  mean (Real) = 27.8825487457743\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# check the results\n",
    "! ogrinfo -al -geom=NO geodata/LST/polygons_zonal_month1.gpkg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0...10...20...30...40...50...60...70...80...90...100 - done.\n",
      "\t- 'VirtualXPath'\t[XML Path Language - XPath]\n",
      "\t- 'VirtualXPath'\t[XML Path Language - XPath]\n",
      "processing layer polygons_zonal\n",
      "0...10...20...30...40...50...60...70...80...90...100 - done.\n"
     ]
    }
   ],
   "source": [
    "# create a vrt and use it for the extraction \n",
    "!  -separate geodata/LST/LST_MOYDmax_allmonth.vrt geodata/LST/LST_MOYDmax_month?.tif geodata/LST/LST_MOYDmax_month??.tif\n",
    "! rm -f geodata/LST/polygons_zonal_allmonth.gpkg\n",
    "! pkextractogr  -f \"GPKG\" -srcnodata -9999 -r mean \\\n",
    "-i geodata/LST/LST_MOYDmax_allmonth.vrt -s  geodata/LST/polygons_zonal.gpkg  -o geodata/LST/polygons_zonal_allmonth.gpkg\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\t- 'VirtualXPath'\t[XML Path Language - XPath]\n",
      "INFO: Open of `geodata/LST/polygons_zonal_allmonth.gpkg'\n",
      "      using driver `GPKG' successful.\n",
      "\n",
      "Layer name: polygons_zonal\n",
      "Geometry: Polygon\n",
      "Feature Count: 3\n",
      "Extent: (28.159091, 0.336869) - (33.182828, 3.912879)\n",
      "Layer SRS WKT:\n",
      "GEOGCRS[\"WGS 84\",\n",
      "    DATUM[\"World Geodetic System 1984\",\n",
      "        ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n",
      "            LENGTHUNIT[\"metre\",1]]],\n",
      "    PRIMEM[\"Greenwich\",0,\n",
      "        ANGLEUNIT[\"degree\",0.0174532925199433]],\n",
      "    CS[ellipsoidal,2],\n",
      "        AXIS[\"geodetic latitude (Lat)\",north,\n",
      "            ORDER[1],\n",
      "            ANGLEUNIT[\"degree\",0.0174532925199433]],\n",
      "        AXIS[\"geodetic longitude (Lon)\",east,\n",
      "            ORDER[2],\n",
      "            ANGLEUNIT[\"degree\",0.0174532925199433]],\n",
      "    USAGE[\n",
      "        SCOPE[\"Horizontal component of 3D system.\"],\n",
      "        AREA[\"World.\"],\n",
      "        BBOX[-90,-180,90,180]],\n",
      "    ID[\"EPSG\",4326]]\n",
      "Data axis to CRS axis mapping: 2,1\n",
      "FID Column = fid\n",
      "Geometry Column = geom\n",
      "ID: Integer (0.0)\n",
      "b0: Real (0.0)\n",
      "b1: Real (0.0)\n",
      "b2: Real (0.0)\n",
      "b3: Real (0.0)\n",
      "b4: Real (0.0)\n",
      "b5: Real (0.0)\n",
      "b6: Real (0.0)\n",
      "b7: Real (0.0)\n",
      "b8: Real (0.0)\n",
      "b9: Real (0.0)\n",
      "b10: Real (0.0)\n",
      "b11: Real (0.0)\n",
      "OGRFeature(polygons_zonal):1\n",
      "  ID (Integer) = 2\n",
      "  b0 (Real) = 32.2190872990031\n",
      "  b1 (Real) = 34.2584177655869\n",
      "  b2 (Real) = 33.1255212012574\n",
      "  b3 (Real) = 29.7562015823403\n",
      "  b4 (Real) = 27.3043843524003\n",
      "  b5 (Real) = 26.0022204431177\n",
      "  b6 (Real) = 25.9803342467279\n",
      "  b7 (Real) = 26.0892614446006\n",
      "  b8 (Real) = 26.5121187622381\n",
      "  b9 (Real) = 26.9326012949042\n",
      "  b10 (Real) = 27.8735073791314\n",
      "  b11 (Real) = 29.8276078033663\n",
      "\n",
      "OGRFeature(polygons_zonal):2\n",
      "  ID (Integer) = 1\n",
      "  b0 (Real) = 37.5875141960014\n",
      "  b1 (Real) = 41.1516169535228\n",
      "  b2 (Real) = 38.9246022332857\n",
      "  b3 (Real) = 32.0599748302381\n",
      "  b4 (Real) = 28.703743676277\n",
      "  b5 (Real) = 27.9588523079267\n",
      "  b6 (Real) = 27.9912325730149\n",
      "  b7 (Real) = 28.0036972196164\n",
      "  b8 (Real) = 29.0760217569295\n",
      "  b9 (Real) = 28.6392086292018\n",
      "  b10 (Real) = 29.1736903942795\n",
      "  b11 (Real) = 32.6712335134043\n",
      "\n",
      "OGRFeature(polygons_zonal):3\n",
      "  ID (Integer) = 3\n",
      "  b0 (Real) = 27.8825487457743\n",
      "  b1 (Real) = 29.5533600442245\n",
      "  b2 (Real) = 29.5085259041395\n",
      "  b3 (Real) = 28.1868469851987\n",
      "  b4 (Real) = 26.6657154694727\n",
      "  b5 (Real) = 25.5369423150765\n",
      "  b6 (Real) = 25.4000528157838\n",
      "  b7 (Real) = 25.7769232047144\n",
      "  b8 (Real) = 26.3030576657987\n",
      "  b9 (Real) = 26.699380307982\n",
      "  b10 (Real) = 26.3889540332556\n",
      "  b11 (Real) = 26.7058148356776\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# check the results\n",
    "! ogrinfo -al -geom=NO geodata/LST/polygons_zonal_allmonth.gpkg"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This exercise can be also solved by building a for loop and extracting the zonal statistic for each file \n",
    "\n",
    "```\n",
    "for file in geodata/LST/LST_MOYDmax_month?.tif ; do\n",
    "filename=$(basename $file .tif)\n",
    "rm -f geodata/LST/${filename}_zonal.gpkg\n",
    "pkextractogr  -f \"GPKG\" -srcnodata -9999 -r mean \\\n",
    "-i $file -s  geodata/LST/polygons_zonal.gpkg -o geodata/LST/${filename}_zonal.gpkg\n",
    "done \n",
    "```\n",
    "\n",
    "Anyway is not the most efficient way because yuo create several files. Rather with the vrt option you concentrate the information in only one file. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4th excercise\n",
    "30% grading value\n",
    "\n",
    "Using the *geodata/shp/Lepidocolaptes_lacrymiger_allpoints.csv* at each location extract pixel value from *./geodata/cloud/SA_intra.tif  ./geodata/cloud/SA_meanannual.tif ./geodata/dem/SA_elevation_mn_GMTED2010_mn.tif ./geodata/vegetation/SA_tree_mn_percentage_GFC2013.tif*. \n",
    "The final result should be a txt file composed as such:  \n",
    "\n",
    "Lat Long SA_intra SA_meanannual SA_elevation SA_tree  \n",
    "-76.18925 3.98125       9           200           97     3\n",
    ".......  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "lon,lat,scientific_name\r",
      "\r\n",
      "-76.18925,3.98125,Lepidocolaptes_lacrymiger\r",
      "\r\n",
      "-76.18406,3.93442,Lepidocolaptes_lacrymiger\r",
      "\r\n",
      "-74.30256,4.60675,Lepidocolaptes_lacrymiger\r",
      "\r\n",
      "-74.30256,4.60675,Lepidocolaptes_lacrymiger\r",
      "\r\n",
      "-76.10394,4.74631,Lepidocolaptes_lacrymiger\r",
      "\r\n",
      "-76.13861,4.74536,Lepidocolaptes_lacrymiger\r",
      "\r\n",
      "-76.13919,4.7445,Lepidocolaptes_lacrymiger\r",
      "\r\n",
      "-76.10431,4.72342,Lepidocolaptes_lacrymiger\r",
      "\r\n",
      "-78.70159,-0.06506,Lepidocolaptes_lacrymiger\r",
      "\r\n"
     ]
    }
   ],
   "source": [
    "! head  geodata/shp/Lepidocolaptes_lacrymiger_allpoints.csv "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-76.18925 3.98125\r\n",
      "-76.18406 3.93442\r\n",
      "-74.30256 4.60675\r\n",
      "-74.30256 4.60675\r\n",
      "-76.10394 4.74631\r\n",
      "-76.13861 4.74536\r\n",
      "-76.13919 4.7445\r\n",
      "-76.10431 4.72342\r\n",
      "-78.70159 -0.06506\r\n",
      "-77.89315 -0.46052\r\n"
     ]
    }
   ],
   "source": [
    "# select only the lon and lat information\n",
    "! awk -F , '{ if(NR>1) print $1,$2}'  geodata/shp/Lepidocolaptes_lacrymiger_allpoints.csv > geodata/shp/Lepidocolaptes_lacrymiger_allpoints_x_y.txt\n",
    "! head geodata/shp/Lepidocolaptes_lacrymiger_allpoints_x_y.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At this point you two options: first using the vrt as input and the second one perform a for loop."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**First option**\n",
    "The first option is the most efficent way. Evrithing is stored in the memory and awk take care of the reformating"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0...10...20...30...40...50...60...70...80...90...100 - done.\n"
     ]
    }
   ],
   "source": [
    "# build a vrt having one tif for each band of the vrt using the -separate option\n",
    "! gdalbuildvrt -overwrite -separate  ./geodata/all_tif.vrt   ./geodata/cloud/SA_intra.tif ./geodata/cloud/SA_meanannual.tif ./geodata/dem/SA_elevation_mn_GMTED2010_mn.tif ./geodata/vegetation/SA_tree_mn_percentage_GFC2013.tif"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "223 8889 1253 6119\n",
      "217 9337 1593 8371\n",
      "85 9762 2574 7015\n",
      "85 9762 2574 7015\n",
      "392 8132 1705 8939\n",
      "305 8571 1744 7383\n",
      "305 8571 1744 7383\n",
      "373 8599 1750 8104\n",
      "191 9701 1669 9513\n",
      "275 9267 1835 4400\n"
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "# extract at each point and twist from 1 column to 4 column using the option awk 'ORS=NR%4?FS:RS'. \n",
    "# The number 4 identifies the number of bands in the vrt \n",
    "gdallocationinfo -geoloc -wgs84 -valonly ./geodata/all_tif.vrt < geodata/shp/Lepidocolaptes_lacrymiger_allpoints_x_y.txt \\\n",
    "| awk 'ORS=NR%4?FS:RS' >  geodata/shp/Lepidocolaptes_lacrymiger_allpoints_extract.txt\n",
    "head geodata/shp/Lepidocolaptes_lacrymiger_allpoints_extract.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Long Lat SA_intra SA_meanannual SA_elevation SA_tree\r\n",
      "-76.18925 3.98125 223 8889 1253 6119\r\n",
      "-76.18406 3.93442 217 9337 1593 8371\r\n",
      "-74.30256 4.60675 85 9762 2574 7015\r\n",
      "-74.30256 4.60675 85 9762 2574 7015\r\n",
      "-76.10394 4.74631 392 8132 1705 8939\r\n",
      "-76.13861 4.74536 305 8571 1744 7383\r\n",
      "-76.13919 4.7445 305 8571 1744 7383\r\n",
      "-76.10431 4.72342 373 8599 1750 8104\r\n",
      "-78.70159 -0.06506 191 9701 1669 9513\r\n"
     ]
    }
   ],
   "source": [
    "! echo \"Long Lat SA_intra SA_meanannual SA_elevation SA_tree\" > geodata/shp/Lepidocolaptes_lacrymiger_allpoints_x_y_extract.txt\n",
    "! paste -d \" \" geodata/shp/Lepidocolaptes_lacrymiger_allpoints_x_y.txt geodata/shp/Lepidocolaptes_lacrymiger_allpoints_extract.txt >> geodata/shp/Lepidocolaptes_lacrymiger_allpoints_x_y_extract.txt\n",
    "! head geodata/shp/Lepidocolaptes_lacrymiger_allpoints_x_y_extract.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Second option**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "for file in ./geodata/cloud/SA_intra.tif ./geodata/cloud/SA_meanannual.tif ./geodata/dem/SA_elevation_mn_GMTED2010_mn.tif ./geodata/vegetation/SA_tree_mn_percentage_GFC2013.tif ; do\n",
    "filename=$(basename $file .tif)\n",
    "gdallocationinfo -geoloc -wgs84 -valonly $file < geodata/shp/Lepidocolaptes_lacrymiger_allpoints_x_y.txt > geodata/shp/$filename.txt\n",
    "done"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At this point it is important to check if the files have the same number of lines and if there are not empity lines"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "223\n",
      "217\n",
      "85\n",
      "85\n",
      "392\n",
      "305\n",
      "305\n",
      "373\n",
      "191\n",
      "275\n",
      " 3438 geodata/shp/SA_intra.txt\n",
      " 3438 geodata/shp/SA_meanannual.txt\n",
      " 3438 geodata/shp/SA_elevation_mn_GMTED2010_mn.txt\n",
      " 3438 geodata/shp/SA_tree_mn_percentage_GFC2013.txt\n",
      "13752 total\n"
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "head geodata/shp/SA_intra.txt\n",
    "wc -l  geodata/shp/{SA_intra,SA_meanannual,SA_elevation_mn_GMTED2010_mn,SA_tree_mn_percentage_GFC2013}.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3438\n",
      "3438\n",
      "3438\n",
      "3438\n"
     ]
    }
   ],
   "source": [
    "! awk '{  if(NF==1) print  NF }' geodata/shp/SA_intra.txt | wc -l \n",
    "! awk '{  if(NF==1) print  NF }' geodata/shp/SA_meanannual.txt | wc -l \n",
    "! awk '{  if(NF==1) print  NF }' geodata/shp/SA_elevation_mn_GMTED2010_mn.txt | wc -l \n",
    "! awk '{  if(NF==1) print  NF }' geodata/shp/SA_tree_mn_percentage_GFC2013.txt | wc -l "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Long Lat SA_intra SA_meanannual SA_elevation SA_tree\n",
      "-76.18925 3.98125 223 8889 1253 6119\n",
      "-76.18406 3.93442 217 9337 1593 8371\n",
      "-74.30256 4.60675 85 9762 2574 7015\n",
      "-74.30256 4.60675 85 9762 2574 7015\n",
      "-76.10394 4.74631 392 8132 1705 8939\n",
      "-76.13861 4.74536 305 8571 1744 7383\n",
      "-76.13919 4.7445 305 8571 1744 7383\n",
      "-76.10431 4.72342 373 8599 1750 8104\n",
      "-78.70159 -0.06506 191 9701 1669 9513\n"
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "echo \"Long Lat SA_intra SA_meanannual SA_elevation SA_tree\" > geodata/shp/Lepidocolaptes_lacrymiger_allpoints_x_y_extract.txt\n",
    "paste -d \" \" geodata/shp/Lepidocolaptes_lacrymiger_allpoints_x_y.txt geodata/shp/{SA_intra,SA_meanannual,SA_elevation_mn_GMTED2010_mn,SA_tree_mn_percentage_GFC2013}.txt >> geodata/shp/Lepidocolaptes_lacrymiger_allpoints_x_y_extract.txt\n",
    "head geodata/shp/Lepidocolaptes_lacrymiger_allpoints_x_y_extract.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Third solution**\n",
    "Also this solution is very effective. The gdallocationinfo is computed on the flight and is kept in the memory and merged with the others using the paste command.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Long Lat SA_intra SA_meanannual SA_elevation SA_tree\n",
      "-76.18925 3.98125 223 8889 1253 6119\n",
      "-76.18406 3.93442 217 9337 1593 8371\n",
      "-74.30256 4.60675 85 9762 2574 7015\n",
      "-74.30256 4.60675 85 9762 2574 7015\n",
      "-76.10394 4.74631 392 8132 1705 8939\n",
      "-76.13861 4.74536 305 8571 1744 7383\n",
      "-76.13919 4.7445 305 8571 1744 7383\n",
      "-76.10431 4.72342 373 8599 1750 8104\n",
      "-78.70159 -0.06506 191 9701 1669 9513\n"
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "echo \"Long Lat SA_intra SA_meanannual SA_elevation SA_tree\" > geodata/shp/Lepidocolaptes_lacrymiger_allpoints_x_y_extract.txt\n",
    "paste -d \" \" geodata/shp/Lepidocolaptes_lacrymiger_allpoints_x_y.txt <(gdallocationinfo -geoloc -wgs84 -valonly ./geodata/cloud/SA_intra.tif < geodata/shp/Lepidocolaptes_lacrymiger_allpoints_x_y.txt) \\\n",
    "<(gdallocationinfo -geoloc -wgs84 -valonly ./geodata/cloud/SA_meanannual.tif < geodata/shp/Lepidocolaptes_lacrymiger_allpoints_x_y.txt) \\\n",
    "<(gdallocationinfo -geoloc -wgs84 -valonly ./geodata/dem/SA_elevation_mn_GMTED2010_mn.tif < geodata/shp/Lepidocolaptes_lacrymiger_allpoints_x_y.txt) \\\n",
    "<(gdallocationinfo -geoloc -wgs84 -valonly ./geodata/vegetation/SA_tree_mn_percentage_GFC2013.tif < geodata/shp/Lepidocolaptes_lacrymiger_allpoints_x_y.txt) \\\n",
    ">> geodata/shp/Lepidocolaptes_lacrymiger_allpoints_x_y_extract.txt\n",
    "head geodata/shp/Lepidocolaptes_lacrymiger_allpoints_x_y_extract.txt "
   ]
  }
 ],
 "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.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}