[postgis-users] Determining if a point is within a polygon (oddfailures)
Gregory S. Williamson
gsw at globexplorer.com
Thu Sep 15 15:02:26 PDT 2005
I'll give this a try -- my earliest attempts actually were using b-boxes but for various reasons I ended up with a brute force approach.
ps is anyone has any ideas about why my single point at a time failed I'd be interested ...
From: postgis-users-bounces at postgis.refractions.net on behalf of Graeme Leeming
Sent: Thu 9/15/2005 2:44 PM
To: PostGIS Users Discussion
Subject: Re: [postgis-users] Determining if a point is within a polygon (oddfailures)
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
- Then create another chunky shape, a multi-hole donut polygon, that
fully excludes all US areas (Continental, Hawaii, Alaska, any other US
- 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
- 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.
Refractions Research Inc.
Gregory S. Williamson wrote:
>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 |
> "req_pnts_x_ndx" btree (y)
> "req_pnts_y_ndx" btree (x)
> "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:
>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);
> SRID=4326;POINT(-73.863730174705 41.185996824525 0.303293077013384)
>cmnd_reqs=# SELECT * FROM check_a_point('01010000A0E61000004A3BED5A477752C0533F73BECE974440600BB85D2769D33F'::geometry);
>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!
>ps -- this is my function:
>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;
> us_geom := (SELECT the_geom FROM us_shape WHERE fips_cntry = ''US'');
> pnt_is_us := within(pnt_2_test, us_geom);
>' LANGUAGE 'plpgsql';
>(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
postgis-users mailing list
postgis-users at postgis.refractions.net
More information about the postgis-users