PostGIS testing if point intersects scaled polygon sometimes gets wrong answer

Question

I need an algorithm where I can check if a point falls within a polygon within a small tolerance (meaning the point is completely within the polygon or sufficiently close to the border). To implement this, I was planning on using the st_scale function to expand the polygon by a factor, then check the point with st_within. Note that I can't use st_buffer because st_buffer expects an absolute radius of the buffer but I need a relative radius (ex. Determine if the point is within the polygon or within the polygon that has been scaled up by 5%).

To accomplish this, I've written the following tester code. I take a polygon, scale it up by 5%, and then test to see if the original polygon's centroid falls within the scaled-up polygon:

SELECT ST_within(
  st_centroid(polygon),
  st_scale(
    polygon,
    'POINT(1.05 1.05)'::geometry,
    st_centroid(polygon)
));

However, sometimes this function fails for no apparent reason (meaning it says the polygon's centroid is not within the scaled-up polygon) and I have no idea why. This feels like it's buggy, but I can't imagine I've found a bug in postGIS.

Passing testcases:

-- drawing a 1x1 box at the origin of the x-y plane
SELECT ST_Within(
  st_centroid('MULTIPOLYGON(((0 0, 1 0, 1 1, 0 1, 0 0)))'::geometry),
  st_scale(
    'MULTIPOLYGON(((0 0, 1 0, 1 1, 0 1, 0 0)))'::geometry,
    'POINT(1.05 1.05)'::geometry,
    st_centroid('MULTIPOLYGON(((0 0, 1 0, 1 1, 0 1, 0 0)))'::geometry))
);

-- box in the first quadrant at an origin of (10, 10)
SELECT ST_Within(
  st_centroid('MULTIPOLYGON(((10 10, 11 10, 11 11, 10 11, 10 10)))'::geometry),
  st_scale(
    'MULTIPOLYGON(((10 10, 11 10, 11 11, 10 11, 10 10)))'::geometry,
    'POINT(1.05 1.05)'::geometry,
    st_centroid('MULTIPOLYGON(((10 10, 11 10, 11 11, 10 11, 10 10)))'::geometry))
);

Failing testcase:

-- 1x1 box with the bottom left corner at (98, 28)
SELECT ST_within(
  st_centroid('MULTIPOLYGON(((98 28,99 28,99 29,98 29,98 28)))'::geometry),
  st_scale(
    'MULTIPOLYGON(((98 28,99 28,99 29,98 29,98 28)))'::geometry,
    'POINT(1.05 1.05)'::geometry,
    st_centroid('MULTIPOLYGON(((98 28,99 28,99 29,98 29,98 28)))'::geometry))
);

Interestingly, though, casting the scaled shape to WKT and then back to a geometry DOES work.

-- 1x1 box with the bottom left corner at (98, 28)
SELECT ST_within(
  st_centroid('MULTIPOLYGON(((98 28,99 28,99 29,98 29,98 28)))'::geometry),
  st_geomfromewkt(st_asewkt(st_scale(
    'MULTIPOLYGON(((98 28,99 28,99 29,98 29,98 28)))'::geometry,
    'POINT(1.05 1.05)'::geometry,
    st_centroid('MULTIPOLYGON(((98 28,99 28,99 29,98 29,98 28)))'::geometry))))
);

So, yeah, what is going on here???? Have I found a bug in PostGIS?

Answer 1

I think you've found a bug in older versions of postgis. I used docker and some servers to test with several versions of postgres and postgis. With version 3.1.4 the code works properly.

In the past I've stumbled in some corner cases were postgis gave unexpected results. Those were related to GEOS, that according to the manual is the library that implements the ST_Within function, but in this case the GEOS library is the same in all tests.

From what your working query that transforms to ewkt and back shows, this should be related with an internal operation in postgis that got corrected. I tried but did not find a bug that I could relate to this.

I post two of the tests that I've done:

postgres 12, postgis 3.0.0

The bug is present:

test=# select postgis_full_version();
                                                                                              postgis_full_version                                                                                              
----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
 POSTGIS="3.0.0 r17983" [EXTENSION] PGSQL="120" GEOS="3.7.1-CAPI-1.11.1 27a5e771" PROJ="Rel. 5.2.0, September 15th, 2018" LIBXML="2.9.4" LIBJSON="0.12.1" LIBPROTOBUF="1.3.1" WAGYU="0.4.3 (Internal)" TOPOLOGY
(1 row)

test=# SELECT ST_within(
  st_centroid('MULTIPOLYGON(((98 28,99 28,99 29,98 29,98 28)))'::geometry),
  st_scale(
    'MULTIPOLYGON(((98 28,99 28,99 29,98 29,98 28)))'::geometry,
    'POINT(1.05 1.05)'::geometry,
    st_centroid('MULTIPOLYGON(((98 28,99 28,99 29,98 29,98 28)))'::geometry))
);
 st_within 
-----------
 f
(1 row)

postgres 13, postgis 3.1.4

The code works as expected:

test=# select postgis_full_version();
                                                                                              postgis_full_version                                                                                               
-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
 POSTGIS="3.1.4 ded6c34" [EXTENSION] PGSQL="130" GEOS="3.7.1-CAPI-1.11.1 27a5e771" PROJ="Rel. 5.2.0, September 15th, 2018" LIBXML="2.9.4" LIBJSON="0.12.1" LIBPROTOBUF="1.3.1" WAGYU="0.5.0 (Internal)" TOPOLOGY
(1 row)

test=# SELECT ST_within(
  st_centroid('MULTIPOLYGON(((98 28,99 28,99 29,98 29,98 28)))'::geometry),
  st_scale(
    'MULTIPOLYGON(((98 28,99 28,99 29,98 29,98 28)))'::geometry,
    'POINT(1.05 1.05)'::geometry,
    st_centroid('MULTIPOLYGON(((98 28,99 28,99 29,98 29,98 28)))'::geometry))
);
 st_within 
-----------
 t
(1 row)