5

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.

Description of problem

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.

  • 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 '13 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. – Paul Ramsey Sep 5 '13 at 18:01
4

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 '13 at 17:13
8

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. – travis Sep 4 '13 at 20:21
  • 3
    But certainly your answer is more full, my was kinda cheap. – Jakub Kania Sep 6 '13 at 5:01

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.