8

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?

Improve this question
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. Commented Jun 3, 2014 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 Commented Jun 3, 2014 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. Commented Jun 3, 2014 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). Commented Jun 3, 2014 at 10:10

3 Answers 3

Reset to default
6

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

Improve this answer
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! Commented Jun 10, 2014 at 12:09
2

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;
Improve this answer
0

Now the qgis plugin Closest_point do this job.

Improve this answer

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.