7

I have a dataset osmMajor containing road data (from OSM; 670,780 features) and another dataset points containing some points of interest (500,000 points). For each feature in my point layer, I want to find the coordinate of the closest point located on the road network.

enter image description here

I have two options to do this: either ArcMap's near tool (Analysis tools/Proximity/Near), or Spatialite. I tried using ArcMap, but it seems a very long task to perform (I had to stop the process after 3 hours). For this reason, I was thinking about using Spatialite to do this, but it seems that the process is less straightforward than what I thought.

How should I write my Spatialite query to get what I want?

For Spatialite, I think that I should follow these steps:

  1. For each point feature, calculate the distance between the point and the closest road
  2. Around each point, create a buffer using the distance to the closest feature
  3. Get the intersection between the buffer and the road data (it should be a point)
  4. Get the coordinates (x,y) of the resulting point

Trying to implement these 4 steps on a sample dataset with only 10 points, and I encountered several issues.

For point (1.), I can find the overall smallest distance between the points and the roads, but I cannot have it for each point:

SELECT test.ID, min(Distance(test.Geometry, osmMajor.Geometry)) AS distance
FROM test, osmMajor
WHERE osmMajor.ROWID IN (
SELECT ROWID 
FROM SpatialIndex
WHERE f_table_name = 'osmMajor'
AND search_frame = test.Geometry)
GROUP BY test.ID

I get a distance for only three points instead of the 10 that are in the test table. I suspect it has something to do with the way the spatial index works, because the query only returns the results for points 2, 3 and 4, which are the closest to the road network. enter image description here

Can anyone give me advice on the best way to get the coordinate of the point on a feature which is closest to another point?

4
  • Two steps should be enough with Spatialite: 1) find the nearest way with min(Distance... 2) Find the closest point from the nearest way with ClosestPoint. Your Spatialite must be compiled with GEOS for ClosestPoint, see www.gaia-gis.it/gaia-sins/spatialite-sql-latest.html. Your spatial index query does nor work as you suppose. Drop out the subquery and test with small amount of features first until you get correct results. It will be slow but using spatial index in Spatialite is a bit tricky. Once the query works, add subquery and check that the result remains good.
    – user30184
    Jun 3 '14 at 9:03
  • Thank you! I didn't know about the ClosestPoint() function, and it does exactly what I want. I just tried it and it works almost perfectly, and very quickly with the spatial index on. First, for one of my points, I obtain 2 closest points instead of one. Second, I still have the issue with the spatial index (I get the results for only 3 out of 10 points): i.imgur.com/tHdRc5Q.png
    – Damien
    Jun 3 '14 at 9:18
  • When you use AND search_frame = test.Geometry your frame is a point. It finds only ways which has such a BBOX that the reference "test.Geometry" happens to be inside that BBOX. You must build bigger search_frame around the points for example with BuildCircleMbr but that brings a problem. If the search frame is too small you won't still find hits always, if it is unnecessary big the query will not be selective and fast. You can decide that you are not interested if the road is farther than x km from the point, or increase the frame iteratively until you find something.
    – user30184
    Jun 3 '14 at 9:30
  • Thanks for the explanation. However, when I set a BBOX using BuildMbrCircle or even Buffer(), I obtain very weird results, such as the green point in this picture i.imgur.com/znNRc4m.png When I do not set a box, I only have results for 3 points, but they are correct (the star). In addition, I am surprised that the ClosestPoint() function returns more than one row for each point. If I remove the "GROUP BY test.ID" clause and without BBOX, I have two results for the point with ID=2 (and even more results if a set a BBOX).
    – Damien
    Jun 3 '14 at 10:10
5

The logic of the Spatialite process will be:

1) Find the road that is closest to the reference point with

select 1.ID,2.ID,Min(Distance(1.geometry,2.geometry);

2) Find the closest point of the closest road that was just identified

select ClosestPoint(1.geometry,2.geometry) from 1,2 where 1.ID=x and 2.ID=y;

3) Tune the subquery to make the search from the spatial b-tree index of Spatialite so that it does not drop correct hits from the resultset but it still gives speed-up for the query.

Edit: Test procedure and results to prove that it is possible

It is often better to do the job step-by-step by using temporary tables. I made a line layer with 1602 features and point layer with 294 points and made a simple test.

First step is to create a table that select the closest road from each point and saves the result into a new table "mindist". Point IDs and line IDs are saved into the table for the next step. This is the slow part and spatial index would make the query faster.

create table mindist as
select "p"."ogc_fid"  pid, 
"l"."ogc_fid"  lid, 
min(Distance("p"."geometry","l"."geometry")) as distance
from "poi" p, "lin" l
group by pid
order by pid;

The next step is to go through the temporary table and select for each row the closest point from the geometry of the road to the point geometry. Points and lines are found with the IDs and because they are primary keys in the main tables the query will be fast.

create table closest_point as
select m.rowid as pid, lid, 
closestpoint(l.geometry,p.geometry)as geometry
from mindist m, lin l, poi p
where m.pid=p.ogc_fid
and m.lid=l.ogc_fid
group by m.rowid;

This image shows the results. Lines are the original line layer (motorways), red dots are the original points (villages) and stars show the closest point on a motorway from the villages. Results look visually good to me.

closest points on motorways from village centers

1
  • Thank you for your very detailled answer. I tested your queries and they work well. For my specific situation, I encountered memory errors due to the huge amount of data. I had to find a workaround using ArcGis batch processes and managed to get it one. Thanks again for your help!
    – Damien
    Jun 10 '14 at 12:09
1

Since version 4.4.0 SpatiaLite supports a KNN virtual table index for nearest neighbor problems. Here is a query the finds the nearest line in a linestring table to each point in a point table.

SELECT k.* FROM knn k, points p
WHERE f_table_name = 'linestrings' 
AND ref_geometry = p.geometry
AND max_items = 1;
0

Now the qgis plugin Closest_point do this job.

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.