Use PKTOOLS for raster/vector operations - osgeo

Setting working directory for the OSGeoLive

PKTOOLS 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)

Open the bash terminal, migrate in the directory, and open the jupter-notebook.ipynb

cd /media/sf_LVM_shared/my_SE_data/exercise
jupyter-lab pktools_osgeo.ipynb

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.

Start to use PKTOOLS

Create a mask

Create a mask by manipulating a text file.

[1]:

%%bash
gdal_translate -of  XYZ   geodata/vegetation/ETmean08-11_crop.tif geodata/vegetation/ETmean08-11_crop.txt
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
gdal_translate -ot Byte geodata/vegetation/ETmean08-11_crop_trh.txt geodata/vegetation/ETmean08-11_crop_trh.tif

Input file size is 240, 240
0...10...20...30...40...50...60...70...80...90...100 - done.
Input file size is 240, 240
0...10...20...30...40...50...60...70...80...90...100 - done.

The same operation can be done more efficient using pkgetmask. We can create 3 masks using different thresholds values.

[2]:

%%bash
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
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
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

0...10...20...30...40...50...60...70...80...90...100 - done.
0...10...20...30...40...50...60...70...80...90...100 - done.
0...10...20...30...40...50...60...70...80...90...100 - done.

Set a mask to other image

Use the prepared mask and apply to the image.

[3]:

%%bash
pksetmask -co COMPRESS=DEFLATE -co ZLEVEL=9 \
-m geodata/vegetation/ETmean08-11_01_trhA.tif  -msknodata 1 -nodata  -9 \
-m geodata/vegetation/ETmean08-11_01_trhB.tif  -msknodata 1 -nodata  -5 \
-m geodata/vegetation/ETmean08-11_01_trhC.tif  -msknodata 1 -nodata -10 \
-i geodata/vegetation/ETmean08-11.tif -o geodata/vegetation/ETmean08-11_01_msk.tif

0...10...20...30...40...50...60...70...80...90...100 - done.

Composite images

create a mask to apply during the composite

[4]:

%%bash
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

0...10...20...30...40...50...60...70...80...90...100 - done.

Calculate mean and standard deviation with several images

[7]:

%%bash
pkcomposite $(for file in geodata/LST/LST_MOYDmax_month??.tif geodata/LST/LST_MOYDmax_month?.tif; do echo -i $file; done) \
-m geodata/LST/LST_MOYDmax_month1-msk.tif -msknodata 0 -cr mean   -dstnodata 0 \
-co  COMPRESS=LZW -co ZLEVEL=9 -o geodata/LST/LST_MOYDmax_monthMean.tif

pkcomposite $(for file in geodata/LST/LST_MOYDmax_month??.tif geodata/LST/LST_MOYDmax_month?.tif; do echo -i $file; done) \
-m geodata/LST/LST_MOYDmax_month1-msk.tif -msknodata 0 -cr stdev   -dstnodata -1 \
-co  COMPRESS=LZW -co ZLEVEL=9 -o geodata/LST/LST_MOYDmax_monthStdev.tif

0...10...20...30...40...50...60...70...80...90...100 - done.
0...10...20...30...40...50...60...70...80...90...100 - done.

An alternative way is to use pkstatprofile

[8]:

%%bash
# Create a multiband vrt
gdalbuildvrt -overwrite -separate geodata/LST/LST_MOYDmax_month.vrt geodata/LST/LST_MOYDmax_month?.tif geodata/LST/LST_MOYDmax_month??.tif
# Calculate mean and standard deviation
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

0...10...20...30...40...50...60...70...80...90...100 - done.
0...10...20...30...40...50...60...70...80...90...100 - done.

Filter/Aggregate images

[9]:

%%bash
# mean aggregation
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
# mean circular moving window
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

opening output image geodata/LST/LST_MOYDmax_monthMean_aggreate10mean.tif with 1 bands
0...10...20...30...40...50...60...70...80...90...100 - done.
opening output image geodata/LST/LST_MOYDmax_monthMean_circular11mean.tif with 1 bands
0...10...20...30...40...50...60...70...80...90...100 - done.

Images statistic and histogram

[10]:

%%bash
pkstat -hist  -src_min 0  -i    geodata/temperature/ug_bio_3.tif > geodata/temperature/ug_bio_3.hist
head geodata/temperature/ug_bio_3.hist

[11]:

%%bash
pkstat -hist  -nbin  20 -src_min 0  -i    geodata/vegetation/GPPstdev08-11.tif

0.1055006266 2058
0.3165018797 3246
0.5275031328 25726
0.7385043859 54955
0.9495056391 43453
1.160506892 81369
1.371508145 53690
1.582509398 10960
1.793510652 33823
2.004511905 42521
2.215513158 13744
2.426514411 4183
2.637515664 1539
2.848516917 633
3.05951817 278
3.270519423 202
3.481520677 252
3.69252193 230
3.903523183 33
4.114524436 1

Images reclassification

[12]:

%%bash
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
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

0...10...20...30...40...50...60...70...80...90...100 - done.

Zonal statistic (polygon extraction)

[13]:

%%bash
rm -f geodata/shp/polygons_stat.*
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
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

# we can also create a csv that can be manipulate later on with awk
rm  -f geodata/shp/polygons_stat.csv
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

processing layer polygons
0...10...20...30...40...50...60...70...80...90...100 - done.
processing layer polygons
0...10...20...30...40...50...60...70...80...90...100 - done.
processing layer polygons
0...10...20...30...40...50...60...70...80...90...100 - done.

Zonal statistic (point extraction)

[14]:

%%bash
# at point location
rm -f geodata/shp/point_stat.csv
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
# at point location + 1 pixel around
rm -f geodata/shp/point-buf_stat.csv
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

processing layer presence
0...10...20...30...40...50...60...70...80...90...100 - done.
processing layer presence
0...10...20...30...40...50...60...70...80...90...100 - done.

Remove all the output

[15]:

%%bash
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