Unexpected results from ST_Intersects(LINESTRING, POLYGON)

Question

I am having a bit of difficultly getting accurate results from a table of mine that contains geometric data. I'm currently using PostGIS 2.1.3, PostgreSQL 9.3.4 and GEOS 3.3.8, though I've also reproduced the issue on early versions.

I have a table, Track, which maintains a LINESTRING summarizing its various locations through time. I need to build the line in real time. When a new Track is added, I create a line from 2 points, then subsequently update the line when new information becomes available. To handle that, I've been using the ST_AddPoint function. I have noticed that I don't get all of the expected tracks when I issue an ST_Intersects query. Lastly, I have tried the same approach with ST_MakeLine which does work.

I have managed to boil the issue down to a very simple test case which shows the failure. I start with a clean spatially enabled database ( postgis.sql and spatial_ref_sys.sql loaded ) .

create table track
(
    id serial
);

select AddGeometryColumn('track', 'track_region', 0, 'LINESTRING', 2);

create index trk_rgn_idx on track using GIST ( track_region  );

insert into track (track_region) values(ST_GeomFromText('LINESTRING(600.6 129.877, 585.785 170.165 )'));

update track set track_region = ST_AddPoint( track_region, ST_GeomFromText('POINT(572.741 211.551)'));
update track set track_region = ST_AddPoint( track_region, ST_GeomFromText('POINT(558.694 251.969)'));
update track set track_region = ST_AddPoint( track_region, ST_GeomFromText('POINT(544.725 292.666)'));
update track set track_region = ST_AddPoint( track_region, ST_GeomFromText('POINT(531.029 333.661)'));
update track set track_region = ST_AddPoint( track_region, ST_GeomFromText('POINT(518.303 375.281)'));
update track set track_region = ST_AddPoint( track_region, ST_GeomFromText('POINT(506.037 416.302)'));
update track set track_region = ST_AddPoint( track_region, ST_GeomFromText('POINT(492.891 455.689)'));
update track set track_region = ST_AddPoint( track_region, ST_GeomFromText('POINT(480.456 491.382)'));
update track set track_region = ST_AddPoint( track_region, ST_GeomFromText('POINT(469.236 522.899)'));

select id from track where( ST_Intersects(track_region, 'POLYGON((100 400,500 400,500 1200,100 1200,100 400))'));

Failed return:

 id
----
(0 rows)

When I attempt the ST_MakeLine approach, below, I get the expect result.

drop table track;

create table track
(
    id serial
);

select AddGeometryColumn('track', 'track_region', 0, 'LINESTRING', 2);

create index trk_rgn_idx on track using GIST ( track_region  );

insert into track (track_region) values(ST_GeomFromText('LINESTRING(600.6 129.877, 585.785 170.165 )'));

update track set track_region = ST_MakeLine( track_region, ST_GeomFromText('POINT(572.741 211.551)'));
update track set track_region = ST_MakeLine( track_region, ST_GeomFromText('POINT(558.694 251.969)'));
update track set track_region = ST_MakeLine( track_region, ST_GeomFromText('POINT(544.725 292.666)'));
update track set track_region = ST_MakeLine( track_region, ST_GeomFromText('POINT(531.029 333.661)'));
update track set track_region = ST_MakeLine( track_region, ST_GeomFromText('POINT(518.303 375.281)'));
update track set track_region = ST_MakeLine( track_region, ST_GeomFromText('POINT(506.037 416.302)'));
update track set track_region = ST_MakeLine( track_region, ST_GeomFromText('POINT(492.891 455.689)'));
update track set track_region = ST_MakeLine( track_region, ST_GeomFromText('POINT(480.456 491.382)'));
update track set track_region = ST_MakeLine( track_region, ST_GeomFromText('POINT(469.236 522.899)'));

select id from track where( ST_Intersects(track_region, 'POLYGON((100 400,500 400,500 1200,100 1200,100 400))'));

Which returns the expected.

 id
----
  1
(1 row)

Can anyone figure what I'm doing wrong? My guess is that it's something with how I'm using ST_AddPoint to create the line because if I just insert the line all at once everything works fine. Also though, if I select the line after it's been added with ST_AddPoint, then use it in a query like so, it works.

select ST_Intersects('LINESTRING(600.6 129.877,585.785 170.165,572.741 211.551,558.694 251.969,544.725 292.666,531.029 333.661,518.303 375.281,506.037 416.302,492.891 455.689,480.456 491.382,469.236 522.899)', 'POLYGON((100 400,500 400,500 1200,100 1200,100 400))' );

 id
----
  1
(1 row)

Thanks!

Answer 1

You have found a bug with ST_Point, which has been simplified here as #2845.

I would consider using a simple method of creating a LineString geometry, rather than repetitively updating the table. ST_MakeLine has an example using an array of point geometries, which I think is a good start. You can also use array_agg to help build these point arrays, if you have tables of point geometries.

Answer 2

In addition to what Mike T has said about you having hit a bug, you can "fix" your query by wrapping the ST_Intersects with ST_MakeValid or by buffering by a tiny amount.

Both

select id from track where ST_Intersects(ST_MakeValid(track_region),
   ST_GeomfromText('POLYGON((100 400,500 400,500 1200,100 1200,100 400))'));

and

 select id from track where ST_Intersects(ST_Buffer(track_region,0.0000001),
   ST_GeomfromText('POLYGON((100 400,500 400,500 1200,100 1200,100 400))'));

return your intersected id, though the second one is admitedly hacky. If you run ST_IsValid(track_region) you get true, but ST_MakeValid "solves" the query problem anyway, even though the geometry is considered to be valid. ST_Buffer(geom, 0) is a hack for fixing node order, which does not work here without a small buffer, but is a useful thing to know.

As an aside, your query above:

select id from track where( ST_Intersects(track_region, 'POLYGON((100 400,500 400,500 1200,100 1200,100 400))'));

is invalid, it should be:

select id from track where ST_Intersects(track_region, ST_GeomFromText('POLYGON((100 400,500 400,500 1200,100 1200,100 400))'));