[postgis-users] Getting TopologyExections when trying to node linestrings to create an overlay
Rémi Cura
remi.cura at gmail.com
Wed Feb 18 04:55:13 PST 2015
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.cura at 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 <bladeoflight16 at gmail.com>:
>
>> On Mon, Feb 9, 2015 at 6:25 AM, Sandro Santilli <strk at 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.cura at 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.cura at 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 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.
>>
>>
>>> 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)
>>>
>>
>> 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.
>>
>> Thanks again to everyone who took a look.
>>
>> _______________________________________________
>> postgis-users mailing list
>> postgis-users at lists.osgeo.org
>> http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
>>
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20150218/b6d72e68/attachment.html>
More information about the postgis-users
mailing list