Tuesday, 4 December 2018

qgis - Find closest point on a MulitLineString for each point of a MultiPoint


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 ?



Answer



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.


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...