st_azimuth geography calculation result mismatch version 3.1.1

[gislars]

I'm doing calculation of directions and with postgis version 3.1.1. I get different or wrong results. It depends on the direction.

I created a test case:

SELECT ST_Azimuth('POINT(0 0)'::geography, 'POINT(0 1)'::geography) AS a_north,
       ST_Azimuth('POINT(0 0)'::geography, 'POINT(1 1)'::geography) AS a_northeast,
       ST_Azimuth('POINT(0 0)'::geography, 'POINT(1 0)'::geography) AS a_east,
       ST_Azimuth('POINT(0 1)'::geography, 'POINT(0 0)'::geography) AS a_north_reverse,
       ST_Azimuth('POINT(1 1)'::geography, 'POINT(0 0)'::geography) AS a_northeast_reverse,
       ST_Azimuth('POINT(1 0)'::geography, 'POINT(0 0)'::geography) AS a_east_reverse
union all
select
       ST_Azimuth('POINT(0 0)'::geometry, 'POINT(0 1)'::geometry) AS a_north,
       ST_Azimuth('POINT(0 0)'::geometry, 'POINT(1 1)'::geometry) AS a_northeast,
       ST_Azimuth('POINT(0 0)'::geometry, 'POINT(1 0)'::geometry) AS a_east,
       ST_Azimuth('POINT(0 1)'::geometry, 'POINT(0 0)'::geometry) AS a_north_reverse,
       ST_Azimuth('POINT(1 1)'::geometry, 'POINT(0 0)'::geometry) AS a_northeast_reverse,
       ST_Azimuth('POINT(1 0)'::geometry, 'POINT(0 0)'::geometry) AS a_east_reverse 

the result:

a_north	a_northeast	a_east	a_north_reverse	a_northeast_reverse	a_east_reverse
3.141592653589793	0.788680084525966	1.5707963267948966	3.141592653589793	5.494352906159104	4.71238898038469
0	0.7853981633974483	1.5707963267948966	3.141592653589793	3.9269908169872423	4.71238898038469

the same result but converted to degrees:

a_north	a_northeast	a_east	a_north_reverse	a_northeast_reverse	a_east_reverse
180	45.18804022935888	90	180	314.80323267835513	270
0	45	90	180	225.00000000000006	270

The first row is displaying the results for geography type, the second for geometry type.

These are the results for postgis version 2.5.2

a_north	a_northeast	a_east	a_north_reverse	a_northeast_reverse	a_east_reverse
0	45.18804022935888	90	180	225.19676732164487	-90
0	45	90	180	225.00000000000006	-90

The value of a_northeast_reverse and a_north value type geography isn't correct, it should be 0 (a_north) or 225 (a_northeast_reverse).

The other values are ok. The difference for a_east_reverse is 360 degree/2pi, the change was introduced in #4718

Am I doing something wrong here? I'm not seeing what the problem might be. I couldn't figure out if it is actual a postgis problem or proj problem. Could someone verify this or help me to get the correct values?

[logi]

I ran into the same problem and traced it back to [4f1fecf2ebc7/git] for fixing #4718.

The final line in lwgeom_azumith_spheroid:

    /* Do the direction calculation */
    az = spheroid_direction(&g1, &g2, spheroid);
    /* Ensure result is positive */
    return az > 0 ? az : M_PI - az;

should read:

    return az >= 0 ? az : 2*M_PI + az;

Subtracting az from π returns a different direction and not a different representation of the direction modulus 2π as intended. The regression test that was added in the same change very luckily hit on the one negative direction where 2π+az = π-az which is az=-π/2 or due West.

Secondly, the comparison needs to be az >= 0 rather than az > 0 so that North is returned as 0 and not 2π (or π with the other bug still present) ​as documented

I am not set up for compiling postgis or testing the results so I'm not submitting a PR with the above. However, in addition to the one-line code change, I propose updating regress/core/tickets.sql with

SELECT '#4888',
        round(degrees(ST_Azimuth('POINT(0 0)'::geography,'POINT(0 0)'::geography))::numeric,3) as same,
        round(degrees(ST_Azimuth('POINT(0 0)'::geography,'POINT(0 1)'::geography))::numeric,3) as N,
        round(degrees(ST_Azimuth('POINT(0 0)'::geography,'POINT(1 1)'::geography))::numeric,3) as NE,
        round(degrees(ST_Azimuth('POINT(0 0)'::geography,'POINT(1 0)'::geography))::numeric,3) as E,
        round(degrees(ST_Azimuth('POINT(0 0)'::geography,'POINT(1 -1)'::geography))::numeric,3) as SE,
        round(degrees(ST_Azimuth('POINT(0 0)'::geography,'POINT(0 -1)'::geography))::numeric,3) as S,
        round(degrees(ST_Azimuth('POINT(0 0)'::geography,'POINT(-1 -1)'::geography))::numeric,3) as SW,
        round(degrees(ST_Azimuth('POINT(0 0)'::geography,'POINT(-1 0)'::geography))::numeric,3) as W,
        round(degrees(ST_Azimuth('POINT(0 0)'::geography,'POINT(-1 1)'::geography))::numeric,3) as NW;

and regress/core/tickets_expected with

#4888|NULL|0.000|45.188|90.000|134.812|180.000|225.118|270.000|314.812

[logi]

I might suggest that this should be more than a medium-priority ticket since it makes the ST_Azimuth function return an entirely wrong value over half its range.

[logi]

OK, this is actually fixed on master. I'm looking at this code which seems to be a bit more sophisticated than my line distinguishing -0 from +0.

    /* Do the direction calculation */
    az = spheroid_direction(&g1, &g2, spheroid);
    /* Ensure result is positive */
    return az < -0 ? 2*M_PI + az : az;
    // return az;

[logi]

See #4840 so this issue seems to be a duplicate of that and the fix is in and has been applied to both master and stable-3.1.

[gislars]

Thank you for your reply and digging in to it. Somehow I missed #4840 but good thing it already got fixed. As a duplicate I guess this issue is solved an can be closed.

[logi]

After building and installing the stable-3.1 branch from ​https://git.osgeo.org/gitea/postgis/postgis.git I get all the expected behavior and our project's unit tests pass again. We'll be testing that today and then deploying deploying to production until 3.1.2 is released with the fix included.

We accidentally updated to 3.1 in production after testing with 3.0 so for a couple of months we've been returning incorrect directions from one of our APIs...

[pramsey]

Milestone renamed

[robe]

In prep for 3.1.3 release

[robe]

Sounds like from what the reporter is saying this is already fixed and was fixed in 3.1.2 so I'm just going to push it back to 3.1.2