[postgis-users] Help with spatial query (bug??)

Brent Wood pcreso at pcreso.com
Tue Apr 5 22:38:55 PDT 2005


--- Carl Anderson <carl.anderson at vadose.org> wrote:

Thanks Carl.

Apologies for the length of this post, but I can't see why what I think is
effectively the same query returns two different results. 


My problem turned out to be that a linestring with 2 points does not intersect
with (or relate to) the polygon containing it if the startpoint = endpoint. So
some 10000 records with 0 length lines disappeared from the join. 

I'm not sure if this constitutes a bug or not (I'm using v0.8 & can't change
that easily)

All the linestrings in this data subset have 2 points.
Each of the two points is the same, so the lines are zero length.
Each cell is a lat/long square, so the bounding box is the actual polygon.
(In the main dataset, where non-zero length lines are present, the queries
appear to work correctly for non-zero length linestrings)

Query one:
psql $DB -c "select t.event_key,
                    c.cell_id,
                    x(startpoint(t.the_geom)),
                    y(startpoint(t.the_geom)),
                    c.the_geom 
             from chat_cells3 c,
                  test_junk t
             where c.the_geom && t.the_geom;"

The output shows lines on the boundary of cells, with all the cells being
identified which touch or contain the line (point). A line(point) on the vertex
of four neighbouring cells gets 4 records returned.


 event_key | cell_id |     x     |    y     |                                  
     the_geom
-----------+---------+-----------+----------+----------------------------------------------------------------------------------------
  13683355 |    3553 |     175.5 |    -44.2 | SRID=4326;POLYGON((175.45
-44.25,175.45 -44.2,175.5 -44.2,175.5 -44.25,175.45 -44.25))
  13683355 |    3554 |     175.5 |    -44.2 | SRID=4326;POLYGON((175.45
-44.2,175.45 -44.15,175.5 -44.15,175.5 -44.2,175.45 -44.2))
  13683355 |    3625 |     175.5 |    -44.2 | SRID=4326;POLYGON((175.5
-44.25,175.5 -44.2,175.55 -44.2,175.55 -44.25,175.5 -44.25))
  13683355 |    3626 |     175.5 |    -44.2 | SRID=4326;POLYGON((175.5
-44.2,175.5 -44.15,175.55 -44.15,175.55 -44.2,175.5 -44.2))
  13677419 |    8186 | 178.68333 | -42.9833 | SRID=4326;POLYGON((178.65
-43,178.65 -42.95,178.7 -42.95,178.7 -43,178.65 -43))
  12583686 |    2910 | 175.01667 | -43.9833 | SRID=4326;POLYGON((175 -44,175
-43.95,175.05 -43.95,175.05 -44,175 -44))
  12578235 |    3198 | 175.23333 |   -43.95 | SRID=4326;POLYGON((175.2
-44,175.2 -43.95,175.25 -43.95,175.25 -44,175.2 -44))
  12578235 |    3199 | 175.23333 |   -43.95 | SRID=4326;POLYGON((175.2
-43.95,175.2 -43.9,175.25 -43.9,175.25 -43.95,175.2 -43.95))
  12497906 |   12481 |    181.65 |    -44.2 | SRID=4326;POLYGON((181.65
-44.25,181.65 -44.2,181.7 -44.2,181.7 -44.25,181.65 -44.25))
  12497906 |   12482 |    181.65 |    -44.2 | SRID=4326;POLYGON((181.65
-44.2,181.65 -44.15,181.7 -44.15,181.7 -44.2,181.65 -44.2))
  12497906 |   12409 |    181.65 |    -44.2 | SRID=4326;POLYGON((181.6
-44.25,181.6 -44.2,181.65 -44.2,181.65 -44.25,181.6 -44.25))
  12497906 |   12410 |    181.65 |    -44.2 | SRID=4326;POLYGON((181.6
-44.2,181.6 -44.15,181.65 -44.15,181.65 -44.2,181.6 -44.2))
  12497919 |    3270 | 175.26667 |   -43.95 | SRID=4326;POLYGON((175.25
-44,175.25 -43.95,175.3 -43.95,175.3 -44,175.25 -44))
  12497919 |    3271 | 175.26667 |   -43.95 | SRID=4326;POLYGON((175.25
-43.95,175.25 -43.9,175.3 -43.9,175.3 -43.95,175.25 -43.95))
  12497936 |    2841 | 174.98333 | -43.8333 | SRID=4326;POLYGON((174.95
-43.85,174.95 -43.8,175 -43.8,175 -43.85,174.95 -43.85))
  12497956 |    4856 | 176.38333 | -43.8666 | SRID=4326;POLYGON((176.35
-43.9,176.35 -43.85,176.4 -43.85,176.4 -43.9,176.35 -43.9))
(16 rows)

the second query uses not disjoint, which with this data, I understand should
give the same result. In fact only 2 records are returned. Same if I use
intersects, relate, add in an "or touches" etc.

psql $DB -c "select event_key,
                    c.cell_id,
                    x(startpoint(t.the_geom)),
                    y(startpoint(t.the_geom)),
                    c.the_geom 
             from chat_cells3 c,
                  test_junk t
             where c.the_geom && t.the_geom
               and not disjoint(t.the_geom, c.the_geom);"



 event_key | cell_id |   x    |   y   |                                      
the_geom
-----------+---------+--------+-------+---------------------------------------------------------------------------------------
  13683355 |    3626 |  175.5 | -44.2 | SRID=4326;POLYGON((175.5 -44.2,175.5
-44.15,175.55 -44.15,175.55 -44.2,175.5 -44.2))
  12497906 |   12482 | 181.65 | -44.2 | SRID=4326;POLYGON((181.65 -44.2,181.65
-44.15,181.7 -44.15,181.7 -44.2,181.65 -44.2))
(2 rows)


I can't see why I get two different results???



I guess the simple fix is to move the (or add a new) endpoint some tiny
distance from the startpoint to create a line of non-zero length.

Is there a (relatively) simple way to do this in PostGIS v0.8?


Thanks again,

  Brent Wood





> Brent Wood wrote:
> 
> >I'm trying to make sense of the various spatial operators. I need to join a
> >table of cells to a table of lines, to create a table where there is a
> record
> >for each case where a line passes entirely or partially through a cell.
> >
> >Looking at the OGC specs was a bit of help, but can anyone explain the
> >difference between the intersects, overlaps & crosses functions:
> >
> >ie:
> >
> >select 
> >.... 
> >where
> >  lines.geom && cell.geom and
> ><one of>
> >  intersects(lines.geom, cell.geom)
> >  overlaps(lines.geom, cell.geom)
> >  crosses(lines.geom, cells.geom)
> >
> >I figure I should be using crosses for this query, but don't undertand the
> >distinction between these.  
> >  
> >
> L -   references the Line geometry
> A -  references the cell, Polygon (area) geometry
> 
> intersects (L,A)  -  line and polygon (cell) touch or overlay in any way
> overlaps (L,A) - is not defined for a line vs polygon  test
> crosses (L,A) - line is part in and part out of the polygon, but not 
> wholly inside
> within (L,A) - the line is wholly inside the polygon
> contains (A,L) - the line is wholly inside the polygon
> 
> you may want 
>     intersects (lines.geom, cells.geom)
> or if you want to ignore the touches only relationship
>    ( crosses (lines.geom, cells.geom) or within(lines.geom, cells.geom) )
> or
>    ( intersects (lines.geom, cells.geom) and not touches(lines.geom, 
> cells.geom) )
> 
> 
> C.
> 
> >
> >Thanks,
> >
> >  Brent Wood
> >_______________________________________________
> >postgis-users mailing list
> >postgis-users at postgis.refractions.net
> >http://postgis.refractions.net/mailman/listinfo/postgis-users
> >  
> >
> 
> _______________________________________________
> postgis-users mailing list
> postgis-users at postgis.refractions.net
> http://postgis.refractions.net/mailman/listinfo/postgis-users
> 



More information about the postgis-users mailing list