[geos-devel] Find self-intersections
Frederik Ramm
frederik at remote.org
Sat Nov 28 18:09:27 EST 2009
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"
More information about the geos-devel
mailing list