I do (it works):
from osgeo import ogr
inputds = ogr.Open(r'c:\temp.shp')
inputlyr = inputds.GetLayer()
for i,feature in enumerate(inputlyr):
print i, feature.GetGeometryRef()
Output:
0 POLYGON ((4.286549707602339 0.526315789473684,3.797114870553876 -2.563854154275787,2.376719651351817 -5.351536733451042,0.164402230527078 -7.563854154275783,-2.623280348648173 -8.984249373477848,-5.713450292397642 -9.473684210526315,-8.803620236147113 -8.984249373477859,-11.591302815322369 -7.563854154275806,-13.803620236147115 -5.351536733451073,-15.224015455349186 -2.563854154275825,-15.71345029239766 0.526315789473643,-15.22401545534921 3.616485733223115,-13.803620236147164 6.404168312398374,-11.591302815322436 8.616485733223128,-8.803620236147188 10.036880952425204,-5.713450292397716 10.526315789473685,-2.623280348648239 10.036880952425237,0.164402230527027 8.616485733223191,2.37671965135178 6.40416831239846,3.797114870553858 3.616485733223211,4.286549707602339 0.526315789473684))
1 POLYGON ((3.894736842105264 1.16374269005848,3.4053020050568 -1.926427253690991,1.984906785854742 -4.714109832866246,-0.227410634969997 -6.926427253690988,-3.015093214145249 -8.346822472893052,-6.105263157894718 -8.836257309941519,-9.195433101644188 -8.346822472893063,-11.983115680819445 -6.926427253691011,-14.195433101644191 -4.714109832866277,-15.61582832084626 -1.926427253691029,-16.105263157894736 1.163742690058438,-15.615828320846287 4.253912633807911,-14.195433101644239 7.04159521298317,-11.983115680819513 9.253912633807921,-9.195433101644262 10.674307853009999,-6.105263157894791 11.163742690058481,-3.015093214145315 10.674307853010031,-0.227410634970049 9.253912633807985,1.984906785854705 7.041595212983256,3.405302005056782 4.253912633808007,3.894736842105264 1.16374269005848))
...
but when I try to do:
print inputlyr[0].GetGeometryRef()
QGIS crash - Crash dumped. minidump written to C...
My goal is not to print the geometry.. I want to Union geometry
poly1 = inputlyr[0].GetGeometryRef()
poly2 = inputlyr[1].GetGeometryRef()
union = poly1.Union(poly2)
print union.ExportToWkt()
But I receives the same error. I test this on 2 comp (QGIS 2.6 and QGIS 2.2)
I do something wrong?
Answer
This is a GDAL/OGR Python gotcha. So when you use print inputlyr[0].GetGeometryRef()
, the reference to the indexed feature has been cleaned up before geometry functions can be used. This is expected behaviour, and there is no roadmap to change it. The issue is that the GDAL/OGR Python bindings are over a decade old, and are tightly coupled with the underlying C++. (Modern projects like fiona and shapely attempt to make this sort of geometry processing more Pythonic.)
Here's the fix. To use two geometry objects, references to the dataset and layer and features must be kept "alive" to prevent them from being dereferenced. Try this:
# Get features and geometries
feat1 = inputlyr.GetFeature(0)
poly1 = feat1.GetGeometryRef()
feat2 = inputlyr.GetFeature(1)
poly2 = feat2.GetGeometryRef()
# Process geometries
union = poly1.Union(poly2)
print union.ExportToWkt()
# Now you can dereference feat1, etc.
del feat1, poly1
del feat2, poly2
No comments:
Post a Comment