[postgis-users] How best to do this?

[Stephen Woodbridge]

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