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@geos.refractions.net
http://geos.refractions.net/mailman/listinfo/geos-devel
_______________________________________________
geos-devel mailing list
geos-devel@geos.refractions.net
http://geos.refractions.net/mailman/listinfo/geos-devel
