Error: Invalid number of points in LinearRing

[tiiipponen]

Why this sql query returns error:

SELECT ST_Intersects(ST_GeomFromText('CURVEPOLYGON(
  (25495445.625 6671632.625, 25495445.625 6671711.375, 25495555.375 6671711.375, 25495555.375 6671632.625, 25495445.625 6671632.625),
  COMPOUNDCURVE(
    CIRCULARSTRING(25495368.0441 6671726.9312,25495368.3959388 6671726.93601515,25495368.7478 6671726.9333),
    (25495368.7478 6671726.9333,25495368.0441 6671726.9312)))'),ST_MakeEnvelope(25495443.625, 6671631.625, 25495556.375, 6671712.375))

ERROR: First argument geometry could not be converted to GEOS: IllegalArgumentException: Invalid number of points in LinearRing found 3 - must be 0 or >= 4
SQL state: XX000
Context: SQL function "st_intersects" statement 1

If I put hole geometry to boundary geometry it doesn't give error:

select ST_Intersects(ST_GeomFromText('CURVEPOLYGON(
  COMPOUNDCURVE(
    CIRCULARSTRING(25495368.0441 6671726.9312,25495368.3959388 6671726.93601515,25495368.7478 6671726.9333),
    (25495368.7478 6671726.9333,25495368.0441 6671726.9312)))'),ST_MakeEnvelope(25495443.625, 6671631.625, 25495556.375, 6671712.375))

st_intersects
boolean

f

[tiiipponen]

I tried to look this a bit more just by looking code, but couldn't find anything clear.
I noticed, that even ST_Intersects does not return error for hole part, it actually returns false, which is not right answer. I noticed also, that validation shows problem for hole too, so I created simpler example:

select ST_IsValidDetail(
  ST_GeomFromText(
    'CURVEPOLYGON(
      COMPOUNDCURVE(
        CIRCULARSTRING(
          1.0441 2.9312,1.3959388 2.93601515,1.7478 2.9333
        ), (
          1.7478 2.9333,1.0441 2.9312
        )
      )
    )'
  )
);

Result:

"(f,"IllegalArgumentException: Invalid number of points in LinearRing found 3 - must be 0 or >= 4",)"

I tried to look code and only thing that makes me suspicious is this:
liblwgeom/lwgeom_geos.c ->
-> LWGEOM2GEOS(const LWGEOM *lwgeom, int autofix)
-> LWGEOM *lwgeom_stroked = lwgeom_stroke(lwgeom, 32);
-> liblwgeom/lwstroke.c -> lwgeom_stroke(const LWGEOM *geom, uint32_t perQuad)

I don't know the stroking value perQuad, but if it is default PostGIS value 32, it leads to arc length 0.78 m, which is more than length of actual arc 0.71m. That would create straight line as a stroking result of that arc and that can give the error message "Invalid number of points...".
I adjusted example coordinates a little so that arc was longer and error message disappeared.

If this is the reason, then this problem is not so easy to fix. We can have very short arcs in our data and using same absolute stroking length means troubles.
I would think more, if things could do completely without stroking.
At least I would use relative stroking compared to arc length.
As a stroking method, I like most the tolerance method. Tolerance means maximum difference between original arc and stroked line. For short arcs the resulting vertice count comes more important.
I created following code with Python for you to copy (and test!) freely:
User can give parameters accuracy "tolerance" or number of "vertices" in resulting stroked line.

    # Radius (r), centroid (cx,cy) and angles (a1 and a3) are now
    # calculated to arc geometry (p1, p2, p3).
    # Then we start stroking the geometry

    angleArc = a3 - a1
    toleranceVertices = 0
    verticesCount = 0
    StrokingTolerance = 0.0
    
    # If user gave the stroking tolerance, we calculate maximum angle of one stroke sector.
    # from equation: cos(angle / 2) = (radius - tolerance) / radius
    # We calculate also minimum count of vertices for whole arc that satisfies tolerance.

    if tolerance > 0.0:
        StrokingTolerance = tolerance
        if StrokingTolerance > 2 * r:
            StrokingTolerance = 2 * r
        angleSectorFromTolerance = 2 * math.acos((r - StrokingTolerance) / r)
        toleranceVertices = int(abs(angleArc / angleSectorFromTolerance))

    # If user gave exact amount of vertices to use,
    # we use it or we use vertice count calculated from tolerance,
    # which ever is bigger so meaning more accurate.

    if vertices > toleranceVertices:
        verticesCount = vertices
    else:
        verticesCount = toleranceVertices

    # Now calculate exact angle of one stroke sector from count of vertices.

    angleSector = angleArc / (verticesCount + 1)

    # Calculate vertice points and create line from them.

    pointList = [(p1x, p1y)]

    for i in range(1, verticesCount + 1):
        currentAngle = a1 + i * angleSector
        newPointX = cx + r * math.cos(currentAngle)
        newPointY = cy + r * math.sin(currentAngle)
        pointList.append((newPointX, newPointY))
        
    pointList.append((p3x, p3y))

Ok. Lots of text and actual error is maybe located somewhere else.
I wish good luck for solving this, because function ST_Intersects is widely used in QGIS, Openlayers, etc. and this error results missing geometries and error messages in data that holds arcs.

[robe]

strk you think your work in #3772 can accommodate this issue.

It's still a problem in after I test with trunk code, but since you are in that area anyway, you might have some thoughts on this.

anyrate I'm pushing this to 2.4 since 2.4 lwstroke is different from 2.3

[strk]

Tricky problem to fix. It's about collapsed polygon due to reduced precision when going from curve to lines. Augmenting the number of segments should fix it. The default is 32, which collapses. Use 64 to get it right:

strk=# select ST_IsValid(ST_CurveToLine(
  ST_GeomFromText(
    'CURVEPOLYGON(
      COMPOUNDCURVE(
        CIRCULARSTRING(
          1.0441 2.9312,1.3959388 2.93601515,1.7478 2.9333
        ), (
          1.7478 2.9333,1.0441 2.9312
        )
      )
    )'
  )
, 32));
NOTICE:  IllegalArgumentException: Invalid number of points in LinearRing found 3 - must be 0 or >= 4
 st_isvalid 
------------
 f
(1 row)

strk=# select ST_IsValid(ST_CurveToLine(
  ST_GeomFromText(
    'CURVEPOLYGON(
      COMPOUNDCURVE(
        CIRCULARSTRING(
          1.0441 2.9312,1.3959388 2.93601515,1.7478 2.9333
        ), (
          1.7478 2.9333,1.0441 2.9312
        )
      )
    )'
  )
, 64));
 st_isvalid 
------------
 t
(1 row)

[strk]

The curved part of the ring is almost a line. Default linearization using 32 segments-per-quadrant end up approximating it with a single straight line.

I guess the fix here would be to enforce a minimum of 2 segments for any non-closed curve and 3 segments for a circle. Dunno if that should be the default or a parametric configuration.

The work done in #3772 could help in that it allows specifying linearization constraints in different ways and allows specifying a flag, but there's no specific flag to ensure at least 2 segments are used per-curve.

[strk]

To make it clear: this ticket is about _curves_, which are a relatively "recent" and still not fully supported feature. Other organizations (big and small) have funded advancement in the support for curves, are you willing to join the effort ?

Patches are also a great way to contribute.

Ensuring each arc is defined by at least 2 segments (reduce max error and max angle, augmenting segments-per-quadrant) would seem a good plan.

[tiiipponen]

Yes, you speak well. We are going to pay, so we can get this fixed. Thanks for your efforts.

[strk]

Now that different tolerance types are supported, I want to report that the subject geometry will be turned invalid with a max-deviation stroking tolerance of 1e-2 units but not with max-deviation stroking tolerance of 1e-3 units. This gives an indication of how close to "invalid" the input is in the first place (a movement of 0.04 units turns it to invalid, one of 0.03 units does not). I still think it makes sense to add a flag to ensure each arc (or maybe each arc-or-quadrant - to generalize it for the closed curve case) has at least 2 segments.

[tiiipponen]

Hi. Nice to hear that this is taking steps forward. I agree the tolerance concept with minimum of two (or five and maximum of ? (10 000?)) segments can be good fix here. I looked from one of our data, what is the smallest tolerance to keep arcs with two segments and here are the results in ascending order:

order_id     tolerance
----------------------
1            4.509183781920001e-8
2            2.0990961502320715e-7
...
20           8.002188010891587e-6
...
61           0.0001082151904370221
...
98           0.0011843443344314863
...
145          0.012285292649195867
...
251          0.1034409626306001
...

That is real life example and "minimum of 2 segments" parameter is needed.

The question then rises: what is the default tolerance for linearization? In Oracle we have used tolerance 0.0005 in metadata "user_sdo_geom_metadata". I have not seen place to put tolerance for table-data in PostGIS. At least it is not in view "public.geometry_columns". I can imagine that not everybody are satisfied with highest possible tolerance, because of its impact to efficiency. On the other hand too big tolerance gives more erroneus results from functions like ST_Intersects() etc. It would be good if users could select tolerance as they want in table level or somehow. I am willing to hear your thoughts about this.

I am also willing to hear your thoughts about PostGIS's future development for arc-geometry handling. Are there plans to do topology- queries etc. without need to linearizing arcs? I can see some similar base development done to PostGIS code already.

I also want to add that we are willing to test this development with our data.

[strk]

The question then rises: what is the default tolerance for linearization?

Default is 32 segments per quadrant (which means single segment if the arc is 2.81250 degrees or less)

The above default is documented in ​http://postgis.net/docs/manual-2.4/ST_CurveToLine.html

[tiiipponen]

Sorry for unclear question. I know the default tolerance is 32 segments per quadrant, which is too parse and which is the reason for all problems in this ticket.

The real question is: if we use this tolerance-parameter you propose inside PostGIS with some default tolerance value, then what would be good value? And because different users want to use different tolerance, how can we make it possible for them to set that parameter in table-level?

Some example calculations about tolerance values:

From equation:
cos(angle_sector / 2) = (radius - tolerance) / radius
=>
tolerance = radius - radius * cos(angle_sector / 2)

If radius = 1000m and angle_sector = (PII/2) / 32 = 0.049rad
=>
tolerance = 30cm and length of arc is 49m
That is just too big tolerance and specially:
all shorter arcs with that kind of radius,
would be straight one segment line after linearization.

If we want tolerance to be 0.0005m then angle_sector would be
angle_sector = 2 * acos((radius - tolerance) / radius)
=>
angle_sector = 0.002rad = 785 segments per quadrant

If we put 10 000 segments per quadrant for arc with radius of 1000m
=>
tolerance = 0.000003m which is quite enough for practical use I think.

Because arc is just mathematical entity to express line and necessarily not any real structure, people can have arcs with radius for example from 1mm to 100km. With 100km radius and 10 000 segments per quadrant, we get tolerance = 0.0003m

So, as a result, either we want possibility to put "segments per quadrant" parameter something between 1000 - 10 000 or we want possibility to put "tolerance" parameter something between 0.0001 - 0.001. The later one is normal and better way to parametrize linearizations in GIS-applications and databases.

P.S I fixed my real data tolerance values in my last message.

[strk]

New ST_CurveToLine signature which landed in 2.4.0 already supports specifying a tolerance in terms of max error (in source units). That tolerance type will result in radius dependent "segments-per-quadrant". Does that solve your problem ?

[tiiipponen]

The minimum solution is to put minimum output segments to two. And why not put it to five segments, as there was in code sample I send to you. Five segments looks better arc-geometry than two.

This later discussion started from your comment 4 weeks ago and concerns accuracy and is somehow linked to this update request. 32 segments per quadrant doesn't give accurate results from spatial queries.

It is good that ST_CurveToLine has developed that way and user can give his/her own tolerance as parameter.

But this ticket concerns on function ST_Intersects() and some other similar functions, where PostGIS is doing linearization internally and temporarily when sending geometry to GEOS library function, who doesn't understand arc-geometries. Many of these functions returns geometries with original arcs, so I can not linearize them beforehand with function ST_CurveToLine.

So how to put proper tolerance for internal functions?

For example in Oracle it can be given in the geometry metadata definition for a table or as an input parameter to certain functions.

PostGIS doesn't have either of these concepts. Now only in function ST_CurveToLine.

So how to give that tolerance for these internal handlings? I am interested to hear your ideas.

If somebody asks it from me, I have some proposals:

1. New tolerance attribute to view "public.geometry_columns".
   
2. Some other place to put this default tolerance parameter for example in database level.
   
3. Tolerance parameter to all functions, who sends its geometries to GEOS.
   
4. Default tolerance to 0.001 or at least 0.01 if efficiency is problem. If it is not, then why not to put default tolerance to 0.0005 or even 0.0001.
   
5. Default tolerance to 10 000 segments per quadrant for all internal linearizations.

[strk]

Tolerance parameter is also being discussed elsewhere, might be mergeable as a concept (check wiki).

A GUC would be another possibility.

Retaining topology could be a non-optional behavior (would mean avoiding to collapse a ring into a flat line, as it was the subject of this specific ticket).

[tiiipponen]

Yes, there are lots of discussions about tolerance in PostGIS Wiki. It means, that users consider it important.

"PostGIS Grand Unified Custom Variables" (GUCs) are just fine for me because all our data has same tolerance anyway. Some other users may disagree and want some more. Anyway the concept of tolerance will be huge step forward for PostGIS.

And yes. Retaining topology by for example keeping some minimum amount of segments in linearization, is reasonable solution here in this ticket.

[tiiipponen]

Here is my proposal for code lines (6 rows) in function lwarc_linearize() in file lwstroke.c.
I didn't test it, but afterwards I am ready to test this fix well with our database data.

static int
lwarc_linearize(POINTARRAY *to,
                 const POINT4D *p1, const POINT4D *p2, const POINT4D *p3,
                 double tol, LW_LINEARIZE_TOLERANCE_TYPE tolerance_type,
                 int flags)
{
	POINT2D center;

...    

row 235:
	/* Angles of each point that defines the arc section */
	a1 = atan2(p1->y - center.y, p1->x - center.x);
	a2 = atan2(p2->y - center.y, p2->x - center.x);
	a3 = atan2(p3->y - center.y, p3->x - center.x);

	LWDEBUGF(2, "lwarc_linearize A1:%g (%g) A2:%g (%g) A3:%g (%g)",
		a1, a1*180/M_PI, a2, a2*180/M_PI, a3, a3*180/M_PI);

	/* Calculate total arc angle, in radians */
	double angle = clockwise ? a1 - a3 : a3 - a1;
	if ( angle < 0 ) angle += M_PI * 2;

	/* Setting minimum of 6 segments if needed */
	if ( (angle / increment) < 6.0 )
	{
		LWDEBUGF(2, "lwarc_linearize: setting minimum of 6 segments");
		increment = angle / 6.0;
	}

	if ( flags & LW_LINEARIZE_FLAG_SYMMETRIC )
	{{
		LWDEBUGF(2, "lwarc_linearize SYMMETRIC requested - total angle %g deg",
			         angle * 180 / M_PI);
		if ( flags & LW_LINEARIZE_FLAG_RETAIN_ANGLE )
		{{
			/* Number of steps */
			int steps = trunc(angle / increment);
... 

[strk]

Why 6 ? Isn't 2 always enough to avoid flattening a curve ? Note that the more segments you force the more likely it is to hit robustness issues in case of a very short arc.

[strk]

To correct myself: 2 is not enough if the arc is 360 degrees (while it is in all other cases). For 360 degrees you need 3.

[tiiipponen]

If you seriously want 2 and 3, then ok. I don't understand what robustness issues you are afraid of. Circle and arc with 6 segments looks like circle and arc and is more accurate and I think it doesn't cost anything. But ok, can do what you think is best.

[strk]

2 and 3 would make sense as being a "minimum". 6 would be just arbitrary. Accuracy for an arc really depends on the arc size.

6 segments on a 1e-15 degrees arc could be too much (even 2 segments could be too much in that case).

If you want more segments in a circle we could set not a minimum number of segments per arc but per arc-or-quadrant. So if the minimum is 2 you would have at least 8 segments for a ricle (2*4)

[tiiipponen]

I actually mean with accuracy also that in many cases we use topological queries like point in polygon and in many applications the point is created near edge of the polygon.
For example with simple arithmetics we can see that
equilateral triangle inside circle covers 20.7% of the circle and
equilateral hexagon inside circle covers 82.7% of the circle.
By putting more segments, we help "point in polygon"-queries to get more reliable results. And the same thing applies to polygons with short arcs too.
Theoretically speaking...
This is maybe not problem in real world, but if there is no extra cost, then why not.

And as you say, very short arcs needs same special handling if they have six or two segments. And in our case, so short arcs are so rare, that they can be edited away by hand if needed.

Putting this "per arc-or-quadrant" can be ok.

But as I said, you have freedom to do, what you think is easiest and best way to do. It is ok to me.

[tiiipponen]

I suggest a compromise that you put minimum of three segments to all arcs and then you don't need to worry, if it is circle or not.

[tiiipponen]

No progress in this. Even this is very important and easy to fix. What a shame.

[dbaston]

PR at ​https://github.com/postgis/postgis/pull/320

[dbaston]

In 16998:

Impose min number of segments per arc during linearization

This commit sets a fixed minimum of two segments per non-circular arc and three
segments per circular arc.

Funded by City of Helsinki.

Fixes #3719
Closes ​https://github.com/postgis/postgis/pull/320

[Raul Marin]

In 17677:

lwarc_linearize: Iterate over the number of segments

Closes ​https://github.com/postgis/postgis/pull/458
Closes #4407
References #3719