Dane,
Couple problems with your data ...
ogrinfo says type='Polygon' and count=246, but I only actually
extracted 245 features from it. Most are actually MultiPolygons, not
Polygons.
from osgeo import ogr
ds = ogr.Open('/Users/seang/data/springmeyer/tm_world_buffer_zero.shp')
world = ds.GetLayer(0)
borders = []
from shapely.wkb import loads
while 1:
f = world.GetNextFeature()
if f is None: break
g = f.geometry()
borders.append(loads(g.ExportToWkb()))
assert len(borders) == 245
a = [g for g in borders if g.geom_type == 'Polygon']
assert len(a) == 100
b = reduce(lambda x,y: x+y, [list(g.geoms) for g in borders if
g.geom_type == 'MultiPolygon'])
assert len(b) == 3238
Cascaded union on the original polygon members of the shapefile works:
from shapely.ops import cascaded_union
ua = cascaded_union(a)
assert ua.geom_type == 'MultiPolygon'
ua.area == 2254.6654667682537
Union of the exploded multis works:
ub = cascaded_union(b)
assert ub.geom_type == 'MultiPolygon'
assert ub.area == 19073.829960315219
Cascaded union of (a+b) failed. Not sure what's up with that.
Sean
On Jan 28, 2009, at 12:57 PM, Dane Springmeyer wrote:
> Hi Sean,
>
> Thanks for those details. The shapefile(s) surely do have invalid
> polygons, at the least.
>
> Focusing on Bjorn's I tried to running the script again after calling
> buffer(0) on the invalid polys to clean them up, but I still hit the
> error.
>
> So, I've posted that shapefile for further testing here:
> http://dbsgeo.com/tmp/
>
> Thanks!
>
> Dane
>
>
> On Jan 28, 2009, at 8:58 AM, Sean Gillies wrote:
>
>> I suspect that Bjorn's shapefile has some features with null
>> geometries. Schuyler's original one did and I removed 5 or 6 of them
>> when packaging it to go with PrimaGIS.
>>
>> Other caveats about the not-ready-for-prime-time cascaded_union:
>>
>> GEOS fails on invalid geometries such as "bowtie" polygons. These are
>> easy to create in a GIS, but rather expensive to check, too expensive
>> to test .is_valid for every item in a batch. Every geometry passed in
>> to cascaded_union should be valid. Expect an exception if even a
>> single one is not.
>>
>> GEOSUnionCascaded might only work on multipolygons, not on other
>> collections. Your watershed data isn't a collection of multipolygons,
>> is it? If so, you might need to explode them and pass a sequence of
>> polgyons to cascaded_union.
>>
>> Thanks for trying it out, Dane. I'm happy to help you use it.
>>
>> Cheers,
>> Sean
>>
>> On Jan 27, 2009, at 10:36 PM, Dane Springmeyer wrote:
>>
>>> Hey Sean,
>>>
>>> I just rebuilt geos from Paul's 3.1.0rc1.tar.bz2 and was excited to
>>> play around with the cascaded union with a big watershed dataset
>>> that's always been a bear.
>>>
>>> I am on mac 10.5.6 and used your script from:
>>> http://sgillies.net/blog/870/a-more-perfect-union-continued
>>>
>>> It seems to hit the last feature and then crashes python:
>>>
>>> [....]
>>> cxx type MultiPolygon
>>> Assertion failed: (!"should never be reached"), function itemsTree,
>>> file AbstractSTRtree.cpp, line 358.
>>> Abort trap
>>>
>>> Then I tried on Bjorn's border shapefile (v3) (TM_WORLD_BORDERS-0.3)
>>> and also on the last feature it fails, but with this exception:
>>>
>>> Traceback (most recent call last):
>>> File "<stdin>", line 1, in <module>
>>> File "/Library/Python/2.5/site-packages/Shapely-1.1.0-py2.5.egg/
>>> shapely/ops.py", line 48, in cascaded_union
>>> return geom_factory(lgeos.GEOSUnionCascaded(collection))
>>> File "/Library/Python/2.5/site-packages/Shapely-1.1.0-py2.5.egg/
>>> shapely/geometry/base.py", line 35, in geom_factory
>>> raise ValueError, "No Shapely geometry can be created from this
>>> null value"
>>> ValueError: No Shapely geometry can be created from this null value
>>>
>>> Seems odd huh? Let me know if you have any ideas for further
>>> testing.
>>>
>>> Thanks!
>>>
>>> Dane
>>>
>>>
>>> On Jan 26, 2009, at 11:11 AM, Sean Gillies wrote:
>>>
>>>> So far so good.
>>>>
>>>> http://sgillies.net/blog/869/efficient-batch-operations-for-shapely
>>>>
>>>> Anyone want to join me as a downstream user and tester of GEOS 3.1?
>>>>
>>>> --
>>>> Sean Gillies
>>>> [email protected]
>>>> http://sgillies.net
>>>>
>>>> _______________________________________________
>>>> Community mailing list
>>>> [email protected]
>>>> http://lists.gispython.org/mailman/listinfo/community
>>>
>>> _______________________________________________
>>> Community mailing list
>>> [email protected]
>>> http://lists.gispython.org/mailman/listinfo/community
>>
>> --
>> Sean Gillies
>> [email protected]
>> http://sgillies.net
>>
>> _______________________________________________
>> Community mailing list
>> [email protected]
>> http://lists.gispython.org/mailman/listinfo/community
>
> _______________________________________________
> Community mailing list
> [email protected]
> http://lists.gispython.org/mailman/listinfo/community
_______________________________________________
Community mailing list
[email protected]
http://lists.gispython.org/mailman/listinfo/community
