Working properly with pyqgqis edit buffer to enable undo commands

Today I wanna share a solution for a problem we faced during the implementation of a tool for our QGIS plugin, DsgTools. The tool in question captures the signal featureAdded to get access to the newly created features for a particular QgsVectorLayer and change their attributes.

But, during our code testing we found a major problem. During the edition, when the undo button was pressed QGIS crashed big time!

So, how can you solve this problem? Surfing the web, lcoandrade and I have come across this: http://gis.stackexchange.com/questions/183423/qgis-crashes-when-doing-a-rollback-after-modifying-values-of-an-user-added-featu

By studing the link above and some other related, we have discovered that our function connected to featureAdded was messing things up, but why? Sure, we have learned some “kludges” or like some might say “McGyver Programming Stunts” to solve the problem, but what is the magic that happens when a user finishes adding a feature in QGIS? We need to understand that in order to code in an elegant way!

First of all, featureAdded() is triggered:


self.iface.actionAddFeature().trigger()

This triggers QgsMapToolAddFeature::addFeature which does the following


QgsFeatureAction *action = new QgsFeatureAction( tr( "add feature" ), *f, vlayer, QString(), -1, this );
bool res = action->addFeature( QgsAttributeMap(), showModal );

Which calls QgsFeatureAction::addFeature


mLayer->beginEditCommand( text() );
mFeatureSaved = mLayer->addFeature( *mFeature );

if ( mFeatureSaved )
mLayer->endEditCommand();
else
mLayer->destroyEditCommand();

So, between:


mLayer->beginEditCommand( text() );

and:


if ( mFeatureSaved )
mLayer->endEditCommand();
else
mLayer->destroyEditCommand();

QgsVectorLayer::addFeature is called


mFeatureSaved = mLayer->addFeature( *mFeature );

Which calls QgsVectorLayerEditBuffer::addFeature. It is only inside this method that the undo stack is updated.


L->undoStack()->push( new QgsVectorLayerUndoCommandAddFeature( this, f ) );

But, until endEditCommand() is called it is not possible to perform undo operations (http://pyqt.sourceforge.net/Docs/PyQt4/qundostack.html#beginMacro). Ok?

So let’s go back to our former coding, keeping in mind that if you decide to perform an Undo operation QGIS will Crash big time (http://qgis-developer.osgeo.narkive.com/5wnziigA/wrapping-changeattributevalue-between-begin-and-end-editcommand#post2).

Ok, on one hand I want to manipulate feature’s attributes when they are added, but on the other hand, if I try to do so, I’ll mess up the undo stack. How should we proceed? What if I only peep while QGIS does its magic and then I barge in and do my own stuff? So I have to have one slot connected to featureAdded signal just to let me know which features should I visit on the editBuffer and another slot connected to editCommandEnded signal, so that I am sure that the undoStack is properly built allowing me to revisit the added features to change its attributes.

Something like this:

myLayer.featureAdded.connect(self.storeFeaturesIds)
mylayer.editCommandEnded.connect(self.updateAttributesAfterAdding)

@pyqtSlot(int)
def storeFeaturesIds(self, featId):
     """
     This method only stores featIds in a class variable
     """
     self.addedFeatures.append(featId)

def updateAttributesAfterAdding(self):
     layer = self.sender() #here I can get the layer that has sent the signal
     while self.addedFeatures:
          featureId = self.addedFeatures.pop()
          #begining the edit command
          layer.beginEditCommand('Your command message here'))
          #Do your stuff here
     layer.endEditCommand()

So there it is, how to access attributes on the fly and not explode your QGIS while doing it!

Advertisements

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!

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?

Have you ever heard of pgModeler?

Someone has given you a logical model with over 500 tables to implement it in PostgreSQL, where to begin? Should I cry first?  Should I just start typing create table statements like there is no tomorrow? No! I present you pgModeler!

pgmodeler

pgModeler is a modelling open software developed by Raphael Araújo e Silva and it allows you to build a database just using a nice graphical interface and after you finish your design, you can just deploy it and there you go! Your database is ready!

For instance, suppose you have a table called road with an text attribute “name”, geometry column “geom” with type MultiLinestring, with epsg 4326. To implement this table you just have to click around and there you go, you have the following table:

new_database

To deploy this, just go to export and choose between .png, sql or direct deploy into postgres. The sql for the table above is:

CREATE TABLE public.road(
id serial NOT NULL,
name text,
geom geometry(MULTILINESTRING, 4326) NOT NULL,
CONSTRAINT road_pk PRIMARY KEY (id)

);

Another great feature of pgModeler is reverse engineering! Just define a database connection and it shows you all tables and relationships.

By the way, I was the guy that was given the huge database to implement and if wasn’t by pgModeler, I think I’d still be crying… =]

Deaggregate geometries with pyqgis

Have you ever needed to explode multi geometry layer into a single geometry layer, using in each new geometry the attributes of the original multi one?  If you were working with FME, for instance, you basically would just use the transformer Deaggregator. Let’s learn how to solve this problem with python and QGIS!

The following code snippet teaches you how to work with QgsVectorLayers, it’s attributes and how to manipulate geometries.

from qgis.core import QgsVectorLayer, QgsFeature, QgsMapLayerRegistry

#fill in your input layer name. In this example, our inputLyrName is input_layer
inputLyrName = 'input_layer'
inputLyr = QgsMapLayerRegistry.instance().mapLayersByName(inputLyrName)[0]

#fill in your output layer name. In this example, our outputLyrName is output_layer
outputName = 'output_layer'
outputLyr = QgsMapLayerRegistry.instance().mapLayersByName(outputLyrName)[0]

#tests type of output: if it is a multi parted geometry or
#a single parted geometry
if outputLyr.wkbType() in [QGis.WKBPoint, QGis.WKBLineString, QGis.WKBPolygon]:
     isMulti = False
else:
     isMulti = True

outputLyr.startEditing()
addList = []
for feat in inputLyr.getFeatures():
     #gets all parts of geometry as an individual single geometry
     parts = feat.geometry().asGeometryCollection()
     #checks if it isMulti, if it is, convert each
     #part in geometryCollection to multi
     if isMulti:
          for part in parts:
               part.convertToMultiType()

     #for each part, get original set of attribute and create a new feat
     #with this set
     for i in range(0,len(parts)):
          #new feature constructor. newFeat has all atributes of feat
          newFeat = QgsFeature(feat)
          #set geometry with part
          newFeat.setGeometry(parts[i])
          #get field id and get defaultValue from provider
          idx = newFeat.fieldNameIndex('id')
          newFeat.setAttribute(idx,provider.defaultValue(idx))
          addList.append(newFeat)
outputLyr.addFeatures(addList,True)
outputLyr.commitChanges()

Hope you guys like it!

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?