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

Carl Anderson carl.anderson at vadose.org
Wed Apr 6 15:04:40 PDT 2005


Brent Wood wrote:

>--- 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. 
>
>  
>
from 099-049 SF-SQL in 2.1.13.3:

Expressed in terms of the DE-9IM:
    a.Disjoint(b) <=> (I(a) n I(b) = 0) or (I(a) n B(b) = 0) or (B(a) n 
I(b) = 0) and (B(a) n B(b) = 0) <=> a.Relate(b,  FF*FF**** )

in my words:
objects are disjoint if
    the interiors (special keyword) of the objects do not intersect
    the interior and boundaries of the objects, in either order do not 
intersect
    the boundaries of the objects do not intersect

For a Linestring the the Boundary is the startpoint and endpoint. Except 
when the startpoint is the same as the endpoint, then the boundary is null
      select boundary('LINESTRING(1 1,1 1)');
              boundary
     ---------------------------
      SRID=-1;MULTIPOINT(EMPTY)

For a Linestring the Interior is the interval (curve) from the 
startpoint to the endpoint

So to my understanding a zero length linestring has a null interior and 
a null boundary
Therefore the disjoint test cannot pass and cannot fail.  A bit like the 
equality operator versus any NULL
select '' = NULL    (answer is NULL)
select '' != NULL   (answer is NULL)
select ( '' = NULL ) is null    (answer is true)

(I looked for but could not find an assertion on Linestrings that 
required the interval to exist.  The document
that would assert that OGC 96-015r1 does not seem to be available anymore).


Why does the BBOX operator && give results that include the 0 length 
linestrings

        select getbbox('LINESTRING(0 0,0 0)'::geometry);
           getbbox
        --------------
         BOX(0 0,0 0)

the && (overlaping) operator for boundingboxes uses left of, right of, 
below, and above tests
and does make use of the interior, boundary, and exterior concepts.

So the two queries are apples and oranges in this case and are testing 
different things.

C.

>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
>>
>>    
>>
>_______________________________________________
>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