ST_Azimuth() results

[Antonio Valanzano]

Can someone tell me why Spatialite is giving results different from PostGIS for the following SELECT ?

SELECT 

       degrees(ST_Azimuth( ST_Point(25, 45),  ST_Point(80, 100))) AS deg_A_B,
       degrees(ST_Azimuth( ST_Point(80, 100), ST_Point(25, 45) )) AS deg_B_A; 

-- SpatiaLite 5.1.0

12.826354 244.496474

-- PostGIS 3.5.2

45 225.00000000000006

Antonio

[Jukka Rahkonen]

Hi,

From https://www.gaia-gis.it/gaia-sins/spatialite-sql-5.1.0.html about ST_Azimuth:

"Starting since 4.1.0 if both points supports long/lat coords the returned Azimuth will be precisely computed on the ellipsoid. "

And from https://postgis.net/docs/ST_Azimuth.html
"Enhanced: 2.0.0 support for geography was introduced"

PostGIS gives the same ellipsoidal result than Spatialite with geography.

SELECT

       degrees(ST_Azimuth( ST_Point(25, 45)::geography,  ST_Point(80, 100)::geography)) AS deg_A_B,
       degrees(ST_Azimuth( ST_Point(80, 100)::geography, ST_Point(25, 45)::geography )) AS deg_B_A;

And Spatialite gives the cartesian answer (45 degrees) when the points are defined to be in some projected coordinate system, for example EPSG:3857.

 
SELECT   degrees(ST_Azimuth( MakePoint(25, 45, 3857),  MakePoint(80, 100, 3857))) AS deg_A_B,

[a.fu...@lqt.it]

On Sun, 23 Feb 2025 04:46:04 -0800 (PST), Jukka Rahkonen wrote:
> Hi,
>
> From https://www.gaia-gis.it/gaia-sins/spatialite-sql-5.1.0.html
> about
> ST_Azimuth:
> "Starting since 4.1.0 if both points supports long/lat coords the
> returned Azimuth will be precisely computed on the ellipsoid. "
>
> And from https://postgis.net/docs/ST_Azimuth.html
> "Enhanced: 2.0.0 support for geography was introduced"
>
> PostGIS gives the same ellipsoidal result than Spatialite with
> geography.
>

> ----------------------------- <snip> ----------------------------

>
> And Spatialite gives the cartesian answer (45 degrees) when the
> points
> are defined to be in some projected coordinate system, for example
> EPSG:3857.
>  

Hi Jukka,

congratulations, you hit the mark.

I was just investigating the problem when your email arrived, and
until that moment I had come to understand two things based on my
very first checks on both Linux and Windows:

A) the result of SpatiaLite's ST_Azimuth() depends on the SRID;
sometimes they are identical to those reported for PostGIS,
but other times they are completely different.
B) even more surprisingly, the results change depending on the
connected DB.

After a more careful check on the source code, here is the
complete explanation.

1. There are two different ways to calculate the azimuth;
on the plane or on the ellipsoid.

2. Therefore the very first action of ST_Azimuth() is to
look for the geodesic parameters of the SRID defined
in the proj4text column of the spatial_ref_sys table.

3. At this point it evaluates the geodesic string to
understand whether it's a planar (projected) or
geographic type SRS, and in this second case it
also retrieves the lengths of the two axes of
the ellipsoid.
Note: if it cannot find the geodesic definition of the
SRID a planar system is assumed.

4. Finally everything is ready to calculate the azimuth
with the most appropriate method.

Drawing the final conclusions:

- In the absence of any explicit indication the default SRID=0
(Undefined - Geographic Long/Lat) will be assumed

- To get the same "planar" results as PostGIS you must be
careful to explicitly declare the SRID=-1
(Undefined - Cartesian)

- But if you use an unknown SRID, or if the DB metadata
are missing or only partial, then you still fall into
the "planar" case.

- Naturally by always using explicit SRID attributions
on a correctly initialized DB, then there will never
be any problems.
Blindly relying on the defaults, on the other hand,
can easily create great confusion. Be careful.

bye Sandro

[Antonio Valanzano]

Dear Jukka and Sandro,

thanks for the detailed explanation.

The problem was created by me using ST_Point() which  does not support the SRID argument.

However there is a lack of information for ST_Point() because it is not clear what type of coordinate are those expressed by this function.

As Sandro explained, the default SRID=0 (Undefined - Geographic Long/Lat) is assumed but never declared (at least I did not find this information).

So if someone wants to have control of the results of ST_Azimuth he has to use ST_MakePoint() and to choose the correct SRID

SRID = 0 per geographic calculation

SRID = -1 per cartesian calculation

or a valid EPSG code

Antonio

--
You received this message because you are subscribed to the Google Groups "SpatiaLite Users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to spatialite-use...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/spatialite-users/499fdbcb491b772ddfb12f42b0ca9865%40lqt.it.

[a.fu...@lqt.it]

On Sun, 23 Feb 2025 18:01:31 +0100, Antonio Valanzano wrote:
> However there is a lack of information for ST_Point() because it is
> not clear what type of coordinate are those expressed by this
> function.
>

Hi Antonio,

The documentation about ST_Point() clearly states that it's simply
an alias name for MakePoint(), except that it doesn't accept the
SRID argument.

It would therefore seem quite clear that the geometries created by
ST_Point() will always refer to the default SRID universally adopted
by all SpatiaLite functions, which is always SRID=0

IMHO repeating this information for the many functions (and there
are really a lot of them) that create brand new geometries without
deriving them from a previous ones would make consulting the
documentation extremely boring and unnecessarily cumbersome.

> As Sandro explained, the default SRID=0 (Undefined - Geographic
> Long/Lat) is assumed but never declared (at least I did not find this
> information).
>

It was announced when v.4.0 was released on 2012-09-05; a dozen
years have passed and it has always remained unchanged.
I would take it for granted that this is now a fact well known
to any SpatiaLite user.

bye Sandro