{ "cells": [ { "cell_type": "markdown", "metadata": { "id": "PBuSNMIUmJ8r" }, "source": [ "# Use PKTOOLS for raster/vector operations - osgeo\n", "\n", "---\n", "\n", "[Recorded lecture: 33:00 - 2.16.00](https://www.youtube.com/watch?v=atchsCMQf5A)\n", "\n", "---\n", "\n" ] }, { "cell_type": "markdown", "metadata": { "id": "kYp80rx3mJ8y" }, "source": [ "**Setting working directory for the OSGeoLive**" ] }, { "cell_type": "markdown", "metadata": { "colab": { "base_uri": "https://localhost:8080/" }, "id": "iwhEWdc8mJ80", "outputId": "ca7c9db2-f86f-4e73-f30d-9646681353c3" }, "source": [ "[PKTOOLS](http://pktools.nongnu.org/html/index.html) is a suite of utilities written in C++ for image processing with a focus on remote sensing applications. It relies on the Geospatial Data Abstraction Library (GDAL and OGR, http://www.gdal.org) (source http://pktools.nongnu.org/html/index.html)\n", "\n", "Open the bash terminal, migrate in the directory, and open the jupter-notebook.ipynb\n", "\n", " cd /media/sf_LVM_shared/my_SE_data/exercise\n", " jupyter-lab pktools_osgeo.ipynb\n", "\n", "We are going to use the jupyter-notebook as an editor for running bash-pktools commands. The beauty of using pktools in the bash environment is that we can combine it with other bash utilities such awk, sed, and other." ] }, { "cell_type": "markdown", "metadata": { "id": "Al7egxCLmJ81" }, "source": [ "## Start to use PKTOOLS\n", "\n", "**Create a mask**\n", "\n", "Create a mask by manipulating a text file." ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "id": "wb2l2rt5mJ81", "outputId": "4262b7a0-0fd9-4a42-ec86-84457c885e03" }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Input file size is 240, 240\n", "0...10...20...30...40...50...60...70...80...90...100 - done.\n", "Input file size is 240, 240\n", "0...10...20...30...40...50...60...70...80...90...100 - done.\n" ] } ], "source": [ "%%bash\n", "gdal_translate -of XYZ geodata/vegetation/ETmean08-11_crop.tif geodata/vegetation/ETmean08-11_crop.txt\n", "awk '{if ($3>2 && $3<4) {print $1,$2,50 } else {print}}' geodata/vegetation/ETmean08-11_crop.txt > geodata/vegetation/ETmean08-11_crop_trh.txt\n", "gdal_translate -ot Byte geodata/vegetation/ETmean08-11_crop_trh.txt geodata/vegetation/ETmean08-11_crop_trh.tif " ] }, { "cell_type": "markdown", "metadata": { "id": "SzC0xzJ0mJ82" }, "source": [ "The same operation can be done more efficient using pkgetmask. We can create 3 masks using different thresholds values." ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "id": "8oLciLK-mJ82", "outputId": "5f4e5c9b-f2d4-46f2-8f71-4490d5fb6bc9" }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0...10...20...30...40...50...60...70...80...90...100 - done.\n", "0...10...20...30...40...50...60...70...80...90...100 - done.\n", "0...10...20...30...40...50...60...70...80...90...100 - done.\n" ] } ], "source": [ "%%bash\n", "pkgetmask -co COMPRESS=DEFLATE -co ZLEVEL=9 -min 1 -max 2 -data 1 -nodata 0 -ot Byte -i geodata/vegetation/ETmean08-11.tif -o geodata/vegetation/ETmean08-11_01_trhA.tif\n", "pkgetmask -co COMPRESS=DEFLATE -co ZLEVEL=9 -min 5 -max 8 -data 1 -nodata 0 -ot Byte -i geodata/vegetation/ETmean08-11.tif -o geodata/vegetation/ETmean08-11_01_trhB.tif\n", "pkgetmask -co COMPRESS=DEFLATE -co ZLEVEL=9 -min 0 -max 10 -data 0 -nodata 1 -ot Byte -i geodata/vegetation/ETmean08-11.tif -o geodata/vegetation/ETmean08-11_01_trhC.tif" ] }, { "cell_type": "markdown", "metadata": { "id": "UC8VPm5CmJ83" }, "source": [ "**Set a mask to other image**\n", "\n", "Use the prepared mask and apply to the image." ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "id": "rOCybsZamJ83", "outputId": "970b90d4-5153-4b54-8940-dea47cfc8774", "scrolled": true }, "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 \\\n", "-m geodata/vegetation/ETmean08-11_01_trhA.tif -msknodata 1 -nodata -9 \\\n", "-m geodata/vegetation/ETmean08-11_01_trhB.tif -msknodata 1 -nodata -5 \\\n", "-m geodata/vegetation/ETmean08-11_01_trhC.tif -msknodata 1 -nodata -10 \\\n", "-i geodata/vegetation/ETmean08-11.tif -o geodata/vegetation/ETmean08-11_01_msk.tif" ] }, { "cell_type": "markdown", "metadata": { "id": "xWjYbcGimJ84" }, "source": [ "**Composite images**\n", "\n", "create a mask to apply during the composite" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "id": "qMEVO1RjmJ84", "outputId": "99d0d72f-4de1-43a8-d9f6-58fc221a5d30" }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0...10...20...30...40...50...60...70...80...90...100 - done.\n" ] } ], "source": [ "%%bash\n", "pkgetmask -co COMPRESS=DEFLATE -co ZLEVEL=9 -min 0 -max 25 -data 0 -nodata 1 -ot Byte -i geodata/LST/LST_MOYDmax_month1.tif -o geodata/LST/LST_MOYDmax_month1-msk.tif" ] }, { "cell_type": "markdown", "metadata": { "id": "c82Wn5bqmJ84" }, "source": [ "Calculate mean and standard deviation with several images" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "id": "th7g2X8omJ85", "outputId": "58bd6a7d-b1dd-4456-ec73-846df8af7c76" }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0...10...20...30...40...50...60...70...80...90...100 - done.\n", "0...10...20...30...40...50...60...70...80...90...100 - done.\n" ] } ], "source": [ "%%bash\n", "pkcomposite $(for file in geodata/LST/LST_MOYDmax_month??.tif geodata/LST/LST_MOYDmax_month?.tif; do echo -i $file; done) \\\n", "-m geodata/LST/LST_MOYDmax_month1-msk.tif -msknodata 0 -cr mean -dstnodata 0 \\\n", "-co COMPRESS=LZW -co ZLEVEL=9 -o geodata/LST/LST_MOYDmax_monthMean.tif\n", "\n", "pkcomposite $(for file in geodata/LST/LST_MOYDmax_month??.tif geodata/LST/LST_MOYDmax_month?.tif; do echo -i $file; done) \\\n", "-m geodata/LST/LST_MOYDmax_month1-msk.tif -msknodata 0 -cr stdev -dstnodata -1 \\\n", "-co COMPRESS=LZW -co ZLEVEL=9 -o geodata/LST/LST_MOYDmax_monthStdev.tif" ] }, { "cell_type": "markdown", "metadata": { "id": "8Jlxls9pmJ85" }, "source": [ "An alternative way is to use pkstatprofile" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "id": "foJoTUH0mJ85", "outputId": "40b097b4-c4a2-4ead-9a7d-c8726f825f76" }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0...10...20...30...40...50...60...70...80...90...100 - done.\n", "0...10...20...30...40...50...60...70...80...90...100 - done.\n" ] } ], "source": [ "%%bash\n", "# Create a multiband vrt\n", "gdalbuildvrt -overwrite -separate geodata/LST/LST_MOYDmax_month.vrt geodata/LST/LST_MOYDmax_month?.tif geodata/LST/LST_MOYDmax_month??.tif\n", "# Calculate mean and standard deviation\n", "pkstatprofile -co COMPRESS=LZW -nodata 0 -f mean -f stdev -i geodata/LST/LST_MOYDmax_month.vrt -o geodata/LST/LST_MOYDmax_month_mean_stdev.tif" ] }, { "cell_type": "markdown", "metadata": { "id": "Do0OMoGLmJ86" }, "source": [ "**Filter/Aggregate images**" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "id": "RH6nU0xTmJ86", "outputId": "656143a2-ce5f-42d8-9cd9-0d5a41efe287" }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "opening output image geodata/LST/LST_MOYDmax_monthMean_aggreate10mean.tif with 1 bands\n", "0...10...20...30...40...50...60...70...80...90...100 - done.\n", "opening output image geodata/LST/LST_MOYDmax_monthMean_circular11mean.tif with 1 bands\n", "0...10...20...30...40...50...60...70...80...90...100 - done.\n" ] } ], "source": [ "%%bash\n", "# mean aggregation \n", "pkfilter -co COMPRESS=DEFLATE -co ZLEVEL=9 -nodata 0 -ot Float32 -dx 10 -dy 10 -f mean -d 10 -i geodata/LST/LST_MOYDmax_monthMean.tif -o geodata/LST/LST_MOYDmax_monthMean_aggreate10mean.tif\n", "# mean circular moving window\n", "pkfilter -co COMPRESS=DEFLATE -co ZLEVEL=9 -nodata 0 -ot Float32 -dx 11 -dy 11 -f mean -c -i geodata/LST/LST_MOYDmax_monthMean.tif -o geodata/LST/LST_MOYDmax_monthMean_circular11mean.tif" ] }, { "cell_type": "markdown", "metadata": { "id": "vvikECXdmJ87" }, "source": [ "**Images statistic and histogram**" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "id": "PuVVTgQpmJ87", "outputId": "9997f3ee-750e-4191-9456-c58d8107b24c" }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0 0\n", "1 0\n", "2 0\n", "3 0\n", "4 0\n", "5 0\n", "6 0\n", "7 0\n", "8 0\n", "9 0\n" ] } ], "source": [ "%%bash\n", "pkstat -hist -src_min 0 -i geodata/temperature/ug_bio_3.tif > geodata/temperature/ug_bio_3.hist\n", "head geodata/temperature/ug_bio_3.hist" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "id": "6UzuZRORmJ87", "outputId": "f1bda511-4de5-4ade-ffc9-ffcdaa20c0aa" }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.1055006266 2058\n", "0.3165018797 3246\n", "0.5275031328 25726\n", "0.7385043859 54955\n", "0.9495056391 43453\n", "1.160506892 81369\n", "1.371508145 53690\n", "1.582509398 10960\n", "1.793510652 33823\n", "2.004511905 42521\n", "2.215513158 13744\n", "2.426514411 4183\n", "2.637515664 1539\n", "2.848516917 633\n", "3.05951817 278\n", "3.270519423 202\n", "3.481520677 252\n", "3.69252193 230\n", "3.903523183 33\n", "4.114524436 1\n" ] } ], "source": [ "%%bash\n", "pkstat -hist -nbin 20 -src_min 0 -i geodata/vegetation/GPPstdev08-11.tif" ] }, { "cell_type": "markdown", "metadata": { "id": "OYjlbLPnmJ87" }, "source": [ "**Images reclassification**" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "id": "2oFcBEgHmJ88", "outputId": "25dad608-fc04-4585-a3a0-2b56f9e13f59" }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0...10...20...30...40...50...60...70...80...90...100 - done.\n" ] } ], "source": [ "%%bash\n", "pkstat -hist -i geodata/temperature/ug_bio_3.tif | grep -v \" 0\" | awk '{ if ($1<75) {print $1, 0} else {print $1 , 1}}' > geodata/temperature/reclass_ug_bio_3.txt\n", "pkreclass -co COMPRESS=DEFLATE -co ZLEVEL=9 -code geodata/temperature/reclass_ug_bio_3.txt -i geodata/temperature/ug_bio_3.tif -o geodata/temperature/reclass_ug_bio_3.tif" ] }, { "cell_type": "markdown", "metadata": { "id": "UelAeY4JmJ88" }, "source": [ "**Zonal statistic (polygon extraction)**" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "id": "zGYUSaTnmJ88", "outputId": "f651593c-e991-4d5f-d2aa-4304dee10f0c" }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "processing layer polygons\n", "0...10...20...30...40...50...60...70...80...90...100 - done.\n", "processing layer polygons\n", "0...10...20...30...40...50...60...70...80...90...100 - done.\n", "processing layer polygons\n", "0...10...20...30...40...50...60...70...80...90...100 - done.\n" ] } ], "source": [ "%%bash\n", "rm -f geodata/shp/polygons_stat.*\n", "pkextractogr -srcnodata -339999995214436424907732413799364296704 -r mean -r stdev -r min -i geodata/vegetation/GPPmean08-11.tif -s geodata/shp/polygons.sqlite -o geodata/shp/polygons_stat.sqlite\n", "pkextractogr -f \"ESRI Shapefile\" -srcnodata -339999995214436424907732413799364296704 -r mean -r stdev -r min -i geodata/vegetation/GPPmean08-11.tif -s geodata/shp/polygons.sqlite -o geodata/shp/polygons_stat.shp\n", "\n", "# we can also create a csv that can be manipulate later on with awk\n", "rm -f geodata/shp/polygons_stat.csv\n", "pkextractogr -f CSV -srcnodata -339999995214436424907732413799364296704 -r mean -r stdev -r min -i geodata/vegetation/GPPmean08-11.tif -s geodata/shp/polygons.sqlite -o geodata/shp/polygons_stat.csv" ] }, { "cell_type": "markdown", "metadata": { "id": "xo_AY1bamJ89" }, "source": [ "**Zonal statistic (point extraction)**" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "id": "dQQp0FR9mJ89", "outputId": "72cb35dd-25f6-473e-c446-fd5bf6d34bb5", "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "processing layer presence\n", "0...10...20...30...40...50...60...70...80...90...100 - done.\n", "processing layer presence\n", "0...10...20...30...40...50...60...70...80...90...100 - done.\n" ] } ], "source": [ "%%bash \n", "# at point location\n", "rm -f geodata/shp/point_stat.csv\n", "pkextractogr -f CSV -srcnodata -339999995214436424907732413799364296704 -r mean -r stdev -r min -i geodata/vegetation/GPPmean08-11.tif -s geodata/shp/presence.shp -o geodata/shp/point_stat.csv\n", "# at point location + 1 pixel around \n", "rm -f geodata/shp/point-buf_stat.csv\n", "pkextractogr -f CSV -buf 2 -srcnodata -339999995214436424907732413799364296704 -r mean -r stdev -r min -i geodata/vegetation/GPPmean08-11.tif -s geodata/shp/presence.shp -o geodata/shp/point-buf_stat.csv" ] }, { "cell_type": "markdown", "metadata": { "id": "gYgDnke_mJ89" }, "source": [ "**Remove all the output**" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "id": "mWC8G93emJ8-" }, "outputs": [], "source": [ "%%bash\n", "rm -f geodata/vegetation/GPPcv08-11.tif geodata/LST/*_crop.tif geodata/vegetation/ETmean08-11_crop_trh.tif geodata/vegetation/ETmean08-11_crop_trh.txt vegetation/ETmean08-11_crop.txt geodata/vegetation/ETmosaic.vrt geodata/vegetation/ETmosaic.tif geodata/vegetation/stack_??.tif geodata/vegetation/stack.vrt geodata/vegetation/tiles.* geodata/vegetation/ETmean08-11_crop_proximity.tif geodata/vegetation/ETmean08-11_crop_buffer.tif geodata/dem/slope.tif geodata/dem/aspect.tif geodata/dem/tri.tif geodata/dem/tpi.tif geodata/dem/roughness.tif geodata/vegetation/ETmean08-11_01_trh?.tif geodata/LST/LST_MOYDmax_month1-msk.tif geodata/LST/LST_MOYDmax_monthStdev.tif geodata/LST/LST_MOYDmax_monthMean.tif geodata/LST/LST_MOYDmax_month_mean_stdev.tif geodata/LST/LST_MOYDmax_month.vrt geodata/LST/LST_MOYDmax_monthMean_aggreate10mean.tif geodata/LST/LST_MOYDmax_monthMean_circular11mean.tif geodata/temperature/reclass_ug_bio_3.tif geodata/temperature/reclass_ug_bio_3.txt geodata/shp/polygons_stat.csv geodata/shp/point-buf_stat.csv geodata/shp/point_stat.csv geodata/shp/polygons_stat.* geodata/shp/TM_LARGE_BORDERS.shp.* geodata/shp/TM_UGANDA_BORDERS-0.3.* geodata/vegetation/ETmean08-11_crop.txt" ] } ], "metadata": { "colab": { "collapsed_sections": [], "name": "02_pktools.ipynb", "provenance": [], "toc_visible": true }, "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.8.10" } }, "nbformat": 4, "nbformat_minor": 4 }