[geos-devel] no outgoing dirEdge found error
Ben Jubb
benjubb at refractions.net
Wed Sep 19 16:54:12 EDT 2007
There is a way to make this particular operation succeed in GEOS, which
is to use the BinaryOp method instead of difference(). For example, the
following code works (in MSVC8):
#include <string>
#include <iostream>
#include <geos/geom/PrecisionModel.h>
#include <geos/geom/GeometryFactory.h>
#include <geos/geom/Geometry.h>
#include <geos/geom/MultiPolygon.h>
#include <geos/io/WKTReader.h>
#include <geos/geom/BinaryOp.h>
#include <geos/operation/overlay/OverlayOp.h>
using namespace geos;
int main (int argc, char *argv[])
{
using namespace operation::overlay;
typedef std::auto_ptr< geom::Geometry > GeomAutoPtr;
std::string mp1s = "MULTIPOLYGON (((88.0541300495810049
-40.3914336120775417,"
" 88.0013586636740683 -40.3952580652168507,"
" 87.9856522556763849 -40.3125682141688557,"
" 88.0410955427405639 -40.3084853410231361,"
" 88.0447770149226727 -40.3320733430391556,"
" 88.0541300495810049 -40.3914336120775417)),"
" ((87.9569549228361041 -40.1614847724049824,"
" 87.9856522556763849 -40.3125682141688557,"
" 87.9594803530699494 -40.3144955269957563,"
" 87.9484023890450430 -40.2597855714132962,"
" 87.9236385364479673 -40.1638760511710515,"
" 87.9569549228361041 -40.1614847724049824)))";
std::string mp2s = "MULTIPOLYGON (((87.7184147523609425
-40.3314878875577705,"
" 87.7106409229361930 -40.2916861830969211,"
" 87.6363928113999862 -39.9074815079801084,"
" 87.8588765594887349 -39.8933226300382628,"
" 88.3245648150878964 -39.8616404104367064,"
" 88.3252910982550503 -39.8653834273614933,"
" 88.4058355104796334 -40.2774282017416283,"
" 88.4066631207306273 -40.2815647501607259,"
" 87.9574816223991007 -40.3146427145853394,"
" 87.7184147523609425 -40.3314878875577705)))";
try {
geos::geom::PrecisionModel *model =
new
geos::geom::PrecisionModel(geos::geom::PrecisionModel::FLOATING);
geos::geom::GeometryFactory *factory = new
geos::geom::GeometryFactory(model);
geos::io::WKTReader *wkt = new geos::io::WKTReader();
geos::geom::MultiPolygon *mp1 = (geos::geom::MultiPolygon
*)wkt->read(mp1s);
geos::geom::MultiPolygon *mp2 = (geos::geom::MultiPolygon
*)wkt->read(mp2s);
std::cout << "Is valid of one = " << mp1->isValid() << std::endl;
std::cout << "Is valid of two = " << mp2->isValid() << std::endl;
//geos::geom::Geometry *diff = mp1->difference(mp2);
GeomAutoPtr gRealRes = BinaryOp(mp1, mp2,
overlayOp(OverlayOp::opDIFFERENCE));
std::cout << mp1->toString() << std::endl;
std::cout << mp2->toString() << std::endl;
//std::cout << diff->toString() << std::endl;
std::cout << gRealRes.get()->toString() << std::endl;
}
catch (std::exception const &se) {
std::cout << "ERROR - " << se.what() << std::endl;
}
catch (...) {
std::cout << "General error" << std::endl;
}
}
BInaryOp does the same difference op, but when that fails it tries again
after shifting the geometry close to the origin. Presumably this gains
a few bits of precision in the intermediate results, that allows the
operation to succeed. This method is perhaps not as robust tho.
b
Martin Davis wrote:
> This computation works fine in the current version of JTS. This may
> indicate a porting bug in GEOS, or a slightly different choice of
> tolerances.
>
> This is a classic case which causes robustness problems - overlaying
> two geometries with linework which is almost coincident.
>
> The only suggestion I can make is to reduce the precision of your data
> slightly in cases which fail. Presumably your source data isn't
> actually accurate to 16 decimal places? It could be a while before
> we can look at the GEOS code for this.
>
> Stuart C Sides wrote:
>>
>> Hi GEOS devels,
>>
>> We are writing a C++ application which is looking for areas of
>> overlap between 1000's of polygons. The app is being
>> developed on Linux with g++ 4.1.0.
>>
>> We first tried the stable release of GEOS 2.2.3 and received this
>> error after processing several polygons:
>>
>> TopologyException: found non-noded intersection between 87.9595
>> -40.3145, 87.9484 -40.2598 and 87.9857 -40.3126, 87.9575 -40.3146
>> 87.9595 -40.3145
>>
>> While reviewing the list archives we found references to major code
>> changes that might fix some of these errors, so we
>> upgraded to RC3 and eventually to an SVN version from 2007-09-18. We
>> got through a lot more polygons without errors,
>> but finally received this error:
>>
>> TopologyException: no outgoing dirEdge found 87.957 -40.1615
>>
>> Note: this coordinate is from the first segment of the second polygon
>> of the first multi-polygon below.
>>
>> In the original code, before I converted the multipolygons to strings
>> to create the example, we received this error:
>>
>> TopologyException: no outgoing dirEdge found 87.7184 -40.3315
>>
>> Note: this coordinate is from the first segment of the second
>> multi-polygon below:
>>
>> I've boiled the error down to just a few lines of code, and would
>> appreciate any suggestions. The error is thrown at the
>> difference.
>>
>> Thanks
>> Stuart
>>
>>
>>
>> #include <string>
>> #include <iostream>
>>
>> #include <geos/geom/PrecisionModel.h>
>> #include <geos/geom/GeometryFactory.h>
>> #include <geos/geom/Geometry.h>
>> #include <geos/geom/MultiPolygon.h>
>> #include <geos/io/WKTReader.h>
>>
>> int main (int argc, char *argv[])
>> {
>> std::string mp1s = "MULTIPOLYGON (((88.0541300495810049
>> -40.3914336120775417,"
>> " 88.0013586636740683 -40.3952580652168507,"
>> " 87.9856522556763849 -40.3125682141688557,"
>> " 88.0410955427405639 -40.3084853410231361,"
>> " 88.0447770149226727 -40.3320733430391556,"
>> " 88.0541300495810049 -40.3914336120775417)),"
>> " ((87.9569549228361041 -40.1614847724049824,"
>> " 87.9856522556763849 -40.3125682141688557,"
>> " 87.9594803530699494 -40.3144955269957563,"
>> " 87.9484023890450430 -40.2597855714132962,"
>> " 87.9236385364479673 -40.1638760511710515,"
>> " 87.9569549228361041 -40.1614847724049824)))";
>>
>> std::string mp2s = "MULTIPOLYGON (((87.7184147523609425
>> -40.3314878875577705,"
>> " 87.7106409229361930 -40.2916861830969211,"
>> " 87.6363928113999862 -39.9074815079801084,"
>> " 87.8588765594887349 -39.8933226300382628,"
>> " 88.3245648150878964 -39.8616404104367064,"
>> " 88.3252910982550503 -39.8653834273614933,"
>> " 88.4058355104796334 -40.2774282017416283,"
>> " 88.4066631207306273 -40.2815647501607259,"
>> " 87.9574816223991007 -40.3146427145853394,"
>> " 87.7184147523609425 -40.3314878875577705)))";
>>
>> try {
>> geos::geom::PrecisionModel *model =
>> new
>> geos::geom::PrecisionModel(geos::geom::PrecisionModel::FLOATING);
>> geos::geom::GeometryFactory *factory = new
>> geos::geom::GeometryFactory(model);
>>
>> geos::io::WKTReader *wkt = new geos::io::WKTReader();
>> geos::geom::MultiPolygon *mp1 = (geos::geom::MultiPolygon
>> *)wkt->read(mp1s);
>> geos::geom::MultiPolygon *mp2 = (geos::geom::MultiPolygon
>> *)wkt->read(mp2s);
>>
>> std::cout << "Is valid of one = " << mp1->isValid() << std::endl;
>> std::cout << "Is valid of two = " << mp2->isValid() << std::endl;
>>
>> geos::geom::Geometry *diff = mp1->difference(mp2);
>>
>> std::cout << mp1->toString() << std::endl;
>> std::cout << mp2->toString() << std::endl;
>> std::cout << diff->toString() << std::endl;
>>
>> }
>> catch (std::exception const &se) {
>> std::cout << "ERROR - " << se.what() << std::endl;
>> }
>> catch (...) {
>> std::cout << "General error" << std::endl;
>> }
>>
>> }
>> ------------------------------------------------------------------------
>>
>> _______________________________________________
>> geos-devel mailing list
>> geos-devel at geos.refractions.net
>> http://geos.refractions.net/mailman/listinfo/geos-devel
>>
>
More information about the geos-devel
mailing list