Получить широту / долготу с учетом текущей точки, расстояния и пеленга
Учитывая существующую точку в широте / долготе, расстояние в (в км) и азимут (в градусах, переведенных в радианы), я хотел бы рассчитать новую широту / длину. Этот сайт появляется снова и снова, но я просто не могу заставить формулу работать на меня.
Формулы, взятые по вышеуказанной ссылке:
lat2 = asin(sin(lat1)*cos(d/R) + cos(lat1)*sin(d/R)*cos(θ))
lon2 = lon1 + atan2(sin(θ)*sin(d/R)*cos(lat1), cos(d/R)−sin(lat1)*sin(lat2))
Приведенная выше формула для MSExcel, где
asin = arc sin()
d = distance (in any unit)
R = Radius of the earth (in the same unit as above)
and hence d/r = is the angular distance (in radians)
atan2(a,b) = arc tan(b/a)
θ is the bearing (in radians, clockwise from north);
Вот код, который я получил на Python.
import math
R = 6378.1 #Radius of the Earth
brng = 1.57 #Bearing is 90 degrees converted to radians.
d = 15 #Distance in km
#lat2 52.20444 - the lat result I'm hoping for
#lon2 0.36056 - the long result I'm hoping for.
lat1 = 52.20472 * (math.pi * 180) #Current lat point converted to radians
lon1 = 0.14056 * (math.pi * 180) #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))
print(lat2)
print(lon2)
я получил
lat2 = 0.472492248844
lon2 = 79.4821662373
15 ответов
Необходимо преобразовать ответы из радианов обратно в градусы. Рабочий код ниже:
import math
R = 6378.1 #Radius of the Earth
brng = 1.57 #Bearing is 90 degrees converted to radians.
d = 15 #Distance in km
#lat2 52.20444 - the lat result I'm hoping for
#lon2 0.36056 - the long result I'm hoping for.
lat1 = math.radians(52.20472) #Current lat point converted to radians
lon1 = math.radians(0.14056) #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)
print(lat2)
print(lon2)
Библиотека Geopy поддерживает это:
import geopy
from geopy.distance import VincentyDistance
# given: lat1, lon1, b = bearing in degrees, d = distance in kilometers
origin = geopy.Point(lat1, lon1)
destination = VincentyDistance(kilometers=d).destination(origin, b)
lat2, lon2 = destination.latitude, destination.longitude
Найдено через /questions/18427015/vyichislenie-koordinatyi-gps-s-uchetom-tochki-azimuta-i-rasstoyaniya/18427023#18427023
Этот вопрос известен как прямая проблема в изучении геодезии.
Это действительно очень популярный вопрос, который постоянно вызывает путаницу. Причина в том, что большинство людей ищут простой и прямой ответ. Но их нет, потому что большинство людей, задающих этот вопрос, не предоставляют достаточно информации просто потому, что не знают, что:
- Земля не идеальная сфера, так как она сплюснута / сжата полюсами
- Из-за (1) земля не имеет постоянного радиуса,
R, Смотрите здесь. - Земля не совсем гладкая (колебания высоты) и т. Д.
- Из-за движения тектонических плит положение географической точки может изменяться на несколько миллиметров (как минимум) каждый год.
Поэтому в различных геометрических моделях используется много разных допущений, которые применяются по-разному, в зависимости от необходимой точности. Поэтому, чтобы ответить на вопрос, нужно подумать, с какой точностью вы хотели бы получить свой результат.
Некоторые примеры:
- Я просто ищу приблизительное местоположение на ближайшие несколько километров для небольших (< 100 км) расстояний в
latitudesмежду0-70 degN | S.(Земля - плоская модель.) - Я хочу получить ответ, который был бы хорош в любой точке земного шара, но с точностью до нескольких метров.
- Я хочу супер точное позиционирование, которое действует вплоть до атомных масштабов
nanometers[Нм]. - Я хочу получить ответы, которые очень быстро и легко рассчитываются и не требуют больших вычислительных ресурсов.
Таким образом, вы можете иметь много вариантов выбора алгоритма. Кроме того, каждый язык программирования имеет свою собственную реализацию или "пакет", умноженный на количество моделей и конкретные потребности разработчиков моделей. Для всех практических целей здесь стоит игнорировать любой другой язык, кроме javascript, поскольку по своей природе он очень похож на псевдокод. Таким образом, он может быть легко преобразован в любой другой язык с минимальными изменениями.
Тогда основными моделями являются:
Euclidian/Flat earth model: хорошо для очень коротких расстояний до ~10 кмSpherical model: подходит для больших продольных расстояний, но с небольшой разницей в ширине. Популярная модель:- Haversine: точность метра в [км] масштабах, очень простой код.
Ellipsoidal models: Наиболее точный на любом широте и долготе и расстоянии, но все еще числовое приближение, которое зависит от того, какая точность вам нужна. Некоторые популярные модели:- Ламберт: ~10 метров с точностью более 1000 км.
- Пол Д. Томас: приближение Андойера-Ламберта
- Vincenty: миллиметровая точность и вычислительная эффективность
- Керни: нанометровая точность
Рекомендации:
- https://en.wikipedia.org/wiki/Reference_ellipsoid
- https://en.wikipedia.org/wiki/Haversine_formula
- https://en.wikipedia.org/wiki/Earth_ellipsoid
- https://en.wikipedia.org/wiki/Geodesics_on_an_ellipsoid
- https://en.wikipedia.org/wiki/Vincenty%27s_formulae
- https://geographiclib.sourceforge.io/scripts/geod-calc.html
Может быть немного поздно для ответа, но после тестирования других ответов, кажется, что они не работают правильно. Вот код PHP, который мы используем для нашей системы. Работаем по всем направлениям.
Код PHP:
lat1 = широта начальной точки в градусах
long1 = долгота начальной точки в градусах
d = расстояние в км
угол = азимут в градусах
function get_gps_distance($lat1,$long1,$d,$angle)
{
# Earth Radious in KM
$R = 6378.14;
# Degree to Radian
$latitude1 = $lat1 * (M_PI/180);
$longitude1 = $long1 * (M_PI/180);
$brng = $angle * (M_PI/180);
$latitude2 = asin(sin($latitude1)*cos($d/$R) + cos($latitude1)*sin($d/$R)*cos($brng));
$longitude2 = $longitude1 + atan2(sin($brng)*sin($d/$R)*cos($latitude1),cos($d/$R)-sin($latitude1)*sin($latitude2));
# back to degrees
$latitude2 = $latitude2 * (180/M_PI);
$longitude2 = $longitude2 * (180/M_PI);
# 6 decimal for Leaflet and other system compatibility
$lat2 = round ($latitude2,6);
$long2 = round ($longitude2,6);
// Push in array and get back
$tab[0] = $lat2;
$tab[1] = $long2;
return $tab;
}
Я перенес ответ от Брэда на ванильный ответ JS, без зависимости от карт Bing
https://jsfiddle.net/kodisha/8a3hcjtd/
// ----------------------------------------
// Calculate new Lat/Lng from original points
// on a distance and bearing (angle)
// ----------------------------------------
let llFromDistance = function(latitude, longitude, distance, bearing) {
// taken from: https://stackru.com/a/46410871/13549
// distance in KM, bearing in degrees
const R = 6378.1; // Radius of the Earth
const brng = bearing * Math.PI / 180; // Convert bearing to radian
let lat = latitude * Math.PI / 180; // Current coords to radians
let lon = longitude * Math.PI / 180;
// Do the math magic
lat = Math.asin(Math.sin(lat) * Math.cos(distance / R) + Math.cos(lat) * Math.sin(distance / R) * Math.cos(brng));
lon += Math.atan2(Math.sin(brng) * Math.sin(distance / R) * Math.cos(lat), Math.cos(distance / R) - Math.sin(lat) * Math.sin(lat));
// Coords back to degrees and return
return [(lat * 180 / Math.PI), (lon * 180 / Math.PI)];
}
let pointsOnMapCircle = function(latitude, longitude, distance, numPoints) {
const points = [];
for (let i = 0; i <= numPoints - 1; i++) {
const bearing = Math.round((360 / numPoints) * i);
console.log(bearing, i);
const newPoints = llFromDistance(latitude, longitude, distance, bearing);
points.push(newPoints);
}
return points;
}
const points = pointsOnMapCircle(41.890242042122836, 12.492358982563019, 0.2, 8);
let geoJSON = {
"type": "FeatureCollection",
"features": []
};
points.forEach((p) => {
geoJSON.features.push({
"type": "Feature",
"properties": {},
"geometry": {
"type": "Point",
"coordinates": [
p[1],
p[0]
]
}
});
});
document.getElementById('res').innerHTML = JSON.stringify(geoJSON, true, 2);
Кроме того, я добавил экспорт geoJSON, чтобы вы могли просто вставить полученный geoJSON в http://geojson.io/ и сразу увидеть результаты.
Результат: 
Быстрый способ с помощью геопии
from geopy import distance
#distance.distance(unit=15).destination((lat,lon),bering)
#Exemples
distance.distance(nautical=15).destination((-24,-42),90)
distance.distance(miles=15).destination((-24,-42),90)
distance.distance(kilometers=15).destination((-24,-42),90)
lon1 и lat1 в градусах
brng = азимут в радианах
d = расстояние в км
R = радиус Земли в км
lat2 = math.degrees((d/R) * math.cos(brng)) + lat1
long2 = math.degrees((d/(R*math.sin(math.radians(lat2)))) * math.sin(brng)) + long1
Я реализовал ваш алгоритм и мой в PHP и протестировал его. Эта версия работала примерно в 50% случаев. Полученные результаты были идентичны, поэтому, похоже, математически эквивалентны.
Я не тестировал код Python выше, поэтому могут быть синтаксические ошибки.
Я портировал Python на Javascript. Это возвращает Bing Maps Location объект, вы можете изменить на что угодно.
getLocationXDistanceFromLocation: function(latitude, longitude, distance, bearing) {
// distance in KM, bearing in degrees
var R = 6378.1, // Radius of the Earth
brng = Math.radians(bearing) // Convert bearing to radian
lat = Math.radians(latitude), // Current coords to radians
lon = Math.radians(longitude);
// Do the math magic
lat = Math.asin(Math.sin(lat) * Math.cos(distance / R) + Math.cos(lat) * Math.sin(distance / R) * Math.cos(brng));
lon += Math.atan2(Math.sin(brng) * Math.sin(distance / R) * Math.cos(lat), Math.cos(distance/R)-Math.sin(lat)*Math.sin(lat));
// Coords back to degrees and return
return new Microsoft.Maps.Location(Math.degrees(lat), Math.degrees(lon));
},
Благодаря @kodisha, вот версия Swift, но с улучшенным и более точным расчетом радиуса Земли:
extension CLLocationCoordinate2D {
func earthRadius() -> CLLocationDistance {
let earthRadiusInMetersAtSeaLevel = 6378137.0
let earthRadiusInMetersAtPole = 6356752.314
let r1 = earthRadiusInMetersAtSeaLevel
let r2 = earthRadiusInMetersAtPole
let beta = latitude
let earthRadiuseAtGivenLatitude = (
( pow(pow(r1, 2) * cos(beta), 2) + pow(pow(r2, 2) * sin(beta), 2) ) /
( pow(r1 * cos(beta), 2) + pow(r2 * sin(beta), 2) )
)
.squareRoot()
return earthRadiuseAtGivenLatitude
}
func locationByAdding(
distance: CLLocationDistance,
bearing: CLLocationDegrees
) -> CLLocationCoordinate2D {
let latitude = self.latitude
let longitude = self.longitude
let earthRadiusInMeters = self.earthRadius()
let brng = bearing.degreesToRadians
var lat = latitude.degreesToRadians
var lon = longitude.degreesToRadians
lat = asin(
sin(lat) * cos(distance / earthRadiusInMeters) +
cos(lat) * sin(distance / earthRadiusInMeters) * cos(brng)
)
lon += atan2(
sin(brng) * sin(distance / earthRadiusInMeters) * cos(lat),
cos(distance / earthRadiusInMeters) - sin(lat) * sin(lat)
)
let newCoordinate = CLLocationCoordinate2D(
latitude: lat.radiansToDegrees,
longitude: lon.radiansToDegrees
)
return newCoordinate
}
}
extension FloatingPoint {
var degreesToRadians: Self { self * .pi / 180 }
var radiansToDegrees: Self { self * 180 / .pi }
}
Также поздно, но для тех, кто может найти это, вы получите более точные результаты, используя библиотеку geographiclib. Ознакомьтесь с описаниями геодезических задач и примерами JavaScript, чтобы легко узнать, как их использовать, чтобы ответить на предметный вопрос, а также для многих других. Реализации на разных языках, включая Python. Намного лучше, чем кодировать свои собственные, если вы заботитесь о точности; лучше, чем VincentyDistance в более ранней рекомендации "использовать библиотеку". Как сказано в документации: "Акцент делается на возвращении точных результатов с ошибками, близкими к округлению (около 5–15 нм)".
Просто поменяйте местами значения в функции atan2(y,x). Не atan2(x,y)!
Я портировал ответ от @David M на java, если кто-то хотел этого... Я получаю немного другой результат: 52.20462299620793, 0.360433887489931
double R = 6378.1; //Radius of the Earth
double brng = 1.57; //Bearing is 90 degrees converted to radians.
double d = 15; //Distance in km
double lat2 = 52.20444; // - the lat result I'm hoping for
double lon2 = 0.36056; // - the long result I'm hoping for.
double lat1 = Math.toRadians(52.20472); //Current lat point converted to radians
double lon1 = Math.toRadians(0.14056); //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.toDegrees(lat2);
lon2 = Math.toDegrees(lon2);
System.out.println(lat2 + ", " + lon2);
Вот версия PHP, основанная на формуляре авиации Эда Уильямса. Модуль обрабатывается немного по-другому в PHP. Это работает для меня.
function get_new_waypoint ( $lat, $lon, $radial, $magvar, $range )
{
// $range in nm.
// $radial is heading to or bearing from
// $magvar for local area.
$range = $range * pi() /(180*60);
$radial = $radial - $magvar ;
if ( $radial < 1 )
{
$radial = 360 + $radial - $magvar;
}
$radial = deg2rad($radial);
$tmp_lat = deg2rad($lat);
$tmp_lon = deg2rad($lon);
$new_lat = asin(sin($tmp_lat)* cos($range) + cos($tmp_lat) * sin($range) * cos($radial));
$new_lat = rad2deg($new_lat);
$new_lon = $tmp_lon - asin(sin($radial) * sin($range)/cos($new_lat))+ pi() % 2 * pi() - pi();
$new_lon = rad2deg($new_lon);
return $new_lat." ".$new_lon;
}
Для тех, кто интересуется решением Java, вот мой код: я заметил, что исходное решение требует некоторых настроек, чтобы вернуть правильное значение долготы, особенно когда точка находится на одном из полюсов. Также иногда требуется операция округления, так как результаты на 0 широте/долготе кажутся немного смещенными от 0. Для небольших расстояний в этом отношении поможет округление.
private static final double EARTH_RADIUS = 6371; // average earth radius
/**
* Returns the coordinates of the point situated at the distance specified, in
* the direction specified. Note that the value is an approximation, not an
* exact result.
*
* @param startPointLatitude
* @param startPointLongitude
* @param distanceInKm
* @param bearing: 0 means moving north, 90 moving east, 180 moving
* south, 270 moving west. Max value 360 min value 0;
* @return new point location
*/
public static LocationDTO getPointAt(double startPointLatitude, double startPointLongitude, double distanceInKm,
double bearing) {
if (Math.abs(startPointLatitude) > 90) {
throw new BadRequestException(ExceptionMessages.INVALID_LATITUDE);
} else if (Math.abs(startPointLatitude) == 90) {
startPointLatitude = 89.99999 * Math.signum(startPointLatitude); // we have to do this conversion else the formula doesnt return the correct longitude value
}
if (Math.abs(startPointLongitude) > 180) {
throw new BadRequestException(ExceptionMessages.INVALID_LONGITUDE);
}
double angularDistance = distanceInKm / EARTH_RADIUS;
bearing = deg2rad(bearing);
startPointLatitude = deg2rad(startPointLatitude);
startPointLongitude = deg2rad(startPointLongitude);
double latitude = Math.asin(Math.sin(startPointLatitude) * Math.cos(angularDistance)
+ Math.cos(startPointLatitude) * Math.sin(angularDistance) * Math.cos(bearing));
double longitude = startPointLongitude
+ Math.atan2(Math.sin(bearing) * Math.sin(angularDistance) * Math.cos(startPointLatitude),
Math.cos(angularDistance) - Math.sin(startPointLatitude) * Math.sin(latitude));
longitude = (rad2deg(longitude) + 540) % 360 - 180; // normalize longitude to be in -180 +180 interval
LocationDTO result = new LocationDTO();
result.setLatitude(roundValue(rad2deg(latitude)));
result.setLongitude(roundValue(longitude));
return result;
}
private static double roundValue(double value) {
DecimalFormat df = new DecimalFormat("#.#####");
df.setRoundingMode(RoundingMode.CEILING);
return Double.valueOf(df.format(value));
}
// This function converts decimal degrees to radians
private static double deg2rad(double deg) {
return (deg * Math.PI / 180.0);
}
// This function converts radians to decimal degrees
private static double rad2deg(double rad) {
return (rad * 180.0 / Math.PI);
}
Очень поздно на вечеринку, но вот ответ на R для всех, кому интересно. Единственное изменение, которое я сделал, это то, что я установил радиус в метры, поэтомуdтакже необходимо установить метры.
get_point_at_distance <- function(lon, lat, d, bearing, R = 6378137) {
# lat: initial latitude, in degrees
# lon: initial longitude, in degrees
# d: target distance from initial point (in m)
# bearing: (true) heading in degrees
# R: mean radius of earth (in m)
# Returns new lat/lon coordinate {d} m from initial, in degrees
## convert to radians
lat1 <- lat * (pi/180)
lon1 <- lon * (pi/180)
a <- bearing * (pi/180)
## new position
lat2 <- asin(sin(lat1) * cos(d/R) + cos(lat1) * sin(d/R) * cos(a))
lon2 <- lon1 + atan2(
sin(a) * sin(d/R) * cos(lat1),
cos(d/R) - sin(lat1) * sin(lat2)
)
## convert back to degrees
lat2 <- lat2 * (180/pi)
lon2 <- lon2 * (180/pi)
## return
return(c(lon2, lat2))
}
lat = 52.20472
lon = 0.14056
distance = 15000
bearing = 90
get_point_at_distance(lon = lon, lat = lat, d = distance, bearing = bearing)
# [1] 0.3604322 52.2045157