5

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: The results. 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?

  • 1
    You never say what the big difference you notice is, but I'd be inclined to suspect that, since you're getting close to half the circumference of the earth that the ST_Segmentize version flips around to the shorter great circle path running in the opposite direction. This is expected, given the way the function is written. – Paul Ramsey Dec 2 '15 at 16:13
  • As far as I know, 19000 km is less than half of the smallest earth's great circle circumference so it is not going the opposite direction. The error is in the range 100-150 km, sorry I was not accurate in my statement. – Redtakfeoh Dec 2 '15 at 16:43
  • 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. – Paul Ramsey Dec 5 '15 at 14:18
3

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.

| improve this answer | |
  • Your statement is right but this is not what I am not doing in the manual case. I always project the start point but at incremented distances. That's why I find it strange... I will add a drawing of the two resulting lines tomorrow. – Redtakfeoh Dec 6 '15 at 21:22
  • Updated above. I tried to find a distance deviation, comparing ST_Length of a raw linestring to a segmentized one, and found almost no differences at all. Are the lengths practically the same, but the paths over land the same in your case? What does the 100km refer to: length, or path difference? – Paul Ramsey Dec 7 '15 at 13:57
  • In my case, I calculated the 100km difference this ways: I iterate on all the point of one of the segmented line and coomputed the distance using ST_Distance. Then I take the max of these values. – Redtakfeoh Dec 8 '15 at 7:17
  • How do you know that it is calculated on the sphere? Is it because it says in the doc for ST_Segmentize that: "Distance computation is performed in 2d only." ? – Redtakfeoh Dec 8 '15 at 7:22
  • I tried using ST_Length on my two lines. If I use the spheroid on the auto line I get the expected distance. So your statement is accurate. It does use the spheroid in that case. However, in both cases, if I don't use the spheroid, I get an length error of 3.6km for the manual line and 4km for the auto line. The total expected distance was 19000km. But I believe in my case, the correct path is the manual line. – Redtakfeoh Dec 8 '15 at 7:23

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.