Re: [postgis-users] Getting TopologyExections when trying to node linestrings to create an overlay
On Mon, Feb 9, 2015 at 6:25 AM, Sandro Santilli s...@keybit.net wrote: First of all I confirm it still happens with GEOS=3.5.0dev-CAPI-1.9.0 r4038. Second, I took a look at a random set (geom_set_id=1) and I found it pretty big. That's to say you could probably further reduce the dataset for the ticket. That set contains 109 polygons, I can get the error by attempting to union the boundaries of the first 40 in that set, and I'm sure you can further reduce the input. Thanks for taking a look. I'll work on doing that when I can find the time, but I don't expect that to be a fast process at all. Even just checking for the pairwise case took a decent amount of time to develop the query and took overnight to finish running (even with optimizations like a bounding box intersection test). I don't really have any good heuristic that could narrow down the possibilies for reproducing, so I don't see much option other than to brute force it possibly with some kind of filter. That's why I didn't put more effort into shrinking the input set to begin with. Is a PostGIS database dump an okay format to provide the shapes, or would you prefer something else? I suppose I could dump groups into shapefiles or something like that if it's more convenient. On Mon, Feb 9, 2015 at 7:00 AM, Rémi Cura remi.c...@gmail.com wrote: this is a precision related issue, coordinates are way too big and should be translated. I understand what you mean by that (as in floating point problems due to size), but these coordinates are pretty typical. This is a standard UTM projection, zone 15N in middle America. It's even predefined in PostGIS' list of spatial references. I'm told ESRI had this kind of problem years ago, but they dealt with it as far as I know. While I would choose PostGIS over ESRI any day, this could be viewed by my coworkers as a good argument against using PostGIS; it represents a serious reliability concern since it applies to a broad range of functions and 2D projections. Basically, if you use any projection with coordinates of this size (of which there are a good number), there seems to be no telling when any function will just blow up in your face at random. On Mon, Feb 9, 2015 at 7:04 AM, Rémi Cura remi.c...@gmail.com wrote: Hey, I executed your data, the following command solve the problem (with very recent GEOS for me) (POSTGIS=2.2.0dev r12846 GEOS=3.5.0dev-CAPI-1.9.0 r0 PROJ=Rel. 4.8.0, 6 March 2012 GDAL=GDAL 2.0.0dev, released 2014/04/16 LIBXML=2.8.0 RASTER) [Snip] The change compared to your approach : convert input to table of simple polygons, (no array, no multi). Then translate to improve precision in geos computing Then the union. I don't really understand what you are trying to do, but ist_union seems dangerous and quit ineffective for that . I'm looking at this function call in your code: ST_Union( ST_MakePolygon(ST_ExteriorRing(geom)) ). That call seems to remove all holes and then create the union of all the covered areas (a single multipolygon that covers all areas covered by the originals). That is not what I'm trying to do; I already have an outer boundary that I could use for that purpose. Here's what I'm trying to do. I'm trying to create an overlay (as I said in the first sentence of my original e-mail but didn't elaborate on), in the sense that the term is used by these two articles: http://boundlessgeo.com/2014/10/postgis-training-creating-overlays/ or http://trac.osgeo.org/postgis/wiki/UsersWikiExamplesOverlayTables. What they do is they create a sort of Venn diagram, if you will; they take all the polygons and create a new polygon for each area with a different set of overlappying polygons. Both of them use ST_Union in the same way I am: they take a set of linestrings and union them to create a fully noded linestring, and then they use this linework to create a bunch of new polygons. I'm leveraging ST_Boundary instead of ST_ExteriorRing in my code because I need to preserve holes, but that shouldn't change the results as far as I know. In what way is ST_Union dangerous and quite ineffective for the purpose of noding lingstrings? If that's true, I'm apparently not the only one who has some misconceptions, since both these articles use it the same way. I should also note that I don't show the whole process here; I'm only showing the part where I'm noding the linestrings because that's where this error occurs. Once I have the polygons, I also go back and relate them to the attributed rows in the original tables (of which there are several). That said, I'm trying out the translate and cutting out multi, but I am still using an array in my actual code for a reason. Namely, I need to be able to do this with *different* sets of geometries. These geometries come as the result of selecting polygons from 3 or 4 different tables, each having completely disparate sets of attributes. So basically, I need to be able to use an arbitrary query to get the group of polygons
Re: [postgis-users] Getting TopologyExections when trying to node linestrings to create an overlay
Hey, last chance for you ^^ I successfully was able to remove errors in 2 ways : - first your geometry array contains a lot of duplicated geometry. ( rest of this is based on geom_set_id= 1) I suspect you forgot a where in an inner join, or something like this. Removing the duplicate in your data seems to solve the problem (removing duplicate : DISTINCT ON (geom)) I generated each pair of geometry in this set, none gave the error individually (which lessen the chance of a legit GEOS problem) - second, without removing the duplicates, a combination of tanslate and snaptogrid(0.5) was sufficient to also remove the error. I would recommend that you analyse a little bit your data. For instance, simply deuplicating on geom make your data set going from --15317 geometry in arrays to --9934 Snapping to a 0.5 grid before deduplicating further reduce the data set to 9617 (which might indicates that somewhere in your workflow you have a precision related issue). I never used ARCGIS, but I can bet you that ARCGIS is not precision-safe. The only product that is truly safe is CGAL, which comes with other type of constraints. You could also probably use GRASS safely if cleaning the data with v.clean on import. Anyway, there is no safe tool, only safe way to use it. You could very easily create another SRS that is a translation of your original srs to increase precision. So you simply change your workflow to ST_Transform before computing, then ST_Transform after computing. Lastly, I don't see the interest of using an array of geom here because you don't need to be able to do my_array[N] (no need to access). So you could simply use a geometry collection, so your input is a geom, and not a geom[]. Cheers, Rémi-C 2015-02-18 10:39 GMT+01:00 BladeOfLight16 bladeofligh...@gmail.com: On Mon, Feb 9, 2015 at 6:25 AM, Sandro Santilli s...@keybit.net wrote: First of all I confirm it still happens with GEOS=3.5.0dev-CAPI-1.9.0 r4038. Second, I took a look at a random set (geom_set_id=1) and I found it pretty big. That's to say you could probably further reduce the dataset for the ticket. That set contains 109 polygons, I can get the error by attempting to union the boundaries of the first 40 in that set, and I'm sure you can further reduce the input. Thanks for taking a look. I'll work on doing that when I can find the time, but I don't expect that to be a fast process at all. Even just checking for the pairwise case took a decent amount of time to develop the query and took overnight to finish running (even with optimizations like a bounding box intersection test). I don't really have any good heuristic that could narrow down the possibilies for reproducing, so I don't see much option other than to brute force it possibly with some kind of filter. That's why I didn't put more effort into shrinking the input set to begin with. Is a PostGIS database dump an okay format to provide the shapes, or would you prefer something else? I suppose I could dump groups into shapefiles or something like that if it's more convenient. On Mon, Feb 9, 2015 at 7:00 AM, Rémi Cura remi.c...@gmail.com wrote: this is a precision related issue, coordinates are way too big and should be translated. I understand what you mean by that (as in floating point problems due to size), but these coordinates are pretty typical. This is a standard UTM projection, zone 15N in middle America. It's even predefined in PostGIS' list of spatial references. I'm told ESRI had this kind of problem years ago, but they dealt with it as far as I know. While I would choose PostGIS over ESRI any day, this could be viewed by my coworkers as a good argument against using PostGIS; it represents a serious reliability concern since it applies to a broad range of functions and 2D projections. Basically, if you use any projection with coordinates of this size (of which there are a good number), there seems to be no telling when any function will just blow up in your face at random. On Mon, Feb 9, 2015 at 7:04 AM, Rémi Cura remi.c...@gmail.com wrote: Hey, I executed your data, the following command solve the problem (with very recent GEOS for me) (POSTGIS=2.2.0dev r12846 GEOS=3.5.0dev-CAPI-1.9.0 r0 PROJ=Rel. 4.8.0, 6 March 2012 GDAL=GDAL 2.0.0dev, released 2014/04/16 LIBXML=2.8.0 RASTER) [Snip] The change compared to your approach : convert input to table of simple polygons, (no array, no multi). Then translate to improve precision in geos computing Then the union. I don't really understand what you are trying to do, but ist_union seems dangerous and quit ineffective for that . I'm looking at this function call in your code: ST_Union( ST_MakePolygon(ST_ExteriorRing(geom)) ). That call seems to remove all holes and then create the union of all the covered areas (a single multipolygon that covers all areas covered by the originals). That is not what I'm trying to do; I already have an outer
Re: [postgis-users] Getting TopologyExections when trying to node linestrings to create an overlay
Computing finished successfuly on your data. The change I did is : ST_Polygonize(ST_Node(ST_Collect(geom)) ) It is awfully slow tough (1 hour for all you dataset I think) Cheers, Rémi-C 2015-02-18 13:07 GMT+01:00 Rémi Cura remi.c...@gmail.com: Hey, last chance for you ^^ I successfully was able to remove errors in 2 ways : - first your geometry array contains a lot of duplicated geometry. ( rest of this is based on geom_set_id= 1) I suspect you forgot a where in an inner join, or something like this. Removing the duplicate in your data seems to solve the problem (removing duplicate : DISTINCT ON (geom)) I generated each pair of geometry in this set, none gave the error individually (which lessen the chance of a legit GEOS problem) - second, without removing the duplicates, a combination of tanslate and snaptogrid(0.5) was sufficient to also remove the error. I would recommend that you analyse a little bit your data. For instance, simply deuplicating on geom make your data set going from --15317 geometry in arrays to --9934 Snapping to a 0.5 grid before deduplicating further reduce the data set to 9617 (which might indicates that somewhere in your workflow you have a precision related issue). I never used ARCGIS, but I can bet you that ARCGIS is not precision-safe. The only product that is truly safe is CGAL, which comes with other type of constraints. You could also probably use GRASS safely if cleaning the data with v.clean on import. Anyway, there is no safe tool, only safe way to use it. You could very easily create another SRS that is a translation of your original srs to increase precision. So you simply change your workflow to ST_Transform before computing, then ST_Transform after computing. Lastly, I don't see the interest of using an array of geom here because you don't need to be able to do my_array[N] (no need to access). So you could simply use a geometry collection, so your input is a geom, and not a geom[]. Cheers, Rémi-C 2015-02-18 10:39 GMT+01:00 BladeOfLight16 bladeofligh...@gmail.com: On Mon, Feb 9, 2015 at 6:25 AM, Sandro Santilli s...@keybit.net wrote: First of all I confirm it still happens with GEOS=3.5.0dev-CAPI-1.9.0 r4038. Second, I took a look at a random set (geom_set_id=1) and I found it pretty big. That's to say you could probably further reduce the dataset for the ticket. That set contains 109 polygons, I can get the error by attempting to union the boundaries of the first 40 in that set, and I'm sure you can further reduce the input. Thanks for taking a look. I'll work on doing that when I can find the time, but I don't expect that to be a fast process at all. Even just checking for the pairwise case took a decent amount of time to develop the query and took overnight to finish running (even with optimizations like a bounding box intersection test). I don't really have any good heuristic that could narrow down the possibilies for reproducing, so I don't see much option other than to brute force it possibly with some kind of filter. That's why I didn't put more effort into shrinking the input set to begin with. Is a PostGIS database dump an okay format to provide the shapes, or would you prefer something else? I suppose I could dump groups into shapefiles or something like that if it's more convenient. On Mon, Feb 9, 2015 at 7:00 AM, Rémi Cura remi.c...@gmail.com wrote: this is a precision related issue, coordinates are way too big and should be translated. I understand what you mean by that (as in floating point problems due to size), but these coordinates are pretty typical. This is a standard UTM projection, zone 15N in middle America. It's even predefined in PostGIS' list of spatial references. I'm told ESRI had this kind of problem years ago, but they dealt with it as far as I know. While I would choose PostGIS over ESRI any day, this could be viewed by my coworkers as a good argument against using PostGIS; it represents a serious reliability concern since it applies to a broad range of functions and 2D projections. Basically, if you use any projection with coordinates of this size (of which there are a good number), there seems to be no telling when any function will just blow up in your face at random. On Mon, Feb 9, 2015 at 7:04 AM, Rémi Cura remi.c...@gmail.com wrote: Hey, I executed your data, the following command solve the problem (with very recent GEOS for me) (POSTGIS=2.2.0dev r12846 GEOS=3.5.0dev-CAPI-1.9.0 r0 PROJ=Rel. 4.8.0, 6 March 2012 GDAL=GDAL 2.0.0dev, released 2014/04/16 LIBXML=2.8.0 RASTER) [Snip] The change compared to your approach : convert input to table of simple polygons, (no array, no multi). Then translate to improve precision in geos computing Then the union. I don't really understand what you are trying to do, but ist_union seems dangerous and quit ineffective for that . I'm looking at this function call in
Re: [postgis-users] Getting TopologyExections when trying to node linestrings to create an overlay
On Sun, Feb 8, 2015 at 9:31 PM, John Abraham j...@hbaspecto.com wrote: Well I could say that using PostGIS ST_Intersects with messy data always seems to give me TopologyExceptions. I've had luck with various combinations of ST_SnapToGrid and ST_Buffer(0), but with messy data there always seems to be some weird case that requires manually edits. Sorry to be the bearer of bad news. Frankly, I don't understand why the GEOS library has to throw that error. I would encourage you to isolate particular problems and file bug reports. Improvements in GEOS to eliminate the underlying error(s) would certainly be welcome. That's pretty much what my e-mail does; it even includes a sample database to reproduce the issue. I brought it up here in case some discussion or slimming down needed/could be done before filing a bug report, and to be extra sure that it's a GEOS issue and not a PostGIS once, given the complexity of the query I'm running. ___ postgis-users mailing list postgis-users@lists.osgeo.org http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
Re: [postgis-users] Getting TopologyExections when trying to node linestrings to create an overlay
Hey, I executed your data, the following command solve the problem (with very recent GEOS for me) (POSTGIS=2.2.0dev r12846 GEOS=3.5.0dev-CAPI-1.9.0 r0 PROJ=Rel. 4.8.0, 6 March 2012 GDAL=GDAL 2.0.0dev, released 2014/04/16 LIBXML=2.8.0 RASTER) : DROP TABLE IF EXISTS unique_polygon ; CREATE TABLE unique_polygon AS SELECT geom_set_id, row_number() over() as gid, ST_Translate(dmp.geom, - 385614, - 4795454 ) AS geom FROM error_generating_polygons,unnest(polygons) as geomn, st_dump(geomn) as dmp; CREATE INDEX ON unique_polygon USING GIST(geom) ; DRoP TABLE IF EXISTS unioned_poly ; CREATE TABLE unioned_poly AS SELECT ST_Union( ST_MakePolygon(ST_ExteriorRing(geom)) ) FROM unique_polygon GROUP BY geom_set_id (150 sec) The change compared to your approach : convert input to table of simple polygons, (no array, no multi). Then translate to improve precision in geos computing Then the union. I don't really understand what you are trying to do, but ist_union seems dangerous and quit ineffective for that . Of course reducing the number of useless points before union make it 10 times faster . DRoP TABLE IF EXISTS unioned_poly ; CREATE TABLE unioned_poly AS SELECT ST_Union( ST_Buffer( ST_MakePolygon( ST_ExteriorRing( ST_SImplifyPreserveTopology( geom ,10 ) ) ) ,1 ) ) FROM unique_polygon GROUP BY geom_set_id (17 sec) Cheers, Rémi-C 2015-02-09 13:00 GMT+01:00 Rémi Cura remi.c...@gmail.com: Hey Sandro, this is a precision related issue, coordinates are way too big and should be translated. Cheers, Rémi-C 2015-02-09 12:25 GMT+01:00 Sandro Santilli s...@keybit.net: On Tue, Feb 03, 2015 at 11:28:35AM -0500, BladeOfLight16 wrote: https://drive.google.com/file/d/0B_6I7kRgE8teVUpha2Q4ZlNDMWs/view?usp=sharing . ... DO $$ DECLARE problem_row error_generating_polygons%ROWTYPE; BEGIN FOR problem_row IN (SELECT * FROM error_generating_polygons) LOOP BEGIN PERFORM ST_Union(ST_Boundary(geom)) FROM UNNEST(problem_row.polygons) p (geom); RAISE NOTICE 'geom_set_id % succeeded', problem_row.geom_set_id; EXCEPTION WHEN OTHERS THEN RAISE NOTICE 'Error for geom_set_id % (Code %): %', problem_row.geom_set_id, SQLSTATE, SQLERRM; END; END LOOP; END $$; First of all I confirm it still happens with GEOS=3.5.0dev-CAPI-1.9.0 r4038. Second, I took a look at a random set (geom_set_id=1) and I found it pretty big. That's to say you could probably further reduce the dataset for the ticket. That set contains 109 polygons, I can get the error by attempting to union the boundaries of the first 40 in that set, and I'm sure you can further reduce the input. So my suggestion: 1) file the ticket 2) attach the _smallest_ input that reproduces the problem About ST_IsValid: lines are always valid, so there's no need to test. Most likely this is a robustness issue failing to deal with very close but not equal lines. NOTE: I've tried my reduced input (~40) geoms against the topology builder and it also resulted in errors, until I specified a tolerance of 1e-4. --strk; () Free GIS Flash consultant/developer /\ http://strk.keybit.net/services.html ___ postgis-users mailing list postgis-users@lists.osgeo.org http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users ___ postgis-users mailing list postgis-users@lists.osgeo.org http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
Re: [postgis-users] Getting TopologyExections when trying to node linestrings to create an overlay
Hey Sandro, this is a precision related issue, coordinates are way too big and should be translated. Cheers, Rémi-C 2015-02-09 12:25 GMT+01:00 Sandro Santilli s...@keybit.net: On Tue, Feb 03, 2015 at 11:28:35AM -0500, BladeOfLight16 wrote: https://drive.google.com/file/d/0B_6I7kRgE8teVUpha2Q4ZlNDMWs/view?usp=sharing . ... DO $$ DECLARE problem_row error_generating_polygons%ROWTYPE; BEGIN FOR problem_row IN (SELECT * FROM error_generating_polygons) LOOP BEGIN PERFORM ST_Union(ST_Boundary(geom)) FROM UNNEST(problem_row.polygons) p (geom); RAISE NOTICE 'geom_set_id % succeeded', problem_row.geom_set_id; EXCEPTION WHEN OTHERS THEN RAISE NOTICE 'Error for geom_set_id % (Code %): %', problem_row.geom_set_id, SQLSTATE, SQLERRM; END; END LOOP; END $$; First of all I confirm it still happens with GEOS=3.5.0dev-CAPI-1.9.0 r4038. Second, I took a look at a random set (geom_set_id=1) and I found it pretty big. That's to say you could probably further reduce the dataset for the ticket. That set contains 109 polygons, I can get the error by attempting to union the boundaries of the first 40 in that set, and I'm sure you can further reduce the input. So my suggestion: 1) file the ticket 2) attach the _smallest_ input that reproduces the problem About ST_IsValid: lines are always valid, so there's no need to test. Most likely this is a robustness issue failing to deal with very close but not equal lines. NOTE: I've tried my reduced input (~40) geoms against the topology builder and it also resulted in errors, until I specified a tolerance of 1e-4. --strk; () Free GIS Flash consultant/developer /\ http://strk.keybit.net/services.html ___ postgis-users mailing list postgis-users@lists.osgeo.org http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users ___ postgis-users mailing list postgis-users@lists.osgeo.org http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
Re: [postgis-users] Getting TopologyExections when trying to node linestrings to create an overlay
On Tue, Feb 03, 2015 at 11:28:35AM -0500, BladeOfLight16 wrote: https://drive.google.com/file/d/0B_6I7kRgE8teVUpha2Q4ZlNDMWs/view?usp=sharing. ... DO $$ DECLARE problem_row error_generating_polygons%ROWTYPE; BEGIN FOR problem_row IN (SELECT * FROM error_generating_polygons) LOOP BEGIN PERFORM ST_Union(ST_Boundary(geom)) FROM UNNEST(problem_row.polygons) p (geom); RAISE NOTICE 'geom_set_id % succeeded', problem_row.geom_set_id; EXCEPTION WHEN OTHERS THEN RAISE NOTICE 'Error for geom_set_id % (Code %): %', problem_row.geom_set_id, SQLSTATE, SQLERRM; END; END LOOP; END $$; First of all I confirm it still happens with GEOS=3.5.0dev-CAPI-1.9.0 r4038. Second, I took a look at a random set (geom_set_id=1) and I found it pretty big. That's to say you could probably further reduce the dataset for the ticket. That set contains 109 polygons, I can get the error by attempting to union the boundaries of the first 40 in that set, and I'm sure you can further reduce the input. So my suggestion: 1) file the ticket 2) attach the _smallest_ input that reproduces the problem About ST_IsValid: lines are always valid, so there's no need to test. Most likely this is a robustness issue failing to deal with very close but not equal lines. NOTE: I've tried my reduced input (~40) geoms against the topology builder and it also resulted in errors, until I specified a tolerance of 1e-4. --strk; () Free GIS Flash consultant/developer /\ http://strk.keybit.net/services.html ___ postgis-users mailing list postgis-users@lists.osgeo.org http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
Re: [postgis-users] Getting TopologyExections when trying to node linestrings to create an overlay
On Tue, Feb 3, 2015 at 11:28 AM, BladeOfLight16 bladeofligh...@gmail.com wrote: I'm trying to create a polygon overlay. The basic process is relatively simple: 1) Get the boundaries 2) Union the boundaries to node the linestrings 3) Polygonize the noded outlines 4) Filter out holes using a contains or intersects test. The problem I'm running into is that I'm getting GEOSUnaryUnion: TopologyException: found non-noded intersection between LINESTRING errors from ST_Union. [a lot snipped] I haven't seen any response to this. I was just wondering if anyone else had a chance or intentions to look this over. Granted, it's pretty long and involved (sorry for that), but I thought all the details I included were important. I do know it went through the mailing list; someone on IRC helped me find it in the... I guess it's not the archives; I don't know what it's called. But the online browsing mechanism. Thanks to anyone who's taking a look. ___ postgis-users mailing list postgis-users@lists.osgeo.org http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
Re: [postgis-users] Getting TopologyExections when trying to node linestrings to create an overlay
Il 09/02/2015 02:31, John Abraham ha scritto: Well I could say that using PostGIS ST_Intersects with messy data always seems to give me TopologyExceptions. I've had luck with various combinations of ST_SnapToGrid and ST_Buffer(0), but with messy data there always seems to be some weird case that requires manually edits. Hi John, have you tried ST_MakeValid? All the best. -- Paolo Cavallini - www.faunalia.eu QGIS PostGIS courses: http://www.faunalia.eu/training.html *New course* QGIS for naturalists: http://www.faunalia.eu/en/nat_course.html ___ postgis-users mailing list postgis-users@lists.osgeo.org http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
Re: [postgis-users] Getting TopologyExections when trying to node linestrings to create an overlay
Well I could say that using PostGIS ST_Intersects with messy data always seems to give me TopologyExceptions. I've had luck with various combinations of ST_SnapToGrid and ST_Buffer(0), but with messy data there always seems to be some weird case that requires manually edits. Sorry to be the bearer of bad news. Frankly, I don't understand why the GEOS library has to throw that error. I would encourage you to isolate particular problems and file bug reports. Improvements in GEOS to eliminate the underlying error(s) would certainly be welcome. -- John Abraham On Feb 8, 2015, at 1:43 PM, BladeOfLight16 bladeofligh...@gmail.com wrote: On Tue, Feb 3, 2015 at 11:28 AM, BladeOfLight16 bladeofligh...@gmail.com mailto:bladeofligh...@gmail.com wrote: I'm trying to create a polygon overlay. The basic process is relatively simple: 1) Get the boundaries 2) Union the boundaries to node the linestrings 3) Polygonize the noded outlines 4) Filter out holes using a contains or intersects test. The problem I'm running into is that I'm getting GEOSUnaryUnion: TopologyException: found non-noded intersection between LINESTRING errors from ST_Union. [a lot snipped] I haven't seen any response to this. I was just wondering if anyone else had a chance or intentions to look this over. Granted, it's pretty long and involved (sorry for that), but I thought all the details I included were important. I do know it went through the mailing list; someone on IRC helped me find it in the... I guess it's not the archives; I don't know what it's called. But the online browsing mechanism. Thanks to anyone who's taking a look. ___ postgis-users mailing list postgis-users@lists.osgeo.org http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users ___ postgis-users mailing list postgis-users@lists.osgeo.org http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
[postgis-users] Getting TopologyExections when trying to node linestrings to create an overlay
I'm trying to create a polygon overlay. The basic process is relatively simple: 1) Get the boundaries 2) Union the boundaries to node the linestrings 3) Polygonize the noded outlines 4) Filter out holes using a contains or intersects test. The problem I'm running into is that I'm getting GEOSUnaryUnion: TopologyException: found non-noded intersection between LINESTRING errors from ST_Union. I've taken the liberty of creating a script that will set up a database to reproduce what I'm seeing. I've gzipped it in the interest of size, since it's got a fairly sizable set of polygons. Let me know if this delivery method is a problem; I'll be happy to use other means. You can find it here: https://drive.google.com/file/d/0B_6I7kRgE8teVUpha2Q4ZlNDMWs/view?usp=sharing. (I originally tried to attach it, but the PostGIS mailing list rejected it.) The contents are a plain SQL file. It won't create a database for you, but it will create the PostGIS extension IF NOT EXISTS. Specifically for this problem, it creates and populates this table: CREATE TABLE error_generating_polygons ( geom_set_id SERIAL PRIMARY KEY, outer_boundary geometry NOT NULL, polygons geometry[] NOT NULL, error_code text NOT NULL, error_message text NOT NULL ); It's a little weird, so I think I should provide some explanation. The primary key is just a surrogate. The outer_boundary can be ignored for the purposes of this error; it's there for me for that filtering out holes and chaff after I get the raw overlay back. The polygons column is the most interesting; it contains a set of polygons that will reproduce this error. (I'll get into why I use an array in a minute.) The error_code and error_message are the values of SQLSTATE and SQLERRM I get when noding the boundaries; I've included them so that anyone can compare if they get different results. The reason I have an array of polygons is that this was the simplest method of providing a large group of polygons. The data you're seeing here is not sample data I've made up. This is real world data provided to me by a client, who has produced them primarily through a mixture of GPS, manual processing, and probably some more or less automated processing. In short, it's pretty messy, and these are the actual polygons I need to overlay. As for why I'm giving you these big groups instead of breaking things down into smaller chunks, I did try that, at least somewhat. I found that if I broke this up into pairs of polygons, I could only reproduce the error with a *single* pair out of thousands of combinations. So sorry for the shear number of polygons here, but apparently, the errors come at least partly from some cumulative effect. I have verified that every polygon in there is a valid one: SELECT geom_set_id, ST_IsValidDetail(geom) FROM (SELECT geom_set_id, UNNEST(polygons) geom FROM error_generating_polygons ) x WHERE NOT ST_IsValid(geom); That query gives me an empty set, plus I do plenty of checking before letting a geometry get into the system. Now, to get down to business. Here's the query you can run to see errors: DO $$ DECLARE problem_row error_generating_polygons%ROWTYPE; BEGIN FOR problem_row IN (SELECT * FROM error_generating_polygons) LOOP BEGIN PERFORM ST_Union(ST_Boundary(geom)) FROM UNNEST(problem_row.polygons) p (geom); RAISE NOTICE 'geom_set_id % succeeded', problem_row.geom_set_id; EXCEPTION WHEN OTHERS THEN RAISE NOTICE 'Error for geom_set_id % (Code %): %', problem_row.geom_set_id, SQLSTATE, SQLERRM; END; END LOOP; END $$; I also tried a variant of this: I replaced ST_Bounrary(geom) with ST_Boundary(ST_SnapToGrid(geom, 0.01)). I tested on two separate installations of PostgreSQL/PostGIS. Believe it or not, I'm getting different results from these. Here are the specs and results for each one: 1) Runs on CentOS, installed from CentOS repositories (I believe. I need to double check with IT to be 100% sure). Version Info: PostgreSQL 9.3.5 on x86_64-unknown-linux-gnu, compiled by gcc (GCC) 4.4.7 20120313 (Red Hat 4.4.7-4), 64-bit POSTGIS=2.1.4 r12966 GEOS=3.4.2-CAPI-1.8.2 r3921 PROJ=Rel. 4.8.0, 6 March 2012 GDAL=GDAL 1.9.2, released 2012/10/08 LIBXML=2.7.6 LIBJSON=UNKNOWN RASTER Every row throws an error on this one. This is where the errors in the table come from. Adding ST_SnapToGrid makes it error out on only one row. This one is more similar to my production machine at present. 2) Runs on Debian, installed from PostgreSQL maintained repository. Version info: PostgreSQL 9.4.0 on x86_64-unknown-linux-gnu, compiled by gcc (Debian 4.7.2-5) 4.7.2, 64-bit POSTGIS=2.1.5 r13152 GEOS=3.3.3-CAPI-1.7.4 PROJ=Rel. 4.7.1, 23 September 2009 GDAL=GDAL 1.9.0, released 2011/12/29 LIBXML=2.8.0 LIBJSON=UNKNOWN RASTER This one errored out less; only about 1/3 of the rows error out. Adding ST_SnapToGrid also makes it error out on only one row. An oddity to note: the Debian install has a more up to date PostGIS
