[postgis-users] Determining if a point is within a polygon (odd failures)

[strk at refractions.net]

For Point-in-Polygon try using distance(point, polygon) == 0
rather then within(). It should be very fast.

As for crashing, if the problem persists with GEOS-2.1.4 I'll
take a look at it.

--strk;

On Thu, Sep 15, 2005 at 02:03:52PM -0700, Gregory S. Williamson wrote:
> Dear peoples,
> 
> This is postgres 8.0.3, running on a linux box with 2 CPUs and 2 gigs of RAM. The postGIS version says:
>  POSTGIS="1.0.1" GEOS="2.1.1" PROJ="Rel. 4.4.9, 29 Oct 2004" USE_STATS DBPROC="0.3.0" RELPROC="0.3.0"
> 
> I've got a new assignment to process data on centerpoints of requests for images on our website. I've got a procedure that can process a day's logs down to the queries of interest (say ~350000 per day) into a table:
> cmnd_reqs=# \d cmnd_reqs
>     Table "public.cmnd_reqs"
>   Column  |   Type   | Modifiers
> ----------+----------+-----------
>  y        | numeric  | not null
>  x        | numeric  | not null
>  z        | numeric  | not null
>  req_pnts | geometry |
>  req_date | date     |
> Indexes:
>     "req_pnts_x_ndx" btree (y)
>     "req_pnts_y_ndx" btree (x)
> Check constraints:
>     "enforce_srid_req_pnts" CHECK (srid(req_pnts) = 4326)
>     "enforce_dims_req_pnts" CHECK (ndims(req_pnts) = 3)
>     "enforce_geotype_req_pnts" CHECK (geometrytype(req_pnts) = 'POINT'::text OR req_pnts IS NULL)
> 
> I create the req_points geometry out of the x,y and z numbers and then the fun begins. I want to createe two output tables, one for points outside of the US and one for points in the US.
> 
> Trying a bulk statement such as "CREATE TABLE us_pnts AS SELECT * FROM cmnd_reqs WHERE within(req_pnts, <big ugly geometry for the US>);" ran for a few hours and then died, telling me that postgres had been restarted.
> 
> Rats. Ok. So I made simple function to teest a point and return the result, and a perl script that iteerates through all of the points one at a time and calls the function, and then writes out the data to the relevant destinations, doing a commit every 1000 rows.
> 
> This ran for several hours and was almost done. There were 365010 rows in the cmnd_reqs table and I processed 306946 points from the US and 15054 others before this blew:
> 
> .........................................commit 1000
> NOTICE:  AssertionFailedException: found single null side
> Can't execute sql to test point! ERROR:  GEOS within() threw an error!
> SQL is <SELECT * FROM check_a_point('01010000A0E61000004A3BED5A477752C0533F73BECE974440600BB85D2769D33F'::geometry);> at ./req_proc.pl line 53.
> 
> Curious, I tried examining the point in question:
> cmnd_reqs=# select asEWKT('01010000A0E61000004A3BED5A477752C0533F73BECE974440600BB85D2769D33F'::geometry);
>                                asewkt
> ---------------------------------------------------------------------
>  SRID=4326;POINT(-73.863730174705 41.185996824525 0.303293077013384)
> (1 row)
> 
> cmnd_reqs=# SELECT * FROM check_a_point('01010000A0E61000004A3BED5A477752C0533F73BECE974440600BB85D2769D33F'::geometry);
>  check_a_point
> ---------------
>  t
> (1 row)
> 
> So it works from psql but not from perl ? Or perhaps I've run into some memory issue or ? What, I wonder ?
> 
> Any advice would be most welcome as I am marooned in a land of frutration and angst!
> 
> Thanks
> 
> Greg Williamson
> DBA
> GlobeXplorer LLC
> 
> ps -- this is my function:
> BEGIN;
> CREATE OR REPLACE FUNCTION check_a_point(GEOMETRY)
>         RETURNS BOOLEAN AS '
> DECLARE pnt_2_test ALIAS FOR $1;
>         us_geom GEOMETRY;
>         pnt_is_us BOOLEAN;
> BEGIN
>         us_geom := (SELECT the_geom FROM us_shape WHERE fips_cntry = ''US'');
>         pnt_is_us := within(pnt_2_test, us_geom);
>         RETURN(pnt_is_us);
> END;
> ' LANGUAGE 'plpgsql';
> COMMIT;
> 
> (yes, I know if would be best to not retrieve the US shape every time but the description of the US outline is so long I haven't been able to make it into a constant!)
> 
> 
> _______________________________________________
> postgis-users mailing list
> postgis-users at postgis.refractions.net
> http://postgis.refractions.net/mailman/listinfo/postgis-users