[postgis-users] Grid of points inside Polygon

Pierre Racine Pierre.Racine at sbf.ulaval.ca
Wed Nov 13 08:49:28 PST 2013


You can use ST_PixelAsCentroids(ST_AsRaster()) to generate a regular grid of points inside a polygon.

Pierre

> -----Original Message-----
> From: postgis-users-bounces at lists.osgeo.org [mailto:postgis-users-
> bounces at lists.osgeo.org] On Behalf Of Rémi Cura
> Sent: Wednesday, November 13, 2013 6:41 AM
> To: PostGIS Users Discussion
> Subject: Re: [postgis-users] Grid of points inside Polygon
> 
> This would be an improvement but still non efficient.
> 
> you have 2 possibilities, supposing that what you want is points 20 meter
> spaced in all road_buf  :
> 
> either you compute for each road_buffer the points inside, one road at a
> time ( figuratively ).
> 
> 
> 	This means you write a function which generates points 20-meter
> spaced in the bounding box of the road you are working on, and keep those
> in the real road buffer, and group result points as a multipoints (for
> cosmetic purpose)
> 	.
> 
> 
> 	You would then use it like this :
> 
> 	SELECT road_id, points_insidea_road(the_geom) AS
> my_points_inside_the_road
> 	FROM road_polygons_table
> 
> 
> 
> 	You would have as output a line per road with a multipoint
> containing all the point 20 meter spaced inside the road.
> 
> 
> 
> Or (what you wrote) you generate all points and keep those intersecting
> one or many road.
> 
> 
> The first one is mandatory, because it avoids to manipulate (incredibly) big
> table of all points spaced by 20 meters for UK (around 500 * 10^6 points ! )
> Even with indexes it's not a good idea to use such number of rows.
> 
> 
> That's the first point (write a function working for one road, then use it for
> all road).
> 
> 
> The second point is the way you compute is very inefficient. If your road is
> going from south-West to NorthEast, you will generate a very big number of
> points, and very few will be in the road_buffer. This is problematic as for a
> road of only 20 kms, you may generate as many as 100k points and keep
> only few of them. If you want to process hundreds of ks or roads it will
> become very problematic. Also you would have to generate points each
> time.
> 
> 
> So here is what  I suggest you : change your strategy :
>  instead of generating all point in bounding box and then keeping only those
> in road_buffer,
> generate a line every 20 meters going North south and an line every 20
> meters going East-West , then use the function ST_Intersection to keep only
> part of this lines being inside the road_polygon, then you have the points
> inside road_polygons as the intersections of these EW lines with the SN
> lines.
> 
> 
> It will be very efficient because you can create a table with all the lines going
> East-West and South-North for great britain (about 25k + 50k lines), and
> build index on it (index on geom and on the column saying if it is SN or EW).
> 
> 
> The trick is the number of lines is around 500km * 50 line/km + 1000km *
> 50 line/km , where the number of points is  500km * 50 line/km * 1000km
> * 50 line/km
> 
> 
> Hope it helps,
> 
> 
> Cheers,
> Rémi-C
> 
> 
> 
> 
> 
> 
> 2013/11/13 James David Smith <james.david.smith at gmail.com>
> 
> 
> 	Hey Remi,
> 
> 	Thanks for your reply. So in your mind you think we should have a
> 	database of say 300 polygons, and then we run a command like this
> 	right?
> 
> 	SELECT
> 	ST_Collect(st_setsrid(ST_POINT(x,y),27700))
> 	FROM
> 	generate_series(53320::int, 667380::int,20) as x,
> 	generate_series(7780::int, 1226580::int,20) as y,
> 	road_polygons_table
> 	WHERE
> 	st_intersects(road_polygons_table.the_geom,
> st_setsrid(ST_POINT(x,y),27700))
> 
> 	What do you think?
> 
> 	Thanks
> 
> 	James
> 
> 
> 
> 
> 
> 	On 11 November 2013 14:51, Rémi Cura <remi.cura at gmail.com>
> wrote:
> 	> Hey,
> 	> the whole point on using a sgbds like postgis is using index.
> 	>
> 	> If you have one line you don't use indexes...
> 	>
> 	> So in short, don't make one polygon with a buffer of all the road,
> but a
> 	> table with a line for the buffer for every road, then do you
> computation to
> 	> create grid of points inside of polygons, then union the result of
> points!
> 	>
> 	> And it s always a bad idea to run a function on big data when you
> have not
> 	> tested it fully (including scaling behavior) on small data.
> 	>
> 	>
> 	> Cheers
> 	> Rémi-C
> 	>
> 	>
> 	> 2013/11/11 James David Smith <james.david.smith at gmail.com>
> 	>>
> 	>> Hi all,
> 	>>
> 	>> Would appreciate some advice on the best way to accomplish
> this please.
> 	>>
> 	>> Our situation is that we have a single polygon which has been
> created
> 	>> by buffering all of the major roads in the UK. Projection is
> OSGB36
> 	>> (27700). Obviously it's quite a big polygon.
> 	>>
> 	>> -->  SELECT st_area(geom) FROM roadbufferunion;
> 	>>      st_area
> 	>> ------------------
> 	>>  77228753220.8271
> 	>>
> 	>> What we now want to do is create a regular grid of 20 metre x 20
> metre
> 	>> points instead the polygon area. So we wrote this function
> (based on
> 	>> some googling, apologies for not being able to recall the exact
> person
> 	>> who originally wrote it):
> 	>>
> 	>> CREATE OR REPLACE FUNCTION makegrid(geometry, integer,
> integer)
> 	>> RETURNS geometry AS
> 	>> 'SELECT ST_Collect(st_setsrid(ST_POINT(x,y),$3)) FROM
> 	>>   generate_series(53320::int, 667380::int,$2) as x
> 	>>   ,generate_series(7780::int, 1226580::int,$2) as y
> 	>> where st_intersects($1,st_setsrid(ST_POINT(x,y),$3))'
> 	>> LANGUAGE sql
> 	>>
> 	>> and we then run this by doing the following:
> 	>>
> 	>> SELECT st_x((ST_Dump(makegrid(geom, 20, 27700))).geom) as x,
> 	>> st_y((ST_Dump(makegrid(geom, 20, 27700))).geom) as y INTO
> grid_points
> 	>> from roadbufferunion;
> 	>>
> 	>> However after over 2 days of the query running on a pretty
> powerful
> 	>> linux cluster, we still have no result.  I'm not sure if it is
> 	>> actually running or not to be honest.
> 	>>
> 	>> Does the query look right?
> 	>> Any ideas how we can make it run quicker?
> 	>>
> 	>> Thanks
> 	>>
> 	>> James
> 	>> _______________________________________________
> 	>> postgis-users mailing list
> 	>> postgis-users at lists.osgeo.org
> 	>> http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
> 	>
> 	>
> 	>
> 	> _______________________________________________
> 	> postgis-users mailing list
> 	> postgis-users at lists.osgeo.org
> 	> http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
> 	_______________________________________________
> 	postgis-users mailing list
> 	postgis-users at lists.osgeo.org
> 	http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
> 
> 



More information about the postgis-users mailing list