Facinating Sean, and awesome to see exploding geometries in python. cool
Seems like a null feature in their could be the culprit. I'll do more
digging...
Thanks!
Dane
On Jan 28, 2009, at 1:23 PM, Sean Gillies wrote:
> 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
_______________________________________________
Community mailing list
[email protected]
http://lists.gispython.org/mailman/listinfo/community
