PostGIS geometry intersection test results different for Point and Linestring

Question

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?

Answer 1

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