Tag Archives: raster

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?

Advertisements

Using GDAL to get raster extent

Hi there guys!!!

Let’s suppose we want to determine the extent of a raster file and we want to use GDAL and Python. How can we do that?

Let’s start importing GDAL:

from osgeo import gdal

Then we need t0 open the raster file:

gdalSrc = gdal.Open('/foo/bar/filename')

To finish getting what we need, let’s get our affine transform coefficients with the following:

upx, xres, xskew, upy, yskew, yres = gdalSrc.GetGeotransform()

Where upx is the upper pixel x, xres is the pixel width and xskew is the shear in the x direction as we can see in the following picture (https://en.wikipedia.org/wiki/Affine_transformation):

xskew

Let’s go forward. The variable upy is the upper pixel y, yres is the pixel height and yskew is the shear in the y direction as we can see in the following picture (https://en.wikipedia.org/wiki/Affine_transformation):

yskew

Understanting the variables we can use the documentation obtained at http://www.gdal.org/gdal_datamodel.html to get the following relationship:

Xgeo = GT(0) + Xpixel*GT(1) + Yline*GT(2)
Ygeo = GT(3) + Xpixel*GT(4) + Yline*GT(5)

Where GT(2) is the xskew and GT(4) is the yskew.

With all of this together, we can make the following code snippet to accomplish our mission:

from osgeo import gdal, ogr
gdalSrc = gdal.Open('/foo/bar/filename')
upx, xres, xskew, upy, yskew, yres = gdalSrc.GetGeotransform()
cols = gdalSrc.RasterXSize
rows = gdalSrc.RasterYSize

ulx = upx + 0*xres + 0*xskew
uly = upy + 0*yskew + 0*yres

llx = upx + 0*xres + rows*xskew
lly = upy + 0*yskew + rows*yres

lrx = upx + cols*xres + rows*xskew
lry = upy + cols*yskew + rows*yres

urx = upx + cols*xres + 0*xskew
ury = upy + cols*yskew + 0*yres

I hope this can be useful to you guys.