[postgis-tickets] [PostGIS] #4888: st_azimuth geography calculation result mismatch version 3.1.1

[PostGIS]

#4888: st_azimuth geography calculation result mismatch version 3.1.1
----------------------+---------------------------
  Reporter:  gislars  |      Owner:  pramsey
      Type:  defect   |     Status:  new
  Priority:  medium   |  Milestone:  PostGIS 3.1.2
 Component:  postgis  |    Version:  3.1.x
Resolution:           |   Keywords:
----------------------+---------------------------

Comment (by logi):

 I ran into the same problem and traced it back to 4f1fecf2ebc7 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) [[https://postgis.net/docs/ST_Azimuth.html|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 '#4718',
         round(degrees(ST_Azimuth('POINT(77.46412
 37.96999)'::geography,'POINT(77.46409 37.96999)'::geography))::numeric,3),
         round(degrees(ST_Azimuth('POINT(77.46412
 37.96999)'::geography,'POINT(77.46429 37.96999)'::geography))::numeric,3),
         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
 {{{
 #4718|270.000|90.000|NULL|0.000|45.188|90.000|134.812|180.000|225.118|270.000|314.812
 }}}

-- 
Ticket URL: <https://trac.osgeo.org/postgis/ticket/4888#comment:2>
PostGIS <http://trac.osgeo.org/postgis/>
The PostGIS Trac is used for bug, enhancement & task tracking, a user and developer wiki, and a view into the subversion code repository of PostGIS project.