[geos-devel] Find self-intersections

[Frederik Ramm]

Hi,

strk wrote:
> What do you really need to do at the end of the day ?
> Want a simple ring ? A polygon ?

I'll happily discuss every last detail of the application, there's 
nothing secret about it, but I'll try to stick to the relevant stuff!

I'm writing an application that creates multipolygons from OpenStreetMap 
data. OSM has a somewhat clumsy implementation of multipolygons, where 
users simply group a number of ways (=linestrings) and then say: These, 
together, form a forest (or so). Such a group of ways may or may not 
make up a valid multipolygon - sometimes the linestrings have 
self-intersections; sometimes they intersect with others; sometimes 
there are gaps in a ring so that it cannot be closed, etc.

In addition to creating a proper polygon, I also want to emit detailed 
error messages so people can fix broken things. Initially I thought I 
could use the Polygonizer, but there were issues. Firstly it didn't to 
multipolygons. But also, I want to know more about what it does; I want 
to know exactly which linestrings have gone into which rings of the 
polygon (above all, whether a certain linestring has become part of an 
inner or outer ring), and if the Polygonizer cannot do my bidding then I 
want to know exactly what the problem was. (I know I can ask the 
Polygonizer for a list of dangles, cut lines, and other brokenness, but 
there are just too many cases where the Polygonizer will just throw an 
exception and not tell me where exactly the problem was.)

I have now built code that "manually" assembles a number of rings from 
the bunch of linestrings, uses GEOS's "contains" method to sort them 
into outers and inners, then creates a number of polygons from the rings 
and ultimately a multipolygon from the polygons. That way, I can detect 
many problems (e.g. gaps) before even trying to make a polygon; but in 
situations where e.g. two inner rings intersect, I will get an error 
from createPolygon, and to pinpoint that, I currently have to compute 
the intersection of almost everything with everything.

Since writing the original e-Mail, I have tried to find out more about 
how I can use GEOS classes to help me here. The available C++ API docs 
seem to only scratch the surface of what is possible - does everyone 
else simply know what to do, or is there a secret bit of documentation I 
haven't found yet? (Is there perhaps JTS documentation that can be 
applied to GEOS?) I'm now playing around with this:

GeometryGraph graph(0, all_my_ways);
LineIntersector li;
graph.computeSelfNodes(&li, true);
vector<Edge *> *edges = graph.getEdges();
for (unsigned int e=0; e<edges->size(); e++)
{
     for(EdgeIntersectionList::const_iterator i = 
edges->at(e)->eiList.begin(); i!= edges->at(e)->eiList.end(); i++)
     {
         EdgeIntersection *ei=*i;
         // ignore intersections at ends of ways
         if (ei->dist != 0 || (ei->segmentIndex>0 && ei->segmentIndex < 
edges->at(e)->getNumPoints()-1))
         {
             // process ei->coord, it seems to be an
             // interesting intersection
         }
     }
}

This seems to catch all of the cases I'm interested in (at least all of 
the testcases that I built), and it is very fast. Whether it is really 
"correct", I have no clue - especially the bit about only using edge 
intersections with dist != 0 and not at the ends of ways is more or less 
a guess.

It seems that I should be able to, instead of iterating over the edges 
as above, iterate over my input linestrings and use the "findEdge" 
method in GeometryGraph to retrieve the edge; this would allow me to 
detect the exact linestring where an intersection has occurred, right?

Bye
Frederik

-- 
Frederik Ramm  ##  eMail frederik at remote.org  ##  N49°00'09" E008°23'33"