Hi Remi, Apologies. One more thing which I am not sure about. We need the final grid of points that we make to be multiples of 20. For example:
53320, 7780 = Yes 53323, 7780 = No 53320, 7794 = No 53321, 7754 = No That is why in the original function that we wrote, we put the numbers in as below. 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))' We did this because the points in the grid have to 'align' with the values of mutliple 20. If we do the N-S and E-W lines solution that you suggest, I don't think that this will work will it? Thanks James On 13 November 2013 11:49, James David Smith <[email protected]> wrote: > Hi Remi, > > Thanks so much for this detailed response. Your idea about creating > the lines and the only storing where they intersect and are within the > polygons is a great idea. I'm going to give that a go now. > > Thanks again, > > James > > On 13 November 2013 11:41, Rémi Cura <[email protected]> wrote: >> 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 <[email protected]> >>> >>> 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 <[email protected]> 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 <[email protected]> >>> >> >>> >> 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 >>> >> [email protected] >>> >> http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users >>> > >>> > >>> > >>> > _______________________________________________ >>> > postgis-users mailing list >>> > [email protected] >>> > http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users >>> _______________________________________________ >>> postgis-users mailing list >>> [email protected] >>> http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users >> >> >> >> _______________________________________________ >> postgis-users mailing list >> [email protected] >> http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users _______________________________________________ postgis-users mailing list [email protected] http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
