Obliczanie szerokości / długości geograficznej X mil od punktu?


46

Chcę znaleźć punkt szerokości i długości geograficznej na podstawie namiaru, odległości oraz początkowej szerokości i długości geograficznej.

Wydaje się, że jest to przeciwieństwo tego pytania ( odległość między punktami szerokości / długości ).

Zajrzałem już do formuły haverine i sądzę, że jej przybliżenie świata jest prawdopodobnie wystarczająco bliskie.

Zakładam, że muszę rozwiązać formułę haverine dla mojego nieznanego lat / long, czy to prawda? Czy są jakieś dobre strony internetowe, które mówią na ten temat? Wygląda na to, że byłoby to powszechne, ale moje google pojawiły się tylko pytania podobne do powyższego.

To, czego naprawdę szukam, to po prostu wzór na to. Chciałbym dać mu początkową długość / długość, namiar i odległość (mile lub kilometry) i chciałbym z niej wydostać parę długość / długość, która reprezentuje miejsce, w którym by się znalazł, gdyby podróżowali ta trasa.


Szukasz narzędzia, które to robi (np. Pe.dll Esri) lub formuły?
Kirk Kuykendall

Przepraszam, że nie byłem konkretny ... Szukam formuły. Zaktualizuję moje pytanie, aby było bardziej szczegółowe.
Jason Whitehorn

Kilka opracowanych próbek matematycznych znajduje się tutaj <a href=" movable-type.co.uk/scripts/latlong.html"> Oblicz odległość, namiar i inne wartości między punktami szerokości / długości geograficznej </a>, co obejmuje rozwiązanie „Miejsce docelowe punkt podana odległość i namiar od punktu początkowego ".
Stephen Quan


oto strona linkujące do lat / long obliczenia [długość / szerokość geograficzna obliczenia] ( movable-type.co.uk/scripts/latlong.html ) również tę stronę Lat / długie obliczenia istnieje kod + Kalkulator
Abo gaseem

Odpowiedzi:


21

Byłbym ciekawy, jak wyniki tej formuły porównać z peri Esri .

( cytowanie ).

Punkt {lat, lon} to odległość d na promieniu tc od punktu 1, jeżeli:

 lat=asin(sin(lat1)*cos(d)+cos(lat1)*sin(d)*cos(tc))
 IF (cos(lat)=0)
    lon=lon1      // endpoint a pole
 ELSE
    lon=mod(lon1-asin(sin(tc)*sin(d)/cos(lat))+pi,2*pi)-pi
 ENDIF

Algorytm ten jest ograniczony do odległości takich, że dlon <pi / 2, tj. Tych, które rozciągają się wokół mniej niż jednej czwartej obwodu Ziemi na długości geograficznej. Całkowicie ogólny, ale bardziej skomplikowany algorytm jest konieczny, jeśli dozwolone są większe odległości:

 lat =asin(sin(lat1)*cos(d)+cos(lat1)*sin(d)*cos(tc))
 dlon=atan2(sin(tc)*sin(d)*cos(lat1),cos(d)-sin(lat1)*sin(lat))
 lon=mod( lon1-dlon +pi,2*pi )-pi

Oto strona HTML do testowania .


Dziękuję za szybka odpowiedz. Pozwól mi przetrawić niektóre z tych informacji, a wrócę z tobą. Na pozór wygląda to jednak idealnie.
Jason Whitehorn

1
Próbowałem bezpośredniego przypadku przy użyciu pe.dll (właściwie libpe.so na solaris) po pobraniu odległości i azymutu przekazującego ze strony HTML dla 32N, 117W do 40N, 82W. Używając 32N, 117W, d = 3255.056515890041, azi = 64.24498012065699, otrzymałem 40N, 82W (faktycznie 82,00000000064).
mkennedy

3
Niesamowite! Dziękuję bardzo za link do artykułu Aviation Formulary autorstwa Eda Williamsa, nie widziałem tego wcześniej, ale jak dotąd okazał się świetną lekturą. Uwaga dla każdego, kto patrzy na to w przyszłości, dane wejściowe i wyjściowe tej formuły są WSZYSTKIE w radianach, nawet na odległość.
Jason Whitehorn

1
Jaka jest jednostka odległości w tym wzorze?
Karthik Murugan,

1
@KarthikMurugan Ed wprowadza, że jednostki odległości są w radianach wzdłuż wielkiego koła.
Kirk Kuykendall

18

Jeśli jesteś w samolocie, to punkt, który jest r metrów na łożysku ciągu stopni na wschód od północy jest przemieszczany przez R * cos (A) w kierunku północnym a R * sin (a) w kierunku wschodnim. (Te stwierdzenia mniej więcej definiują sinus i cosinus.)

Chociaż nie jesteś w samolocie - pracujesz na powierzchni zakrzywionej elipsoidy, która modeluje powierzchnię Ziemi - każda odległość mniejsza niż kilkaset kilometrów obejmuje tak niewielką część powierzchni, że dla większości praktycznych celów może być uważanym za płaski. Jedyną pozostałą komplikacją jest to, że jeden stopień długości geograficznej nie obejmuje tej samej odległości co stopień szerokości geograficznej. W sferycznym modelu Ziemi jeden stopień długości geograficznej jest tylko cos (szerokość geograficzna) tak długi, jak stopień szerokości geograficznej. (W modelu elipsoidalnym jest to nadal doskonałe przybliżenie, dobre do około 2,5 znaczących liczb).

Wreszcie jeden stopień szerokości geograficznej wynosi około 10 ^ 7/90 = 111 111 metrów. Mamy teraz wszystkie informacje potrzebne do przeliczenia liczników na stopnie:

Przemieszczenie na północ wynosi r * cos (a) / 111111 stopni;

Przemieszczenie na wschód to r * sin (a) / cos (szerokość geograficzna) / 111111 stopni.

Na przykład na szerokości -0,31399 stopni i namiocie a = 30 stopni na wschód od północy możemy obliczyć

cos(a) = cos(30 degrees) = cos(pi/6 radians) = Sqrt(3)/2 = 0.866025.
sin(a) = sin(30 degrees) = sin(pi/6 radians) = 1/2 = 0.5.
cos(latitude) = cos(-0.31399 degrees) = cos(-0.00548016 radian) = 0.999984984.
r = 100 meters.
east displacement = 100 * 0.5 / 0.999984984 / 111111 = 0.000450007 degree.
north displacement = 100 * 0.866025 / 111111 = 0.000779423 degree.

Skąd, począwszy od (-78,4437, -0,31399), nowa lokalizacja to (-78,4437 + 0,00045, -0,31399 + 0,0007794) = (-78,4432, -0,313211).

Bardziej dokładna odpowiedź w nowoczesnym systemie odniesienia ITRF00 to (-78.4433, -0.313207): jest to 0,43 metra od przybliżonej odpowiedzi, co wskazuje na błędy przybliżania o 0,43% w tym przypadku. Aby osiągnąć wyższą dokładność, musisz użyć albo elipsoidalnych wzorów odległości (które są znacznie bardziej skomplikowane), albo rzetelnego rzutu konformalnego o zerowej dywergencji (aby łożysko było prawidłowe).


2
+1 za prawidłowe zrozumienie kontekstu matematycznego (tj. Jego płaszczyzny lokalnej)
Frank Conry

4

Jeśli potrzebujesz rozwiązania JavaScript, rozważ te functionsi to skrzypce :

var gis = {
  /**
  * All coordinates expected EPSG:4326
  * @param {Array} start Expected [lon, lat]
  * @param {Array} end Expected [lon, lat]
  * @return {number} Distance - meter.
  */
  calculateDistance: function(start, end) {
    var lat1 = parseFloat(start[1]),
        lon1 = parseFloat(start[0]),
        lat2 = parseFloat(end[1]),
        lon2 = parseFloat(end[0]);

    return gis.sphericalCosinus(lat1, lon1, lat2, lon2);
  },

  /**
  * All coordinates expected EPSG:4326
  * @param {number} lat1 Start Latitude
  * @param {number} lon1 Start Longitude
  * @param {number} lat2 End Latitude
  * @param {number} lon2 End Longitude
  * @return {number} Distance - meters.
  */
  sphericalCosinus: function(lat1, lon1, lat2, lon2) {
    var radius = 6371e3; // meters
    var dLon = gis.toRad(lon2 - lon1),
        lat1 = gis.toRad(lat1),
        lat2 = gis.toRad(lat2),
        distance = Math.acos(Math.sin(lat1) * Math.sin(lat2) +
            Math.cos(lat1) * Math.cos(lat2) * Math.cos(dLon)) * radius;

    return distance;
  },

  /**
  * @param {Array} coord Expected [lon, lat] EPSG:4326
  * @param {number} bearing Bearing in degrees
  * @param {number} distance Distance in meters
  * @return {Array} Lon-lat coordinate.
  */
  createCoord: function(coord, bearing, distance) {
    /** http://www.movable-type.co.uk/scripts/latlong.html
    * φ is latitude, λ is longitude, 
    * θ is the bearing (clockwise from north), 
    * δ is the angular distance d/R; 
    * d being the distance travelled, R the earth’s radius*
    **/
    var 
      radius = 6371e3, // meters
      δ = Number(distance) / radius, // angular distance in radians
      θ = gis.toRad(Number(bearing));
      φ1 = gis.toRad(coord[1]),
      λ1 = gis.toRad(coord[0]);

    var φ2 = Math.asin(Math.sin(φ1)*Math.cos(δ) + 
      Math.cos(φ1)*Math.sin(δ)*Math.cos(θ));

    var λ2 = λ1 + Math.atan2(Math.sin(θ) * Math.sin(δ)*Math.cos(φ1),
      Math.cos(δ)-Math.sin(φ1)*Math.sin(φ2));
    // normalise to -180..+180°
    λ2 = (λ2 + 3 * Math.PI) % (2 * Math.PI) - Math.PI; 

    return [gis.toDeg(λ2), gis.toDeg(φ2)];
  },
  /**
   * All coordinates expected EPSG:4326
   * @param {Array} start Expected [lon, lat]
   * @param {Array} end Expected [lon, lat]
   * @return {number} Bearing in degrees.
   */
  getBearing: function(start, end){
    var
      startLat = gis.toRad(start[1]),
      startLong = gis.toRad(start[0]),
      endLat = gis.toRad(end[1]),
      endLong = gis.toRad(end[0]),
      dLong = endLong - startLong;

    var dPhi = Math.log(Math.tan(endLat/2.0 + Math.PI/4.0) / 
      Math.tan(startLat/2.0 + Math.PI/4.0));

    if (Math.abs(dLong) > Math.PI) {
      dLong = (dLong > 0.0) ? -(2.0 * Math.PI - dLong) : (2.0 * Math.PI + dLong);
    }

    return (gis.toDeg(Math.atan2(dLong, dPhi)) + 360.0) % 360.0;
  },
  toDeg: function(n) { return n * 180 / Math.PI; },
  toRad: function(n) { return n * Math.PI / 180; }
};

Więc jeśli chcesz obliczyć nową współrzędną, może to wyglądać tak:

var start = [15, 38.70250];
var end = [21.54967, 38.70250];
var total_distance = gis.calculateDistance(start, end); // meters
var percent = 10;
// this can be also meters
var distance = (percent / 100) * total_distance;
var bearing = gis.getBearing(start, end);
var new_coord = gis.createCoord(icon_coord, bearing, distance);

2

Mam to działa w ObjectiveC. Kluczem tutaj jest wiedza, że ​​potrzebujesz punktów łac. / Lng w radianach i musisz zastosować je z powrotem do łac. / Lng po zastosowaniu równania. Wiedz także, że odległość i tc podano w radianach.

Oto oryginalne równanie:

 lat=asin(sin(lat1)*cos(d)+cos(lat1)*sin(d)*cos(tc))
 IF (cos(lat)=0)
    lon=lon1      // endpoint a pole
 ELSE
    lon=mod(lon1-asin(sin(tc)*sin(d)/cos(lat))+pi,2*pi)-pi
 ENDIF

Tutaj jest implementowany w ObjC, gdzie radian jest radianem mierzonym w kierunku przeciwnym do ruchu wskazówek zegara od N (np. PI / 2 to W, PI to S, 2 PI / 3 to E), a odległość w kilometrach.

+ (CLLocationCoordinate2D)displaceLatLng:(CLLocationCoordinate2D)coordinate2D withRadian:(double)radian
                            withDistance:(CGFloat)distance {
  double lat1Radians = coordinate2D.latitude * (M_PI / 180);
  double lon1Radians = coordinate2D.longitude * (M_PI / 180);
  double distanceRadians = distance / 6371;
  double lat = asin(sin(lat1Radians)*cos(distanceRadians)+cos(lat1Radians)*sin(distanceRadians)
      *cos(radian));
  double lon;
  if (cos(lat) == 0) {
    lon = lon1Radians;
  } else {
    lon = fmodf((float) (lon1Radians - asin(sin(radian) * sin(distanceRadians) / cos(lat)) + M_PI),
        (float) (2 * M_PI)) - M_PI;
  }
  double newLat = lat * (180 / M_PI);
  double newLon = lon * (180 / M_PI);
  return CLLocationCoordinate2DMake(newLat, newLon);
}

Szukam rozwiązania, w którym chcę uzyskać 4 lat długości od miejsca, w którym stoję do 50 mil na północ, 50 mil na zachód, 50 mil na wschód i tak dalej ... W powyższej odpowiedzi możesz wyjaśnić, w jaki sposób mogę osiągnąć to?
Rahul Vyas

1

Jeśli jesteś zainteresowany lepszą dokładnością, jest Vincenty . (Link do formularza „bezpośredniego”, którego dokładnie szukasz).

Istnieje wiele istniejących implementacji, jeśli szukasz kodu.

A także pytanie: nie zakładasz, że podróżny utrzymuje to samo nastawienie przez całą podróż, prawda? Jeśli tak, to te metody nie odpowiadają na właściwe pytanie. (Lepiej byłoby rzutować na mercator, rysować linię prostą, a następnie cofać rzutowanie).


Bardzo dobre pytanie, pomimo sformułowania w moim pytaniu, które mogło wskazywać, że obliczam miejsce docelowe dla podróżnego, nie jestem. Dobra racja jednak. Było to głównie po to, abym mógł obliczyć obszar ograniczający (na małe zamówienie, powiedzmy 50 mil) ... Mam nadzieję, że to ma sens.
Jason Whitehorn,

gis.stackexchange.com/questions/3264/… miał to samo pytanie (konstruowanie obszarów granicznych z punktu i odległości)
Dan S.

0

Oto rozwiązanie Python:

    def displace(self, theta, distance):
    """
    Displace a LatLng theta degrees counterclockwise and some
    meters in that direction.
    Notes:
        http://www.movable-type.co.uk/scripts/latlong.html
        0 DEGREES IS THE VERTICAL Y AXIS! IMPORTANT!
    Args:
        theta:    A number in degrees.
        distance: A number in meters.
    Returns:
        A new LatLng.
    """
    theta = np.float32(theta)

    delta = np.divide(np.float32(distance), np.float32(E_RADIUS))

    def to_radians(theta):
        return np.divide(np.dot(theta, np.pi), np.float32(180.0))

    def to_degrees(theta):
        return np.divide(np.dot(theta, np.float32(180.0)), np.pi)

    theta = to_radians(theta)
    lat1 = to_radians(self.lat)
    lng1 = to_radians(self.lng)

    lat2 = np.arcsin( np.sin(lat1) * np.cos(delta) +
                      np.cos(lat1) * np.sin(delta) * np.cos(theta) )

    lng2 = lng1 + np.arctan2( np.sin(theta) * np.sin(delta) * np.cos(lat1),
                              np.cos(delta) - np.sin(lat1) * np.sin(lat2))

    lng2 = (lng2 + 3 * np.pi) % (2 * np.pi) - np.pi

    return LatLng(to_degrees(lat2), to_degrees(lng2))

-2

Korzystam z podejścia opisanego poniżej przy określaniu następnej współrzędnej, biorąc pod uwagę namiar i odległość od poprzedniej współrzędnej. Mam problem z dokładnością z innym podejściem, które czytam z Internetu.

Używam tego do określania obszaru ziemi, który jest wielokątem, i wykreślam go w Google Earth. Tytuł lądu ma zapisane namiary i odległości: „NorthOrSouth x stopnie i minuty EastOrWest, z metrów do punktu n;”.

Rozpoczynając od podanych współrzędnych punktu odniesienia, najpierw obliczam odległość na jeden stopień szerokości geograficznej i jeden stopień długości geograficznej, stosując wzór haverine, ponieważ zmienia się on w zależności od lokalizacji. Następnie określam następną współrzędną ze wzoru sinus i cosinus trygonometrii.

Poniżej znajduje się javascript:

var mapCenter = new google.maps.LatLng(referencePointLatitude, referencePointLongitude); //the ref point lat and lon must be given, usually a land mark (BLLM)
var latDiv = latDiv(mapCenter); //distance per one degree latitude in this location
var lngDiv = lngDiv(mapCenter); //distance per one degree longitude in this location
var LatLng2 = NextCoord(PrevCoord,NorthOrSouth,x,y,EastOrWest,z); //next coordinate given the bearing and distance from previous coordinate
var Lat2 = LatLng2.lat(); //next coord latitude in degrees
var Lng2 = LatLng2.lng(); //next coord longitude in degrees
var polygon=[p1,p2,...,pn-1,pn,p1]; //p1,p2,etc. are the coordinates of points of the polygon, i.e. the land Title. Be sure to close the polygon to the point of beginning p1
var area = Area(polygon); //area of the polygon in sq.m.
function NextCoord(PrevCoord,NorthOrSouth,x,y,EastOrWest,z) {
  var angle = ( x + ( y / 60 ) ) * Math.PI / 180;
  var a = 1;
  var b = 1;
  if (NorthOrSouth == 'South') { a = -1; }
  if (EastOrWest == 'West') { b = -1; }
  var nextLat = PrevCoord.lat() +  ( a * z * Math.cos(angle) / latDiv );
  var nextLng = PrevCoord.lng() +  ( b * z * Math.sin(angle) / lngDiv );
  var nextCoord = new google.maps.LatLng(nextLat, nextLng);
  return nextCoord;
}
function latDiv(mapCenter) {
  var p1 = new google.maps.LatLng(mapCenter.lat()-0.5, mapCenter.lng());
  var p2 = new google.maps.LatLng(mapCenter.lat()+0.5, mapCenter.lng());
  return dist(p1,p2);
}
function lngDiv(mapCenter) {
  var p3 = new google.maps.LatLng(mapCenter.lat(), mapCenter.lng()-0.5);
  var p4 = new google.maps.LatLng(mapCenter.lat(), mapCenter.lng()+0.5);
  return dist(p3,p4);
}
function dist(pt1, pt2) {
    var dLat  = ( pt2.lat() - pt1.lat() ) * Math.PI / 180;
    var dLng = ( pt2.lng() - pt1.lng() ) * Math.PI / 180;
    var a = Math.sin(dLat/2) * Math.sin(dLat/2) +                 
            Math.cos(rad(pt1.lat())) * Math.cos(rad(pt2.lat())) *
            Math.sin(dLng/2) * Math.sin(dLng/2);
    var R = 6372800; //earth's radius
    var distance = R * 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1-a));
    return distance;
}
function Area(polygon) {
  var xPts=[];
  for (i=1; i&lt;polygon.length; i++) {
    xPts[i-1] = ( polygon[i].lat() - polygon[0].lat() ) * latDiv;
  }
  var yPts=[];
  for (i=1; i&lt;polygon.length; i++) {
    yPts[i-1] = ( polygon[i].lng() - polygon[0].lng() ) * lngDiv;
  }
  var area = 0;
  j = polygon.length-2;
  for (i=0; i&lt;polygon.length-1; i++) {
    area = area +  ( xPts[j] + xPts[i] ) * ( yPts[j] - yPts[i] );
    j = i;
  }
  area = Math.abs(area/2);
  return area;
}

1
Wzór kartezjański dla obszaru wielokąta, który próbujesz tutaj zastosować, nie ma zastosowania do odległości i kątów obliczonych na zakrzywionej powierzchni (np. Sferoidy). Ta formuła powoduje dodatkowy błąd, wykorzystując szerokość i długość geograficzną, jakby były współrzędnymi kartezjańskimi. Jedynymi okolicznościami, w których można rozważyć jego zastosowanie, byłyby dokładnie te (w przypadku bardzo małych wielokątów), w których formuła haverine i tak jest niepotrzebna. Ogólnie wydaje się, że ten kod działa o wiele za ciężko, aby nie zyskać.
whuber
Korzystając z naszej strony potwierdzasz, że przeczytałeś(-aś) i rozumiesz nasze zasady używania plików cookie i zasady ochrony prywatności.
Licensed under cc by-sa 3.0 with attribution required.