Elevation Profiles and Flightlines with PostGIS
A community member on the postgis-users mailing list had a question recently:
I have a table of elevation points, and I would like to figure out an elevation profile for a flightline running through those points. How?
This question is a nice showcase of some of my favorite spatial tools with indexing, point to point distance queries on a sphere, and nearest neighbor queries. I thought it would make a great post. The original question author was nice enough to share his elevation data, so I can show you some pictures and an approach.
Elevation Points every 2KM
The data set is 1.5M elevation points in Europe, one point for every 2km. When plotted on a map you can see this is quite dense.
The table structure is simple.
Column | Type | Collation | Nullable | Default
-----------+------------------+-----------+----------+---------
latlng | geography | | |
lat | double precision | | |
lng | double precision | | |
elevation | real | | |
All we need is a spatial index on the latlng column to make our nearest
neighbor searches fast.
CREATE INDEX elevation_latlng_x ON elevation USING GIST (latlng)
Flightlines
Flightlines are different from overland lines because they are not only "as the crow flies", but "as the crow flies over a sphere". Over longer distances, and at higher latitudes, the effect is more visible, but great circle routes are not quite the same as point-point routes.
For example, a straight line between London and Warsaw can be constructed like this:
SELECT ST_MakeLine(ST_MakePoint(-0.061,51.482), -- London
ST_MakePoint(21.252,52.150)) -- Warsaw
But a straight line (shortest path over a plane) is not what we want, we want a
great circle (shortest path over a sphere). So, we flip the line into the
geography type, which uses great circles to interpolate between vertices, and
then ST_Segmentize the line, to
add new vertices between the endpoints.
SELECT ST_Segmentize(geography(
ST_MakeLine(
ST_MakePoint(-0.061,51.482),
ST_MakePoint(21.252,52.150))
), 4000)
The red line is the great circle, the blue line is the original cartesian line.
Nearest to a Point
For the final query, we will construct a flight path between Madrid and Guadalajara, two cities close enough together that we can see the actual elevation points. We'll put a vertex every 4KM along that path, and find the closest elevation point to each vertex.
To find the closest point, we use the <->
"index ordered distance"
operator. When placed in the ORDER BY clause of a query <-> not only
calculates the distance between the left and right operands, it also returns the
distances starting from the minimum possible distance.
So, the query for the nearest elevation point to Madrid (-3.6996,40.4009) is:
SELECT
e.elevation,
ST_AsText(e.latlng, 3) AS coords,
geography(ST_MakePoint(-3.6996,40.4009)) <-> e.latlng AS dist
FROM elevation e
ORDER BY dist
LIMIT 1
Which returns (in 2ms, because of the index assist):
elevation | coords | dist
-----------+----------------------+--------------------
650.8142 | POINT(-3.707 40.409) | 1045.4185642910293
Nearest to Flight Path
The index-assisted nearest neighbor operator <-> is an odd one, as it
requires a literal value on one side of the operator. It cannot take in two
variable arguments, so how do we join up our flight path, which will have more
than one point, to the elevation grid?
The answer is to use a LATERAL JOIN to control the iteration and supply the flight path points, one at a time, to the nearest neighbor calculation.
Taken in chunks, the query first builds the flight path.
WITH path AS (
SELECT ST_Segmentize(geography(
ST_MakeLine(
ST_MakePoint(-3.6996,40.4009), -- Madrid
ST_MakePoint(-3.1675,40.6582)) -- Guadalajara
), 4000) AS geog
)
Then it takes apart the path into one record for each vertex using ST_DumpPoints
points AS (
SELECT (ST_DumpPoints(geog::geometry)).geom FROM path
)
Finally, it does a lateral join to the elevation points, finding the nearest elevation for each path point.
SELECT
d.elevation,
round(d.dist) AS distance,
ST_AsText(d.latlng, 3) AS elevation_pt,
ST_AsText(points.geom, 3) AS path_pt
FROM points
CROSS JOIN LATERAL (
SELECT
e.elevation, e.latlng,
points.geom::geography <-> e.latlng AS dist
FROM elevation e
ORDER BY dist
LIMIT 1
) d
The final result take about 14ms to run, and finds an elevation for every point on the flightline.
elevation | distance | elevation_pt | path_pt
-----------+----------+----------------------+----------------------
650.8142 | 1045 | POINT(-3.707 40.409) | POINT(-3.7 40.401)
630.8268 | 563 | POINT(-3.66 40.415) | POINT(-3.666 40.417)
701.6845 | 728 | POINT(-3.641 40.436) | POINT(-3.633 40.433)
632.26556 | 882 | POINT(-3.595 40.443) | POINT(-3.6 40.449)
611.845 | 733 | POINT(-3.575 40.464) | POINT(-3.567 40.465)
607.7102 | 701 | POINT(-3.533 40.488) | POINT(-3.534 40.482)
644.5825 | 1057 | POINT(-3.509 40.491) | POINT(-3.501 40.498)
637.69604 | 166 | POINT(-3.467 40.515) | POINT(-3.467 40.514)
666.1685 | 1335 | POINT(-3.447 40.536) | POINT(-3.434 40.53)
673.436 | 373 | POINT(-3.401 40.543) | POINT(-3.401 40.546)
717.8294 | 938 | POINT(-3.358 40.567) | POINT(-3.368 40.562)
696.11383 | 910 | POINT(-3.335 40.57) | POINT(-3.334 40.578)
709.38855 | 723 | POINT(-3.292 40.594) | POINT(-3.301 40.594)
716.1025 | 697 | POINT(-3.273 40.615) | POINT(-3.268 40.61)
655.39514 | 862 | POINT(-3.226 40.621) | POINT(-3.234 40.626)
659.81757 | 501 | POINT(-3.207 40.642) | POINT(-3.201 40.642)
654.524 | 943 | POINT(-3.164 40.666) | POINT(-3.168 40.658)
To summarize - my recommended approach to calculating elevation profiles and flightlines is:
- create a spatial index to make the nearest neighbor searches fast
- query between lines using the geography type so the shortest path is over a sphere not a straight line
- use the index ordered distance operator to find the shortest distance between two points
- use the nearest neighbor operator with a lateral join to elevation points.
