# Finding points along a path

I am unsure of the best approach to this problem. I can think of some ways of doing it which I will list below but I am looking for the best and most efficient way of doing this.

# Problem

• Given a variable number of points, create a "route" or line between the points.
• Create a bounding box (of 50km[arbitrary]) around the route.
• Return all points within this bounding box
• Points must be returned in order along the route (A->D)

Note: The radius does not have to be circles. A,B,C,D can use rectangular bounding boxes. Data

• Table with points (x,y - technically lat/lon)
• Input points to create a route
• Distance from route to include points

We are using PostgreSQL (8.4) and PostGIS (1.3.6) and Python.

This is the solution I came up with. Thoughts? Ideas? Input?

# Possible Solution

• Create three polygon "tubes" (A->B, B->C, C->D)
• Filter points to only those in the polygons.
• Subtract B->C from A->B so the polygons do not overlap (no duplicates)
• Subtract C->D from B->C so the polygons do not overlap (no duplicates)
• For each polygon
• Calculate vector (A->B, B->C, C->D)
• Find furthest point opposite of vector in the polygon.
• Calculate distance from furthest point to each site in polygon.
• Return points in order of distance (smallest -> largest)
• Join all the points together (A->B points, B->C points, C->D points)

# Solution

My approach was naive, after taking some time and looking at the PostGIS methods I came up with this single SQL call. Note its specific to our database but might help someone in the future. Im sure it could be improved as well.

``````    SELECT aerodrome_id
FROM
geo_aerodromecustom_view,
(SELECT ST_Buffer(
ST_Transform(
(SELECT ST_MakeLine(geom)
FROM
(SELECT geom
FROM geo_aerodromecustom_view
WHERE aerodrome_id IN ('CYOW','CYVR')
ORDER BY CASE aerodrome_id WHEN 'CYOW' THEN 1 WHEN 'CYVR' THEN 2 ELSE 5 END
) sites
),
900913
),
50000) as geo
) as line
WHERE ST_Intersects(ST_Transform(geom, 900913), line.geo)
ORDER BY ST_Line_Locate_Point(
ST_Transform(
(SELECT ST_MakeLine(geom)
FROM
(SELECT geom
FROM geo_aerodromecustom_view
WHERE aerodrome_id IN ('CYOW','CYVR')
ORDER BY CASE aerodrome_id WHEN 'CYOW' THEN 1 WHEN 'CYVR' THEN 2 ELSE 5 END
) sites
),
900913
),
ST_Transform(geom, 900913)
)
``````

I had to project the data because I am using PostGIS 1.3.4 which doesnt support Geography type.

Basically what I am doing is I am using ST_MakeLine and a query to locate "aerodromes" and return their geometry.

I had to order them (using the CASE directive) so that the line would be connected in the right order.

I then project and buffer this line to create a Polygon that I can then use to see what other aerodromes intersect with the buffered polygon.

Using the unbuffered line (route) and the call ST_Line_Locate_Line I then order the discovered aerodromes as they appear appear along the path.

• All but the last part of your question are answered in (multiple) threads on this site: there are many GPS-oriented questions telling how to make routes out of points; you want to buffer the route (a "bounding box" is altogether a different thing), and you want to select points by polygons (a "point-in-polygon" problem or "spatial selection"). The remaining question boils down to how to sort a set of points by position along a given oriented polyline. Sep 4, 2013 at 18:27
• Replace the pattern ST_Intersects(A, ST_Buffer(B, D)) with ST_DWithin(A, B, D) for a vastly more efficient query. Sep 5, 2013 at 18:01

You should be able to run a SQL statement using the PostGIS function ST_LineLocatePoint to select the relative distance along the route. You will need a single route, a filter for the max distance from the route, and an order by clause.

This is untested, but should get you started:

``````Select    points.pointID
from      points,
rountes
where     rountes.routeID = :routeid
and     ST_Intersects(points.pointShape, ST_Buffer(route.rounteShape, 50000))
order by  ST_LineLocatePoint(route.rounteShape, points.pointShape);
``````

You could filter by ST_Distance(points.pointShape, route.routeShape) <= 50000, but I'm not sure if that is going to use the spatial index.

• Never mind, looks like Jukub beat me to it. Sep 4, 2013 at 20:21
• But certainly your answer is more full, my was kinda cheap. Sep 6, 2013 at 5:01

You need two postgis functions `ST_Buffer` and `ST_Line_Locate_Point`.

• These were the calls I was looking for. Thank you! Sep 5, 2013 at 17:13