[postgis-users] How best to do this?

Paul Ramsey pramsey at refractions.net
Thu Feb 1 21:14:22 PST 2007


GEOS 3 does include a number of fixes that will help difference and  
union.

P

On 1-Feb-07, at 9:01 PM, Stephen Woodbridge wrote:

> OK, one step forward and two steps backward ...
>
> I was about to get rid of my gemetrycollections with:
>
> update newtable set the_geom =
>  (select geomunion(geom1) from
>    (select (dump(the_geom)).geom as geom1) as foo
>     where geometrytype(geom1)='POLYGON'
>  )
>  where geometrytype(the_geom)='GEOMETRYCOLLECTION';
>
> This worked very nicely once I figured out that this was the way to  
> solve the problem. So now each geometry collection is converted  
> into a polygon or multipolygon and all the linestring and point  
> artifacts of the intersection are removed.
>
> select distinct geometrytype(the_geom) from soils2;
> POLYGON
> MULTIPOLYGON
>
> select * from soils where not isvalid(the_geom);
> reports 0 rows
>
> insert into newtable (..., the_geom)
>   select A....,
>          difference(A.the_geom, B.the_geom) as the_geom
>   from A, B
>   where A.the_geom && B.the_geom and overlaps(A.the_geom, B.the_geom);
>
> NOTICE:  TopologyException: no outgoing dirEdge found  
> (287030,892621,892621)
>
> ERROR:  GEOS difference() threw an error!
>
> Now what do I do??? If I have to upgrade it might be easiest to  
> upgrade my WinXP laptop if there is an easy way to do that, then  
> dump and load the data to that.
>
> select postgis_full_version();
> "POSTGIS="1.1.5" GEOS="2.2.3-CAPI-1.1.1" PROJ="Rel. 4.4.9, 29 Oct  
> 2004" USE_STATS (procs from 1.1.1 need upgrade)"
>
> select version();
> "PostgreSQL 8.0.8 on i386-portbld-freebsd5.3, compiled by GCC cc  
> (GCC) 3.4.2 [FreeBSD] 20040728"
>
> -Steve
>
> Stephen Woodbridge wrote:
>> Hi all,
>> Sorry the subject line is not more useful, but let me explain.
>> 1) I have a set of polgons A and another set of polygon B
>> 2) I need to a new set of polygons C that are the union of
>>    A intersect B and A difference B
>> So basically B partially covers A and I want to split polygons in  
>> a into the covered piece(s) and that that is not covered.
>> On the surface this is very straight forward, but it isn't ...
>> insert into newtable (....)
>>    select ..., intersection(A.the_geom, B.the_geom) as the_geom
>>      from A, B
>>      where A.the_geom && B.the_geom and
>>        intersects(A.the_geom, B.the_geom);
>> This works, but I get a bunch of GEOMETRYCOLLECTION objects in the  
>> results and out of 1022 collection I have 2780 additional POLYGONs.
>> So I could do something like:
>> insert into newtable (....)
>>     select * from
>>       (select ..., (dump(the_geom)).geom as the_geom from newtable
>>          where geometrytype(the_geom)='GEOMETRYCOLLECTION' ) as foo
>>     where geometrytype(the_geom)='POLYGON';
>> delete from newtable where geometrytype(the_geom) 
>> ='GEOMETRYCOLLECTION';
>> So now comes the tricky part how do you do the difference? I was  
>> thinking I could take A difference newtable but because I have  
>> multiple little pieces instead of a single piece, I think I am  
>> getting:
>>   (A - Part1), (A - Part2), (A - Part3) instead of
>>   (A - Part1 - Part2 - Part3)
>> So I think I need to union the polygon pieces from each row that  
>> is dumped above. Is that as simple as tossing part of the dump  
>> into a geomunion()? How would that look?
>> -Steve W
>> _______________________________________________
>> postgis-users mailing list
>> postgis-users at postgis.refractions.net
>> http://postgis.refractions.net/mailman/listinfo/postgis-users
>
> _______________________________________________
> postgis-users mailing list
> postgis-users at postgis.refractions.net
> http://postgis.refractions.net/mailman/listinfo/postgis-users




More information about the postgis-users mailing list