postgresql 13.1: precision of spatial operations

Hi Vadim,

On 29.11.2022 19:39, Вадим Самохин wrote:

Hi there,
I have polygons in a table and I'm fetching those that contain a
specific point. The problem is that when I'm checking against a point
which is really close (~5 meters) but is not contained within a
polygon, it gets fetched nevertheless.

Steps to reproduce:
1.
create table zones (
    zone_id int,
    zone_polygon polygon,
    description text
);
create index zones__zone_polygon on zones using gist(zone_polygon
poly_ops);

2. insert into zones (zone_polygon) values
('(37.6040241,55.7609641),(37.6240129,55.7519367),(37.6215344,55.7536616),(37.6172064,55.7559509),(37.6126178,55.7584013),(37.6088694,55.7622611),(37.60747,55.7633072),(37.6040241,55.7609641)');

3. Here is this polygon on a map:
https://www.keene.edu/campus/maps/tool/?coordinates=37.6040241%2C%2055.7609641%0A37.6240129%2C%2055.7519367%0A37.6215344%2C%2055.7536616%0A37.6172064%2C%2055.7559509%0A37.6126178%2C%2055.7584013%0A37.6088694%2C%2055.7622611%0A37.6074700%2C%2055.7633072%0A37.6040241%2C%2055.7609641

4. Check whether a point is contained within a polygon: select
count(1) from zones where zone_polygon @>
'(37.617635,55.755814)'::polygon;
 count
-------
     1
(1 row)

5. But actually it's not (sorry, couldn't find a way to represent this
point on the same map. Use street view, it's more convenient to see
that):
https://www.keene.edu/campus/maps/tool/?coordinates=37.617635%2C%2055.755814

First of all, you could reduce your large example to a smaller one:
select
'(37.6220129,55.7519367),(37.6215344,55.7536616),(37.6172064,55.7559509),(37.6220129,55.7519367)'::polygon
@> '(37.617900,55.755814)'::point;

6. Just in case, here are the images.
The polygon:
image.png
And the point:
image.png

gnuplot could be used for easier visualization.

Am I doing anything wrong? Any idea how to fix that?

The problem is that the geometric comparison operations are "fuzzy", see
here
https://github.com/postgres/postgres/blob/master/src/include/utils/geo_decls.h
I can recommend you either multiply your coordinates by 10000 or rebuild
Postgres with EPSILON = 0.

Regards,
Ivan

Thank you so much Ivan, it worked!

ср, 30 нояб. 2022 г. в 00:22, Ivan Panchenko <i.panchenko@postgrespro.ru>:

On 11/30/2022 9:48 AM, Вадим Самохин wrote:

Thank you so much Ivan, it worked!

Can you give any more detail on which approach you took, for the sake of
future followers?

ср, 30 нояб. 2022 г., 20:51 Brad White <b55white@gmail.com>:

On 11/30/2022 9:48 AM, Вадим Самохин wrote:

Thank you so much Ivan, it worked!

Can you give any more detail on which approach you took, for the sake of
future followers?

Sure, I multiplied all points' coordinates by a factor of 10^6.

Here are the steps to demonstrate that the solution Ivan gave worked:

1. create table zones (
zone_id int,
zone_polygon polygon,
description text
);
create index zones__zone_polygon on zones using gist(zone_polygon poly_ops);

2. insert into zones (zone_polygon) values
('(37622012.9,55751936.7),(37621534.4,55753661.6),(37617206.4,55755950.9),(37622012.9,55751936.7)');

3. select count(1) from zones where zone_polygon @>
'(37617635,55755814)'::polygon;
count
-------
0
(1 row)