Re: [postgis-users] Grid of points inside Polygon

[James David Smith]

Hi Remi / all,

I am just coming back to this project. I have made some progress. But
would appreciate advice on how best to proceed now. What I have now is
a table with two columns, 'st_x' and 'st_y'. They are the 20 metre by
20 metre grid points that are within 1km of main roads in the UK. The
table has 499,677,666 rows!  Because of the method that I used to
create this however, I have alot of duplication of points. So I think
I need to do something like:

CREATE TABLE grid_no_duplications AS (SELECT DISTINCT(x, y) FROM grid);

However I thought that first I could put an index on it?

CREATE INDEX grid_x_y ON grid (st_x, st_y);

Is this a good approach?

Cheers

James

On 15 November 2013 16:39, Rémi Cura remi.c...@gmail.com wrote:
 good =)

 It is better to do the buffer computation first,
 because this way you are sure youwon't compute the same buffer many times.

 Something like (not so sure without trying it)
 UPDATE ukmajrdbuffer SET geom = ST_Buffer(geom,1000);

 You may have an error if you hav e a fixed geom data type, then you could
 just add a new column like

 ALTER TABLE ukmajrdbuffer ADD COLUMN the_geom geometry ;
 UPDATE ukmajrdbuffer SET the_geom = ST_Buffer(geom,1000);


 Cheers,

 Rémi-C



 2013/11/15 James David Smith james.david.sm...@gmail.com

 Hey Remi,

 I've actually managed to get the file 'ukmajorroads' already and have
 loaded it into my database.  There are 395356 rows in the database. There is
 a field called 'geom' and I have built an index on it as below;

 CREATE INDEX ukrds_index ON ukrds USING GIST (geom);

 So should I now run the same query, but do the buffering of the roads 'on
 the fly' maybe?  For example as below?


 CREATE TABLE lines_for_each_road AS
 WITH all_lines AS (
 SELECT *
 FROM all_lines
 ),
 cutted_lines AS ( --we cut the line to keep only part of lines inside
 road_buffer
 SELECT ST_Intersection(all_lines.the_geom,st_buffer(ukmajrdbuffer.geom,
 1000)) as lines_cutted, direction
 FROM ukmajrdbuffer, all_lines
 WHERE ST_Intersects(st_buffer(ukmajrdbuffer.geom, 1000),
 all_lines.the_geom)=TRUE
 ),
 cutted_lines_SN AS ( --this is the cutted lines which going from South to
 North
 SELECT *
 FROM cutted_lines
 WHERE direction = 'SN'
 ),
 cutted_lines_EW AS ( --this is the cutted lines going from East toWest
 SELECT *
 FROM cutted_lines
 WHERE direction = 'EW'
 ),
 points AS ( -- we take the intersection of EW  lines with SN lines , that
 is the points on the grid.
 SELECT ST_Intersection(clSN.lines_cutted, clEW.lines_cutted) AS point
 FROM cutted_lines_SN as clSN, cutted_lines_EW AS clEW
 WHERE ST_Intersects(clSN.lines_cutted, clEW.lines_cutted)=TRUE --no point
 ot compute an intersection if lines don't intersect
 )
 SELECT row_number() over() AS id  , point
 FROM points ;


 On 15 November 2013 15:42, Rémi Cura remi.c...@gmail.com wrote:

 Still is a shame :
 with proper data we should have some result with about one hour i guess.

 Cheers,
 Rémi C


 2013/11/15 James David Smith james.david.sm...@gmail.com

 Remi,

 Ok. Cool. I've set the below query running. On Monday I will also
 attempt to get the original road lines file too. Let's see if we have
 a result on Monday.
 --
 CREATE TABLE lines_for_each_road AS
 WITH all_lines AS (
 SELECT *
 FROM all_lines
 ),
 cutted_lines AS ( --we cut the line to keep only part of lines inside
 road_buffer
 SELECT ST_Intersection(all_lines.the_geom,ukmajrdbuffer.geom) as
 lines_cutted, direction
 FROM ukmajrdbuffer, all_lines
 WHERE ST_Intersects(ukmajrdbuffer.geom, all_lines.the_geom)=TRUE
 ),
 cutted_lines_SN AS ( --this is the cutted lines which going from South
 to North
 SELECT *
 FROM cutted_lines
 WHERE direction = 'SN'
 ),
 cutted_lines_EW AS ( --this is the cutted lines going from East toWest
 SELECT *
 FROM cutted_lines
 WHERE direction = 'EW'
 ),
 points AS ( -- we take the intersection of EW  lines with SN lines ,
 that is the points on the grid.
 SELECT ST_Intersection(clSN.lines_cutted, clEW.lines_cutted) AS point
 FROM cutted_lines_SN as clSN, cutted_lines_EW AS clEW
 WHERE ST_Intersects(clSN.lines_cutted, clEW.lines_cutted)=TRUE --no
 point ot compute an intersection if lines don't intersect
 )
 SELECT row_number() over() AS id  , point
 FROM points ;
 --
 Thanks

 James

 On 15 November 2013 15:37, Rémi Cura remi.c...@gmail.com wrote:
 
  You should try to get the original road file, because i will be very
  easy to
  reapply some buffer and filtering to keep only major road.
  Yet, why not let it run during the week end?
 
  Cheers,
 
  Rémi-C
 
 
  2013/11/15 James David Smith james.david.sm...@gmail.com
 
  Hey Remi,
 
  Do you think before I try running the big query you have just sent
  me,
  that I should go back and try to get the original file of uk roads? I
  mean the very original file that has not had any buffers applied or
  any merging done.
 
  Or shall we just go for it and see if it has finished when I come
  into

Re: [postgis-users] Grid of points inside Polygon

[James David Smith]

Hi Remi,

Thanks for your continued guidance with this big data problem.  I
wasn't in the office yesterday so am looking at it again now. I've
just ran the below below query running, and will let you know how many
rows there are once it's finished.

CREATE TABLE lines_for_each_road AS
WITH all_lines AS (
SELECT *
FROM all_lines
)
SELECT row_number() over() AS id--,*
FROM ukmajrdbuffer, all_lines
WHERE ST_Intersects(ukmajrdbuffer.geom, all_lines.the_geom)=TRUE;

In the meantime however, I have noticed something else that is a
worry. I did this:

SELECT st_isvalid(geom) FROM ukmajrdbuffer;

NOTICE:  Ring Self-intersection at or near point 367541.09074241668
311890.25257437676
NOTICE:  Ring Self-intersection at or near point 209901.5000186
706890.5000373
NOTICE:  Ring Self-intersection at or near point 444961.46434461698
161217.66612910852
NOTICE:  Ring Self-intersection at or near point 559818.64931676909
103089.14606037363
NOTICE:  Ring Self-intersection at or near point 420114.82424666081
298567.38722311892
 st_isvalid

 f
 f
 t
 f
 f
 f
 t
 t
(8 rows)

I guess that I should try to fix these geometries first before we
attempt to continue any further right?

Thanks

James







On 13 November 2013 18:04, Rémi Cura remi.c...@gmail.com wrote:
 All lines should be in one table, you can add a column to differentiate
 between line going SN and EW,
 this is what the query I wrote does.

 Not changing your code much, it gives for the line table :

 CREATE TABLE all_lines tablespace jamestabs AS
 ( SELECT st_setsrid(st_makeline(st_makepoint(x, 7780), st_makepoint(x,
 1226580)),27700) as the_geom, 'SN' AS direction

 FROM generate_series(53320::int, 667380::int, 20) as x)
 UNION
 ( SELECT st_setsrid(st_makeline(st_makepoint(53320, y),
 st_makepoint(667380, y)),27700) as the_geom, 'EW' AS direction

 FROM generate_series(7780::int, 1226580::int, 20) as y)

 then building index on it.


 -- Merge the road buffer polygons into one big polygon
 -- Should maybe consider leaving them as they are, but am merging for now.
 SELECT st_union(geom) as geom INTO union_roads FROM ukmajrdbuffer;
 You must not do it !
 As I told you blabla index blabla.

 instead, create an index on table ukmajrdbuffer, and use it without
 unioning.

 Then you query is not going to work at all ... (or *very very* slowly)
 FROM east_to_west, north_to_south, buff_union_roads : it potentially
 generates 1000*50*1000*50 * number_of_road rows to filter, very bad idea,
 you'll get the same problem as with points...

 (please see some example on internet to understand WITH, it is simple ,
 useful, and makes query much more clear)

 First you could try :

 CREATE TABLE lines_for_each_road tablespace jamestabs AS
 WITH all_lines AS ( --this subquerry gets all the lines

 SELECT *
 FROM all_lines

 ),
 SELECT row_number() over() AS id--,*
 FROM ukmajrdbuffer, all_lines
 WHERE ST_Intersects(ukmajrdbuffer.the_geom, all_lines.the_geom)=TRUE;

 This will give you (if you uncomment the --,*, ), for every road, all the
 lines intersecting the road
 example :
 road_1 | line_2_EW
 road_1 | line_3_EW
 road_1 | line_78_SN
 road_1 | line_7_EW
 road_2 ...

 How long is it (don't uncomment plz), and how many rows do you get (select
 count(*) from lines_for_each_road )?

 If it is not too long, we will go on and cut the lines so that for every
 road, we keep the part of the lines that are inside the road_buffer.
 (I have to leave now, tomorrow)

 Cheers,
 Rémi-C





 2013/11/13 James David Smith james.david.sm...@gmail.com

 Hey Remi,

 I don't quite get what the query you gave does to be honest. But
 working with some of the ideas and notions you gave me I have put the
 below together. The final bit of the code is the key query. The tables
 of NS and EW lines were done very quickly. Nice. But now doing the
 'intersects within polygon' is taking some time so I'll see how it
 gets on when I get in to work tomorrow morning:

 DROP TABLE IF EXISTS north_to_south;
 DROP TABLE IF EXISTS east_to_west;

 -- Make a table of north to south lines
 CREATE TABLE north_to_south tablespace jamestabs AS
 ( SELECT st_setsrid(st_makeline(st_makepoint(x, 7780), st_makepoint(x,
 1226580)),27700) as the_geom
 FROM generate_series(53320::int, 667380::int, 20) as x);

 -- Add an index
 CREATE INDEX north_to_south_geom_index ON north_to_south USING GIST (
 the_geom );

 -- Make a table of east to west lines
 CREATE TABLE east_to_west tablespace jamestabs AS
 ( SELECT st_setsrid(st_makeline(st_makepoint(53320, y),
 st_makepoint(667380, y)),27700) as the_geom
 FROM generate_series(7780::int, 1226580::int, 20) as y);

 -- Add an index
 CREATE INDEX east_to_west_geom_index ON east_to_west USING GIST ( the_geom
 );

 -- Import the Road Buffer file using the pg_loader interface

 -- Update the geometry field as the loader didn't seem to do this.
 SELECT UpdateGeometrySRID('ukmajrdbuffer','geom',27700);

 -- Merge the road buffer polygons into one big polygon
 -- Should maybe 

Re: [postgis-users] Grid of points inside Polygon

[Rémi Cura]

Yep, maybe something like

UPDATE ukmajrdbuffer SET the_geom = ST_MakeValid(the_geom)
WHERE ST_IsValid(the_geom) = FALSE


Cheers,
Rémi-C


2013/11/15 James David Smith james.david.sm...@gmail.com

 Hi Remi,

 Thanks for your continued guidance with this big data problem.  I
 wasn't in the office yesterday so am looking at it again now. I've
 just ran the below below query running, and will let you know how many
 rows there are once it's finished.

 CREATE TABLE lines_for_each_road AS
 WITH all_lines AS (
 SELECT *
 FROM all_lines
 )
 SELECT row_number() over() AS id--,*
 FROM ukmajrdbuffer, all_lines
 WHERE ST_Intersects(ukmajrdbuffer.geom, all_lines.the_geom)=TRUE;

 In the meantime however, I have noticed something else that is a
 worry. I did this:

 SELECT st_isvalid(geom) FROM ukmajrdbuffer;

 NOTICE:  Ring Self-intersection at or near point 367541.09074241668
 311890.25257437676
 NOTICE:  Ring Self-intersection at or near point 209901.5000186
 706890.5000373
 NOTICE:  Ring Self-intersection at or near point 444961.46434461698
 161217.66612910852
 NOTICE:  Ring Self-intersection at or near point 559818.64931676909
 103089.14606037363
 NOTICE:  Ring Self-intersection at or near point 420114.82424666081
 298567.38722311892
  st_isvalid
 
  f
  f
  t
  f
  f
  f
  t
  t
 (8 rows)

 I guess that I should try to fix these geometries first before we
 attempt to continue any further right?

 Thanks

 James







 On 13 November 2013 18:04, Rémi Cura remi.c...@gmail.com wrote:
  All lines should be in one table, you can add a column to differentiate
  between line going SN and EW,
  this is what the query I wrote does.
 
  Not changing your code much, it gives for the line table :
 
  CREATE TABLE all_lines tablespace jamestabs AS
  ( SELECT st_setsrid(st_makeline(st_makepoint(x, 7780), st_makepoint(x,
  1226580)),27700) as the_geom, 'SN' AS direction
 
  FROM generate_series(53320::int, 667380::int, 20) as x)
  UNION
  ( SELECT st_setsrid(st_makeline(st_makepoint(53320, y),
  st_makepoint(667380, y)),27700) as the_geom, 'EW' AS direction
 
  FROM generate_series(7780::int, 1226580::int, 20) as y)
 
  then building index on it.
 
 
  -- Merge the road buffer polygons into one big polygon
  -- Should maybe consider leaving them as they are, but am merging for
 now.
  SELECT st_union(geom) as geom INTO union_roads FROM ukmajrdbuffer;
  You must not do it !
  As I told you blabla index blabla.
 
  instead, create an index on table ukmajrdbuffer, and use it without
  unioning.
 
  Then you query is not going to work at all ... (or *very very* slowly)
  FROM east_to_west, north_to_south, buff_union_roads : it potentially
  generates 1000*50*1000*50 * number_of_road rows to filter, very bad idea,
  you'll get the same problem as with points...
 
  (please see some example on internet to understand WITH, it is simple ,
  useful, and makes query much more clear)
 
  First you could try :
 
  CREATE TABLE lines_for_each_road tablespace jamestabs AS
  WITH all_lines AS ( --this subquerry gets all the lines
 
  SELECT *
  FROM all_lines
 
  ),
  SELECT row_number() over() AS id--,*
  FROM ukmajrdbuffer, all_lines
  WHERE ST_Intersects(ukmajrdbuffer.the_geom, all_lines.the_geom)=TRUE;
 
  This will give you (if you uncomment the --,*, ), for every road, all
 the
  lines intersecting the road
  example :
  road_1 | line_2_EW
  road_1 | line_3_EW
  road_1 | line_78_SN
  road_1 | line_7_EW
  road_2 ...
 
  How long is it (don't uncomment plz), and how many rows do you get
 (select
  count(*) from lines_for_each_road )?
 
  If it is not too long, we will go on and cut the lines so that for
 every
  road, we keep the part of the lines that are inside the road_buffer.
  (I have to leave now, tomorrow)
 
  Cheers,
  Rémi-C
 
 
 
 
 
  2013/11/13 James David Smith james.david.sm...@gmail.com
 
  Hey Remi,
 
  I don't quite get what the query you gave does to be honest. But
  working with some of the ideas and notions you gave me I have put the
  below together. The final bit of the code is the key query. The tables
  of NS and EW lines were done very quickly. Nice. But now doing the
  'intersects within polygon' is taking some time so I'll see how it
  gets on when I get in to work tomorrow morning:
 
  DROP TABLE IF EXISTS north_to_south;
  DROP TABLE IF EXISTS east_to_west;
 
  -- Make a table of north to south lines
  CREATE TABLE north_to_south tablespace jamestabs AS
  ( SELECT st_setsrid(st_makeline(st_makepoint(x, 7780), st_makepoint(x,
  1226580)),27700) as the_geom
  FROM generate_series(53320::int, 667380::int, 20) as x);
 
  -- Add an index
  CREATE INDEX north_to_south_geom_index ON north_to_south USING GIST (
  the_geom );
 
  -- Make a table of east to west lines
  CREATE TABLE east_to_west tablespace jamestabs AS
  ( SELECT st_setsrid(st_makeline(st_makepoint(53320, y),
  st_makepoint(667380, y)),27700) as the_geom
  FROM generate_series(7780::int, 1226580::int, 20) as y);
 
  -- Add an 

Re: [postgis-users] Grid of points inside Polygon

[Sandro Santilli]

On Fri, Nov 15, 2013 at 11:50:42AM +0100, Rémi Cura wrote:
 Yep, maybe something like
 
 UPDATE ukmajrdbuffer SET the_geom = ST_MakeValid(the_geom)
 WHERE ST_IsValid(the_geom) = FALSE

ST_MakeValid internally checks for ST_IsValid, so no need
to add the condition (which would run the test twice).

--strk;
___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

Re: [postgis-users] Grid of points inside Polygon

[James David Smith]

Thanks both. Geometries now fixed.

The query 'CREATE TABLE lines_for_each_road' has now been set
running. Will report back when it's done. I suspect it may take a
while!

James

On 15 November 2013 11:03, Sandro Santilli s...@keybit.net wrote:
 On Fri, Nov 15, 2013 at 11:50:42AM +0100, Rémi Cura wrote:
 Yep, maybe something like

 UPDATE ukmajrdbuffer SET the_geom = ST_MakeValid(the_geom)
 WHERE ST_IsValid(the_geom) = FALSE

 ST_MakeValid internally checks for ST_IsValid, so no need
 to add the condition (which would run the test twice).

 --strk;
 ___
 postgis-users mailing list
 postgis-users@lists.osgeo.org
 http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

Re: [postgis-users] Grid of points inside Polygon

[Rémi Cura]

It should be finished by now,
can you check you have geom indexes on :
ukmajrdbuffer.geom and  all_lines.the_geom


How many geoms do you have in ukmajrdbuffer?

Cheers,
Rémi-C


2013/11/15 Rémi Cura remi.c...@gmail.com

 Hey Sandro,

 Thanks for this, it is at least twice faster =)

 Cheers,
 Rémi-C




 2013/11/15 James David Smith james.david.sm...@gmail.com

 Thanks both. Geometries now fixed.

 The query 'CREATE TABLE lines_for_each_road' has now been set
 running. Will report back when it's done. I suspect it may take a
 while!

 James

 On 15 November 2013 11:03, Sandro Santilli s...@keybit.net wrote:
  On Fri, Nov 15, 2013 at 11:50:42AM +0100, Rémi Cura wrote:
  Yep, maybe something like
 
  UPDATE ukmajrdbuffer SET the_geom = ST_MakeValid(the_geom)
  WHERE ST_IsValid(the_geom) = FALSE
 
  ST_MakeValid internally checks for ST_IsValid, so no need
  to add the condition (which would run the test twice).
 
  --strk;
  ___
  postgis-users mailing list
  postgis-users@lists.osgeo.org
  http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
 ___
 postgis-users mailing list
 postgis-users@lists.osgeo.org
 http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users



___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

Re: [postgis-users] Grid of points inside Polygon

[James David Smith]

Hey.

Yes, it's done. Was just getting some lunch! :-)

select count(*) from lines_for_each_road
Result = 187033

I have also just ran 'VACUUM ANALYZE' on the tables
'lines_for_each_road' as well as the table 'all_lines'

I also can confirm that I have ran the following commands:

CREATE INDEX all_lines_index ON all_lines USING GIST ( the_geom )
CREATE INDEX ukmajrdbuffer_index ON ukmajrdbuffer USING GIST (geom);

Should we now uncomment this line from the previous query?

  SELECT row_number() over() AS id--,*  

Thanks again Remi,

James



On 15 November 2013 13:14, Rémi Cura remi.c...@gmail.com wrote:
 Also if you do have indexes,
 can you run a VACUUM ANALYZE, so that the indexes will be used?

 Cheers,

 Rémi-C


 2013/11/15 Rémi Cura remi.c...@gmail.com

 It should be finished by now,
 can you check you have geom indexes on :
 ukmajrdbuffer.geom and  all_lines.the_geom


 How many geoms do you have in ukmajrdbuffer?

 Cheers,
 Rémi-C


 2013/11/15 Rémi Cura remi.c...@gmail.com

 Hey Sandro,

 Thanks for this, it is at least twice faster =)

 Cheers,
 Rémi-C




 2013/11/15 James David Smith james.david.sm...@gmail.com

 Thanks both. Geometries now fixed.

 The query 'CREATE TABLE lines_for_each_road' has now been set
 running. Will report back when it's done. I suspect it may take a
 while!

 James

 On 15 November 2013 11:03, Sandro Santilli s...@keybit.net wrote:
  On Fri, Nov 15, 2013 at 11:50:42AM +0100, Rémi Cura wrote:
  Yep, maybe something like
 
  UPDATE ukmajrdbuffer SET the_geom = ST_MakeValid(the_geom)
  WHERE ST_IsValid(the_geom) = FALSE
 
  ST_MakeValid internally checks for ST_IsValid, so no need
  to add the condition (which would run the test twice).
 
  --strk;
  ___
  postgis-users mailing list
  postgis-users@lists.osgeo.org
  http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
 ___
 postgis-users mailing list
 postgis-users@lists.osgeo.org
 http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users





 ___
 postgis-users mailing list
 postgis-users@lists.osgeo.org
 http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

Re: [postgis-users] Grid of points inside Polygon

[James David Smith]

Hey Remi,

I'll do a few checks and get back to you. Maybe I did something wrong
because I set this query going on my local PostgreSQL installation but
also on our Linux Cluster machine which is much more powerful. And it
finished on my local installation BEFORE the cluster. So maybe I did
something wrong on the local query.

The query that has finished took about 1 hour.
The query on our cluster is still running.

select count(*) from ukmajrdbuffer = 8
This is because before I was given the data the roads buffer had
already been dissolved unfortunately.

James




On 15 November 2013 13:43, Rémi Cura remi.c...@gmail.com wrote:
 Good !

 Something is strange with the result,
 you get only 190k lines intersecting road buffer,
 it is very few , I expected at least 10times this!
 How many time took this computing by the way?

 How many road buffer geom do you have ? (select count(*) from ukmajrdbuffer
 ).

 Cheers,
 Rémi-C


 2013/11/15 James David Smith james.david.sm...@gmail.com

 Hey.

 Yes, it's done. Was just getting some lunch! :-)

 select count(*) from lines_for_each_road
 Result = 187033

 I have also just ran 'VACUUM ANALYZE' on the tables
 'lines_for_each_road' as well as the table 'all_lines'

 I also can confirm that I have ran the following commands:

 CREATE INDEX all_lines_index ON all_lines USING GIST ( the_geom )
 CREATE INDEX ukmajrdbuffer_index ON ukmajrdbuffer USING GIST (geom);

 Should we now uncomment this line from the previous query?

   SELECT row_number() over() AS id--,*  

 Thanks again Remi,

 James



 On 15 November 2013 13:14, Rémi Cura remi.c...@gmail.com wrote:
  Also if you do have indexes,
  can you run a VACUUM ANALYZE, so that the indexes will be used?
 
  Cheers,
 
  Rémi-C
 
 
  2013/11/15 Rémi Cura remi.c...@gmail.com
 
  It should be finished by now,
  can you check you have geom indexes on :
  ukmajrdbuffer.geom and  all_lines.the_geom
 
 
  How many geoms do you have in ukmajrdbuffer?
 
  Cheers,
  Rémi-C
 
 
  2013/11/15 Rémi Cura remi.c...@gmail.com
 
  Hey Sandro,
 
  Thanks for this, it is at least twice faster =)
 
  Cheers,
  Rémi-C
 
 
 
 
  2013/11/15 James David Smith james.david.sm...@gmail.com
 
  Thanks both. Geometries now fixed.
 
  The query 'CREATE TABLE lines_for_each_road' has now been set
  running. Will report back when it's done. I suspect it may take a
  while!
 
  James
 
  On 15 November 2013 11:03, Sandro Santilli s...@keybit.net wrote:
   On Fri, Nov 15, 2013 at 11:50:42AM +0100, Rémi Cura wrote:
   Yep, maybe something like
  
   UPDATE ukmajrdbuffer SET the_geom = ST_MakeValid(the_geom)
   WHERE ST_IsValid(the_geom) = FALSE
  
   ST_MakeValid internally checks for ST_IsValid, so no need
   to add the condition (which would run the test twice).
  
   --strk;
   ___
   postgis-users mailing list
   postgis-users@lists.osgeo.org
   http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
  ___
  postgis-users mailing list
  postgis-users@lists.osgeo.org
  http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
 
 
 
 
 
  ___
  postgis-users mailing list
  postgis-users@lists.osgeo.org
  http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
 ___
 postgis-users mailing list
 postgis-users@lists.osgeo.org
 http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users



 ___
 postgis-users mailing list
 postgis-users@lists.osgeo.org
 http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

Re: [postgis-users] Grid of points inside Polygon

[Rémi Cura]

Outch,
you have only 8 roads in GB? =)

This is not going to go fast till you don't have some dozens k's of roads.

I'm guessing you are not at will to send those roads?


The next step will be (when you'll be sure everything is OK)


CREATE TABLE lines_for_each_road AS
WITH all_lines AS (
SELECT *
FROM all_lines
),
cutted_lines AS ( --we cut the line to keep only part of lines inside
road_buffer
SELECT ST_Intersection(all_lines.the_geom,ukmajrdbuffer.geom) as
lines_cutted, direction
FROM ukmajrdbuffer, all_lines
WHERE ST_Intersects(ukmajrdbuffer.geom, all_lines.the_geom)=TRUE
),
cutted_lines_SN AS ( --this is the cutted lines which going from South to
North
SELECT *
FROM cutted_lines
WHERE direction = 'SN'
),
cutted_lines_EW AS ( --this is the cutted lines going from East toWest
SELECT *
FROM cutted_lines
WHERE direction = 'EW'
),
points AS ( -- we take the intersection of EW  lines with SN lines , that
is the points on the grid.
SELECT ST_Intersection(clSN.lines_cutted, clEW.lines_cutted) AS point
FROM cutted_lines_SN as clSN, cutted_lines_EW AS clEW
WHERE ST_Intersects(clSN.lines_cutted, clEW.lines_cutted)=TRUE --no point
ot compute an intersection if lines don't intersect
)
SELECT row_number() over() AS id  , point
FROM points ;





Cheers,

Rémi-C



2013/11/15 James David Smith james.david.sm...@gmail.com

 Hey Remi,

 I'll do a few checks and get back to you. Maybe I did something wrong
 because I set this query going on my local PostgreSQL installation but
 also on our Linux Cluster machine which is much more powerful. And it
 finished on my local installation BEFORE the cluster. So maybe I did
 something wrong on the local query.

 The query that has finished took about 1 hour.
 The query on our cluster is still running.

 select count(*) from ukmajrdbuffer = 8
 This is because before I was given the data the roads buffer had
 already been dissolved unfortunately.

 James




 On 15 November 2013 13:43, Rémi Cura remi.c...@gmail.com wrote:
  Good !
 
  Something is strange with the result,
  you get only 190k lines intersecting road buffer,
  it is very few , I expected at least 10times this!
  How many time took this computing by the way?
 
  How many road buffer geom do you have ? (select count(*) from
 ukmajrdbuffer
  ).
 
  Cheers,
  Rémi-C
 
 
  2013/11/15 James David Smith james.david.sm...@gmail.com
 
  Hey.
 
  Yes, it's done. Was just getting some lunch! :-)
 
  select count(*) from lines_for_each_road
  Result = 187033
 
  I have also just ran 'VACUUM ANALYZE' on the tables
  'lines_for_each_road' as well as the table 'all_lines'
 
  I also can confirm that I have ran the following commands:
 
  CREATE INDEX all_lines_index ON all_lines USING GIST ( the_geom )
  CREATE INDEX ukmajrdbuffer_index ON ukmajrdbuffer USING GIST (geom);
 
  Should we now uncomment this line from the previous query?
 
SELECT row_number() over() AS id--,*  
 
  Thanks again Remi,
 
  James
 
 
 
  On 15 November 2013 13:14, Rémi Cura remi.c...@gmail.com wrote:
   Also if you do have indexes,
   can you run a VACUUM ANALYZE, so that the indexes will be used?
  
   Cheers,
  
   Rémi-C
  
  
   2013/11/15 Rémi Cura remi.c...@gmail.com
  
   It should be finished by now,
   can you check you have geom indexes on :
   ukmajrdbuffer.geom and  all_lines.the_geom
  
  
   How many geoms do you have in ukmajrdbuffer?
  
   Cheers,
   Rémi-C
  
  
   2013/11/15 Rémi Cura remi.c...@gmail.com
  
   Hey Sandro,
  
   Thanks for this, it is at least twice faster =)
  
   Cheers,
   Rémi-C
  
  
  
  
   2013/11/15 James David Smith james.david.sm...@gmail.com
  
   Thanks both. Geometries now fixed.
  
   The query 'CREATE TABLE lines_for_each_road' has now been set
   running. Will report back when it's done. I suspect it may take a
   while!
  
   James
  
   On 15 November 2013 11:03, Sandro Santilli s...@keybit.net
 wrote:
On Fri, Nov 15, 2013 at 11:50:42AM +0100, Rémi Cura wrote:
Yep, maybe something like
   
UPDATE ukmajrdbuffer SET the_geom = ST_MakeValid(the_geom)
WHERE ST_IsValid(the_geom) = FALSE
   
ST_MakeValid internally checks for ST_IsValid, so no need
to add the condition (which would run the test twice).
   
--strk;
___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
   ___
   postgis-users mailing list
   postgis-users@lists.osgeo.org
   http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
  
  
  
  
  
   ___
   postgis-users mailing list
   postgis-users@lists.osgeo.org
   http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
  ___
  postgis-users mailing list
  postgis-users@lists.osgeo.org
  http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
 
 
 
  

Re: [postgis-users] Grid of points inside Polygon

[Rémi Cura]

You should try to get the original road file, because i will be very easy
to reapply some buffer and filtering to keep only major road.
Yet, why not let it run during the week end?

Cheers,

Rémi-C


2013/11/15 James David Smith james.david.sm...@gmail.com

 Hey Remi,

 Do you think before I try running the big query you have just sent me,
 that I should go back and try to get the original file of uk roads? I
 mean the very original file that has not had any buffers applied or
 any merging done.

 Or shall we just go for it and see if it has finished when I come into
 work on Monday?! haha.

 Thanks

 James

 On 15 November 2013 14:01, Rémi Cura remi.c...@gmail.com wrote:
  Outch,
  you have only 8 roads in GB? =)
 
  This is not going to go fast till you don't have some dozens k's of
 roads.
 
  I'm guessing you are not at will to send those roads?
 
 
  The next step will be (when you'll be sure everything is OK)
 
 
  CREATE TABLE lines_for_each_road AS
  WITH all_lines AS (
  SELECT *
  FROM all_lines
  ),
  cutted_lines AS ( --we cut the line to keep only part of lines inside
  road_buffer
  SELECT ST_Intersection(all_lines.the_geom,ukmajrdbuffer.geom) as
  lines_cutted, direction
  FROM ukmajrdbuffer, all_lines
  WHERE ST_Intersects(ukmajrdbuffer.geom, all_lines.the_geom)=TRUE
  ),
  cutted_lines_SN AS ( --this is the cutted lines which going from South to
  North
  SELECT *
  FROM cutted_lines
  WHERE direction = 'SN'
  ),
  cutted_lines_EW AS ( --this is the cutted lines going from East toWest
  SELECT *
  FROM cutted_lines
  WHERE direction = 'EW'
  ),
  points AS ( -- we take the intersection of EW  lines with SN lines ,
 that is
  the points on the grid.
  SELECT ST_Intersection(clSN.lines_cutted, clEW.lines_cutted) AS point
  FROM cutted_lines_SN as clSN, cutted_lines_EW AS clEW
  WHERE ST_Intersects(clSN.lines_cutted, clEW.lines_cutted)=TRUE --no
 point ot
  compute an intersection if lines don't intersect
  )
  SELECT row_number() over() AS id  , point
  FROM points ;
 
 
 
 
 
  Cheers,
 
  Rémi-C
 
 
 
  2013/11/15 James David Smith james.david.sm...@gmail.com
 
  Hey Remi,
 
  I'll do a few checks and get back to you. Maybe I did something wrong
  because I set this query going on my local PostgreSQL installation but
  also on our Linux Cluster machine which is much more powerful. And it
  finished on my local installation BEFORE the cluster. So maybe I did
  something wrong on the local query.
 
  The query that has finished took about 1 hour.
  The query on our cluster is still running.
 
  select count(*) from ukmajrdbuffer = 8
  This is because before I was given the data the roads buffer had
  already been dissolved unfortunately.
 
  James
 
 
 
 
  On 15 November 2013 13:43, Rémi Cura remi.c...@gmail.com wrote:
   Good !
  
   Something is strange with the result,
   you get only 190k lines intersecting road buffer,
   it is very few , I expected at least 10times this!
   How many time took this computing by the way?
  
   How many road buffer geom do you have ? (select count(*) from
   ukmajrdbuffer
   ).
  
   Cheers,
   Rémi-C
  
  
   2013/11/15 James David Smith james.david.sm...@gmail.com
  
   Hey.
  
   Yes, it's done. Was just getting some lunch! :-)
  
   select count(*) from lines_for_each_road
   Result = 187033
  
   I have also just ran 'VACUUM ANALYZE' on the tables
   'lines_for_each_road' as well as the table 'all_lines'
  
   I also can confirm that I have ran the following commands:
  
   CREATE INDEX all_lines_index ON all_lines USING GIST ( the_geom )
   CREATE INDEX ukmajrdbuffer_index ON ukmajrdbuffer USING GIST (geom);
  
   Should we now uncomment this line from the previous query?
  
 SELECT row_number() over() AS id--,*  
  
   Thanks again Remi,
  
   James
  
  
  
   On 15 November 2013 13:14, Rémi Cura remi.c...@gmail.com wrote:
Also if you do have indexes,
can you run a VACUUM ANALYZE, so that the indexes will be used?
   
Cheers,
   
Rémi-C
   
   
2013/11/15 Rémi Cura remi.c...@gmail.com
   
It should be finished by now,
can you check you have geom indexes on :
ukmajrdbuffer.geom and  all_lines.the_geom
   
   
How many geoms do you have in ukmajrdbuffer?
   
Cheers,
Rémi-C
   
   
2013/11/15 Rémi Cura remi.c...@gmail.com
   
Hey Sandro,
   
Thanks for this, it is at least twice faster =)
   
Cheers,
Rémi-C
   
   
   
   
2013/11/15 James David Smith james.david.sm...@gmail.com
   
Thanks both. Geometries now fixed.
   
The query 'CREATE TABLE lines_for_each_road' has now been
 set
running. Will report back when it's done. I suspect it may take
 a
while!
   
James
   
On 15 November 2013 11:03, Sandro Santilli s...@keybit.net
wrote:
 On Fri, Nov 15, 2013 at 11:50:42AM +0100, Rémi Cura wrote:
 Yep, maybe something like

 UPDATE ukmajrdbuffer SET the_geom = ST_MakeValid(the_geom)
 WHERE ST_IsValid(the_geom) = FALSE

 

Re: [postgis-users] Grid of points inside Polygon

[James David Smith]

Remi,

Ok. Cool. I've set the below query running. On Monday I will also
attempt to get the original road lines file too. Let's see if we have
a result on Monday.
--
CREATE TABLE lines_for_each_road AS
WITH all_lines AS (
SELECT *
FROM all_lines
),
cutted_lines AS ( --we cut the line to keep only part of lines inside
road_buffer
SELECT ST_Intersection(all_lines.the_geom,ukmajrdbuffer.geom) as
lines_cutted, direction
FROM ukmajrdbuffer, all_lines
WHERE ST_Intersects(ukmajrdbuffer.geom, all_lines.the_geom)=TRUE
),
cutted_lines_SN AS ( --this is the cutted lines which going from South to North
SELECT *
FROM cutted_lines
WHERE direction = 'SN'
),
cutted_lines_EW AS ( --this is the cutted lines going from East toWest
SELECT *
FROM cutted_lines
WHERE direction = 'EW'
),
points AS ( -- we take the intersection of EW  lines with SN lines ,
that is the points on the grid.
SELECT ST_Intersection(clSN.lines_cutted, clEW.lines_cutted) AS point
FROM cutted_lines_SN as clSN, cutted_lines_EW AS clEW
WHERE ST_Intersects(clSN.lines_cutted, clEW.lines_cutted)=TRUE --no
point ot compute an intersection if lines don't intersect
)
SELECT row_number() over() AS id  , point
FROM points ;
--
Thanks

James

On 15 November 2013 15:37, Rémi Cura remi.c...@gmail.com wrote:

 You should try to get the original road file, because i will be very easy to
 reapply some buffer and filtering to keep only major road.
 Yet, why not let it run during the week end?

 Cheers,

 Rémi-C


 2013/11/15 James David Smith james.david.sm...@gmail.com

 Hey Remi,

 Do you think before I try running the big query you have just sent me,
 that I should go back and try to get the original file of uk roads? I
 mean the very original file that has not had any buffers applied or
 any merging done.

 Or shall we just go for it and see if it has finished when I come into
 work on Monday?! haha.

 Thanks

 James

 On 15 November 2013 14:01, Rémi Cura remi.c...@gmail.com wrote:
  Outch,
  you have only 8 roads in GB? =)
 
  This is not going to go fast till you don't have some dozens k's of
  roads.
 
  I'm guessing you are not at will to send those roads?
 
 
  The next step will be (when you'll be sure everything is OK)
 
 
  CREATE TABLE lines_for_each_road AS
  WITH all_lines AS (
  SELECT *
  FROM all_lines
  ),
  cutted_lines AS ( --we cut the line to keep only part of lines inside
  road_buffer
  SELECT ST_Intersection(all_lines.the_geom,ukmajrdbuffer.geom) as
  lines_cutted, direction
  FROM ukmajrdbuffer, all_lines
  WHERE ST_Intersects(ukmajrdbuffer.geom, all_lines.the_geom)=TRUE
  ),
  cutted_lines_SN AS ( --this is the cutted lines which going from South
  to
  North
  SELECT *
  FROM cutted_lines
  WHERE direction = 'SN'
  ),
  cutted_lines_EW AS ( --this is the cutted lines going from East toWest
  SELECT *
  FROM cutted_lines
  WHERE direction = 'EW'
  ),
  points AS ( -- we take the intersection of EW  lines with SN lines ,
  that is
  the points on the grid.
  SELECT ST_Intersection(clSN.lines_cutted, clEW.lines_cutted) AS point
  FROM cutted_lines_SN as clSN, cutted_lines_EW AS clEW
  WHERE ST_Intersects(clSN.lines_cutted, clEW.lines_cutted)=TRUE --no
  point ot
  compute an intersection if lines don't intersect
  )
  SELECT row_number() over() AS id  , point
  FROM points ;
 
 
 
 
 
  Cheers,
 
  Rémi-C
 
 
 
  2013/11/15 James David Smith james.david.sm...@gmail.com
 
  Hey Remi,
 
  I'll do a few checks and get back to you. Maybe I did something wrong
  because I set this query going on my local PostgreSQL installation but
  also on our Linux Cluster machine which is much more powerful. And it
  finished on my local installation BEFORE the cluster. So maybe I did
  something wrong on the local query.
 
  The query that has finished took about 1 hour.
  The query on our cluster is still running.
 
  select count(*) from ukmajrdbuffer = 8
  This is because before I was given the data the roads buffer had
  already been dissolved unfortunately.
 
  James
 
 
 
 
  On 15 November 2013 13:43, Rémi Cura remi.c...@gmail.com wrote:
   Good !
  
   Something is strange with the result,
   you get only 190k lines intersecting road buffer,
   it is very few , I expected at least 10times this!
   How many time took this computing by the way?
  
   How many road buffer geom do you have ? (select count(*) from
   ukmajrdbuffer
   ).
  
   Cheers,
   Rémi-C
  
  
   2013/11/15 James David Smith james.david.sm...@gmail.com
  
   Hey.
  
   Yes, it's done. Was just getting some lunch! :-)
  
   select count(*) from lines_for_each_road
   Result = 187033
  
   I have also just ran 'VACUUM ANALYZE' on the tables
   'lines_for_each_road' as well as the table 'all_lines'
  
   I also can confirm that I have ran the following commands:
  
   CREATE INDEX all_lines_index ON all_lines USING GIST ( the_geom )
   CREATE INDEX ukmajrdbuffer_index ON ukmajrdbuffer USING GIST (geom);
  
   Should we 

Re: [postgis-users] Grid of points inside Polygon

[James David Smith]

Hey Remi,

I've actually managed to get the file 'ukmajorroads' already and have
loaded it into my database.  There are 395356 rows in the database. There
is a field called 'geom' and I have built an index on it as below;

CREATE INDEX ukrds_index ON ukrds USING GIST (geom);

So should I now run the same query, but do the buffering of the roads 'on
the fly' maybe?  For example as below?

CREATE TABLE lines_for_each_road AS
WITH all_lines AS (
SELECT *
FROM all_lines
),
cutted_lines AS ( --we cut the line to keep only part of lines inside
road_buffer
SELECT ST_Intersection(all_lines.the_geom,st_buffer(ukmajrdbuffer.geom,
1000)) as lines_cutted, direction
FROM ukmajrdbuffer, all_lines
WHERE ST_Intersects(st_buffer(ukmajrdbuffer.geom, 1000),
all_lines.the_geom)=TRUE
),
cutted_lines_SN AS ( --this is the cutted lines which going from South to
North
SELECT *
FROM cutted_lines
WHERE direction = 'SN'
),
cutted_lines_EW AS ( --this is the cutted lines going from East toWest
SELECT *
FROM cutted_lines
WHERE direction = 'EW'
),
points AS ( -- we take the intersection of EW  lines with SN lines , that
is the points on the grid.
SELECT ST_Intersection(clSN.lines_cutted, clEW.lines_cutted) AS point
FROM cutted_lines_SN as clSN, cutted_lines_EW AS clEW
WHERE ST_Intersects(clSN.lines_cutted, clEW.lines_cutted)=TRUE --no point
ot compute an intersection if lines don't intersect
)
SELECT row_number() over() AS id  , point
FROM points ;


On 15 November 2013 15:42, Rémi Cura remi.c...@gmail.com wrote:

 Still is a shame :
 with proper data we should have some result with about one hour i guess.

 Cheers,
 Rémi C


 2013/11/15 James David Smith james.david.sm...@gmail.com

 Remi,

 Ok. Cool. I've set the below query running. On Monday I will also
 attempt to get the original road lines file too. Let's see if we have
 a result on Monday.
 --
 CREATE TABLE lines_for_each_road AS
 WITH all_lines AS (
 SELECT *
 FROM all_lines
 ),
 cutted_lines AS ( --we cut the line to keep only part of lines inside
 road_buffer
 SELECT ST_Intersection(all_lines.the_geom,ukmajrdbuffer.geom) as
 lines_cutted, direction
 FROM ukmajrdbuffer, all_lines
 WHERE ST_Intersects(ukmajrdbuffer.geom, all_lines.the_geom)=TRUE
 ),
 cutted_lines_SN AS ( --this is the cutted lines which going from South to
 North
 SELECT *
 FROM cutted_lines
 WHERE direction = 'SN'
 ),
 cutted_lines_EW AS ( --this is the cutted lines going from East toWest
 SELECT *
 FROM cutted_lines
 WHERE direction = 'EW'
 ),
 points AS ( -- we take the intersection of EW  lines with SN lines ,
 that is the points on the grid.
 SELECT ST_Intersection(clSN.lines_cutted, clEW.lines_cutted) AS point
 FROM cutted_lines_SN as clSN, cutted_lines_EW AS clEW
 WHERE ST_Intersects(clSN.lines_cutted, clEW.lines_cutted)=TRUE --no
 point ot compute an intersection if lines don't intersect
 )
 SELECT row_number() over() AS id  , point
 FROM points ;
 --
 Thanks

 James

 On 15 November 2013 15:37, Rémi Cura remi.c...@gmail.com wrote:
 
  You should try to get the original road file, because i will be very
 easy to
  reapply some buffer and filtering to keep only major road.
  Yet, why not let it run during the week end?
 
  Cheers,
 
  Rémi-C
 
 
  2013/11/15 James David Smith james.david.sm...@gmail.com
 
  Hey Remi,
 
  Do you think before I try running the big query you have just sent me,
  that I should go back and try to get the original file of uk roads? I
  mean the very original file that has not had any buffers applied or
  any merging done.
 
  Or shall we just go for it and see if it has finished when I come into
  work on Monday?! haha.
 
  Thanks
 
  James
 
  On 15 November 2013 14:01, Rémi Cura remi.c...@gmail.com wrote:
   Outch,
   you have only 8 roads in GB? =)
  
   This is not going to go fast till you don't have some dozens k's of
   roads.
  
   I'm guessing you are not at will to send those roads?
  
  
   The next step will be (when you'll be sure everything is OK)
  
  
   CREATE TABLE lines_for_each_road AS
   WITH all_lines AS (
   SELECT *
   FROM all_lines
   ),
   cutted_lines AS ( --we cut the line to keep only part of lines inside
   road_buffer
   SELECT ST_Intersection(all_lines.the_geom,ukmajrdbuffer.geom) as
   lines_cutted, direction
   FROM ukmajrdbuffer, all_lines
   WHERE ST_Intersects(ukmajrdbuffer.geom, all_lines.the_geom)=TRUE
   ),
   cutted_lines_SN AS ( --this is the cutted lines which going from
 South
   to
   North
   SELECT *
   FROM cutted_lines
   WHERE direction = 'SN'
   ),
   cutted_lines_EW AS ( --this is the cutted lines going from East
 toWest
   SELECT *
   FROM cutted_lines
   WHERE direction = 'EW'
   ),
   points AS ( -- we take the intersection of EW  lines with SN lines ,
   that is
   the points on the grid.
   SELECT ST_Intersection(clSN.lines_cutted, clEW.lines_cutted) AS point
   FROM cutted_lines_SN as clSN, cutted_lines_EW AS clEW
   WHERE 

Re: [postgis-users] Grid of points inside Polygon

[Rémi Cura]

good =)

It is better to do the buffer computation first,
because this way you are sure youwon't compute the same buffer many times.

Something like (not so sure without trying it)
UPDATE ukmajrdbuffer SET geom = ST_Buffer(geom,1000);

You may have an error if you hav e a fixed geom data type, then you could
just add a new column like

ALTER TABLE ukmajrdbuffer ADD COLUMN the_geom geometry ;
UPDATE ukmajrdbuffer SET the_geom = ST_Buffer(geom,1000);


Cheers,

Rémi-C



2013/11/15 James David Smith james.david.sm...@gmail.com

 Hey Remi,

 I've actually managed to get the file 'ukmajorroads' already and have
 loaded it into my database.  There are 395356 rows in the database. There
 is a field called 'geom' and I have built an index on it as below;

 CREATE INDEX ukrds_index ON ukrds USING GIST (geom);

 So should I now run the same query, but do the buffering of the roads 'on
 the fly' maybe?  For example as below?


 CREATE TABLE lines_for_each_road AS
 WITH all_lines AS (
 SELECT *
 FROM all_lines
 ),
 cutted_lines AS ( --we cut the line to keep only part of lines inside
 road_buffer
 SELECT ST_Intersection(all_lines.the_geom,st_buffer(ukmajrdbuffer.geom,
 1000)) as lines_cutted, direction
 FROM ukmajrdbuffer, all_lines
 WHERE ST_Intersects(st_buffer(ukmajrdbuffer.geom, 1000),
 all_lines.the_geom)=TRUE
 ),
 cutted_lines_SN AS ( --this is the cutted lines which going from South to
 North
 SELECT *
 FROM cutted_lines
 WHERE direction = 'SN'
 ),
 cutted_lines_EW AS ( --this is the cutted lines going from East toWest
 SELECT *
 FROM cutted_lines
 WHERE direction = 'EW'
 ),
 points AS ( -- we take the intersection of EW  lines with SN lines , that
 is the points on the grid.
 SELECT ST_Intersection(clSN.lines_cutted, clEW.lines_cutted) AS point
 FROM cutted_lines_SN as clSN, cutted_lines_EW AS clEW
 WHERE ST_Intersects(clSN.lines_cutted, clEW.lines_cutted)=TRUE --no point
 ot compute an intersection if lines don't intersect
 )
 SELECT row_number() over() AS id  , point
 FROM points ;


 On 15 November 2013 15:42, Rémi Cura remi.c...@gmail.com wrote:

 Still is a shame :
 with proper data we should have some result with about one hour i guess.

 Cheers,
 Rémi C


 2013/11/15 James David Smith james.david.sm...@gmail.com

 Remi,

 Ok. Cool. I've set the below query running. On Monday I will also
 attempt to get the original road lines file too. Let's see if we have
 a result on Monday.
 --
 CREATE TABLE lines_for_each_road AS
 WITH all_lines AS (
 SELECT *
 FROM all_lines
 ),
 cutted_lines AS ( --we cut the line to keep only part of lines inside
 road_buffer
 SELECT ST_Intersection(all_lines.the_geom,ukmajrdbuffer.geom) as
 lines_cutted, direction
 FROM ukmajrdbuffer, all_lines
 WHERE ST_Intersects(ukmajrdbuffer.geom, all_lines.the_geom)=TRUE
 ),
 cutted_lines_SN AS ( --this is the cutted lines which going from South
 to North
 SELECT *
 FROM cutted_lines
 WHERE direction = 'SN'
 ),
 cutted_lines_EW AS ( --this is the cutted lines going from East toWest
 SELECT *
 FROM cutted_lines
 WHERE direction = 'EW'
 ),
 points AS ( -- we take the intersection of EW  lines with SN lines ,
 that is the points on the grid.
 SELECT ST_Intersection(clSN.lines_cutted, clEW.lines_cutted) AS point
 FROM cutted_lines_SN as clSN, cutted_lines_EW AS clEW
 WHERE ST_Intersects(clSN.lines_cutted, clEW.lines_cutted)=TRUE --no
 point ot compute an intersection if lines don't intersect
 )
 SELECT row_number() over() AS id  , point
 FROM points ;
 --
 Thanks

 James

 On 15 November 2013 15:37, Rémi Cura remi.c...@gmail.com wrote:
 
  You should try to get the original road file, because i will be very
 easy to
  reapply some buffer and filtering to keep only major road.
  Yet, why not let it run during the week end?
 
  Cheers,
 
  Rémi-C
 
 
  2013/11/15 James David Smith james.david.sm...@gmail.com
 
  Hey Remi,
 
  Do you think before I try running the big query you have just sent me,
  that I should go back and try to get the original file of uk roads? I
  mean the very original file that has not had any buffers applied or
  any merging done.
 
  Or shall we just go for it and see if it has finished when I come into
  work on Monday?! haha.
 
  Thanks
 
  James
 
  On 15 November 2013 14:01, Rémi Cura remi.c...@gmail.com wrote:
   Outch,
   you have only 8 roads in GB? =)
  
   This is not going to go fast till you don't have some dozens k's of
   roads.
  
   I'm guessing you are not at will to send those roads?
  
  
   The next step will be (when you'll be sure everything is OK)
  
  
   CREATE TABLE lines_for_each_road AS
   WITH all_lines AS (
   SELECT *
   FROM all_lines
   ),
   cutted_lines AS ( --we cut the line to keep only part of lines
 inside
   road_buffer
   SELECT ST_Intersection(all_lines.the_geom,ukmajrdbuffer.geom) as
   lines_cutted, direction
   FROM ukmajrdbuffer, all_lines
   WHERE ST_Intersects(ukmajrdbuffer.geom, 

Re: [postgis-users] Grid of points inside Polygon

[James David Smith]

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.c...@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.sm...@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@lists.osgeo.org
 http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users



 ___
 postgis-users mailing list
 postgis-users@lists.osgeo.org
 http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

Re: [postgis-users] Grid of points inside Polygon

[Rémi Cura]

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.sm...@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.c...@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.sm...@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 

Re: [postgis-users] Grid of points inside Polygon

[James David Smith]

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.c...@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.sm...@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.c...@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.sm...@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:
 
  

Re: [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.sm...@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.c...@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.sm...@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.c...@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.sm...@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 

Re: [postgis-users] Grid of points inside Polygon

[Rémi Cura]

Using the lines you exactly do the same as with points.
The lines will be spaced exactly by 20 meters, end and start being exactly
on the 20 meters grid, line being either perfectly veritcla or perfectly
horizontal.
So intersections are guaranteed to be on the 20 meter grid, because this is
the definition of a 20 meter grid =) .


Here is a not optimal way to do it :
For instance if UK is enclosed in a bouding box (xmin,ymin,xmax,ymax),
where xmin-max and ymin-max are rounded to the nearest 20 meters multiple
(ie x or y % 20 = 0)
(Please note that this is not clever because we generate too much SN lines,
GB being only 500 km wide at most, so we would need to filter, but you
don't care, following would take less than a minute anyway, and you need to
do it only once).

WITH minmax AS (
SELECT xmin,ymin,xmax,ymax, bbox_of_england
),
WITH series AS (
SELECT s
FROM generate_series(0,1000*(1000/20)::bigint,20) AS s --casting to bigint
to be safe against going over int limit
),
WITH lines AS (
SELECT ST_MakeLine (ST_MakePoint(xmin+s,ymin),ST_MakePoint(xmin+s,ymax)) AS
lineSN,
 ST_MakeLine(ST_MakePoint(xmin, ymin+s), ST_MakePoint(xmax,ymin+s) AS lineEW
FROM series, minmax
),
unioned_lines AS (

SELECT lineSN AS line, 'SN' as direction
FROM lines,minmax
WHERE ST_Interesects(lineSN,bbox_of_england)= TRUE

UNION

SELECT lineEW AS line, 'EW' AS direction
FFROM lines,minmax
WHERE ST_Interesects(lineEW,bbox_of_england)= TRUE

)
SELECT row_number() over() AS uid, line, direction
FROM unioned_line

then creating table and index on geom and on direction (and possibly on
uid, but better create a primary key then if you use it)

Of course I didn't test this query, but you should be able to use it easily.

Cheers,
Rémi-C




2013/11/13 James David Smith james.david.sm...@gmail.com

 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.sm...@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.c...@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 

Re: [postgis-users] Grid of points inside Polygon

[James David Smith]

Hey Remi,

I don't quite get what the query you gave does to be honest. But
working with some of the ideas and notions you gave me I have put the
below together. The final bit of the code is the key query. The tables
of NS and EW lines were done very quickly. Nice. But now doing the
'intersects within polygon' is taking some time so I'll see how it
gets on when I get in to work tomorrow morning:

DROP TABLE IF EXISTS north_to_south;
DROP TABLE IF EXISTS east_to_west;

-- Make a table of north to south lines
CREATE TABLE north_to_south tablespace jamestabs AS
( SELECT st_setsrid(st_makeline(st_makepoint(x, 7780), st_makepoint(x,
1226580)),27700) as the_geom
FROM generate_series(53320::int, 667380::int, 20) as x);

-- Add an index
CREATE INDEX north_to_south_geom_index ON north_to_south USING GIST (
the_geom );

-- Make a table of east to west lines
CREATE TABLE east_to_west tablespace jamestabs AS
( SELECT st_setsrid(st_makeline(st_makepoint(53320, y),
st_makepoint(667380, y)),27700) as the_geom
FROM generate_series(7780::int, 1226580::int, 20) as y);

-- Add an index
CREATE INDEX east_to_west_geom_index ON east_to_west USING GIST ( the_geom );

-- Import the Road Buffer file using the pg_loader interface

-- Update the geometry field as the loader didn't seem to do this.
SELECT UpdateGeometrySRID('ukmajrdbuffer','geom',27700);

-- Merge the road buffer polygons into one big polygon
-- Should maybe consider leaving them as they are, but am merging for now.
SELECT st_union(geom) as geom INTO union_roads FROM ukmajrdbuffer;

-- The unioned roads file isn't valid geom. So I simplify it a little
to make it valid.
-- We simplify with a 70 metres factor. Geom is now valid.
SELECT ST_SimplifyPreserveTopology(geom, 70) AS geom INTO
buff_union_roads FROM union_roads;

-- Now calculate the points where the lines intersect and also where
they are within the polygon.
CREATE TABLE gridpoints tablespace jamestabs
AS (SELECT ST_Intersection(east_to_west.the_geom, north_to_south.the_geom)
FROM east_to_west, north_to_south, buff_union_roads
WHERE st_within(buff_union_roads.geom,
ST_Intersection(east_to_west.the_geom, north_to_south.the_geom)) =
TRUE);

Thanks

James




On 13 November 2013 13:33, Rémi Cura remi.c...@gmail.com wrote:
 Using the lines you exactly do the same as with points.
 The lines will be spaced exactly by 20 meters, end and start being exactly
 on the 20 meters grid, line being either perfectly veritcla or perfectly
 horizontal.
 So intersections are guaranteed to be on the 20 meter grid, because this is
 the definition of a 20 meter grid =) .


 Here is a not optimal way to do it :
 For instance if UK is enclosed in a bouding box (xmin,ymin,xmax,ymax), where
 xmin-max and ymin-max are rounded to the nearest 20 meters multiple (ie x or
 y % 20 = 0)
 (Please note that this is not clever because we generate too much SN lines,
 GB being only 500 km wide at most, so we would need to filter, but you don't
 care, following would take less than a minute anyway, and you need to do it
 only once).

 WITH minmax AS (
 SELECT xmin,ymin,xmax,ymax, bbox_of_england
 ),
 WITH series AS (
 SELECT s
 FROM generate_series(0,1000*(1000/20)::bigint,20) AS s --casting to bigint
 to be safe against going over int limit
 ),
 WITH lines AS (
 SELECT ST_MakeLine (ST_MakePoint(xmin+s,ymin),ST_MakePoint(xmin+s,ymax)) AS
 lineSN,
  ST_MakeLine(ST_MakePoint(xmin, ymin+s), ST_MakePoint(xmax,ymin+s) AS lineEW
 FROM series, minmax
 ),
 unioned_lines AS (

 SELECT lineSN AS line, 'SN' as direction
 FROM lines,minmax
 WHERE ST_Interesects(lineSN,bbox_of_england)= TRUE

 UNION

 SELECT lineEW AS line, 'EW' AS direction
 FFROM lines,minmax
 WHERE ST_Interesects(lineEW,bbox_of_england)= TRUE

 )
 SELECT row_number() over() AS uid, line, direction
 FROM unioned_line

 then creating table and index on geom and on direction (and possibly on uid,
 but better create a primary key then if you use it)

 Of course I didn't test this query, but you should be able to use it easily.

 Cheers,
 Rémi-C




 2013/11/13 James David Smith james.david.sm...@gmail.com

 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.sm...@gmail.com wrote:
  Hi Remi,
 
  Thanks so much for this detailed response. Your idea about creating
  

Re: [postgis-users] Grid of points inside Polygon

[Rémi Cura]

All lines should be in one table, you can add a column to differentiate
between line going SN and EW,
this is what the query I wrote does.

Not changing your code much, it gives for the line table :

CREATE TABLE all_lines tablespace jamestabs AS
( SELECT st_setsrid(st_makeline(st_makepoint(x, 7780), st_makepoint(x,
1226580)),27700) as the_geom, 'SN' AS direction
FROM generate_series(53320::int, 667380::int, 20) as x)
UNION
( SELECT st_setsrid(st_makeline(st_makepoint(53320, y),
st_makepoint(667380, y)),27700) as the_geom, 'EW' AS direction
FROM generate_series(7780::int, 1226580::int, 20) as y)

then building index on it.


-- Merge the road buffer polygons into one big polygon
-- Should maybe consider leaving them as they are, but am merging for now.
SELECT st_union(geom) as geom INTO union_roads FROM ukmajrdbuffer;
*You must not do it ! *
As I told you blabla index blabla.

instead, create an index on table ukmajrdbuffer, and use it without
unioning.

Then you query is not going to work at all ... (or *very very* slowly)
FROM east_to_west, north_to_south, buff_union_roads : it potentially
generates 1000*50*1000*50 * number_of_road rows to filter, very bad idea,
you'll get the same problem as with points...

(please see some example on internet to understand WITH, it is simple
, useful, and makes query much more clear)

First you could try :

CREATE TABLE lines_for_each_road tablespace jamestabs AS
WITH all_lines AS ( --this subquerry gets all the lines

SELECT *
FROM all_lines

),
SELECT row_number() over() AS id--,*
FROM ukmajrdbuffer, all_lines
WHERE ST_Intersects(ukmajrdbuffer.the_geom, all_lines.the_geom)=TRUE;

This will give you (if you uncomment the --,*, ), for every road, all the
lines intersecting the road
example :
road_1 | line_2_EW
road_1 | line_3_EW
road_1 | line_78_SN
road_1 | line_7_EW
road_2 ...

*How long is it (don't uncomment plz), and how many rows do you get (select
count(*) from lines_for_each_road )?*

If it is not too long, we will go on and cut the lines so that for every
road, we keep the part of the lines that are inside the road_buffer.
(I have to leave now, tomorrow)

Cheers,
Rémi-C





2013/11/13 James David Smith james.david.sm...@gmail.com

 Hey Remi,

 I don't quite get what the query you gave does to be honest. But
 working with some of the ideas and notions you gave me I have put the
 below together. The final bit of the code is the key query. The tables
 of NS and EW lines were done very quickly. Nice. But now doing the
 'intersects within polygon' is taking some time so I'll see how it
 gets on when I get in to work tomorrow morning:

 DROP TABLE IF EXISTS north_to_south;
 DROP TABLE IF EXISTS east_to_west;

 -- Make a table of north to south lines
 CREATE TABLE north_to_south tablespace jamestabs AS
 ( SELECT st_setsrid(st_makeline(st_makepoint(x, 7780), st_makepoint(x,
 1226580)),27700) as the_geom
 FROM generate_series(53320::int, 667380::int, 20) as x);

 -- Add an index
 CREATE INDEX north_to_south_geom_index ON north_to_south USING GIST (
 the_geom );

 -- Make a table of east to west lines
 CREATE TABLE east_to_west tablespace jamestabs AS
 ( SELECT st_setsrid(st_makeline(st_makepoint(53320, y),
 st_makepoint(667380, y)),27700) as the_geom
 FROM generate_series(7780::int, 1226580::int, 20) as y);

 -- Add an index
 CREATE INDEX east_to_west_geom_index ON east_to_west USING GIST ( the_geom
 );

 -- Import the Road Buffer file using the pg_loader interface

 -- Update the geometry field as the loader didn't seem to do this.
 SELECT UpdateGeometrySRID('ukmajrdbuffer','geom',27700);

 -- Merge the road buffer polygons into one big polygon
 -- Should maybe consider leaving them as they are, but am merging for now.
 SELECT st_union(geom) as geom INTO union_roads FROM ukmajrdbuffer;

 -- The unioned roads file isn't valid geom. So I simplify it a little
 to make it valid.
 -- We simplify with a 70 metres factor. Geom is now valid.
 SELECT ST_SimplifyPreserveTopology(geom, 70) AS geom INTO
 buff_union_roads FROM union_roads;

 -- Now calculate the points where the lines intersect and also where
 they are within the polygon.
 CREATE TABLE gridpoints tablespace jamestabs
 AS (SELECT ST_Intersection(east_to_west.the_geom, north_to_south.the_geom)
 FROM east_to_west, north_to_south, buff_union_roads
 WHERE st_within(buff_union_roads.geom,
 ST_Intersection(east_to_west.the_geom, north_to_south.the_geom)) =
 TRUE);

 Thanks

 James




 On 13 November 2013 13:33, Rémi Cura remi.c...@gmail.com wrote:
  Using the lines you exactly do the same as with points.
  The lines will be spaced exactly by 20 meters, end and start being
 exactly
  on the 20 meters grid, line being either perfectly veritcla or perfectly
  horizontal.
  So intersections are guaranteed to be on the 20 meter grid, because this
 is
  the definition of a 20 meter grid =) .
 
 
  Here is a not optimal way to do it :
  For instance if UK is enclosed in a bouding box (xmin,ymin,xmax,ymax),
 where
  xmin-max 

Re: [postgis-users] Grid of points inside Polygon

[Rémi Cura]

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.sm...@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@lists.osgeo.org
 http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users