Friday 27 February 2015

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




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




  • 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

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.



Answer



You need two postgis functions ST_Buffer and ST_Line_Locate_Point.


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...