Getting raster histogram in QGIS using SEXTANTE and R

S.S. Rebelious

9 years ago

[This article was first published on Misanthrope's Thoughts, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

The issue with the broken histogram creation tool in QGIS annoyed me far too long. Sometimes you just need a quick glance on the histogram of a raster just to make a decision on how to process it or just to assess distribution of classes. But as you know, QGIS histogram tool creates broken plots: they are just useless (don’t know if something have changed in Master).

Fortunately there is the SEXTANTE – extremely useful plugin that provides access to hundreds of geoprocessing algorithms and allows to create your own very quickly. SEXTANTE grants access to R functions, and one of the out-of-the-box example scripts, surprise-surprise, is about raster histogram creation. “Hooray! The problem is solved!” – one may cry out. Unfortunately it’s not. “Raster histogram” tool is somewhat suitable only for single-band rasters (its algorithm is as simple as hist(as.matrix(raster))). Obviously we need to create our own tool. Luckily it is extremely simple if one is already familiar with R.

To create your own R tool for QGIS go to Analisys menu and activate SEXTANTE toolbox. Then go to R-scripts -> Tools and double-click ‘Create new R script’ and just copy and paste this code (an explanatoin will follow the description):

##layer=raster
##Tools=group
##dens_or_hist=string
##showplots
library('rpanel')
library('rasterVis')
str <- dens_or_hist
if (str !='dens' & str != 'hist'){
rp.messagebox('you must enter "dens" or "hist"', title = 'oops!')
} else {
    im <- stack(layer)
    if (str == 'dens') {
    densityplot(im)
    } else if (str == 'hist') {
    histogram(im)
  }
}

Now save the script. It will appear among other R scripts. Launch it and you will see this window:

Script window

In ‘layer‘ drop-down menu choose needed raster. In ‘dens or hist‘ field type dens or hist to identify which representation of the raster cell data you prefer: density plot or a histogram. [ It would be great to know if it is possible to embed these options in drop-down menu as well (help me if you know how). ] In ‘R plots’ field set the destination of the graph to be saved (optional).

If you chose density plot you will get something like this:

Density plot

If you choose histogram, you will get something like this:

Histogram

And if you fail to provide correct value in ‘dens or hist‘ field you will get this:

Fail window

How it works.

##layer=raster provides us with the drop-down menu to select needed raster and it will be addresed as layer in our script.

##Tools=group will assign our script to the Tools group of R scripts in SEXTANTE.

##dens_or_hist=string provides us with the string input field dens or hist.

##showplots will show the plot in Results window after execution.

library(‘rpanel’) this library is needed to show pop-up window with the error when the user fails to type 4 letters correctly.

library(‘rasterVis’) this library works on top of raster package and it is responsible for plotting our graphs.

Then we define condition to pop up error window via rp.messagebox if the plot type is misspelled or missing.

Finally, we transform our raster to the RasterStack object to be able to run either densityplot() or histogram() function.

Final notes

In a blink of an eye we were able to create quite handy tool to replace a long broken one. But there is the downside. This script won’t be able to handle huge rasters – I wasn’t able to process my IKONOS imagery where sides of the rasters exceeds 10000 pixels. Also I would not recommend to use in on a hyperspectral imagery with the dozens of bands because the process will run extremely sloooow (it will take hours and tons of RAM) and the result will be obviously unreadable.

To leave a comment for the author, please follow the link and comment on their blog: Misanthrope's Thoughts.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.