Use PKTOOLS for raster/vector operations - colab

Setting working directory for the Google Collaboratory

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)

[1]:
from google.colab import drive
drive.mount('/gdrive')
Mounted at /gdrive

List the mounted gdrive. If you have folders and files in you gdrive you should be able to see them

[2]:
! ls /gdrive/MyDrive
 auto
 Big_blue
'Colab Notebooks'
'Dataset for global stream network and MOSQLAND.gdoc'
'Discharge Tasks.gdoc'
 Elena_intro.gdoc
 GeoCompCourse
 GeoCompYaleMay2020
 grants
'ID PASSPORT'
 LandscapeGenetic
 Lawsuit
 Margosa
 panniers-trailers.pdf
 RDOC_CMS
 review
 Sofia_intro.gdoc
 Stream_Variables_dataset_selected.gsheet
 Taino
 yale

Download the SE_data folder yousing git

First remove the directory if already exist, then enter in the directory.

[3]:
%rm -rf /gdrive/MyDrive/SE_data
%cd /gdrive/MyDrive
%ls
/gdrive/MyDrive
 auto/
 Big_blue/
'Colab Notebooks'/
'Dataset for global stream network and MOSQLAND.gdoc'
'Discharge Tasks.gdoc'
 Elena_intro.gdoc
 GeoCompCourse/
 GeoCompYaleMay2020/
 grants/
'ID PASSPORT'/
 LandscapeGenetic/
 Lawsuit/
 Margosa/
 panniers-trailers.pdf
 RDOC_CMS/
 review/
 Sofia_intro.gdoc
 Stream_Variables_dataset_selected.gsheet
 Taino/
 yale/

Then perform the downloading

[4]:
!git clone https://github.com/selvaje/SE_data.git
Cloning into 'SE_data'...
remote: Enumerating objects: 195, done.
remote: Counting objects: 100% (195/195), done.
remote: Compressing objects: 100% (155/155), done.
remote: Total 262 (delta 22), reused 190 (delta 20), pack-reused 67
Receiving objects: 100% (262/262), 78.99 MiB | 17.46 MiB/s, done.
Resolving deltas: 100% (25/25), done.
Checking out files: 100% (162/162), done.

Check if the SE_data folder is full

[5]:
! ls -l SE_data/*
-rw------- 1 root root   58 Jan 19 22:45 SE_data/README.md

SE_data/exercise:
total 187
-rw------- 1 root root   8591 Jan 19 22:45 00_Setting_Colab_for_for_Spatial_Ecology_course.ipynb
-rw------- 1 root root 153327 Jan 19 22:45 01_gdal.ipynb
-rw------- 1 root root  24271 Jan 19 22:45 02_pktools.ipynb
drwx------ 7 root root   4096 Jan 19 22:45 geodata

SE_data/pktools_local:
total 15
-rw------- 1 root root 11141 Jan 19 22:45 00_pktools_gdrive_install.ipynb
drwx------ 5 root root  4096 Jan 19 22:45 usr

Install pktools

pktools is not present in Colab but is possible to install via the “apt install” comand. Pay attention that if you re-open another colab you have to re-install pktools

[7]:
! apt install pktools
Reading package lists... Done
Building dependency tree
Reading state information... Done
The following additional packages will be installed:
  libalgorithms1 libfann2 libfileclasses1 libgsl23 libgslcblas0
  libimageclasses1 liblas3 liblasclasses1
Suggested packages:
  libfann-dev libfann-doc gsl-ref-psdoc | gsl-doc-pdf | gsl-doc-info
  | gsl-ref-html libgeotiff-epsg
The following NEW packages will be installed:
  libalgorithms1 libfann2 libfileclasses1 libgsl23 libgslcblas0
  libimageclasses1 liblas3 liblasclasses1 pktools
0 upgraded, 9 newly installed, 0 to remove and 12 not upgraded.
Need to get 2,347 kB of archives.
After this operation, 9,702 kB of additional disk space will be used.
Get:1 http://archive.ubuntu.com/ubuntu bionic/universe amd64 libgslcblas0 amd64 2.4+dfsg-6 [79.7 kB]
Get:2 http://archive.ubuntu.com/ubuntu bionic/universe amd64 libgsl23 amd64 2.4+dfsg-6 [823 kB]
Get:3 http://archive.ubuntu.com/ubuntu bionic/universe amd64 libimageclasses1 amd64 2.6.7.3+ds-1 [57.9 kB]
Get:4 http://archive.ubuntu.com/ubuntu bionic/universe amd64 libalgorithms1 amd64 2.6.7.3+ds-1 [130 kB]
Get:5 http://archive.ubuntu.com/ubuntu bionic/universe amd64 libfann2 amd64 2.2.0+ds-3 [64.5 kB]
Get:6 http://archive.ubuntu.com/ubuntu bionic/universe amd64 libfileclasses1 amd64 2.6.7.3+ds-1 [12.7 kB]
Get:7 http://archive.ubuntu.com/ubuntu bionic/universe amd64 liblas3 amd64 1.8.1-6build1 [270 kB]
Get:8 http://archive.ubuntu.com/ubuntu bionic/universe amd64 liblasclasses1 amd64 2.6.7.3+ds-1 [15.2 kB]
Get:9 http://archive.ubuntu.com/ubuntu bionic/universe amd64 pktools amd64 2.6.7.3+ds-1 [894 kB]
Fetched 2,347 kB in 2s (1,435 kB/s)
Selecting previously unselected package libgslcblas0:amd64.
(Reading database ... 146364 files and directories currently installed.)
Preparing to unpack .../0-libgslcblas0_2.4+dfsg-6_amd64.deb ...
Unpacking libgslcblas0:amd64 (2.4+dfsg-6) ...
Selecting previously unselected package libgsl23:amd64.
Preparing to unpack .../1-libgsl23_2.4+dfsg-6_amd64.deb ...
Unpacking libgsl23:amd64 (2.4+dfsg-6) ...
Selecting previously unselected package libimageclasses1.
Preparing to unpack .../2-libimageclasses1_2.6.7.3+ds-1_amd64.deb ...
Unpacking libimageclasses1 (2.6.7.3+ds-1) ...
Selecting previously unselected package libalgorithms1.
Preparing to unpack .../3-libalgorithms1_2.6.7.3+ds-1_amd64.deb ...
Unpacking libalgorithms1 (2.6.7.3+ds-1) ...
Selecting previously unselected package libfann2:amd64.
Preparing to unpack .../4-libfann2_2.2.0+ds-3_amd64.deb ...
Unpacking libfann2:amd64 (2.2.0+ds-3) ...
Selecting previously unselected package libfileclasses1.
Preparing to unpack .../5-libfileclasses1_2.6.7.3+ds-1_amd64.deb ...
Unpacking libfileclasses1 (2.6.7.3+ds-1) ...
Selecting previously unselected package liblas3.
Preparing to unpack .../6-liblas3_1.8.1-6build1_amd64.deb ...
Unpacking liblas3 (1.8.1-6build1) ...
Selecting previously unselected package liblasclasses1.
Preparing to unpack .../7-liblasclasses1_2.6.7.3+ds-1_amd64.deb ...
Unpacking liblasclasses1 (2.6.7.3+ds-1) ...
Selecting previously unselected package pktools.
Preparing to unpack .../8-pktools_2.6.7.3+ds-1_amd64.deb ...
Unpacking pktools (2.6.7.3+ds-1) ...
Setting up liblas3 (1.8.1-6build1) ...
Setting up liblasclasses1 (2.6.7.3+ds-1) ...
Setting up libfann2:amd64 (2.2.0+ds-3) ...
Setting up libfileclasses1 (2.6.7.3+ds-1) ...
Setting up libgslcblas0:amd64 (2.4+dfsg-6) ...
Setting up libgsl23:amd64 (2.4+dfsg-6) ...
Setting up libimageclasses1 (2.6.7.3+ds-1) ...
Setting up libalgorithms1 (2.6.7.3+ds-1) ...
Setting up pktools (2.6.7.3+ds-1) ...
Processing triggers for libc-bin (2.27-3ubuntu1.3) ...
/sbin/ldconfig.real: /usr/local/lib/python3.6/dist-packages/ideep4py/lib/libmkldnn.so.0 is not a symbolic link

Processing triggers for man-db (2.8.3-2ubuntu0.1) ...

Test if pktools is working fine by running a pktools command and getting the help

[8]:
! pkgetmask --help
   -i       --input                Input image file
   -o       --output               Output mask file
   -min     --min                  Values smaller than min threshold(s) are masked as invalid. Use one threshold for each band
   -max     --max                  Values greater than max threshold(s) are masked as invalid. Use one threshold for each band
   -data    --data                 value(s) for valid pixels: between min and max (default: 1)
   -nodata   --nodata               value(s) for invalid pixels: not between min and max (default: 0)
   -b       --band                 band(s) used for mask (default: 0)
   -p       --operator             Operator: [AND,OR]. (default: OR)
   -ot      --otype                Data type for output image ({Byte/Int16/UInt16/UInt32/Int32/Float32/Float64/CInt16/CInt32/CFloat32/CFloat64}). Empty string: inherit type from input image (default: Byte)
   -of      --oformat              Output image format (see also gdal_translate). (default: GTiff)
   -co      --co                   Creation option for output file. Multiple options can be specified.
   -ct      --ct                   color table (file with 5 columns: id R G B ALFA (0: transparent, 255: solid)

Usage: pkgetmask -i input -o output

short option -h shows basic options only, use long option --help to show all options

Start to use PKTOOLS

Create a mask

Create a mask with by manipulate a text file.

[ ]:
%%bash
cd /gdrive/MyDrive/SE_data
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.

[ ]:
%%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.
[ ]:
%%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
[ ]:
%%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

[ ]:
%%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

[ ]:
%%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

[ ]:
%%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

[ ]:
%%bash
pkstat -hist  -src_min 0  -i    geodata/temperature/ug_bio_3.tif
0 0
1 0
2 0
3 0
4 0
5 0
6 0
7 0
8 0
9 0
10 0
11 0
12 0
13 0
14 0
15 0
16 0
17 0
18 0
19 0
20 0
21 0
22 0
23 0
24 0
25 0
26 0
27 0
28 0
29 0
30 0
31 0
32 0
33 0
34 0
35 0
36 0
37 0
38 0
39 0
40 0
41 0
42 0
43 0
44 0
45 0
46 0
47 0
48 0
49 0
50 0
51 0
52 0
53 0
54 0
55 0
56 0
57 0
58 0
59 0
60 0
61 2
62 9
63 31
64 138
65 235
66 288
67 481
68 589
69 803
70 1989
71 3471
72 4490
73 5513
74 6938
75 12332
76 21295
77 23643
78 20937
79 24687
80 36864
81 32858
82 35070
83 33486
84 30120
85 25312
86 19138
87 11203
88 5892
89 5799
90 4372
91 5354
92 2456
93 168
[ ]:
%%bash
pkstat -hist  -nbin  100 -src_min 0  -i    geodata/vegetation/GPPstdev08-11.tif
0.02110012531 616
0.06330037594 0
0.1055006266 30
0.1477008772 624
0.1899011278 650
0.2321013784 555
0.2743016291 589
0.3165018797 602
0.3587021303 634
0.4009023809 695
0.4431026316 761
0.4853028822 1072
0.5275031328 2287
0.5697033834 5055
0.6119036341 9187
0.6541038847 12469
0.6963041353 12988
0.7385043859 11639
0.7807046366 10078
0.8229048872 8994
0.8651051378 7597
0.9073053885 6994
0.9495056391 6978
0.9917058897 8092
1.03390614 9392
1.076106391 10288
1.118306642 11877
1.160506892 13784
1.202707143 15860
1.244907393 17719
1.287107644 18117
1.329307895 16990
1.371508145 14917
1.413708396 11927
1.455908647 8212
1.498108897 4854
1.540309148 2430
1.582509398 1295
1.624709649 1298
1.6669099 1862
1.70911015 2812
1.751310401 3990
1.793510652 5286
1.835710902 6235
1.877911153 6864
1.920111403 7441
1.962311654 7356
2.004511905 7345
2.046712155 8170
2.088912406 9272
2.131112657 9107
2.173312907 7296
2.215513158 4944
2.257713408 3490
2.299913659 2446
2.34211391 1688
2.38431416 1415
2.426514411 1180
2.468714662 995
2.510914912 860
2.553115163 706
2.595315413 564
2.637515664 441
2.679715915 347
2.721916165 318
2.764116416 269
2.806316667 256
2.848516917 217
2.890717168 170
2.932917418 115
2.975117669 112
3.01731792 107
3.05951817 86
3.101718421 63
3.143918672 44
3.186118922 61
3.228319173 53
3.270519423 44
3.312719674 45
3.354919925 41
3.397120175 31
3.439320426 37
3.481520677 42
3.523720927 38
3.565921178 44
3.608121428 55
3.650321679 54
3.69252193 51
3.73472218 46
3.776922431 62
3.819122682 45
3.861322932 49
3.903523183 34
3.945723433 17
3.987923684 15
4.030123935 6
4.072324185 9
4.114524436 1
4.156724687 0
4.198924937 1

Images reclassification

[ ]:
%%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)

[ ]:
%%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)

[ ]:
%%bash
# at point location
rm 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 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

[ ]:
%%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/polygons.sqlite 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