Witam,
Od dłuższego czasu poszukuje informacji jak zamienić 2 punkty współrzędnych geograficznych na odległość między nimi.
Myślę, że przyjęcie, że ziemia jest kulą było by mało dokładne, a nie widziałem wzorów z zastosowaniem elipsoidy.
0
0
Mam nadzieję, że Ci pomoże
Trygonometria sferyczna
0
Jesli chcesz miec bardzo dokladne wyniki to - wg mnie oczywiscie - musialbys uzyc calek (tj zcalkowac po wspolrzednych od pierwszego do drugiego punktu wzor na droge na okreslonej szerokosci/dlugosci geograficznej). Nie wiem czy daloby sie analitycznie jakos to rozwiazac bo jesli nie to pozostaje calkowanie numeryczne (woooolne i troche klepania jest).
Duzo lepiej jest zastosowac przyblizenie, ze ziemia jest kula :)
0
Analitycznie by się nie dało, nawet odległości między biegunami (=połowa obwodu elipsy) nie da się wyliczyć analitycznie. Pojawi się całka od zera do pi z sqrt(a2cos2(t)+b2cos2(t)) która jest nieobliczalna analitycznie dla a różnego od b
0
Ja liczę tak:
function Dist(Lat1,Lon1,Lat2,Lon2:Double):Double;
const dz=12756.274;//średnica Ziemi na równiku [km]
var a,b:Double;
begin
a:=(Lon2-Lon1)*Cos(Lat1*pi/180);
b:=(Lat2-Lat1);
Result:=Sqrt(a*a+b*b)*pi*dz/360;//[km]
end;
Współrzędne geogr. w stopniach
0
@pelsta, to jest na elipsoidzie czy na sferze ?
0
myślałem o tym, żeby zrobić to na kuli co będzie myślę wystarczające dla małych odległości, jednak przy wiekszych wystąpią duże przekłamania
0
Wg. mnie na sferze.
Wzór ten używam do obliczania przebytej odległości na podstawie tzw. trackloga - zapis trasy podróży, generowanego przez urządzenia do nawigacji (konkretnie Garmin Nuvi 200).
0
Dla sfery będzie dużo prościej. Jak chcesz oszacować błąd, to oblicz odległość między biegunami:
dla sfery = pi*promień
dla elipsoidy = całka którą napisałem wcześniej, gdzie a = promień biegunowy, b = promień równikowy
0
Jeśli zdecydujesz się na elipsoidę, to napiszę ci odpowiednie całki, które sobie porozwiązujesz numerycznie.
pozdrawiam
0
@adf88 mylisz się, tych całek nie można obliczyć analitycznie.
// po wysłaniu powyższego na serwer zobaczyłem, że post poprzednika się ulotniłł
0
no niestety dopiero w sobotę będę mógł kontynuować pisanie kodu także zrobię dla sfery i sprawdzę niedokładność w obliczaniu odległości.
0
@adf88, o ile dobrze zapamiętałem twój "znikający post", to napisałeś że najkrótszą drogą łącząca dwa punkty na elipsoidzie obrotowej E jest łuk elipsy będącej przekrojem elipsoidy E i płaszczyzny przechodzącej przez środek elipsoidy E. Wyczytałeś to gdzieś, czy tak ci się zdawało (przez analogię ze sferą) ? Próbowałem policzyć i wyszło mi, że opisana wyżej elipsa, nie jest geodezyjną, nie może więc być najkrótszą drogą. Pytam ponieważ nie mam pełnego zaufania do własnych rachunków.
pozdrawiam
0
Heh, no właśnie skasowałem swój post bo stwierdziłem, że to nie musi być prawda. Dowodu jednak nie mam.
A jak to liczyłeś ?
0
Elipsa będąca przekrojem elipsoidy obrotowej o "promieniu" równikowym b i "promieniu" biegunowym a oraz płaszczyzny o równaniu z=my ma równanie parametryczne a(t) = (bcos(t),((ba)/sqrt(m2+a2))sin(t),m((ba)/sqrt(m2+a2))*sin(t))
Jest to opis naturalny, tzn. wektor a'(t) ma stałą długość lub inaczej wektory a'(t) oraz a''(t) są prostopadłe.
Zatem warunek geodezyjności (składowa wektora a(t) prostopadła do a'(t) jest prostopadła do powierzchni) przyjmuje prostszą postać: wektor a(t) jest prostopadły do powierzchni. A nie jest. Nawet dla zwykłej elipsy na płaszczyźnie nie jest prawdą, że wektor wodzący punktu na elipsie jest prostopadły do elipsy (z wyjątkiem wierzchołków elipsy).
0
hmm na moje oko jezeli uproscisz ziemie do sfery, to mozesz wyznaczyc dwa wektory normalne powierzchni stycznej do sfery w punktach a i b, policzyc kat miedzy nimi i na tej podstawie obliczyc dlugosc luku o promieniu rownym prominiowi ziemi :P
0
@cepa masz całkowicie rację, ja też bym tak liczył. Autor wątku chce jednak "liczyć dokładniej" przyjmując, że kula ziemska jest elipsoidą obrotową.
0
Sprawdziłem różne sposoby liczenia odległości punktów na kuli ziemskiej, dla uproszczenia przyjąłem że punkty A i B leżą na jednym południku - są wyznaczone przez podanie szerokości geograficznej (0-90).
Stosowałem pięć sposobów:
- Punkty A i B leżą na kuli o promieniu Rr=6378,245 km (promień równikowy)
- Punkty A i B leżą na kuli o promieniu Rb=6356,863 km (promień biegunowy)
- Punkty A i B leżą na kuli o promieniu (Rr+Rb)/2
- Kula ziemska jest elipsoidą obrotową, z równania elipsoidy i szerokości geograficznych wynikają
odległości r1 i r2 punktów A i B od środka kuli ziemskiej. Odległość między A i B liczę zakładając, że
punkty te leżą na kuli o promieniu (r1+r2)/2 - Liczę długość łuku elipsy korzystając z definicji długości krzywej: jest to kres górny długości łamanych
których wierzchołki leżą na krzywej (zaczynam od łamanej jednoodcinkowej, w kolejnym kroku
zwiększam dwukrotnie ilość segmentów łamanej, kończę gdy różnica między kolejnymi wynikami jest
< 0,00001.
Przykładowe wyniki:
00,000 90,000 10018,924 09985,337 10002,130 10002,130 10002,137
00,000 01,000 00111,321 00110,948 00111,135 00111,321 00111,321
89,000 90,000 00111,321 00110,948 00111,135 00110,948 00110,948
30,000 60,000 03339,641 03328,446 03334,043 03334,033 03334,048
30,000 45,000 01669,821 01664,223 01667,022 01667,715 01667,740
05,000 45,000 04452,855 04437,928 04445,391 04449,057 04449,812
Morał: Spośród sposobów nie liczących długości łuku elipsy zdecydowanie najlepszy jest sposób 4).
Należy zatem liczyć tak (zmodyfikowany wzór Pelsty):
function Dist(Lat1,Lon1,Lat2,Lon2:Double):Double;
const dr=12756.490;//średnica Ziemi na równiku [km], inna niż u Pelsty - taką znalazłem (nie mierzyłem osobiście)
const db=12717,776;//średnica biegunowa [km]
var a,b,d1,d2,tg,dz:Double;
begin
if (Lat1==90) then d1:=db
else
begin
tg:=Tan(Lat1);
d1:=dr*db*Sqrt(1+tg*tg)/Sqrt(db*db+dr*dr*tg*tg)
end;
if (Lat2==90) then d2:=db
else
begin
tg:=Tan(Lat2);
d2:=dr*db*Sqrt(1+tg*tg)/Sqrt(db*db+dr*dr*tg*tg)
end;
dz:=(d1+d2)/2;
a:=(Lon2-Lon1)*Cos(Lat1*pi/180);
b:=(Lat2-Lat1);
Result:=Sqrt(a*a+b*b)*pi*dz/360;//[km]
end;
0
bogdans napisał(a)
wzór Pelsty
Ten wzór znalazłem gdzieś w Internecie, nie przypisuję sobie do niego żadnych praw!
Średnicę Ziemi znalazłem w Wikipedii. Nie wiem na ile jest to dobra wartość.
W Delphi należały by napisać:
...
if (Lat1=90) then d1:=db;
...
if (Lat2=90) then d2:=db;
...
W wolnym czasie przetestuję Twoje rozwiązanie. Myślę, że przy małych odległościach nie będzie żadnych różnic.
0
Będą, jeżeli oba punkty leżą blisko biegunów.
Wzory i twierdzenia nie zawsze wiązane są z nazwiskiem rzeczywistego odkrywcy. Spróbuję wspomniany wzór wylansować jako wzór Pelsty. ;-)
Pisząc wzór zapomniałem, że poniżej równika jest półkula południowa, w konsekwencji szerokość geograficzna może być ujemna. Aby wzór objął półkulę południową trzeba wprowadzić takie zmiany:
if (Abs(Lat1)==90) ...
if (Abs(Lat2)==90) ...
tg:=Abs(Tan(Lat1));
tg:=Abs(Tan(Lat2));
0
Tutaj http://www.movable-type.co.uk/scripts/latlong.html rozgryzają podobny problem.
A w Twojej funkcji są błędy formalne.
0
Jakie ?
Funkcja była pisana w javie i w niej kompilowana. Na forum zmodyfikowałem twoją funkcję dlatego użyłem Delphi. Nie chciało mi się uruchamiać Delphi, a zakładałem że każdy sobie poradzi jeśli zrobię jakiś błąd formalny.
0
To, co zauważyłem:
const db=12717,776;*średnica biegunowa [km] -> const db=12717.776;*średnica biegunowa [km]
if (Lat1==90) then d1:=db -> if (Lat1=90) then d1:=db
tg:=Tan(Lat1); -> Lat1 trzeba przeliczyć na radiany (w Javie nie było potrzeby?)
if (Lat2==90) then d2:=db -> if (Lat2=90) then d2:=db
tg:=Tan(Lat2); -> Lat2 przeliczyć na radiany
A generalnie, to nie rozumiem Twoich przeliczeń.
0
Dzięki.
W Javie też była konieczność przeliczenia stopni na radiany.
Przeliczenie bierze się stąd, iż odległość punktu A o danej szerokości geograficznej szer od środka kuli ziemskiej zależy od tej szerokości. W szczególach:
równanie elipsoidy ziemskiej, to (x/Rr)2+(y/Rr)2+(z/Rb)2=1,
znajduję punkt przecięcia P tej elipsoidy z prostą o równaniu z=x*tg(szer) i przyjmuję, że odległość punktu A od środka jest identyczna z odległością punktu P od środka. Stąd się biorą wzory przeliczające.
Jak pisałem wcześniej zrobiłem program, który liczył odległość (dla punktów leżących na jednym południku) twoim wzorem przy czterech różnych sposobach wyznaczenia "średnicy" kuli ziemskiej, jednocześnie liczył przyjmując, że kula ziemska jest elipsoidą (długość łuku elipsoidy liczyłem jako kres górny łamanych wpisanych w ten łuk).
Przyjęcie, że "średnica" jest średnicą biegunową, średnicą równikową lub ich średnią arytmetyczną prowadziło du dużego błędu. Błąd zależał bardzo mocno od położenia punktów, np. gdy punkty leżały blisko równika, to drugi i trzeci sposób były dobre, pierwszy nie. Natomiast czwarty nieco skomplikowany sposób dawał dużą dokładność przy każdym położeniu punktów.
0
Uruchomiłem Delphi i skompilowałem funkcje:
function radius(szer:Double):Double;
const Rr=6378.245;
const Rb=6356.863;
var tg: Double;
begin
if (Abs(szer)=90) then Result:=Rb
else
begin
tg:=Tan(GradToRad(Abs(szer)));
Result:=Rr*Rb*Sqrt(1+tg*tg)/Sqrt(Rb*Rb+Rr*Rr*tg*tg);
end;
end;
function distance(szer1,dlug1,szer2,dlug2,promien:Double):Double;
var a,b:Double;
begin
a:=(dlug2-dlug1)/Cos(szer1*pi/180);
b:=szer2-szer1;
Result:=Sqrt(a*a+b*b)*pi*promien/180;
end;
//Edit Podczas naprawiania roweru (zbawienny wpływ pracy fizycznej na pracę umysłu) olśniło mnie że opisany wyżej sposób jest dobry dla punktów leżących na jednej półkuli - dla takich punktów był testowany. Dla punktów leżących na różnych półkulach jest gorzej, jeśli np. będziemy liczyć tą metodą odległość biegunów, to dostaniemy pi*Rb a odległość jest równa połowie obwodu elipsy o półosiach Rr i Rb. Wzory podam później, idea jest taka by dla danych punktów A i B na sferze leżących na różnych półkulach dołączyć odpowiedni punkt pośredni S leżący na równiku i liczyć tak: d(A,B)=d(A,S)+d(S,B).
Ze słownika idiomów: Być raz na wozie, raz pod wozem: To work as a car mechanic
0
Oto kod źródłowy programu obliczającego odległośc na podstawie współrzędnych geograficznych
napisany na podstawie dwóch algorytmów zamieszczonych w Jean Meeus "Astronomical Algorithms"
http://chomikuj.pl/fcvreqnynw/DISTAN~1.C
#include<stdio.h>
#include<conio.h>
#include<math.h>
#define sqr(a) ((a)*(a))
#define a 6378.137
#define f ((21.385)/(6378.137))
#define Rm 6371
double deg2dms(double deg);
double dms2deg(double dms);
double rad2deg(double rad);
double deg2rad(double deg);
double dms2rad(double dms);
double distance01(double,double,double,double);
double distance02(double,double,double,double);
int main(){
char ch;
double phi[2],theta[2];
clrscr();
do{
printf("Podaj wspolrzedne geograficzne punktu 1 \n");
scanf("%lf %lf",&phi[0],&theta[0]);
printf("Podaj wspolrzedne geograficzne punktu 2 \n");
scanf("%lf %lf",&phi[1],&theta[1]);
printf("Przyblizona odleglosc miedzy punktami to : %.10lf\n",distance01(phi[0],theta[0],phi[1],theta[1]));
printf("Przyblizona odleglosc miedzy punktami to : %.10lf\n",distance02(phi[0],theta[0],phi[1],theta[1]));
ch=getch();
}
while(ch!=27);
return 0;
}
double rad2deg(double rad){
return 180*rad/M_PI;
}
double deg2rad(double deg){
return deg*M_PI/180;
}
double deg2dms(double deg){
double m,s;
m=((deg-(int)deg)*60);
s=((m-(int)m)*60);
return ((int)(deg)+m/100.0+s/10000.0);
}
double dms2deg(double dms){
double d,m,s;
d=(int)(dms);
m=(int)(100*(dms-d));
s=100*((100*(dms-d)-m));
return (int)(d)+m/60.0+s/3600.0;
}
double dms2rad(double dms){
double min,sec,r;
min=100*(dms-(int)(dms));
sec=100*(min-(int)(min));
r=(int)dms+(int)(min)/60.0+sec/3600.0;
return M_PI*r/180.0;
}
double distance01(double phi0,double theta0,double phi1, double theta1){
double phi[2],theta[2];
double d,q;
phi[0]=deg2rad(dms2deg(phi0));
phi[1]=deg2rad(dms2deg(phi1));
theta[0]=deg2rad(dms2deg(theta0));
theta[1]=deg2rad(dms2deg(theta1));
q=rad2deg(acos(sin(phi[0])*sin(phi[1])+cos(phi[0])*cos(phi[1])*cos(theta[0]-theta[1])));
d=(2*M_PI*Rm*q)/360;
return d;
}
double distance02(double phi1,double L1,double phi2,double L2){
double F,G,lambda,S,C,omega,R,D,s,H1,H2;
F=(dms2rad(phi1)+dms2rad(phi2))/2.0;
G=(dms2rad(phi1)-dms2rad(phi2))/2.0;
lambda=(dms2rad(L1)-dms2rad(L2))/2.0;
S=sqr(sin(G))*sqr(cos(lambda))+sqr(cos(F))*sqr(sin(lambda));
C=sqr(cos(G))*sqr(cos(lambda))+sqr(sin(F))*sqr(sin(lambda));
if(C!=0.0) omega=atan(sqrt(S/C));
else omega=M_PI/2.0;
R=sqrt(S*C)/omega;
D=2.0*omega*a;
H1=((3.0*R-1.0)/(2.0*C));
H2=((3.0*R+1.0)/(2.0*C));
s=D*(1+f*H1*sqr(sin(F))*sqr(cos(G))-f*H2*sqr(cos(F))*sqr(sin(G)));
return s;
}
0
D = R * 2 * asin(sqrt((sin((B1-B2)/2))2 + cos(B1)cos (B2)(sin((L1-L2)/2))2))
D: szukana odległość
B1, B2: szerokości geogr. punktów
L1, L2: długości geogr. punktów
promień Ziemi; Ziemia nie jest kulą, więc nie istnieje coś takiego jak jeden „jedynie słuszny” promień Ziemi — stosuje się m.in. następujące wartości:
Rn = 6366,710 km (definicja mili morskiej)
Rf = 6371,000 km (dawny standard FAI)
Rmax = 6378,137 km (maksymalny w WGS-84)
Rmin = 6356,752 km (minimalny w WGS-84)
w celu wyeliminowania wpływu niekulistości Ziemi można zastosować wzór:
R = Rmax-(Rmax-Rmin)*sin((B1+B2)/2)
Macie Jajogłowi :)
0
Kanciastogłowy (tzn. @oaieuy) niczego nowego nie napisałeś. Natomiast zaproponowana przez Ciebie poprawka
R = Rmax-(Rmax-Rmin)*sin((B1+B2)/2)
jest bez sensu. Przy liczeniu odległości między biegunami (niekulistość Ziemi ma wtedy największy wpływ) otrzymujemy, że R = Rmax.
0
@oaieuy: Nekromancja? Wątek ma prawie okrągły rok!
0
oaieuy napisał(a)
D = R * 2 * asin(sqrt((sin((B1-B2)/2))2 + cos(B1)cos (B2)(sin((L1-L2)/2))2))
D: szukana odległość
B1, B2: szerokości geogr. punktów
L1, L2: długości geogr. punktów
promień Ziemi; Ziemia nie jest kulą, więc nie istnieje coś takiego jak jeden „jedynie słuszny” promień Ziemi — stosuje się m.in. następujące wartości:
Rn = 6366,710 km (definicja mili morskiej)
Rf = 6371,000 km (dawny standard FAI)
Rmax = 6378,137 km (maksymalny w WGS-84)
Rmin = 6356,752 km (minimalny w WGS-84)
w celu wyeliminowania wpływu niekulistości Ziemi można zastosować wzór:
R = Rmax-(Rmax-Rmin)*sin((B1+B2)/2)Macie Jajogłowi :)
Ja wziąłem te wzorki z książki Jeana Meeusa Astronomical Algorithms
Poza tym napisałem że to są wartości przybliżone