Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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?

share|improve this question
    
Can you post the WKT of the multipolygon? –  underdark Apr 26 '12 at 20:52
    
Sure - I just used the st_astext function - I'm not sure why it didn't say SRID=4326; at the front, but it is of the same type. –  Petriborg Apr 26 '12 at 21:00
1  
To get the SRID, use ST_AsEWKT (extended well-known text) –  Mike T Apr 26 '12 at 22:21
add comment

1 Answer 1

up vote 5 down vote accepted

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
share|improve this answer
    
Unfortunately we're using 1.5 - is there an alternative method? –  Petriborg Apr 26 '12 at 21:25
1  
For 1.5, you can try casting from geometry::box2d::geometry, as updated above. –  Mike T Apr 26 '12 at 22:19
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.