Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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.

share|improve this question
    
Is your polygon valid? –  johanvdw Oct 19 '12 at 13:58
    
please tell us the contents of the first message box ("the first says something with QgisVectorGeometry"). When I print my polygon geometry I get <qgis.core.QgsGeometry object at 0x0D1894F8>. –  andytilia Oct 22 '12 at 14:07
    
    
Thanks nickves, but I'm afraid my problem starts earlier, as I'm not even able to get the vertex points out of my polygon geometry without qgis being closing because of segmentation fault… Is this a QGIS bug? Should I write to qgis-dev? –  Dominik Froehlich Oct 24 '12 at 8:29
add comment

3 Answers

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()))
share|improve this answer
    
As thepoly is only one geometry of the layer, that won't work: AttributeError: 'QgsGeometry' object has no attribute 'vectorType'. The Layer itself is of type polygon. –  Dominik Froehlich Oct 19 '12 at 14:01
    
@DominikFroehlich, oops! Yes, indeed. I fixed the code snippet in the answer. Please try again. Sorry about that. –  andytilia Oct 19 '12 at 15:02
add comment
up vote 3 down vote accepted

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!

share|improve this answer
add comment

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)

share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.