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

[Gregory S. Williamson]

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!)