Slope, curvature, aspect with RQIS 3

Question

I'm trying to calculate morfometrics index provided by SAGA using RQGIS 3 with QGIS 3.4 version.

I`ve followed tutorial examples, but i dont have been successful. The script run, but don't produce any results. Any suggestion?

Here my script:

find_algorithms(search_term = "slope",name_only = TRUE)

get_usage("saga:slopeaspectcurvature")

params2 = get_args_man(alg = "saga:slopeaspectcurvature")

wd = "E:/artigo_declividade/mdesrtm/saga" #workdirectory

params2$ELEVATION = RQGIS3::dem
params2$METHOD = 6
params2$UNIT_SLOPE = 1
params2$UNIT_ASPECT = 1
params2$SLOPE = file.path(wd, "slope.tif")
params2$ASPECT = file.path( wd, "aspect.tif")
params2$C_GENE = file.path( wd, "generalCurvature.tif")
params2$C_PLAN = file.path( wd, "planCurvature.tif")
params2$C_PROF = file.path( wd, "profCurvature.tif")
params2$C_TANG = file.path( wd, "tangencialCurvature.tif")
params2$C_LONG = file.path( wd, "longitudinalCurvature.tif")
params2$C_CROS = file.path( wd, "crossSectionalCurvature.tif")
params2$C_MINI = file.path( wd, "minimalCurvature.tif")
params2$C_MAXI = file.path( wd, "maximalCurvature.tif")
params2$C_TOTA = file.path( wd, "totalCurvature.tif")
params2$C_ROTO = file.path( wd, "flowLineCurvature.tif")


run_qgis(alg = "saga:slopeaspectcurvature",
         params = params2)

Results at Console:

$`ASPECT`
[1] "E:/artigo_declividade/mdesrtm/saga/aspect.tif"

$C_CROS
[1] "E:/artigo_declividade/mdesrtm/saga/crossSectionalCurvature.tif"

$C_GENE
[1] "E:/artigo_declividade/mdesrtm/saga/generalCurvature.tif"

$C_LONG
[1] "E:/artigo_declividade/mdesrtm/saga/longitudinalCurvature.tif"

$C_MAXI
[1] "E:/artigo_declividade/mdesrtm/saga/maximalCurvature.tif"

.
.
.

EDITED:

Im trying this. Its works until "Terrain Convexity". After that show the same behavior (morphometric functions run but with empty results).

devtools::install_github("jannes-m/RQGIS3")
require(devtools)
require(raster)
require(rgdal)
require(RQGIS3)


qgis_env <-set_env(root ='C:/Program Files/QGIS 3.4')
set_env()
open_app()


dem = RQGIS3::dem #file
wd = "E:/a" #output

morfometria = function(dem, wd){ #function that calculate many morfometric index

#aspect

  params2 = get_args_man(alg = "saga:slopeaspectcurvature")
  params2$ELEVATION = dem
  params2$METHOD = 6
  params2$UNIT_SLOPE = 1
  params2$UNIT_ASPECT = 1
  params2$SLOPE = file.path(wd, "slope.sdat")
  params2$ASPECT = file.path( wd, "aspect.sdat")
  params2$C_GENE = file.path( wd, "generalCurvature.sdat")
  params2$C_PLAN = file.path( wd, "planCurvature.sdat")
  params2$C_PROF = file.path( wd, "profCurvature.sdat")
  params2$C_TANG = file.path( wd, "tangencialCurvature.sdat")
  params2$C_LONG = file.path( wd, "longitudinalCurvature.sdat")
  params2$C_CROS = file.path( wd, "crossSectionalCurvature.sdat")
  params2$C_MINI = file.path( wd, "minimalCurvature.sdat")
  params2$C_MAXI = file.path( wd, "maximalCurvature.sdat")
  params2$C_TOTA = file.path( wd, "totalCurvature.sdat")
  params2$C_ROTO = file.path( wd, "flowLineCurvature.sdat")

  run_qgis(alg = "saga:slopeaspectcurvature",
           params = params2)


  # Valley and Ridge Detection

  paramValley = get_args_man(alg = "saga:valleyandridgedetectiontophatapproach")

  paramValley$DEM = dem
  paramValley$VALLEY = file.path( wd, "valleyDepth.sdat")
  paramValley$HILL = file.path( wd, "hillHeight.sdat")
  paramValley$VALLEY_IDX = file.path( wd, "valleyIDX.sdat")
  paramValley$HILL_IDX = file.path( wd, "hillIDX.sdat")
  paramValley$SLOPE_IDX = file.path( wd, "slopeIDX.sdat")

  run_qgis(alg = "saga:valleyandridgedetectiontophatapproach",
           params = paramValley)


  #Terrain Surface Convexity

  paramTerrainSurfaceConvexity = get_args_man(alg = "saga:terrainsurfaceconvexity")
  paramTerrainSurfaceConvexity$DEM = dem
  paramTerrainSurfaceConvexity$CONVEXITY = file.path( wd, "TSConvexity.sdat")
  paramTerrainSurfaceConvexity$METHOD = 1
  paramTerrainSurfaceConvexity$DW_WEIGHTING = 3


  run_qgis(alg = "saga:terrainsurfaceconvexity",
           params = paramTerrainSurfaceConvexity)  


  #Terrain Surface Texture

  paramTerrainSurfaceTexture = get_args_man(alg = "saga:terrainsurfacetexture")
  paramTerrainSurfaceTexture$DEM = dem
  paramTerrainSurfaceTexture$TEXTURE = file.path( wd, "TSTexture")
  paramTerrainSurfaceTexture$METHOD = 1
  paramTerrainSurfaceTexture$DW_WEIGHTING = 3


  run_qgis(alg = "saga:terrainsurfacetexture",
           params = paramTerrainSurfaceTexture)  


  #Terrain Surface Classification

  paramTSC = get_args_man(alg = "saga:terrainsurfaceclassificationiwahashiandpike")
  paramTSC$DEM = dem
  paramTSC$LANDFORMS = file.path(wd, "Landforms1.sdat")

  run_qgis(alg = "saga:terrainsurfaceclassificationiwahashiandpike",
           params = paramTSC)  

  #TPI Based Landform Classification 

  paramTPILand = get_args_man(alg = "saga:tpibasedlandformclassification")
  paramTPILand$DEM = dem
  paramTPILand$LANDFORMS = file.path( wd, "Landforms2.sdat")

  run_qgis(alg = "saga:tpibasedlandformclassification",
           params = paramTPILand) 

  # Topographic Position Index (TPI)

  paramTPI = get_args_man(alg = "saga:topographicpositionindextpi")
  paramTPI$DEM = dem
  paramTPI$TPI = file.path(wd, "TPI.sdat")

  run_qgis(alg = "saga:topographicpositionindextpi",
           params = paramTPI) 


  # Vector Ruggedness Index (VRI) 

  paramVector = get_args_man(alg = "saga:vectorruggednessmeasurevrm")
  paramVector$DEM = dem
  paramVector$VRM = file.path(wd, "vectorRuggednessIndex.sdat")

  run_qgis(alg = "saga:vectorruggednessmeasurevrm",
           params = paramVector)

  # Terrain Ruggedness Index (TRI)

  paramReggedness = get_args_man(alg = "saga:terrainruggednessindextri")

  paramReggedness$DEM = dem

  paramReggedness$TRI = file.path(wd, "terrainRuggednessIndex.sdat")

  run_qgis(alg = "saga:terrainruggednessindextri",
           params = paramReggedness)

  # Mass Balance Index

  paramMass = get_args_man(alg = "saga:massbalanceindex")
  paramMass$DEM = dem
  paramMass$MBI = file.path(wd, "massBalanceIndex.sdat")

  run_qgis(alg = "saga:massbalanceindex",
           params = paramMass)

  # Downslope Distance Gradient 

  paramGradient = get_args_man(alg = "saga:downslopedistancegradient")
  paramGradient$DEM = dem
  paramGradient$GRADIENT = file.path(wd, "gradient.sdat")
  paramGradient$DIFFERENCE = file.path(wd, "difference.sdat")

  run_qgis(alg = "saga:downslopedistancegradient",
           params = paramGradient)

  # Multiresolution Index of Valley Bottom Flatness (MRVBF)

  paramMultires = get_args_man(alg = "saga:multiresolutionindexofvalleybottomflatnessmrvbf")
  paramMultires$DEM = dem
  paramMultires$MRVBF = file.path(wd, "MRVBF.sdat")
  paramMultires$MRRTF = file.path(wd, "MRRTF.sdat")

  run_qgis(alg = "saga:multiresolutionindexofvalleybottomflatnessmrvbf",
           params = paramMultires)


  # Morphometric.Protection.Index

  paramMorphometric = get_args_man(alg = "saga:morphometricprotectionindex")
  paramMorphometric$DEM = dem
  paramMorphometric$RADIUS = 1000
  paramMorphometric$PROTECTION = file.path(wd, "morphoProtectionIndex.sdat")

  run_qgis(alg = "saga:morphometricprotectionindex",
           params = paramMorphometric)

  #Real Surface Area

  paramRealSA = get_args_man(alg = "saga:realsurfacearea")
  paramRealSA$DEM = dem
  paramRealSA$AREA = file.path(wd, "realSurfaceArea.sdat")

  run_qgis(alg = "saga:realsurfacearea",
           params = paramRealSA)


  #Curvature classification

  paramCurvature = get_args_man(alg = "saga:curvatureclassification")
  paramCurvature$DEM = dem
  paramCurvature$CLASS = file.path(wd, "curvatureClassification.sdat")

  run_qgis(alg = "saga:curvatureclassification",
           params = paramCurvature)

  #Surface.Specific.Points

  paramSSP = get_args_man(alg = "saga:surfacespecificpoints")
  paramSSP$ELEVATION = dem
  paramSSP$METHOD = 1
  paramSSP$RESULT = file.path(wd, "surfaceSpecificPoints.sdat")

  run_qgis(alg = "saga:surfacespecificpoints",
           params = paramSSP)

  # convergence index

  paramConvergence = get_args_man(alg = "saga:convergenceindex")
  paramConvergence$ELEVATION = dem
  paramConvergence$NEIGHBOURS = 1
  paramConvergence$RESULT = file.path(wd, "convergenceIndex.sdat") 

  run_qgis(alg = "saga:convergenceindex",
           params = paramConvergence)


  # wetness index

  paramWetness = get_args_man(alg = "saga:sagawetnessindex")
  paramWetness$DEM = dem
  paramWetness$TWI = file.path(wd, "topographicWetIndex.sdat") 
  paramWetness$AREA = file.path(wd, "AREAtopographicWetIndex.sdat") 
  paramWetness$SLOPE = file.path(wd, "SLOPEtopographicWetIndex.sdat") 
  paramWetness$AREA_MOD = file.path(wd, "AREAMODtopographicWetIndex.sdat") 

  run_qgis(alg = "saga:sagawetnessindex",
           params = paramWetness) 



}


morfometria(dem = dem, wd = wd)

Where are you installing RQGIS3 from - it doesn't appear to be on CRAN. — Spacedman, Commented Apr 2, 2019 at 21:48

Spacedman · Accepted Answer · 2019-04-02 22:07:14Z

For me, run_qgis creates a bundle of .sdat files in the working directory:

> list.files(wd,pattern="*.sdat$")
 [1] "aspect.sdat"                  "crossSectionalCurvature.sdat"
 [3] "flowLineCurvature.sdat"       "generalCurvature.sdat"       
 [5] "longitudinalCurvature.sdat"   "maximalCurvature.sdat"       
 [7] "minimalCurvature.sdat"        "planCurvature.sdat"          
 [9] "profCurvature.sdat"           "slope.sdat"                  
[11] "tangencialCurvature.sdat"     "totalCurvature.sdat"

These can be read in using the raster package's raster function:

> aspect = raster(file.path(wd,"aspect.sdat"))
> slope = raster(file.path(wd,"slope.sdat"))
> plot(aspect)
>

I don't know why you don't get tif files created or why run_qgis returns nothing and only prints the things it prints... I couldn't get the load_output option to work...

Worked perfectly! I advise you to convert .sdat files to .tif because it has a faster processing time, even after loading it in R. — Arthur, Commented Apr 3, 2019 at 11:40

Spacedman · Accepted Answer · 2019-05-07 16:32:29Z

The following code works including the load_output option. Just out of curiosity, are you running RQGIS3 from within RStudio? Because in my case running RQGIS3::open_app() crashes the RStudio R session. This was also confirmed by several other users (see the github issue tracker of RQGIS3). Since RQGIS3 works when run from the CLI, I am not really sure what to do about it... BTW running RQGIS3 in a rocker-docker container, i.e. in a docker container running RStudio, also works.

# attach libraries
library("RQGIS3")
#> Loading required package: reticulate
library("raster")
#> Loading required package: sp
# report QGIS3 session info
qgis_session_info()
#> Assuming that your root path is '/usr'!
#> $gdal
#> [1] "2.3.2"
#> 
#> $grass7
#> [1] "7.6.0"
#> 
#> $qgis_version
#> [1] "3.4.6-Madeira"
#> 
#> $saga
#> [1] "2.3.1"
# retrieve and specify parameters of a specific geoalgorithm
alg = "saga:slopeaspectcurvature"
params = get_args_man(alg)
#> Choosing default values for following parameters:
#> METHOD: 0
#> UNIT_SLOPE: 0
#> UNIT_ASPECT: 0
#> See get_options('saga:slopeaspectcurvature') for all available options.
params$ELEVATION = RQGIS3::dem
params$METHOD = 6
params$UNIT_SLOPE = 1
params$UNIT_ASPECT = 1
wd = tempdir()
params$SLOPE = file.path(wd, "slope.sdat")
params$ASPECT = file.path(wd, "aspect.sdat")
params$C_GENE = file.path(wd, "generalCurvature.sdat")
params$C_PLAN = file.path(wd, "planCurvature.sdat")
params$C_PROF = file.path(wd, "profCurvature.sdat")
params$C_TANG = file.path(wd, "tangencialCurvature.sdat")
params$C_LONG = file.path(wd, "longitudinalCurvature.sdat")
params$C_CROS = file.path(wd, "crossSectionalCurvature.sdat")
params$C_MINI = file.path(wd, "minimalCurvature.sdat")
params$C_MAXI = file.path(wd, "maximalCurvature.sdat")
params$C_TOTA = file.path(wd, "totalCurvature.sdat")
params$C_ROTO = file.path(wd, "flowLineCurvature.sdat")
out = run_qgis(alg, params = params, load_output = TRUE)
#> $ASPECT
#> [1] "/tmp/RtmpyQkN4o/aspect.sdat"
#> 
#> $C_CROS
#> [1] "/tmp/RtmpyQkN4o/crossSectionalCurvature.sdat"
#> 
#> $C_GENE
#> [1] "/tmp/RtmpyQkN4o/generalCurvature.sdat"
#> 
#> $C_LONG
#> [1] "/tmp/RtmpyQkN4o/longitudinalCurvature.sdat"
#> 
#> $C_MAXI
#> [1] "/tmp/RtmpyQkN4o/maximalCurvature.sdat"
#> 
#> $C_MINI
#> [1] "/tmp/RtmpyQkN4o/minimalCurvature.sdat"
#> 
#> $C_PLAN
#> [1] "/tmp/RtmpyQkN4o/planCurvature.sdat"
#> 
#> $C_PROF
#> [1] "/tmp/RtmpyQkN4o/profCurvature.sdat"
#> 
#> $C_ROTO
#> [1] "/tmp/RtmpyQkN4o/flowLineCurvature.sdat"
#> 
#> $C_TANG
#> [1] "/tmp/RtmpyQkN4o/tangencialCurvature.sdat"
#> 
#> $C_TOTA
#> [1] "/tmp/RtmpyQkN4o/totalCurvature.sdat"
#> 
#> $SLOPE
#> [1] "/tmp/RtmpyQkN4o/slope.sdat"
plot(stack(out))

Im using R-studio. I got some erros when loading RQGIS3, but nothing that crashes it. Another problem that im having is that im executing many morphometric index, one after other. After execute 5/6 indexs, happens what i said in the topic. Results doesnt appear. I tryed in RQGIS with QGIS 2.18, but looks the same. I tested it on CML too, but got the same problem. Any ideia? — Arthur, Commented Apr 4, 2019 at 1:56
I have edited my reply to include all morphometric indices and it works just fine. Yeah, my bad, you are under Windows, where RQGIS3 doesn't crash RStudio. — Jannes, Commented Apr 4, 2019 at 8:27
Dear Jannes, thanks for the anwser. I did it like you. But when i tryed many morfometric index, like edited question, it goes wrong. Can you run my edited question to check if im doing somethin wrong? Thank for the help. — Arthur, Commented Apr 15, 2019 at 14:11

Stack Exchange Network

Slope, curvature, aspect with RQIS 3

2 Answers 2

Your Answer

Not the answer you're looking for? Browse other questions tagged
qgis
r
rqgis
or ask your own question.

Hot Network Questions

Slope, curvature, aspect with RQIS 3

2 Answers 2

Your Answer

Sign up or log in

Post as a guest

Not the answer you're looking for? Browse other questions tagged qgisrrqgis or ask your own question.

Related

Hot Network Questions

Not the answer you're looking for? Browse other questions tagged
qgis
r
rqgis
or ask your own question.