Thread: BUG #4967: polygon @> point incorrect for points close to shared borders

The following bug has been logged online:

Bug reference:      4967
Logged by:          Paul Matthews
Email address:      plm@netspace.net.au
PostgreSQL version: 8.3.7
Operating system:   Linux Open Suse 11.0 + 11.1
Description:        polygon @> point incorrect for points close to shared
borders
Details:

Have two polygons, both with many vertex, sharing a common edge. Several
thousand points where then tested to see which of the polygons the points
fell into using the "polygon @> point" operator.

A number of points close to the common border claimed they fell into both of
the polygons.

A quick perl+DBI+GD application was developed to plot both the polygons, the
polygon boundaries, as well as the points that thought they where in both.

This showed points close, but still several pixels away from the shared
border, thinking they where in both.

Guidance on this matter would be appreciated.

"Paul Matthews" <plm@netspace.net.au> writes:
> A number of points close to the common border claimed they fell into both of
> the polygons.

How close is "close"?  There's some pretty arbitrary fuzzy-comparisons
logic in the geometric datatypes ... see FPeq() and friends.  That might
be doing it to you.

            regards, tom lane

Tom Lane wrote:

  "Paul Matthews" <plm@netspace.net.au> writes:


    A number of points close to the common border claimed they fell into both of
the polygons.



How close is "close"?  There's some pretty arbitrary fuzzy-comparisons
logic in the geometric datatypes ... see FPeq() and friends.  That might
be doing it to you.

            regards, tom lane


I'll try to figure out how "relatively" close tonight, this stuff is
sub-metre resolution GPS data. The attached picture shows the two
polygons, the shared border, a road in this case, and the houses that
think they are on both sides of the road. Houses and other features are
located with latitude+longitude.

Last night I plugged in the polygon contains point code from
http://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html.
This resolved to houses correctly. If it helps in anyway please see the
attached. No use of fuzziness. Opaque yes, fuzzy no. <span
 class="moz-smiley-s1"> :-)  . Use in any way you
see fit.

#include "postgres.h"
#include "utils/geo_decls.h"
#include "fmgr.h"

#ifdef PG_MODULE_MAGIC
PG_MODULE_MAGIC;
#endif

PG_FUNCTION_INFO_V1(kontains);

Datum
kontains(PG_FUNCTION_ARGS)
{
  POLYGON* polygon;
  Point*   point;
  int      isin;
 
  polygon = PG_GETARG_POLYGON_P(0);
  point   = PG_GETARG_POINT_P(1);
  isin    = contains( polygon->npts, polygon->p, point );
 
  PG_RETURN_BOOL(isin);
}

int contains( int nvert, Point* vertex, Point* test )
{
  int i, j, c = 0;
  for( i=0, j=nvert-1; i<nvert; j=i++ ) {
    if( ((vertex[i].y>test->y) != (vertex[j].y>test->y))
&&
     (test->x < (vertex[j].x-vertex[i].x) *
(test->y-vertex[i].y) /
         (vertex[j].y-vertex[i].y) + vertex[i].x) )
      c = !c;
  }
  return c;
}

Paul Matthews <plm@netspace.net.au> writes:
>> How close is "close"?  There's some pretty arbitrary fuzzy-comparisons
>> logic in the geometric datatypes ... see FPeq() and friends.  That might
>> be doing it to you.

> I'll try to figure out how "relatively" close tonight, this stuff is
> sub-metre resolution GPS data. The attached picture shows the two
> polygons, the shared border, a road in this case, and the houses that
> think they are on both sides of the road. Houses and other features are
> located with latitude+longitude.<br>

Hmm ... just out of curiosity, why aren't you using PostGIS for this?
Our built-in geometric types are mostly an academic proof-of-concept,
they aren't industrial strength IMHO.

            regards, tom lane