User Tools

Site Tools


wiki:gdalpktools_kenya

Differences

This shows you the differences between two versions of the page.

Link to this comparison view

Both sides previous revision Previous revision
Previous revision
wiki:gdalpktools_kenya [2019/09/26 18:58]
wiki:gdalpktools_kenya [2020/05/12 11:28] (current)
new_ost4sem_admin [PKTOOLS Commands]
Line 1: Line 1:
 +====== Use GDAL/OGR and PKTOOLS for raster/​vector operations ======
 +
 +==== Enter in the AWS Ubuntu Machine ​ ====
 +\\
 +
 +On **Linux or Mac**:
 +  ​
 +  wget -N http://​www.spatial-ecology.net/​ost4sem/​exercise/​id_rsa_class ​ -O id_rsa_class
 +  # Mac user can use this
 +  # curl -o id_rsa_class ​  ​http://​www.spatial-ecology.net/​ost4sem/​exercise/​id_rsa_class
 +  chmod 600 id_rsa_class
 +  # login using the ssh. Change user???? in accordance to your user ID 
 +  ssh -i  id_rsa_class -X -Y user??​@52.138.103.28 ​
 +
 +
 +On **Windows**:​\\
 +
 +Download [[ https://​www.puttygen.com/​ | PuTTYgen ]] - copy the key in a text file and save as id_rsa_class - and import the [[  http://​www.spatial-ecology.net/​ost4sem/​exercise/​id_rsa_class | id_rsa_class ]] file.\\
 +Use [[ https://​putty.org | PuTTY ]]  to ssh to uuser??​@52.138.103.28.
 +
 +==== In the Ubuntu VM ====
 +Data Directory: $HOME/​ost4sem/​exercise/​geodata\_SDM/​ \\
 +
 +  ##### Download the data ONLY if your following the exercise in your machine. ​
 +  ##### I you are using the Linux Virtuan Machine you do not need to dowload the data. 
 +  mkdir -p $HOME/​ost4sem/​exercise
 +  cd $HOME/​ost4sem/​exercise
 +  wget -N http://​spatial-ecology.net/​ost4sem/​exercise/​geodata_SDM.tar.gz
 +  tar xvzf geodata_SDM.tar.gz
 +  cd $HOME/​ost4sem/​exercise/​geodata_SDM
 +  wget -N http://​spatial-ecology.net/​ost4sem/​exercise/​geodata_SDM/​gdal-ogr-pktools.sh
 +  emacs -nw  $HOME/​ost4sem/​exercise/​geodata_SDM/​gdal-ogr-pktools.sh
 +  ## gedit $HOME/​ost4sem/​exercise/​geodata_SDM/​gdal-ogr-pktools.sh &
 +
 +===== GDAL Commands =====
 +
 +=== Change the directory, and explore the files ===
 +
 +<code bash>
 +cd $HOME/​ost4sem/​exercise/​geodata_SDM
 +ls */​*.tif ​ # list all the files with the extension .tif
 +gdalinfo vegetation/​ETmean08-11.tif
 +
 +# Visualize the image
 +
 +openev vegetation/​ETmean08-11.tif
 +</​code>​
 +
 +Reply to the following questions:​\\
 +What is the pixel size?\\
 +How do you get min and max pixel values?\\
 +
 +=== Understand data type ===
 +  ​
 +                        Ranges of GDAL data types        Image
 +  GDAL data type        ​minimum ​ maximum ​       Size 
 +  Byte       ​0 ​     255        39M
 +  UInt16  ​    ​0 ​ 65,​535 ​       78M
 +  Int16, CInt16        ​-32,​768 ​ 32,​767 ​       78M
 +  UInt32  ​    ​0 ​   4,​294,​967,​295 ​       155M  ​
 +  Int32, CInt32 -2,​147,​483,​648 ​   2,​147,​483,​647 ​       155M
 +  Float32, CFloat32 ​     -3.4E38 ​          ​3.4E38 ​       155M
 +  Float64, CFloat64 ​   -1.79E308 ​        ​1.79E308 ​       309M
 +
 +
 +=== Understand NoData Value ===
 +
 +<code bash>
 +gdalinfo -mm vegetation/​ETmean08-11.tif
 +gdal_edit.py -a_nodata -9999 vegetation/​ETmean08-11.tif
 +gdalinfo -mm vegetation/​ETmean08-11.tif
 +gdal_edit.py -a_nodata -339999995214436424907732413799364296704.000 ​   vegetation/​ETmean08-11.tif
 +</​code>​
 +
 +  ​
 +=== Calculate minimum and maximum values for all the images ===
 +<code bash>
 +for file in vegetation/​*.tif ; do 
 +gdalinfo ​ -mm  $file  | grep Computed
 +done
 +</​code>​
 +
 +=== Exercise: ===
 +Use one of the gdal commands and try to clip the vegetation/​ETmean08-11.tif image using a square box (numbers are in degree):\\
 +xmin : 32\\
 +xmax : 34\\
 +ymin : 1\\
 +ymax : 3\\
 +\\
 +
 +   ​gdal............................ ​ vegetation/​ETmean08-11.tif vegetation/​ETmean08-11_crop.tif
 +
 +
 +=== Create a Coefficient of variation ​ ===
 +
 +GPPmean08-11.tif temporal mean of Gross Primary Productivity \\
 +GPPstdev08-11.tif temporal standard deviation of Gross Primary Productivity \\
 +
 +<code bash>
 +gdal_calc.py --type=Float32 --overwrite ​ -A vegetation/​GPPstdev08-11.tif -B  vegetation/​GPPmean08-11.tif \
 +   ​--calc="​( A.astype(float) / ( B.astype(float) + 0.000000001 ) )" --outfile=vegetation/​GPPcv08-11.tif
 +</​code>​
 +
 +=== Change pixel resolution ===
 +<code bash>
 +
 +# looping trough the images
 +for file in LST/​LST_MOYDmax_month?​.tif LST/​LST_MOYDmax_month??​.tif;​ do 
 +  filename=$(basename $file .tif )
 +  gdalwarp -overwrite -te 29 -1 35.0000000048 4.000000004 -tr 0.008333333340000 0.008333333340000 -co COMPRESS=DEFLATE -co ZLEVEL=9 $file LST/​${filename}_crop.tif ​
 +done
 +</​code>​
 +
 +=== Get value at one pixel/line image-location ===
 +
 +<code bash>
 +# looping trough the images
 +for file in LST/​LST_MOYDmax_month?​.tif LST/​LST_MOYDmax_month??​.tif;​ do 
 +   ​gdallocationinfo -valonly $file  20 20  ​
 +done
 +</​code>​
 +
 +=== Get value at lat/​long ​ image-location ===
 +
 +<code bash>
 +# looping trough the images
 +for file in LST/​LST_MOYDmax_month?​.tif LST/​LST_MOYDmax_month??​.tif;​ do 
 +   ​gdallocationinfo -valonly -geoloc $file  33 2  ​
 +done
 +</​code>​
 +
 +=== Get value at multiple lat/​long ​ image-location === 
 +
 +<code bash>
 +# create the lat long file 
 +echo 32.5 2.5 > LST/x_y.txt
 +echo 31.1 2.1 >> LST/x_y.txt
 +# looping trough the images
 +for file in LST/​LST_MOYDmax_month?​.tif LST/​LST_MOYDmax_month??​.tif;​ do 
 +   ​gdallocationinfo -valonly -geoloc $file < LST/​x_y.txt  ​
 +   echo ""​
 +done
 +</​code>​
 +
 +=== From Image to text and from txt to image ===
 +
 +<code bash>
 +gdal_translate -of  XYZ   ​vegetation/​ETmean08-11_crop.tif vegetation/​ETmean08-11_crop.txt
 +awk '{if ($3>2 && $3<4) {print $1,$2,50 } else {print}}'​ vegetation/​ETmean08-11_crop.txt > vegetation/​ETmean08-11_crop_trh.txt
 +gdal_translate -ot Byte vegetation/​ETmean08-11_crop_trh.txt vegetation/​ETmean08-11_crop_trh.tif ​
 +</​code>​
 +
 +=== The use of VRT to stack images and tiling ===
 +
 +<code bash>
 +# stack the tif
 +gdalbuildvrt -overwrite -separate ​  ​vegetation/​stack.vrt ​  ​vegetation/​ETmean08-11.tif vegetation/​ETstdev08-11.tif
 +gdalinfo vegetation/​stack.vrt
 +
 +# split in tiles 
 +gdal_translate -srcwin 0     0 360 300 vegetation/​stack.vrt vegetation/​stack_UL.tif
 +gdal_translate -srcwin 0   300 360 300 vegetation/​stack.vrt vegetation/​stack_UR.tif
 +gdal_translate -srcwin 360   0 360 300 vegetation/​stack.vrt vegetation/​stack_LL.tif
 +gdal_translate -srcwin 360 300 360 300 vegetation/​stack.vrt vegetation/​stack_LR.tif
 +  ​
 +# recompose the image 
 +gdalbuildvrt -overwrite vegetation/​ETmosaic.vrt vegetation/​stack_UL.tif vegetation/​stack_UR.tif vegetation/​stack_LL.tif vegetation/​stack_LR.tif
 +gdal_translate -co COMPRESS=DEFLATE -co ZLEVEL=9 vegetation/​ETmosaic.vrt vegetation/​ETmosaic.tif
 +</​code>​
 +
 +=== Create shp border tiles  ===
 +<code bash>
 +rm -f vegetation/​tiles.*
 +gdaltindex ​ vegetation/​tiles.shp ​ vegetation/​stack_??​.tif
 +</​code>​
 +
 +=== Surface distance/​buffer ===
 +<code bash>
 +# Continues distance surface
 +gdal_proximity.py -co COMPRESS=DEFLATE -co ZLEVEL=9 -values 0 -distunits PIXEL  -ot UInt32 ​  ​vegetation/​ETmean08-11_crop_trh.tif vegetation/​ETmean08-11_crop_proximity.tif
 +
 +# Fix buffer surface
 +gdal_proximity.py -fixed-buf-val 4 -maxdist 4 -nodata 10 -co COMPRESS=DEFLATE -co ZLEVEL=9 -values 0 -distunits PIXEL -ot Byte vegetation/​ETmean08-11_crop_trh.tif vegetation/​ETmean08-11_crop_buffer.tif
 +</​code>​
 +
 +=== Topography variables ===
 +<code bash>
 +# calculate ​ slope 
 +gdaldem slope -s 111120 -co COMPRESS=DEFLATE -co ZLEVEL=9 dem/​GMTED2010.tif dem/​slope.tif ​
 + 
 +# calculate ​ apect
 +gdaldem aspect -zero_for_flat -co COMPRESS=DEFLATE -co ZLEVEL=9 dem/​GMTED2010.tif dem/​aspect.tif ​
 + 
 +# calculate ​ Terrain Ruggedness Index TRI  ​
 +gdaldem TRI -co COMPRESS=DEFLATE -co ZLEVEL=9 dem/​GMTED2010.tif dem/​tri.tif ​
 + 
 +# calculate ​ Topographic Position Index TPI
 +gdaldem TPI -co COMPRESS=DEFLATE -co ZLEVEL=9 dem/​GMTED2010.tif dem/​tpi.tif ​
 + 
 +# calculate ​ Roughness ​  
 +gdaldem roughness -co COMPRESS=DEFLATE -co ZLEVEL=9 dem/​GMTED2010.tif dem/​roughness.tif ​
 +</​code>​
 +
 +
 +
 +
 +===== OGR Commands =====
 +
 +=== Select a polygon and create a new shape file === 
 +
 +<code bash>
 +ogrinfo -al -geom=NO ​  ​shp/​TM_WORLD_BORDERS.shp
 +# base on an attribute
 +rm -f shp/​TM_UGANDA_BORDERS-0.3.*
 +ogr2ogr ​ -overwrite ​ -f "ESRI Shapefile" ​ -where "NAME = '​Uganda'"​ shp/​TM_UGANDA_BORDERS-0.3.shp shp/​TM_WORLD_BORDERS.shp
 +
 +# base on dimension of the polygons
 +rm -f shp/​TM_LARGE_BORDERS.*
 +ogr2ogr ​ -overwrite ​ -f "ESRI Shapefile"​ -sql "​SELECT * FROM TM_WORLD_BORDERS WHERE OGR_GEOM_AREA > 100 " shp/​TM_LARGE_BORDERS.shp shp/​TM_WORLD_BORDERS.shp
 +</​code>​
 +
 +=== Add attribute and calculate a new value === 
 +
 +<code bash>
 +rm -f shp/​TM_WORLD_BORDERS_Area100.shp
 +ogr2ogr ​ shp/​TM_WORLD_BORDERS_Area100.shp shp/​TM_WORLD_BORDERS.shp
 +
 +ogrinfo -sql "ALTER TABLE TM_WORLD_BORDERS_Area100 ADD COLUMN Area100 real(12,​0)"​ shp/​TM_WORLD_BORDERS_Area100.shp
 +
 +ogrinfo ​ -dialect SQLite -sql "​UPDATE TM_WORLD_BORDERS_Area100 set Area100 = AREA + 10 " shp/​TM_WORLD_BORDERS_Area100.shp ​
 +</​code>​
 +
 +more example at http://​www.gdal.org/​ogr_sql.html ​
 +
 +
 +===== PKTOOLS Commands =====
 +
 +=== Create a mask === 
 +
 +The same operation that we have done before ​
 +<code bash>
 +gdal_translate -of  XYZ   ​vegetation/​ETmean08-11_crop.tif vegetation/​ETmean08-11_crop.txt
 +awk '{if ($3>2 && $3<4) {print $1,$2,50 } else {print}}'​ vegetation/​ETmean08-11_crop.txt > vegetation/​ETmean08-11_crop_trh.txt
 +gdal_translate -ot Byte vegetation/​ETmean08-11_crop_trh.txt vegetation/​ETmean08-11_crop_trh.tif ​
 +</​code>​
 +
 +can be done more efficient using pkgetmask. We can create 3 masks using different thresholds values. ​
 +
 +<code bash>
 +pkgetmask ​ -co COMPRESS=DEFLATE -co ZLEVEL=9 -min  1  -max  2 -data 1 -nodata 0 -ot Byte  -i  vegetation/​ETmean08-11.tif ​ -o vegetation/​ETmean08-11_01_trhA.tif
 +pkgetmask ​ -co COMPRESS=DEFLATE -co ZLEVEL=9 -min  5  -max  8 -data 1 -nodata 0 -ot Byte  -i  vegetation/​ETmean08-11.tif ​ -o vegetation/​ETmean08-11_01_trhB.tif
 +pkgetmask ​ -co COMPRESS=DEFLATE -co ZLEVEL=9 -min  0  -max  10 -data 0 -nodata 1 -ot Byte  -i  vegetation/​ETmean08-11.tif ​ -o vegetation/​ETmean08-11_01_trhC.tif
 +</​code>​
 +
 +=== Set a mask to other image ===
 +
 +Use the prepared mask and apply to the image.
 +<code bash>
 +pksetmask -co COMPRESS=DEFLATE -co ZLEVEL=9 \
 +-m vegetation/​ETmean08-11_01_trhA.tif ​ -msknodata 1 -nodata ​ -9 \
 +-m vegetation/​ETmean08-11_01_trhB.tif ​ -msknodata 1 -nodata ​ -5 \
 +-m vegetation/​ETmean08-11_01_trhC.tif ​ -msknodata 1 -nodata -10 \
 +-i vegetation/​ETmean08-11.tif -o vegetation/​ETmean08-11_01_msk.tif
 +</​code>​
 +
 +=== Composite images ===
 +
 +<code bash>
 +# create a mask to apply during the composite
 +pkgetmask ​ -co COMPRESS=DEFLATE -co ZLEVEL=9 -min 0 -max 25 -data 0 -nodata 1 -ot Byte -i LST/​LST_MOYDmax_month1.tif -o LST/​LST_MOYDmax_month1-msk.tif
 +</​code>​
 +
 +<code bash>
 +# Calculate mean and standard deviation with several images
 +pkcomposite $(for file in LST/​LST_MOYDmax_month??​.tif LST/​LST_MOYDmax_month?​.tif ​ ; do echo -i $file ;  done ) -m LST/​LST_MOYDmax_month1-msk.tif -msknodata 0 -cr mean   ​-dstnodata 0 -co  COMPRESS=LZW -co ZLEVEL=9 -o LST/​LST_MOYDmax_monthMean.tif
 +pkcomposite $(for file in LST/​LST_MOYDmax_month??​.tif LST/​LST_MOYDmax_month?​.tif ​ ; do echo -i $file ;  done ) -m LST/​LST_MOYDmax_month1-msk.tif -msknodata 0 -cr stdev   ​-dstnodata -1 -co  COMPRESS=LZW -co ZLEVEL=9 -o LST/​LST_MOYDmax_monthStdev.tif
 +</​code>​
 +Pkcomposite is very flexible, it can deal with different extents, resolutions,​ and can be used in combination with gdal comand to make composite of modis immages (see [[ wiki:​modiswget | Download MODIS Land Surface Temperature (*.hdf) from Nasa ftp ]])
 +
 +An alternative way is to use pkstatprofile
 +
 +<code bash>
 +# Create a multiband vrt 
 +gdalbuildvrt -overwrite -separate LST/​LST_MOYDmax_month.vrt LST/​LST_MOYDmax_month?​.tif LST/​LST_MOYDmax_month??​.tif ​
 +# Calculate mean and standard deviation
 +pkstatprofile -co  COMPRESS=LZW -nodata 0 -f mean -f stdev  -i LST/​LST_MOYDmax_month.vrt -o LST/​LST_MOYDmax_month_mean_stdev.tif ​
 +
 +</​code>​
 +
 +=== Filter/​Aggregate images ===
 +
 +<code bash>
 +# mean aggregation ​
 +pkfilter -co COMPRESS=DEFLATE -co ZLEVEL=9 -nodata 0 -ot Float32 -dx 10 -dy 10 -f mean -d 10 -i LST/​LST_MOYDmax_monthMean.tif -o 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 LST/​LST_MOYDmax_monthMean.tif -o LST/​LST_MOYDmax_monthMean_circular11mean.tif
 +</​code>​
 +
 +=== Images statistic and histogram ===
 +
 +<code bash>
 +pkstat -hist  -nbin  100 -src_min 0  -i   ​vegetation/​GPPstdev08-11.tif
 +pkstat -hist  -src_min 0             ​-i ​  ​temperature/​ug_bio_3.tif
 +</​code>​
 +
 +=== Images reclassification ===
 +
 +<code bash>
 +pkstat ​ -hist -i temperature/​ug_bio_3.tif ​ | grep -v " 0" ​ | awk '{ if ($1<75) { print $1 , 0 } else { print $1 , 1 }  }' > temperature/​reclass_ug_bio_3.txt
 +pkreclass -co COMPRESS=DEFLATE -co ZLEVEL=9 -code temperature/​reclass_ug_bio_3.txt -i temperature/​ug_bio_3.tif ​ -o temperature/​reclass_ug_bio_3.tif
 +</​code>​
 +
 +=== Zonal statistic (polygon extraction) ===
 +
 +<code bash>
 +rm -f shp/​polygons_stat.*
 +pkextractogr -srcnodata -339999995214436424907732413799364296704 ​  -r mean -r stdev -r min    -i vegetation/​GPPmean08-11.tif -s  shp/​polygons.sqlite ​   -o   ​shp/​polygons_stat.sqlite
 +
 +pkextractogr -f "ESRI Shapefile"​ -srcnodata -339999995214436424907732413799364296704 ​  -r mean -r stdev -r min -i vegetation/​GPPmean08-11.tif -s  shp/​polygons.sqlite -o   ​shp/​polygons_stat.shp
 +
 +# we can also create a csv that can be manipulate later on with awk
 +rm  -f shp/​polygons_stat.csv
 +pkextractogr -f CSV -srcnodata -339999995214436424907732413799364296704 ​  -r mean -r stdev -r min    -i vegetation/​GPPmean08-11.tif -s  shp/​polygons.sqlite ​   -o   ​shp/​polygons_stat.csv
 +</​code>​
 +
 +=== Zonal statistic (point extraction) ===
 +
 +<code bash>
 +# at point location
 +rm shp/​point_stat.csv
 +pkextractogr -f CSV -srcnodata -339999995214436424907732413799364296704 -r mean -r stdev -r min    -i vegetation/​GPPmean08-11.tif -s  shp/​presence.shp -o   ​shp/​point_stat.csv
 +# at point location + 1 pixel around ​
 +rm shp/​point-buf_stat.csv
 +pkextractogr -f CSV -buf 2 -srcnodata -339999995214436424907732413799364296704 -r mean -r stdev -r min -i vegetation/​GPPmean08-11.tif -s shp/​presence.shp -o shp/​point-buf_stat.csv
 +</​code>​
 +
 +Remove all the output
 +
 +<code bash>
 +rm -f  vegetation/​GPPcv08-11.tif LST/​*_crop.tif vegetation/​ETmean08-11_crop_trh.tif vegetation/​ETmean08-11_crop_trh.txt vegetation/​ETmean08-11_crop.txt vegetation/​ETmosaic.vrt vegetation/​ETmosaic.tif ​ vegetation/​stack_UL.tif vegetation/​stack_UR.tif vegetation/​stack_LL.tif vegetation/​stack_LR.tif vegetation/​stack.vrt vegetation/​tiles.* vegetation/​ETmean08-11_crop_proximity.tif vegetation/​ETmean08-11_crop_buffer.tif ​ dem/​slope.tif dem/​aspect.tif ​ dem/tri.tif dem/tpi.tif dem/​roughness.tif vegetation/​ETmean08-11_01_trh?​.tif LST/​LST_MOYDmax_month1-msk.tif LST/​LST_MOYDmax_monthStdev.tif LST/​LST_MOYDmax_monthMean.tif ​ LST/​LST_MOYDmax_month_mean_stdev.tif LST/​LST_MOYDmax_month.vrt LST/​LST_MOYDmax_monthMean_aggreate10mean.tif LST/​LST_MOYDmax_monthMean_circular11mean.tif ​ temperature/​reclass_ug_bio_3.tif temperature/​reclass_ug_bio_3.txt shp/​polygons_stat.csv shp/​polygons_stat.* shp/​TM_WORLD_BORDERS_Area100.*
 +</​code>​
 +
 +
 +
 +