All posts by lcoandrade

Performing a HSV pan-sharpening in QGIS

Hi there guys!

Today I wanna talk about image pan-sharpening in QGIS.

I don’t want to enter in deeper technical details regarding pan-sharpening algorithms or something like that. I just want to focus on visual perception and image quality, ok?

QGIS in its processing toolbox has the ability to use the Orfeo toolbox (https://www.orfeo-toolbox.org). Orfeo is a great open-source software to perform image analysis. It is really great, believe me! But when talking about pan-sharpening it lacks the ability of perform HSV fusions (https://www.orfeo-toolbox.org/CookBook/Applications/app_Pansharpening.html).

But how can this affect our work?

Lets consider that I want to perform a pan-sharpening using this following images:

Using Orfeo toolbox I can achieve the following results: (Bayes, LMVM and RCS respectively)

Looking at the results made by Orfeo we can choose the RCS output as the best looking pan-sharpened image. In my opinion, it seems to get more details of that part of the land surrounded by water on top of the image, don’t you agree?

Now, the other 2, Bayes and LMVM outputs don’t look quite good in my humble opinion. The Bayes output delivered almost no change and the LMVM output seems a bit blurred.

Considering this, you can try to perform a pan-sharpening operation using HSV fusion, but this is not available neither in QGIS (natively) nor in Orfeo toolbox.

So, having that in mind, my co-workers and I decided to develop a script to perform pan-sharpening using HSV fusion. The result of our work is available as a tool in our QGIS plugin called DSG Tools (https://plugins.qgis.org/plugins/DsgTools/). In another another opportunity I will explain more about DSG Tools.

DSG Tools is a QGIS plugin developed by the Geographic Service at Brazilian Army. It was developed to make it possible to acquire vector data in QGIS in accordance to the Brazilian Geospatial Vector Data Structure, but besides that, DSG Tools provides a considerable toolbox that allow among other things, perform pan-sharpening using HSV fusion.

The script we made can be installed into QGIS’ toolbox and it will be available in a group called DSG as we can see next:

dsg_toolbox_group

Well, starting the script and setting its parameters (remember to superimpose the images prior the fusion – Orfeo has a tool for this) we can achieve the following result:

dsgtools_hsv_fusion

In my opinion, the pan-sharpened image generated by our HSV fusion script seems very close to the real terrain. Let’s now compare Orfeo’s RCS output and our HSV fusion output side by side:

Which one would you choose? I would choose the HSV fusion… But remember, this is just my opinion…

See you later! Bye!

Advertisements

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 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.

Detecting duplicated geometries in a PostGIS table

Hi there!

Today I want to make a small post and talk about another annoying problem: Duplicated Geometries!

How can you identify which geometries are duplicated in your table? In some cases we have duplicated geometries with different attributes, which one is the correct one?

To deal with this, a simple query can help you in this matter.

Concepts to keep:

  • PARTITION BY: Is a window function that receives a column, this means it will partition the column in groups with equal values (this is what we want for our “geom” column: duplicates!)
  • ROW_NUMBER(): gives a row number, inside a partition, counting from 1. This means that any number greater than 1 is a duplicate.

Here it is:

select * from (
SELECT id, ROW_NUMBER() OVER(PARTITION BY geom ORDER BY id asc) AS Row,
geom FROM ONLY your_schema.your_table
) dups where dups.Row > 1

Now you can check the duplicates and choose which ones you want to keep.

That’s it!

Detecting spikes in geometries using angle threshold

Hi there,

We all know that search for errors in geometries can be quite a journey. One of the errors we need to fix is the presence of spikes in our geometries. One way to determine the location of those spikes is to determine the angles and check if them are smaller than a predefined threshold.

Ok them, one way to tackle this problem is to load your data into a PostGIS layer and use the available ST functions.

In this post I’ll show you guys a SQL query to solve this.

Let’s suppose we have geometries with problems like these here:

angle_before

To solve this problems we will use the following ST PostGIS functions:

  • ST_Azimuth
  • ST_PointN
  • ST_NPoints
  • ST_Boundary
  • ST_Dump
  • ST_ForceRHR

For more information about these functions check this: http://postgis.net/docs/manual-1.3/ch06.html

The following query can show us where the problems are located, let’s use a limit of 10 degrees here to determine the spikes:

WITH result AS (SELECT points.id, points.anchor, (degrees
                            (
                                ST_Azimuth(points.anchor, points.pt1) - ST_Azimuth(points.anchor, points.pt2)
                            )::decimal + 360) % 360 as angle
            FROM
            (SELECT
                  ST_PointN(geom, generate_series(1, ST_NPoints(geom)-1)) as pt1,
                  ST_PointN(geom, generate_series(1, ST_NPoints(geom)-1) %  (ST_NPoints(geom)-1)+1) as anchor,
                  ST_PointN(geom, generate_series(2, ST_NPoints(geom)) %  (ST_NPoints(geom)-1)+1) as pt2,
                  linestrings.id as id
                FROM
                  (SELECT id as id, ST_Boundary((ST_Dump(ST_ForceRHR(geom))).geom) as geom
                   FROM only my schema.mylayer -- enter your layer name here
                   ) AS linestrings WHERE ST_NPoints(linestrings.geom) > 2 ) as points)
select distinct id, anchor, angle from result where (result.angle % 360) < 10 or result.angle > (360.0 - (10 % 360.0)) -- the 10 here is our threshold

With the results from this query we can locate our spikes using the “anchor” value returned in the query. Something like this:

angle_later

Very 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?