User Tools

Site Tools


wiki:grass:grasstopovar

A suite of global, cross-scale topographic variables for environmental and biodiversity modeling

Scripting procedure to recreate derived and aggregated topographic variables using Global Multi-resolution Terrain Elevation Data 2010 (GMTED2010) Script written by Giuseppe Amatulli

This example contains the following steps:

  • Download the GMTED2010 digital elevation model (DEM)
  • Calculate the topographic variables using gdal commands
  • Calculate the topographic variables using grass7 commands
  • Aggregate the topographic variables to different resolution using pktools
  • Create a plot in R
These scripts have been tested in Ubuntu Linux environment. Microsoft Windows users: consider installing the OSGEO-Live Virtual machine that has GRASS GIS 7 already installed and follow the installation instruction to install pktools .

Use gdaldem to create topographic variables

 # download the script 
 wget https://raw.githubusercontent.com/selvaje/spatial-ecology-codes/master/wiki/grass/sc1_topography_varGDAL.sh
 bash sc1_topography_varGDAL.sh
wiki/grass/sc1_topography_varGDAL.sh
# get information of the files in the GMTED repository new 
 
gdalinfo  /vsizip/vsicurl/http://edcintl.cr.usgs.gov/downloads/sciweb1/shared/topo/downloads/GMTED/Grid_ZipFiles/md75_grd.zip
 
echo create a VRT with the extent of the study region
 
gdalbuildvrt -overwrite  -te  5.83333333333333  44.833333333333333333333333  10.83333333333333333333333   48.166666666666666 elevation.vrt  /vsizip/vsicurl/http://edcintl.cr.usgs.gov/downloads/sciweb1/shared/topo/downloads/GMTED/Grid_ZipFiles/md75_grd.zip/md75_grd/w00*x.adf
 
gdal_translate -co COMPRESS=DEFLATE -co ZLEVEL=9   elevation.vrt  elevation.tif
 
echo calculate  slope 
gdaldem slope  -s 111120 -co COMPRESS=DEFLATE -co ZLEVEL=9  elevation.tif    slope.tif 
 
echo calculate  apect
gdaldem aspect  -zero_for_flat -co COMPRESS=DEFLATE -co ZLEVEL=9 elevation.tif   aspect.tif 
 
echo calculate  cosine and sine of aspect
gdal_calc.py  --NoDataValue=-9999 --co=COMPRESS=DEFLATE --co=ZLEVEL=9 --overwrite --type=Float32  -A aspect.tif  --calc="(sin(A.astype(float)  * 3.141592 / 180 ))" --outfile  aspectsine.tif 
gdal_calc.py  --NoDataValue=-9999 --co=COMPRESS=DEFLATE --co=ZLEVEL=9 --overwrite --type=Float32  -A aspect.tif  --calc="(cos(A.astype(float)* 3.141592 / 180 ))" --outfile  aspectcosine.tif 
 
echo calculate Eastness  Northness 
gdal_calc.py --NoDataValue=-9999 --co=COMPRESS=DEFLATE --co=ZLEVEL=9 --overwrite --type=Float32   -A slope.tif -B  aspectsine.tif   --calc="((sin(A.astype(float) * 3.141592 / 180 ))*B)" --outfile   eastness.tif 
gdal_calc.py --NoDataValue=-9999 --co=COMPRESS=DEFLATE --co=ZLEVEL=9 --overwrite --type=Float32   -A slope.tif -B aspectcosine.tif  --calc="((sin(A.astype(float) * 3.141592 / 180 ))*B)" --outfile   northness.tif
 
echo  calculate  Terrain Ruggedness Index TRI  
gdaldem TRI -co COMPRESS=DEFLATE -co ZLEVEL=9   elevation.tif   tri.tif 
 
echo  calculate  Topographic Position Index TPI
gdaldem TPI  -co COMPRESS=DEFLATE -co ZLEVEL=9  elevation.tif   tpi.tif 
 
echo  calculate  Roughness   
gdaldem  roughness   -co COMPRESS=DEFLATE -co ZLEVEL=9 elevation.tif  roughness.tif 

Use GRASS to create others topographic variables

 wget https://raw.githubusercontent.com/selvaje/spatial-ecology-codes/master/wiki/grass/sc2_topography_varGRASS.sh

Create the GRASS GIS data base and enter GRASS:

 grass70 -text -c elevation.tif  grass_location
 # run the script
 bash sc2_topography_varGRASS.sh
 
wiki/grass/sc2_topography_varGRASS.sh
 
# Read data into GRASS 
 
r.in.gdal input=elevation.tif    output=elevation
 
g.extension  extension=r.vector.ruggedness
 
echo calculate  Vector Ruggedness Measure
r.vector.ruggedness      elevation=elevation   output=vrm   --o 
 
# export the results 
r.out.gdal -c  createopt="COMPRESS=DEFLATE,ZLEVEL=9" format=GTiff  type=Float32 input=vrm output=vrm.tif 
 
echo calculate curvature variables 
r.slope.aspect elevation=elevation  precision=FCELL  pcurv=pcurv tcurv=tcurv dx=dx  dy=dy  dxx=dxx dyy=dyy  --o
 
# export the results 
 
r.out.gdal -c  createopt="COMPRESS=DEFLATE,ZLEVEL=9" format=GTiff  type=Float32 input=pcurv output=pcurv.tif 
r.out.gdal -c  createopt="COMPRESS=DEFLATE,ZLEVEL=9" format=GTiff  type=Float32 input=tcurv output=tcurv.tif 
 
r.out.gdal -c  createopt="COMPRESS=DEFLATE,ZLEVEL=9" format=GTiff  type=Float32 input=dx output=dx.tif 
r.out.gdal -c  createopt="COMPRESS=DEFLATE,ZLEVEL=9" format=GTiff  type=Float32 input=dxx output=dxx.tif 
 
r.out.gdal -c  createopt="COMPRESS=DEFLATE,ZLEVEL=9" format=GTiff  type=Float32 input=dy output=dy.tif 
r.out.gdal -c  createopt="COMPRESS=DEFLATE,ZLEVEL=9" format=GTiff  type=Float32 input=dyy output=dyy.tif 
 
echo calculate geomorphological landforms 
 
g.extension  extension=r.geomorphon
r.geomorphon   elevation=elevation  forms=geomorphon  search=3 skip=0 flat=1 dist=0 step=0 start=0 --overwrite
 
r.out.gdal -c  createopt="COMPRESS=DEFLATE,ZLEVEL=9" format=GTiff  type=Byte input=geomorphon output=geomorphon.tif 

Additional topographic variables can be create using:

  • r.diversity which calculate Pielou, Renyi, Shannon and Simpson indices.
  • r.texture which generate images with textural features as first order statistic and second order (so-called grey level co-occurrence matrices).

Use PKtools to aggregate categorical topographic variables

Use PKtools to aggregate continues topographic variables

 wget https://raw.githubusercontent.com/selvaje/spatial-ecology-codes/master/wiki/grass/sc3_aggregation_varPKTOOLS.sh
 bash sc3_aggregation_varPKTOOLS.sh
wiki/grass/sc3_aggregation_varPKTOOLS.sh
# aggregate variables 
 
for file in *.tif ; do  
    for km in 1 5 10 50 100 ; do 
	km=$km
	res=$(expr $km \* 4)
	for aggreg in min max mean median stdev ; do 	    
	    filename=$( basename $file .tif) 
	    pkfilter -co COMPRESS=DEFLATE -co ZLEVEL=9 -ot Float32 -of GTiff -nodata -9999 -dx $res -dy $res -f $aggreg  -d $res -i $file  -o ${filename}_${km}KM${aggreg}.tif    
	done 
    done 
done 
 
 
echo majority 
 
for km in 1 5 10 50 100 ; do 
    km=$km
    res=$(expr $km \* 4)
    for aggreg in mode countid ; do
 	if [ $aggreg = "mode" ] ; then naggreg=maj ; fi 
 	if [ $aggreg = "countid" ] ; then naggreg=count ; fi 
 	pkfilter -co COMPRESS=DEFLATE  -co ZLEVEL=9 -ot Byte -nodata 255 -dx $res -dy $res -f $aggreg  -d $res -i   geomorphon.tif   -o geom_${km}KM${naggreg}.tif
    done 
done
 
 
for km in 1 5 10 50 100 ; do 
    km=$km
    res=$(expr $km \* 4)
    for class in $(seq 1 10) ;  do  
	if [ $class  = "1"  ]    ; then class2=flat     ;   fi                                                                                                                              
	if [ $class  = "2"  ]    ; then class2=peak     ;   fi                                                                                                                              
 	if [ $class  = "3"  ]    ; then class2=ridge    ;   fi                                                                                                                              
 	if [ $class  = "4"  ]    ; then class2=shoulder ;   fi                                                                                                                              
 	if [ $class  = "5"  ]    ; then class2=spur     ;   fi                                                                                                                              
 	if [ $class  = "6"  ]    ; then class2=slope    ;   fi                                                                                                                              
 	if [ $class  = "7"  ]    ; then class2=hollow   ;   fi                                                                                                                              
 	if [ $class  = "8"  ]    ; then class2=footslope;   fi                                                                                                                              
 	if [ $class  = "9"  ]    ; then class2=valley   ;   fi                                                                                                                              
 	if [ $class  = "10" ]    ; then class2=pit      ;   fi 
 	pkfilter     -co  COMPRESS=DEFLATE  -co ZLEVEL=9  -ot  Float32 -nodata 255 -dx $res -dy $res -f density -class $class -d $res -i geomorphon.tif -o geom_class${class}_${km}KM.tif
 	gdal_calc.py --co=COMPRESS=DEFLATE --co=ZLEVEL=9  --NoDataValue=0 --type=UInt16 -A geom_class${class}_${km}KM.tif   --calc="( A.astype(float) * 100 )"  --outfile=geom${class2}_${km}KMperc.tif
    done 
done  
 
echo  calculate shannon index
 
for km in 1 5 10 50 100 ; do 
    km=$km
    res=$(expr $km \* 4)
    for class in $(seq 1 10) ;  do  
	pksetmask -msknodata 0 -nodata 255 -m   geom_class${class}_${km}KM.tif  -i  geom_class${class}_${km}KM.tif  -o   geom_class255_${class}_${km}KM.tif 
	gdal_calc.py  --co=COMPRESS=LZW --co=ZLEVEL=9 --NoDataValue=0 --type=Float32 --overwrite -A  geom_class${class}_${km}KM.tif -B geom_class255_${class}_${km}KM.tif   --calc="((log ( B  / 100 )) * ( A / 100 ))"  --outfile=geom_classH_${class}_${km}KM.tif
        gdal_edit.py -a_nodata -1  geom_classH_${class}_${km}KM.tif   	
        rm -f  geom_class255_${class}_${km}KM.tif 
    done     
    gdal_calc.py  --co=COMPRESS=LZW --co=ZLEVEL=9  --type=Float32 --overwrite -A geom_classH_1_${km}KM.tif -B geom_classH_2_${km}KM.tif -C geom_classH_3_${km}KM.tif -D geom_classH_4_${km}KM.tif  -E geom_classH_5_${km}KM.tif -F geom_classH_6_${km}KM.tif -G geom_classH_7_${km}KM.tif -H geom_classH_8_${km}KM.tif -I geom_classH_9_${km}KM.tif -J geom_classH_10_${km}KM.tif --calc="((A + B + C + D + E + F + G + H + I + J) * -1  )"  --outfile=geom_${km}KMsha.tif   	
 
done

Use R to plot a variable

 wget https://raw.githubusercontent.com/selvaje/spatial-ecology-codes/master/wiki/grass/sc4_plot_varR.R   

start R

 R
 source ("sc4_plot_varR.R")
wiki/grass/sc4_plot_varR.R
library(raster)
library(rasterVis)
 
raster=raster("vrm.tif")
 
png("vrm.png",   res=100 , width = 800, height = 800  )
 
max=0.004 ; min=0 ; des="Vector Ruggedness Measure"
 
n=100
at=seq(min,max,length=n)
colR=colorRampPalette(c("blue","green","yellow", "orange" , "red", "brown", "black" ))
 
cols=colR(n)
res=1e7 
par(cex.axis=2, cex.lab=2, cex.main=2, cex.sub=2)
raster[raster>max] <- max
raster[raster<min] <- min
print ( levelplot(raster,col.regions=colR(n),  scales=list(cex=1.5) ,   cuts=99,at=at,colorkey=list(space="right",adj=2 , labels=list( cex=1.5)), panel=panel.levelplot.raster,   margin=T,maxpixels=res,ylab="", xlab=list(paste(des,sep="") , cex=1.5 , space="left" ) ,useRaster=T )  )
 
dev.off()
 
 
raster=raster("geom_1KMsha.tif")
 
png("geom_1KMsha.png",   res=100 , width = 800, height = 800  )
 
min=0.234 ; max=2.220 ; des="Shanon index of geomorphic landform"
n=100
at=seq(min,max,length=n)
colR=colorRampPalette(c("blue","green","yellow", "orange" , "red", "brown", "black" ))
 
cols=colR(n)
res=1e7 
par(cex.axis=2, cex.lab=2, cex.main=2, cex.sub=2)
raster[raster>max] <- max
raster[raster<min] <- min
print ( levelplot(raster,col.regions=colR(n),  scales=list(cex=1.5) ,   cuts=99,at=at,colorkey=list(space="right",adj=2 , labels=list( cex=1.5)), panel=panel.levelplot.raster,   margin=T,maxpixels=res,ylab="", xlab=list(paste(des,sep="") , cex=1.5 , space="left" ) ,useRaster=T ) )
 
dev.off()

wiki/grass/grasstopovar.txt · Last modified: 2018/12/18 21:57 (external edit)