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.


  • 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.

Description of problem


  • 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)


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
      (SELECT ST_Buffer(
              (SELECT ST_MakeLine(geom) 
                  (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
          50000) as geo
      ) as line
    WHERE ST_Intersects(ST_Transform(geom, 900913), line.geo)
    ORDER BY ST_Line_Locate_Point(
          (SELECT ST_MakeLine(geom) 
              (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
      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.

  • 2
    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.
    – whuber
    Sep 4, 2013 at 18:27
  • 1
    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

2 Answers 2


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, 
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.
    – travis
    Sep 4, 2013 at 20:21
  • 3
    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!
    – jreid42
    Sep 5, 2013 at 17:13

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.