2

I'm working on a project to connect buildings to electric power plant via a network. To this end I'm beginning by trying to connect every buildings to the network.

I've importing the shp files to a database via shp2pgsql and now I've got a table of MultiPoints (buildings, with p the number of rows) and another one of MultiLineString (network, with n the number of rows).

I used ST_ClosestPoint and ST_ShortestLine on these tables but I'm facing the same problem with both of them : it returns a table of n x p rows because it computes the closest point (or the shortest line) from each point to each line of the network. My goal is to find the closest point on any line of the network for each building.

Right now the best I've got is this query :

SELECT min(ST_Length(ST_ShortestLine(network.geom, buildings.geom))) as length, id FROM network, building group by id

It returns a table of the length of the shortest line for each building but I can't find any way to get the geometry of the concern line (in order to know the location of the point on the network).

I could do a sub-query to do it point by point, order the n rows by calculing the length of each line and keep the one with the shortest length but the final goal is to use this query on files with over 200 000 buildings and a network of over 50 000 lines so I think I need to find the most efficient way.

Moreover, I did not find any way to find the "closest closest point" using the result of the query ST_ClosestPoints(network, buildings).

Does anyone think of a way of doing this efficiently ?

2

What you need to do is first find the closes geometry, and then find the closest point on that geometry.

This can be done in a few ways:

select st_closestpoint(a_geom, b_geom), a_id, b_id from 
(
  select distinct on (a.id) st_distance(a.geom, b.geom) dist,
  a.geom a_geom, b.geom b_geom, a.id a_id, b.id b_id from a, b order 
  by a.id,dist
) c 

Here we calculate every distance combination, which is slow on big data sets

select st_closestpoint(a_geom, b_geom), a_id, b_id from 
(
  select distinct on (a.gid) st_distance(a.geom, b.geom) dist,
  a.geom a_geom, b.geom b_geom, a.id a_id, b.id b_id from a, b where st_dwithin(a.geom,b.geom,1000) order by a.id,dist
) c 

Here we use st_dwithin to first find combinations that is closer than 1000 meters (or whatever unit your projection is in). ST_Dwiting uses the spatial index so is faster

select st_closestpoint(a.geom, b.geom), a.gid, b.gid 
from a cross join lateral 
  (select * from  b order by a.geom<->b.geom limit 1) b

This is a quite new one. I haven't tested the code above and it can be some typos or bugs. But it is described here Here we use knn calculation which is index based. So we get the gain of using index without guessing about how far away they might be.

  • Thanks for trying to help me. You're right, I think the best way is to find the closest geometry to a point and then find the closest point on that geometry. However, I've been trying to solutions without sucess. I understand what you're trying to do in the second solution and it seems to me like a great solution but I think that the syntax is not right and I do not suceed to correct it. – Florian L. Apr 6 '17 at 14:14
  • Ok, I will take a look and get back – Nicklas Avén Apr 6 '17 at 15:52
  • 1
    @FlorianL. I'm sorry, I guess that was some world record of typos. Lesson learned: never post without double checking and testing. I hope this works better. – Nicklas Avén Apr 6 '17 at 16:06
  • It works ! Thanks. Just to clarify, if you want st_closestpoint to return the closest point on the line, you have to pass the argument in the right order (network, points). So in your example, a would be the network and b the points. But to make the subquery work, the select distinct should be on the points ids (b.id) and it should be order by b.id too. – Florian L. Apr 7 '17 at 8:05
  • I found an extension in QGIS which does exactly what I was looking for. It is called Network and the function is Connect. It connects each point of a Point vector layer to the closest point on a PolyLine vector layer. – Florian L. Apr 11 '17 at 9:49

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.