[postgis-users] Determining if a point is within apolygon (oddfailures)
Gregory S. Williamson
gsw at globexplorer.com
Thu Sep 15 18:39:03 PDT 2005
Stephen --
Thanks for the suggestion (I was the original poster). I tried Graeme's suggestion and it seems to work reasonably quickly and didn't suffer any unexpected failures.
Still not sure why the original error occurred, but I seem to have a work-around.
Thanks, all.
Greg
-----Original Message-----
From: postgis-users-bounces at postgis.refractions.net on behalf of Stephen Woodbridge
Sent: Thu 9/15/2005 5:00 PM
To: PostGIS Users Discussion
Cc:
Subject: Re: [postgis-users] Determining if a point is within apolygon (oddfailures)
Graeme,
Try the PostGIS Simplify(geom, tolerance) function to create a
simplified version of your big US polygon. Try various tolerance values
until you get the simplest geometry that you can live with. Then
manipulate it to simplify some bays and coastal regions.
-Steve
Gregory S. Williamson wrote:
> Graeme --
>
> 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.
>
> Thanks,
>
> Greg
>
> ps is anyone has any ideas about why my single point at a time failed I'd be interested ...
>
>
> -----Original Message-----
> 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
> Cc:
> Subject: Re: [postgis-users] Determining if a point is within a polygon (oddfailures)
> 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
>>
>>
>
>
> _______________________________________________
> postgis-users mailing list
> postgis-users at postgis.refractions.net
> http://postgis.refractions.net/mailman/listinfo/postgis-users
>
>
>
>
>
>
> _______________________________________________
> postgis-users mailing list
> postgis-users at postgis.refractions.net
> http://postgis.refractions.net/mailman/listinfo/postgis-users
>
_______________________________________________
postgis-users mailing list
postgis-users at postgis.refractions.net
http://postgis.refractions.net/mailman/listinfo/postgis-users
!DSPAM:432a09c6323101561020452!
More information about the postgis-users
mailing list