<div dir="ltr"><div>On Mon, Feb 9, 2015 at 6:25 AM, Sandro Santilli <span dir="ltr"><<a href="mailto:strk@keybit.net" target="_blank">strk@keybit.net</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">First of all I confirm it still happens with GEOS="3.5.0dev-CAPI-1.9.0 r4038".<br>
Second, I took a look at a random set (geom_set_id=1) and I found it pretty<br>
big. That's to say you could probably further reduce the dataset for the<br>
ticket. That set contains 109 polygons, I can get the error by attempting<br>
to union the boundaries of the first 40 in that set, and I'm sure you can<br>
further reduce the input.<br></blockquote><div><br></div>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.<br><br></div>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.<br><div><div class="gmail_extra"><br><div class="gmail_quote">On Mon, Feb 9, 2015 at 7:00 AM, Rémi Cura <span dir="ltr"><<a href="mailto:remi.cura@gmail.com" target="_blank">remi.cura@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div>this is a precision related issue, coordinates are way too big and should be translated.</div></div></blockquote><div><br></div><div>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.<br></div></div></div><div class="gmail_extra"><br><div class="gmail_quote">On Mon, Feb 9, 2015 at 7:04 AM, Rémi Cura <span dir="ltr"><<a href="mailto:remi.cura@gmail.com" target="_blank">remi.cura@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div style="font-size:12.8px"><div>Hey,<br></div>I executed your data,<br></div><div style="font-size:12.8px">the following command solve the problem (with very recent GEOS for me)<br>(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)<br></div></div></blockquote><blockquote style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex" class="gmail_quote"><div> </div></blockquote><blockquote style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex" class="gmail_quote"><div>[Snip] <br></div></blockquote><blockquote style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex" class="gmail_quote"><div> </div></blockquote><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div style="font-size:12.8px"></div><div style="font-size:12.8px">The change compared to your approach : convert input to table of simple polygons, (no array, no multi).<br></div><div style="font-size:12.8px">Then translate to improve precision in geos computing<br><div style="font-size:12.8px">Then the union.<br></div><div style="font-size:12.8px">I don't really understand what you are trying to do,<br></div>but ist_union seems dangerous and quit ineffective for that .</div></div></blockquote><div><br></div><div>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.<br><br>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: <a href="http://boundlessgeo.com/2014/10/postgis-training-creating-overlays/">http://boundlessgeo.com/2014/10/postgis-training-creating-overlays/</a> or <a href="http://trac.osgeo.org/postgis/wiki/UsersWikiExamplesOverlayTables">http://trac.osgeo.org/postgis/wiki/UsersWikiExamplesOverlayTables</a>. 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).<br><br></div><div>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 <i>different</i> 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 to be passed into a reusable overlay process. As a result, a function is a natural fit. I opted to use a normal function that takes an array. The only other option I can think of is some kind of aggregate function, which I didn't investigate doing as I'm not sure whether that might be more or less reliable. If you think an aggregate would be better or know of a better way to accomplish that, I'm all ears.<br></div><div style="font-size:12.8px"> <br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div style="font-size:12.8px">Of course reducing the number of useless points before union make it 10 times faster .<br><br>DRoP TABLE IF EXISTS unioned_poly ; <br>CREATE TABLE unioned_poly AS <br>SELECT ST_Union( <br> ST_Buffer(<br> ST_MakePolygon(<br> ST_ExteriorRing(<br> ST_SImplifyPreserveTopology(<br> geom<br> ,10<br> )<br> )<br> )<br> ,1 ) <br> )<br>FROM unique_polygon<br>GROUP BY geom_set_id<br></div><div style="font-size:12.8px">(17 sec)<br></div></div></blockquote><div><br></div><div>10 meters is a lot of precision to lose. That's over a tenth of a football field. I might be willing to simplify by 1 meter at most, but even that's a little big for my taste. My client and users expect to see the same shapes that got input come back out.<br><br></div><div>Thanks again to everyone who took a look.<br></div></div></div></div></div>