Calcular coordenadas dado um rolamento e um raio
-
22-08-2019 - |
Pergunta
Estou tendo problemas de execução da função descrita aqui aqui .
Esta é a minha aplicação Java:
private static double[] pointRadialDistance(double lat1, double lon1,
double radianBearing, double radialDistance) {
double lat = Math.asin(Math.sin(lat1)*Math.cos(radialDistance)+Math.cos(lat1)
*Math.sin(radialDistance)*Math.cos(radianBearing));
double lon;
if(Math.cos(lat) == 0) { // Endpoint a pole
lon=lon1;
}
else {
lon = ((lon1-Math.asin(Math.sin(radianBearing)*Math.sin(radialDistance)/Math.cos(lat))
+Math.PI) % (2*Math.PI)) - Math.PI;
}
return (new double[]{lat, lon});
}
converter o rolamento grau em radianos e converter a distância (km) em uma distância radianos antes de chamar a função -. Modo que não é o problema
No entanto, quando eu coordenadas de entrada, tais como: lat = 49,25705; lon = -123,140259; com um rolamento de 225 (Sudoeste) e uma distância de um km
Recebo esta retornado: lat: -1,0085434360125864 lon: -3,7595299668539504
Sua obviamente não é correto, alguém pode ver o que estou fazendo de errado?
Graças
Solução
Parece que estes são os problemas em seu código:
- Você precisa converter
lat1
elon1
para radianos antes de chamar a sua função. - Você pode estar escalando
radialDistance
incorretamente. - Testing um número de ponto flutuante de igualdade é perigoso. Dois números que são iguais depois poder aritmética exata não ser exatamente igual depois de aritmética de ponto flutuante. Assim
abs(x-y) < threshold
é mais seguro do quex == y
para testar dois números de ponto flutuantex
ey
pela igualdade. - Eu acho que você deseja converter
lat
elon
de radianos para graus.
Aqui está minha implementação do seu código em Python:
#!/usr/bin/env python
from math import asin,cos,pi,sin
rEarth = 6371.01 # Earth's average radius in km
epsilon = 0.000001 # threshold for floating-point equality
def deg2rad(angle):
return angle*pi/180
def rad2deg(angle):
return angle*180/pi
def pointRadialDistance(lat1, lon1, bearing, distance):
"""
Return final coordinates (lat2,lon2) [in degrees] given initial coordinates
(lat1,lon1) [in degrees] and a bearing [in degrees] and distance [in km]
"""
rlat1 = deg2rad(lat1)
rlon1 = deg2rad(lon1)
rbearing = deg2rad(bearing)
rdistance = distance / rEarth # normalize linear distance to radian angle
rlat = asin( sin(rlat1) * cos(rdistance) + cos(rlat1) * sin(rdistance) * cos(rbearing) )
if cos(rlat) == 0 or abs(cos(rlat)) < epsilon: # Endpoint a pole
rlon=rlon1
else:
rlon = ( (rlon1 - asin( sin(rbearing)* sin(rdistance) / cos(rlat) ) + pi ) % (2*pi) ) - pi
lat = rad2deg(rlat)
lon = rad2deg(rlon)
return (lat, lon)
def main():
print "lat1 \t lon1 \t\t bear \t dist \t\t lat2 \t\t lon2"
testcases = []
testcases.append((0,0,0,1))
testcases.append((0,0,90,1))
testcases.append((0,0,0,100))
testcases.append((0,0,90,100))
testcases.append((49.25705,-123.140259,225,1))
testcases.append((49.25705,-123.140259,225,100))
testcases.append((49.25705,-123.140259,225,1000))
for lat1, lon1, bear, dist in testcases:
(lat,lon) = pointRadialDistance(lat1,lon1,bear,dist)
print "%6.2f \t %6.2f \t %4.1f \t %6.1f \t %6.2f \t %6.2f" % (lat1,lon1,bear,dist,lat,lon)
if __name__ == "__main__":
main()
Aqui está a saída:
lat1 lon1 bear dist lat2 lon2
0.00 0.00 0.0 1.0 0.01 0.00
0.00 0.00 90.0 1.0 0.00 -0.01
0.00 0.00 0.0 100.0 0.90 0.00
0.00 0.00 90.0 100.0 0.00 -0.90
49.26 -123.14 225.0 1.0 49.25 -123.13
49.26 -123.14 225.0 100.0 48.62 -122.18
49.26 -123.14 225.0 1000.0 42.55 -114.51
Outras dicas
Eu acho que há um problema no algoritmo fornecido na mensagem 5.
Ele funciona, mas apenas para a latitude, para a longitude não é um problema por causa do sinal.
Os dados falam por si mesmos:
49,26 -123,14 225,0 1,0 49.25 -123,13
Se você começar a partir de -123,14 ° e vá para o oeste você deve ter algo FAR no Ocidente. Aqui vamos nós de volta no EAST (-123,13)!
A fórmula deve inclui algum lugar:
degreeBearing = ((360-degreeBearing)% 360)
Antes de conversão radiano.
Quando eu implementei isso, meus latitudes resultantes foram corretas, mas as longitudes estavam errados. Por exemplo ponto de partida: 36.9460678N 9.434807E, Rolamento 45,03334, Distância 15,0083313 km O resultado foi 37.0412865N 9.315302E Isso é mais a oeste do que o meu ponto de partida, em vez de mais a leste. Na verdade, é como se o rolamento foi 315.03334 graus.
Mais web busca me levou a: http: //www.movable- type.co.uk/scripts/latlong.html O código de longitude é mostrado abaixo (em C # com tudo em radianos)
if ((Math.Cos(rLat2) == 0) || (Math.Abs(Math.Cos(rLat2)) < EPSILON))
{
rLon2 = rLon1;
}
else
{
rLon2 = rLon1 + Math.Atan2(Math.Sin(rBearing) * Math.Sin(rDistance) * Math.Cos(rLat1), Math.Cos(rDistance) - Math.Sin(rLat1) * Math.Sin(rLat2));
}
Isso parece funcionar bem para mim. Espero que seja útil.
Obrigado por seu código python Eu tentei configurá-lo no meu caso de uso onde eu estou tentando encontrar o lon lat de um ponto entre dois outros a uma distância definida a partir do primeiro ponto por isso é bastante similare ao seu código appart que meu rolamento é calculado dinamicamente
startpoint (lat1) lon1 / lat1 = 55,625541, -21,142463
ponto final (Lat2) lon2 / Lat2 = 55,625792, -22,142248
meu resultado deve ser um ponto entre estes dois em lon3 / lat3 unfortunetly I get lon3 / lat3 = 0.0267695450609,0.0223553243666
Eu pensei que isso poderia ser uma diferença na lon lat mas nenhuma quando eu adicionar ou sub-lo não é bom
qualquer conselho seria realmente grande Graças
aqui está a minha aplicação
= 0,001 distância epsilon = 0.000001
cálculo rolamento dinamicamente
y = math.sin(distance) * math.cos(lat2);
x = math.cos(lat1)*math.sin(lat2) - math.sin(lat1)*math.cos(lat2)*math.cos(distance);
bearing = math.atan2(y, x)
cálculo lat3 lon3 dinamicamente
rlat1 = (lat1 * 180) / math.pi
rlon1 = (lon1 * 180) / math.pi
rbearing = (bearing * 180) / math.pi
rdistance = distance / R # normalize linear distance to radian angle
rlat = math.asin( math.sin(rlat1) * math.cos(rdistance) + math.cos(rlat1) * math.sin(rdistance) * math.cos(rbearing) )
if math.cos(rlat) == 0 or abs(math.cos(rlat)) < epsilon: # Endpoint a pole
rlon=rlon1
else:
rlon = ( (rlon1 + math.asin( math.sin(rbearing)* math.sin(rdistance) / math.cos(rlat) ) + math.pi ) % (2*math.pi) ) - math.pi
lat3 = (rlat * math.pi)/ 180
lon3 = (rlon * math.pi)/ 180
Tudo está funcionando como previsto, mas o problema é que um simples cálculo assume a Terra era uma esfera, quando na realidade se aproxima um elipsóide.
A rede de arrasto rápida de seu motor de busca preferido para 'Vincenty Formula' venha a revelar-se útil.