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

[Graeme Leeming]

Greg,

I can't answer why your query is crashing, but I can suggest something 
else to attempt if you need to make it run faster.

What about trying to make the query more efficient by only using the big 
ugly US polygon where really necessary?  BBoxes in PostGIS make spatial 
queries a lot faster and more efficient, so you could try taking 
advantage of that concept by making your own simple bounding polygons.
- First create a chunky polygon that is fully inside the continental US 
and captures as much area inside as possible (even just a rectangle 
would help)
- Then create another chunky shape, a multi-hole donut polygon, that 
fully excludes all US areas (Continental, Hawaii, Alaska, any other US 
islands, etc.).
- Now for the queries: if the point falls in the first polygon, it is in 
the US so add to us_pnts.  Else if the point falls in the second polygon 
it is not in the US, so ignore it.  Both of these queries should be 
really quick if the polygons you created don't have an excessive number 
of points.
- For anything remaining, you're left with a near-border edge case and 
you can just run these against the full US polygon as originally 
planned.  The fewer the cases in this category, the better.

If this works, you could then edit and tune both polygons based on your 
expected request distribution.  The goal is to make the appropriate 
polygons capture as many high frequency areas near the border as 
possible (e.g. New York, L.A., Detroit, Vancouver (for the excluding 
one), etc.).  This will further reduce the need to use the detailed query.

Graeme Leeming
Refractions Research Inc.

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