[postgis-users] Grid of points inside Polygon

[James David Smith]

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
<james.david.smith at gmail.com> 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 <remi.cura at gmail.com> 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 <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
>>
>>
>>
>> _______________________________________________
>> postgis-users mailing list
>> postgis-users at lists.osgeo.org
>> http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users