[GIS] PostGIS geometry intersection test results different for Point and Linestring

intersectionlinestringpostgis

Can anyone please tell me what I'm doing wrong, and how do I avoid the problem? I'll attempt to clarify what I'm doing as things go along.

The geom column is a multipolygon column, all values are SRID=4326.

archive=> select st_intersects( GeomFromewkt('SRID=4326;LINESTRING(-16.899216666666668 12.610816666666667,-16.899216666666668 12.610816666666667)') 
, z.geom) as intersects
from zone z
where z.zone_id=2;
 intersects 
------------
 f
(1 row)

archive=> select   st_intersects( GeomFromewkt('SRID=4326;POINT(-16.899216666666668 12.610816666666667)') 
, z.geom) as intersects
from zone z
where z.zone_id=2;
 intersects 
------------
 t
(1 row)

This is the text representation of the geom:

MULTIPOLYGON(((-18.7273147400488 12.8862644964914,-18.7334078119781 12.3609663554343,-16.5597662148828 12.3796303212717,-16.6302631480867 13.1209501528692,-18.7273147400488 12.8862644964914)))

How can the linestring not intersect with the geometry, yet a single point does? Is this some sort of rounding error, or a problem because the linestring has a length of zero? If the length is the problem, how do I avoid this problem?

Best Answer

Yup, it looks like that is the behaviour from JTS and GEOS. The problem is that your LINESTRING is invalid. If you have PostGIS 2.0, you can use ST_MakeValid(geometry) to fix the LINESTRING to a POINT.

This query verifies your bug, and uses ST_MakeValid as a workaround.

WITH data AS (SELECT 'POLYGON((150 280, 99 215, 190 210, 150 280))'::geometry AS poly,
                     'POINT(170 240)'::geometry AS pt)

SELECT ST_Intersects(poly, pt) AS intersects_poly_pt,
       ST_Intersects(poly, ST_MakeLine(pt, pt)) AS intersects_poly_line,
       ST_IsValid(ST_MakeLine(pt, pt)) AS isvalid_line,
       ST_AsText(ST_MakeValid(ST_MakeLine(pt, pt))) AS valid_geom,
       ST_Intersects(poly, ST_MakeValid(ST_MakeLine(pt, pt))) AS intersects_poly_valid_geom
FROM data;

with results (using psql's \x option):

NOTICE:  Too few points in geometry component at or near point 170 240
-[ RECORD 1 ]--------------+---------------
intersects_poly_pt         | t
intersects_poly_line       | f
isvalid_line               | f
valid_geom                 | POINT(170 240)
intersects_poly_valid_geom | t

If you are using a previous version of PostGIS (pre 2.0), then you can cast from an invalid LINESTRING to a box2d, then back to a GEOMETRY. For a two-vertex LINESTRING with same coordinates, this turns into a POINT. Here is the PostGIS 1.5 version of the above:

WITH data AS (SELECT 'POLYGON((150 280, 99 215, 190 210, 150 280))'::geometry AS poly,
                     'LINESTRING(170 240, 170 240)'::geometry AS line)

SELECT ST_Intersects(poly, line) AS intersects_bad_line,
       ST_IsValid(line) AS isvalid_line,
       ST_AsText(CASE WHEN NOT ST_IsValid(line) THEN line::box2d::geometry
                      ELSE line END) as valid_geom,
       ST_Intersects(poly, CASE WHEN NOT ST_IsValid(line) THEN line::box2d::geometry
                                ELSE line END) AS intersects_poly_valid_geom
FROM data;

with results:

NOTICE:  Too few points in geometry component at or near point 170 240
NOTICE:  Too few points in geometry component at or near point 170 240
NOTICE:  Too few points in geometry component at or near point 170 240
-[ RECORD 1 ]--------------+---------------
intersects_bad_line        | f
isvalid_line               | f
valid_geom                 | POINT(170 240)
intersects_poly_valid_geom | t
Related Question