User Tools

Site Tools


wiki:gdalpktools_kenya

Differences

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

Link to this comparison view

wiki:gdalpktools_kenya [2018/08/17 13:14] (current)
Line 1: Line 1:
 +====== Use GDAL/OGR and PKTOOLS for raster/​vector operations ======
 +
 +\\
 +Data Directory: /​home/​user/​ost4sem/​exercise/​geodata\_SDM/​ \\
 +Script: /​home/​user/​ost4sem/​exercise/​geodata_SDM/​gdal\_pktools.sh \\
 +
 +  cd /​home/​user/​ost4sem/​exercise/​geodata_SDM
 +  wget -N http://​spatial-ecology.net/​ost4sem/​exercise/​geodata_SDM/​gdal-ogr-pktools.sh
 +  rstudio /​home/​user/​ost4sem/​exercise/​geodata_SDM/​gdal-ogr-pktools.sh & 
 +
 +===== GDAL Commands =====
 +
 +=== Change the directory, and explore the files ===
 +
 +<code bash>
 +cd /​home/​user/​ost4sem/​exercise/​geodata_SDM
 +ls */​*.tif ​ # list all the files with the extension .tif
 +gdalinfo vegetation/​ETmean08-11.tif
 +
 +# Visualize the image
 +
 +qgis vegetation/​ETmean08-11.tif
 +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>​
 +
 +=== 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
 +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
 +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>​
 +
 +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>​
 +
 +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.sqlite
 +pkextractogr -srcnodata -339999995214436424907732413799364296704 ​  -r mean -r stdev -r min    -i vegetation/​GPPmean08-11.tif -s  shp/​polygons.sqlite ​   -o   ​shp/​polygons_stat.sqlite
 +
 +# 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>​
  
wiki/gdalpktools_kenya.txt ยท Last modified: 2018/08/17 13:14 (external edit)