How to extract area programmatically?

Question

is there a possibility to extract the area covered by one geometry of a shapefile (polygon)?

I tried with

thepoly.area()

but get "nan" as result.

---

Ok. I think you are right and there is something wrong with thepoly, so I'll need to post some code. I have a buildings vector (polygon) layer containing three polygons (which I created for testing). This layer is called "build" in the code. To extract its geometries, I wrote a little function:

def extract_poly_geometry(self, layer): 
    feature = QgsFeature()
    geom_list = []
    while layer.dataProvider().nextFeature(feature):
        geom_list.append(feature.geometry())
    return geom_list

In the main part of my program it is called by

theINPolyShapeField = self.extract_poly_geometry(build)

Later I iterate over the polygons and try to run some code for each of them:

for a in range(0, (len(theINPolyShapeField) - 1), 1): # going through buildings
   thepoly = theINPolyShapeField[a] # extract building number a
   QMessageBox.information(None, "DEBUG:", str(thepoly))
   # polyarea = float(thepoly.area())
   QMessageBox.information(None, "DEBUG:", str(thepoly.asPolygon()))

I use the QMessageBoxes for debuging. The first says

qgis.core.QgsGeometry object at 0xbbe889c

(thanks for the reminder, andytillia), which I assume is correct. The second one however deliveres an empty list, "[]".

---

I tried to convert to MultyPoint by

QMessageBox.information(None, "DEBUG:", str(thepoly.asMultiPoint()))

to have a look at the points, but this crashes qGIS with segmentation fault…

---

I checked the type of thepoly on the console by

print(thepoly.wkbType())

which returns "3". According to Quantum GIS API Documentation this means "WKBPolygon". So thepoly is a polygon, is it?

---

Hm, i tried thepoly.validateGeometry() on the python console. When I tried it for my first geometry, it returned "[]". A call for the second one crashed QGIS with segmentation fault…

The theINPolyShapeField list should not be the problem, as the issue also occures on the console, where I use geom_list directly.

Yes, exactly, I see three polygons (that I drew before). I just created them using the "New shapefile" tool. I think I selected the right layer as I checked with build.name() and got the right "test_buildings", what is correct.

len(theINPolyShapeField) returnes 3.

Unfortunatelly thepoly.area() does not return anything as QGIS crashes with segmentation fault… Oo

Thanks for the loop. It really looks lot nicer than mine. However I think this is not the problem, as I diden't use the loop on the console and the same issue occures.

Answer 1

First, I would make sure that thepoly is a polygon. You could see what this returns:

geom = thepoly.geometry()
print geom.asPolygon()

If it's a set of vertices, then thepoly is okay and this isn't your answer.

If not, then thepoly isn't really a polygon, and you should try to find out why. If that's hard to do, then post the code here ("many eyes make light work").

---

EDIT after code posted in the question

(1) What does thepoly.validateGeometry() return?

(2) Maybe there's something going on with the theINPolyShapeField list.

2A: Since you've run a line from the python console, can we assume you see three polygons in the QGIS interface? And you're sure it's the same layer that you're pointing to with the code?

2B: What is len(theINPolyShapeField)? You would expect 3, I think.

2C: Does thepoly.area() return a value for any of the features in theINPolyShapeField?

2D: Does a more pythonic loop handle it any differently? Try this adaptation of your code above:

for thepoly in theINPolyShapeField:
   QMessageBox.information(None, "DEBUG:", str(thepoly))
   QMessageBox.information(None, "DEBUG:", str(thepoly.asPolygon()))

Answer 2

As I diden't know what to do I asked the Qgis-developer mailing list and got following response by Martin D., that solved the problem:

The problem in your code is that you are storing QgsGeometry instances which are owned by the QgsFeature instance. So when another feature is stored into your QgsFeature instance (in nextFeature() call) or when feature instance goes out of scope, the geometry instance becomes invalid - you still have a python object, but it contains a dangling pointer, hence the crashes. The solution is to store a copy of the geometry:

geom_list.append( QgsGeometry( feature.geometry() ) )

Hope this will help others too!

Answer 3

When you "return" geometry from a function it still corrupts (returns a garbage). I appended geometry to a list and "return'ed" this list - no geometry!!! the only way to "return" the geometry by a function is to convert it to a list; for example

"return aFeat.geometry().asPolyline()"

(qgis 1.8, python 2.7, windows 7)