Postgis i problem z wyliczeniem współrzędnych przy zadaniu dystansu i kąta

Witajcie,

Nigdy wcześniej nie bawiłem się z postgisem i pewnie z tego wynika moja niewiedza i nieprawidłowe działanie. Mam jednak nadzieje, że mnie oświecicie.

W postgisie jest taka funkcja
https://postgis.net/docs/ST_Project.html

jak widać można za jej pomoca wyliczyć współrzędne punktu od punktu początkowego na zadanym dystansie i kącie. Innymi słowy stojąc w wawie chcecie zobaczyć punkt 100 km na północ od wawy to ST_Project zwróci współrzędne tej lokalizacji.

W teorii super, a w praktyce ... już pokazuję:

SELECT ST_AsText(ST_Project('POINT(52 21)'::geography, 1500, radians(0)));

zwraca

POINT(52.00000000000001 21.013548040821505)

tu od razu wtf przecież 0 stopni to 0 radianów i wg dokuentacji powinno być względem północy?

The azimuth (also known as heading or bearing) is given in radians. It is measured clockwise from true north (azimuth zero). East is azimuth π/2 (90 degrees); south is azimuth π (180 degrees); west is azimuth 3π/2 (270 degrees). Negative azimuth values and values greater than 2π (360 degrees) are supported.

dobra może coś źle rozumiem ale sprawdźmy dystans tego punktu

jak k..a 900 jak podałem 1500 m...

dodam, że im większy dystans tym większy mam rozjazd. Przy 100 km było to prawie 16 km różnicy.

Przypuszczam, że problemem jest wyliczanie tych współrzędnych ponieważ użyłem strony
https://www.fcc.gov/media/radio/find-terminal-coordinates
podałem ten same parametry:

i dostałem:

przy wrzuceniu na mapę wszystko się zgadza:

(te 13 m to dlatgo, że trzeba ręcznie kliknąć na mapę - lub ja nie wiem jak to wyliczyć bez klikania)

Macie pomysł dlaczego postgis mi takie głupoty zwraca?

a wez uzyj https://postgis.net/docs/manual-1.4/ST_Distance_Spheroid.html i powiedz co wyszlo

@masterc: a krzyczy mi, że nie ma takiej funkcjo na ST_Distance_Spheroid

ERROR:  function st_distance_spheroid(geometry, geometry, unknown) does not exist
LINE 2:   ST_Distance_Spheroid(ST_Centroid(the_geom), ST_GeomFromTex...
          ^
HINT:  No function matches the given name and argument types. You might need to add explicit type casts.

to stworz ja na chwile i zobacz czy dziala

CREATE FUNCTION `st_distance_sphere`(`pt1` POINT, `pt2` POINT) RETURNS 
    decimal(10,2)
    BEGIN
    return 6371000 * 2 * ASIN(SQRT(
       POWER(SIN((ST_Y(pt2) - ST_Y(pt1)) * pi()/180 / 2),
       2) + COS(ST_Y(pt1) * pi()/180 ) * COS(ST_Y(pt2) *
       pi()/180) * POWER(SIN((ST_X(pt2) - ST_X(pt1)) *
       pi()/180 / 2), 2) ));
    END

Ogólnie poradziłem sobie nieco inaczej wykorzystując plpython :) Zostawiam dla potomnych

CREATE TYPE sd.tlatlon AS (
  lat VARCHAR,
  long VARCHAR
);

COMMENT ON TYPE sd.tlatlon
IS 'Typ dla LAT i LONG';

CREATE OR REPLACE FUNCTION sd.fcalclatlon (
  abearing double precision,
  adistance double precision,
  alat double precision,
  alon double precision
)
RETURNS sd.tlatlon AS
$body$
import math

R = 6378.1 #Radius of the Earth
brng = abearing #Bearing degrees converted to radians.
d = adistance #Distance in km

lat1 = math.radians(alat) #Current lat point converted to radians
lon1 = math.radians(alon) #Current long point converted to radians

lat2 = math.asin( math.sin(lat1)*math.cos(d/R) +
     math.cos(lat1)*math.sin(d/R)*math.cos(brng))

lon2 = lon1 + math.atan2(math.sin(brng)*math.sin(d/R)*math.cos(lat1),
             math.cos(d/R)-math.sin(lat1)*math.sin(lat2))

lat2 = math.degrees(lat2)
lon2 = math.degrees(lon2)

return [(lat2),(lon2)]
$body$
LANGUAGE 'plpython3u'
VOLATILE
CALLED ON NULL INPUT
SECURITY INVOKER
PARALLEL UNSAFE
COST 100;