ST_Project vs ST_Segmentize Postgis

Question

I noticed some inconsistency on geographies between ST_Project and ST_Segmentize in postgis when dealing with trying to segmentize over very long distances (around 18000-19000 km)

Below are two plpgsql functions that should give about the same result but they don't. I tried it using different versions of Postgis (up to 2.1.8).

For both functions I use a start point lat and long, an azimuth, a total length and a segment length.

1. segmentize_line_manually: I project the starting point at incremented distances (segment length) with the given azimuth using ST_Project.

2. segmentize_line_auto: I project the starting point at the total distance with the given azimuth, then I call ST_Segmentize with the segment length.

Around 18000-19000 km I notice a rather important difference between the two lines. In my opinion, the segmentize_line_manually gives the correct result.

CREATE OR REPLACE FUNCTION segmentize_line_manually(startLat DECIMAL, 

startLong DECIMAL, azimuth DECIMAL, max_distance_meter INT, interval_meter INT)
RETURNS geometry AS
$$
DECLARE
    line_points geometry[];
    startPoint geography := ST_MakePoint(startLong, startLat);
    distance_meter INT := 0;
    i INT := 1;
BEGIN
    WHILE distance_meter < max_distance_meter LOOP
        distance_meter := distance_meter + interval_meter;
        line_points[i] = geometry(ST_Project(startPoint, distance_meter, azimuth));
        i := i + 1;
    END LOOP;

    RETURN ST_MakeLine(line_points);
END;
$$
LANGUAGE 'plpgsql' IMMUTABLE;

CREATE OR REPLACE FUNCTION segmentize_line_auto(startLat DECIMAL, startLong DECIMAL, azimuth DECIMAL, max_distance_meter INT, interval_meter INT)
RETURNS geometry AS
$$
DECLARE
    line_points geometry[];
    azimuth_line geometry;
BEGIN
    line_points[1] := geography(ST_MakePoint(startLong, startLat));
    line_points[2] := ST_Project(line_points[1], max_distance_meter, azimuth);

    azimuth_line := ST_MakeLine(line_points);
    RETURN st_segmentize(geography(azimuth_line), interval_meter);
END;
$$
LANGUAGE 'plpgsql' IMMUTABLE;

I call both methods using the following parameters:

SELECT segmentize_line_****(0, 0, 0.628973508651991, 19000000, 50000);

The results below. In red, the line created using ST_Project, in blue with ST_Segmentize:

Is there a problem with my code, or is it really an issue with ST_Segmentize?

Answer 1

You're building two rather different lines, I think.

Your "manual" line is a series of rhumb lines between points, with an unchanging bearing all the way. Your "auto" line is an orthodrome (great circle) between two points.

I'd actually expect, that if you plot them on a map you'll find they not only differ in length a little, they also go to two very different parts of the world.

In words: the place you get to by holding a constant bearing of 90 for 1000km is not the same as the place you get to by bearing 90 at the start and then travelling 1000km without turning.

In pseudo code, this:

orient(west)
move(1000km)

Is not the same as this:

for i in 1..1000
   orient(west)
   move(1km)

(Unless you start on the equator.)

Update: Per the comment, then the issue is spheroid vs sphere. The ST_Segmentize version will be generating a great circle on a sphere, while the manual process will be working on a spheroid.