Category Archives: grass

Decision tree classifier in QGIS using GRASS

Hi there guys!

Maybe you’ve already needed to perform a decision tree classification in some point in your GIS life. You might have used ENVI’s decision tree to do the job, right? But, what about QGIS? Can we achieve the same results using QGIS processing? Definitely yes!!!

Today I’ll gonna show how this can be simple.

Imagine you want to classify this image:

image1

Using this one to help filter the land cover:

image2

On ENVI, this can be done with a decision tree like this one (one of my co-workers made it):

decision_tree

Kinda easy right? But, what about QGIS? Let’s take a look!

In QGIS we have the powerful processing toolbox. Inside it we have the awesome GRASS to use as much as we want. To accomplish this job we will use r.mapcalculator (it provides a GUI frontend to r.mapcalc – https://grass.osgeo.org/grass73/manuals/r.mapcalc.html).

Ok, we can see that the decision tree uses several conditions to discover what pixels it should classify in a specific class. This conditions are translated to GRASS using this:

if(x,a,b)               a if x not zero, b otherwise

Where “x” is a condition to be followed.

So, to classify water (let’s say pixel value 0) to a class with value 5 and classify everything else to 255 we need to make this expression if(x,5,255), where x is the condition (e.g. image==0).

Quite easy, right?

Moving forward, with more levels (like in the decision tree above), we need to use nested if conditions (i.e if conditions inside another if condition).

So, the next step is to open the processing toolbox in QGIS and search for r.mapcalculator to open the following dialog:

r-mapcalulator

We can choose several images (A, B, and so on…). In my case, I’ll choose Raster Layer A and Raster Layer B (the first two images in this article). Inside the formula I’ll insert my nested if conditions taken from my decision tree, as follows:

(if((A==0),6,if((A>=0 && A<=4),if((B==1),3,4),if((A>=3 && A<=6),if((B==2),1,2),if((A>=6 && A<=12),if((B==3),8,7),if((A>=12),if((B==1),10,9),5))))))

With this set, I just need to to click on “Run” and wait for my classified image. By the way, the process is quite fast thanks to GRASS…

As we can see, the results are identical (another point for QGIS and GRASS for the awesome results they deliver together!!!). Below, on the right side we have the QGIS’ result and on the left side we have ENVI’s result (I made a small subset to make it easier to check both results).

What do you guys think? Pretty nice, right?

Using GRASS with pyqgis to clean up geometries

We all know that GRASS is a great GIS software. Combined with QGIS it is even more great!

Using GRASS from within QGIS is very useful to deal with daily GIS problems. Everyone that works with geospatial data knows how annoying is to clean up geometries full of errors. The manual process demands lots of time and we can always forget something in the end. Do this kind of job automatically is faster and safer.

Let’s se how to do this using pyqgis. Imagine that we have a database layer like this:

error

A good way to clean problems like those shown above and at the same time solve snapping problems is to use the following tools in v.clean.advanced provided by GRASS:

  • break
  • rmsa
  • rmdangle

If you want a description on how those tools work, take a look at: https://grass.osgeo.org/grass73/manuals/v.clean.html

To clean it using grass we can use the following piece of code:


#choosing the algorithm
alg = 'grass7:v.clean.advanced'

#getting the vector layer we want to clean
input = iface.activeLayer()

#setting tools
tools = 'break,rmsa,rmdangle'
threshold = -1
#getting mapcanvas extent (bounding box) supposing we can see our data
e = iface.mapCanvas().extent()
xmax = e.xMaximum()
ymax = e.yMaximum()
xmin = e.xMinimum()
ymin = e.yMinimum()

extent = '{0},{1},{2},{3}'.format(xmin, xmax, ymin, ymax)

#setting parameters: choose them according to your data
snap = 100.0
minArea = 0.001

#running the grass algorithm
ret = processing.runalg(alg, input, tools, threshold, extent, snap, minArea, None, None)

#getting output layer
outputLayer = processing.getObject(ret['output'])
#Adding to registry
QgsMapLayerRegistry.instance().addMapLayer(outputLayer)

#getting error flags
errorLayer = processing.getObject(ret['error'])
#Adding to registry
QgsMapLayerRegistry.instance().addMapLayer(errorLayer) 

After running, we can see results like this:

cleaned

Quite good, right?