如何计算与另一个点相距一定距离的点的纬度?

IT技术 javascript google-maps
2021-01-16 12:55:24

为了在地图上画一个圆,我有一个中心 GLatLng (A) 和一个以米为单位的半径 (r)。

这是一个图表:

           -----------
        --/           \--
      -/                 \-
     /                     \
    /                       \
   /                   r     \
   |            *-------------*
   \             A           / B
    \                       /
     \                     /
      -\                 /-
        --\           /--
           -----------

如何计算位置 B 的 GLatLng?假设r平行于赤道。

使用 GLatLng.distanceFrom() 方法在给定 A 和 B 时获取半径是微不足道的 - 但反过来则不然。似乎我需要做一些更重的数学运算。

6个回答

我们需要一个方法,当给定方位角和从源点行进的距离时返回目标点。幸运的是,Chris Veness 在计算纬度/经度点之间的距离、方位等方面有一个非常好的 JavaScript 实现

以下内容已经过调整以适用于google.maps.LatLng该类:

Number.prototype.toRad = function() {
   return this * Math.PI / 180;
}

Number.prototype.toDeg = function() {
   return this * 180 / Math.PI;
}

google.maps.LatLng.prototype.destinationPoint = function(brng, dist) {
   dist = dist / 6371;  
   brng = brng.toRad();  

   var lat1 = this.lat().toRad(), lon1 = this.lng().toRad();

   var lat2 = Math.asin(Math.sin(lat1) * Math.cos(dist) + 
                        Math.cos(lat1) * Math.sin(dist) * Math.cos(brng));

   var lon2 = lon1 + Math.atan2(Math.sin(brng) * Math.sin(dist) *
                                Math.cos(lat1), 
                                Math.cos(dist) - Math.sin(lat1) *
                                Math.sin(lat2));

   if (isNaN(lat2) || isNaN(lon2)) return null;

   return new google.maps.LatLng(lat2.toDeg(), lon2.toDeg());
}

您只需按如下方式使用它:

var pointA = new google.maps.LatLng(25.48, -71.26); 
var radiusInKm = 10;

var pointB = pointA.destinationPoint(90, radiusInKm);

这是使用Google Maps API v3的完整示例

<!DOCTYPE html>
<html> 
<head> 
   <meta http-equiv="content-type" content="text/html; charset=UTF-8"/> 
   <title>Google Maps Geometry</title> 
   <script src="http://maps.google.com/maps/api/js?sensor=false" 
           type="text/javascript"></script> 
</head> 
<body> 
   <div id="map" style="width: 400px; height: 300px"></div> 

   <script type="text/javascript"> 
      Number.prototype.toRad = function() {
         return this * Math.PI / 180;
      }

      Number.prototype.toDeg = function() {
         return this * 180 / Math.PI;
      }

      google.maps.LatLng.prototype.destinationPoint = function(brng, dist) {
         dist = dist / 6371;  
         brng = brng.toRad();  

         var lat1 = this.lat().toRad(), lon1 = this.lng().toRad();

         var lat2 = Math.asin(Math.sin(lat1) * Math.cos(dist) + 
                              Math.cos(lat1) * Math.sin(dist) * Math.cos(brng));

         var lon2 = lon1 + Math.atan2(Math.sin(brng) * Math.sin(dist) *
                                      Math.cos(lat1), 
                                      Math.cos(dist) - Math.sin(lat1) *
                                      Math.sin(lat2));

         if (isNaN(lat2) || isNaN(lon2)) return null;

         return new google.maps.LatLng(lat2.toDeg(), lon2.toDeg());
      }

      var pointA = new google.maps.LatLng(40.70, -74.00);   // Circle center
      var radius = 10;                                      // 10km

      var mapOpt = { 
         mapTypeId: google.maps.MapTypeId.TERRAIN,
         center: pointA,
         zoom: 10
      };

      var map = new google.maps.Map(document.getElementById("map"), mapOpt);

      // Draw the circle
      new google.maps.Circle({
         center: pointA,
         radius: radius * 1000,       // Convert to meters
         fillColor: '#FF0000',
         fillOpacity: 0.2,
         map: map
      });

      // Show marker at circle center
      new google.maps.Marker({
         position: pointA,
         map: map
      });

      // Show marker at destination point
      new google.maps.Marker({
         position: pointA.destinationPoint(90, radius),
         map: map
      });
   </script> 
</body> 
</html>

截屏:

谷歌地图几何

更新:

作为对下面Paul评论的回复,这就是当圆圈环绕其中一个极点时会发生的情况。

pointA在北极附近绘制,半径为 1,000 公里:

  var pointA = new google.maps.LatLng(85, 0);   // Close to north pole
  var radius = 1000;                            // 1000km

截图pointA.destinationPoint(90, radius)

靠近北极

@Nirmal:不,这取决于您传递给的第一个参数destinationPoint()90 是东,但您可以使用任何方位,从 0 = 北开始,顺时针移动。
2021-03-16 12:55:24
@保罗:有趣的观察。我已经用一个例子更新了答案。Google 地图使用墨卡托投影 ( en.wikipedia.org/wiki/Mercator_projection ),随着比例从赤道到两极增加,在那里变得无限大,它会扭曲大型物体的大小和形状。使用上面的示例,var pointA = new google.maps.LatLng(85, 0);半径为 1000 公里var radius = 1000;似乎有效。圆帽环绕北极,因此在 2D 投影中它不再具有圆形形状。目标点似乎是正确的。
2021-03-20 12:55:24
@Daniel,对于最初的问题,这是一个高质量的答案。请保持好的工作状态!
2021-03-30 12:55:24
这会在两极附近工作吗?即它会产生一些在地图上看起来不是圆形的东西,但在地球上是圆形的吗?
2021-04-11 12:55:24
那么这是否总是给出位于东方的目的地点?
2021-04-12 12:55:24

要计算给定方位的纬度和经度点以及与另一个的距离,您可以使用谷歌的 JavaScript 实现:

var pointA = new google.maps.LatLng(25.48, -71.26); 
var distance = 10; // 10 metres
var bearing = 90; // 90 degrees
var pointB = google.maps.geometry.spherical.computeOffset(pointA, distance, bearing);

有关文档,请参阅https://developers.google.com/maps/documentation/javascript/reference#sphereal

简单,简短,干净:D
2021-03-21 12:55:24

如果您在地球表面上的 2 个纬度/经度点之间的距离之后,那么您可以在此处找到 javascript:

http://www.movable-type.co.uk/scripts/latlong-vincenty.html

这与 android API 中使用的公式相同 android.location.Location::distanceTo

您可以轻松地将代码从 javascript 转换为 java。

如果你想计算给定起点、方位和距离的目的地点,那么你需要这个方法:

http://www.movable-type.co.uk/scripts/latlong-vincenty-direct.html

下面是java中的公式:

public class LatLngUtils {

  /**
   * @param lat1
   *          Initial latitude
   * @param lon1
   *          Initial longitude
   * @param lat2
   *          destination latitude
   * @param lon2
   *          destination longitude
   * @param results
   *          To be populated with the distance, initial bearing and final
   *          bearing
   */

  public static void computeDistanceAndBearing(double lat1, double lon1,
      double lat2, double lon2, double results[]) {
    // Based on http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf
    // using the "Inverse Formula" (section 4)

    int MAXITERS = 20;
    // Convert lat/long to radians
    lat1 *= Math.PI / 180.0;
    lat2 *= Math.PI / 180.0;
    lon1 *= Math.PI / 180.0;
    lon2 *= Math.PI / 180.0;

    double a = 6378137.0; // WGS84 major axis
    double b = 6356752.3142; // WGS84 semi-major axis
    double f = (a - b) / a;
    double aSqMinusBSqOverBSq = (a * a - b * b) / (b * b);

    double L = lon2 - lon1;
    double A = 0.0;
    double U1 = Math.atan((1.0 - f) * Math.tan(lat1));
    double U2 = Math.atan((1.0 - f) * Math.tan(lat2));

    double cosU1 = Math.cos(U1);
    double cosU2 = Math.cos(U2);
    double sinU1 = Math.sin(U1);
    double sinU2 = Math.sin(U2);
    double cosU1cosU2 = cosU1 * cosU2;
    double sinU1sinU2 = sinU1 * sinU2;

    double sigma = 0.0;
    double deltaSigma = 0.0;
    double cosSqAlpha = 0.0;
    double cos2SM = 0.0;
    double cosSigma = 0.0;
    double sinSigma = 0.0;
    double cosLambda = 0.0;
    double sinLambda = 0.0;

    double lambda = L; // initial guess
    for (int iter = 0; iter < MAXITERS; iter++) {
      double lambdaOrig = lambda;
      cosLambda = Math.cos(lambda);
      sinLambda = Math.sin(lambda);
      double t1 = cosU2 * sinLambda;
      double t2 = cosU1 * sinU2 - sinU1 * cosU2 * cosLambda;
      double sinSqSigma = t1 * t1 + t2 * t2; // (14)
      sinSigma = Math.sqrt(sinSqSigma);
      cosSigma = sinU1sinU2 + cosU1cosU2 * cosLambda; // (15)
      sigma = Math.atan2(sinSigma, cosSigma); // (16)
      double sinAlpha = (sinSigma == 0) ? 0.0 : cosU1cosU2 * sinLambda
          / sinSigma; // (17)
      cosSqAlpha = 1.0 - sinAlpha * sinAlpha;
      cos2SM = (cosSqAlpha == 0) ? 0.0 : cosSigma - 2.0 * sinU1sinU2
          / cosSqAlpha; // (18)

      double uSquared = cosSqAlpha * aSqMinusBSqOverBSq; // defn
      A = 1 + (uSquared / 16384.0) * // (3)
          (4096.0 + uSquared * (-768 + uSquared * (320.0 - 175.0 * uSquared)));
      double B = (uSquared / 1024.0) * // (4)
          (256.0 + uSquared * (-128.0 + uSquared * (74.0 - 47.0 * uSquared)));
      double C = (f / 16.0) * cosSqAlpha * (4.0 + f * (4.0 - 3.0 * cosSqAlpha)); // (10)
      double cos2SMSq = cos2SM * cos2SM;
      deltaSigma = B
          * sinSigma
          * // (6)
          (cos2SM + (B / 4.0)
              * (cosSigma * (-1.0 + 2.0 * cos2SMSq) - (B / 6.0) * cos2SM
                  * (-3.0 + 4.0 * sinSigma * sinSigma)
                  * (-3.0 + 4.0 * cos2SMSq)));

      lambda = L
          + (1.0 - C)
          * f
          * sinAlpha
          * (sigma + C * sinSigma
              * (cos2SM + C * cosSigma * (-1.0 + 2.0 * cos2SM * cos2SM))); // (11)

      double delta = (lambda - lambdaOrig) / lambda;
      if (Math.abs(delta) < 1.0e-12) {
        break;
      }
    }

    double distance = (b * A * (sigma - deltaSigma));
    results[0] = distance;
    if (results.length > 1) {
      double initialBearing = Math.atan2(cosU2 * sinLambda, cosU1 * sinU2
          - sinU1 * cosU2 * cosLambda);
      initialBearing *= 180.0 / Math.PI;
      results[1] = initialBearing;
      if (results.length > 2) {
        double finalBearing = Math.atan2(cosU1 * sinLambda, -sinU1 * cosU2
            + cosU1 * sinU2 * cosLambda);
        finalBearing *= 180.0 / Math.PI;
        results[2] = finalBearing;
      }
    }
  }

  /*
   * Vincenty Direct Solution of Geodesics on the Ellipsoid (c) Chris Veness
   * 2005-2012
   * 
   * from: Vincenty direct formula - T Vincenty, "Direct and Inverse Solutions
   * of Geodesics on the Ellipsoid with application of nested equations", Survey
   * Review, vol XXII no 176, 1975 http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf
   */

  /**
   * Calculates destination point and final bearing given given start point,
   * bearing & distance, using Vincenty inverse formula for ellipsoids
   * 
   * @param lat1
   *          start point latitude
   * @param lon1
   *          start point longitude
   * @param brng
   *          initial bearing in decimal degrees
   * @param dist
   *          distance along bearing in metres
   * @returns an array of the desination point coordinates and the final bearing
   */

  public static void computeDestinationAndBearing(double lat1, double lon1,
      double brng, double dist, double results[]) {
    double a = 6378137, b = 6356752.3142, f = 1 / 298.257223563; // WGS-84
                                                                 // ellipsiod
    double s = dist;
    double alpha1 = toRad(brng);
    double sinAlpha1 = Math.sin(alpha1);
    double cosAlpha1 = Math.cos(alpha1);

    double tanU1 = (1 - f) * Math.tan(toRad(lat1));
    double cosU1 = 1 / Math.sqrt((1 + tanU1 * tanU1)), sinU1 = tanU1 * cosU1;
    double sigma1 = Math.atan2(tanU1, cosAlpha1);
    double sinAlpha = cosU1 * sinAlpha1;
    double cosSqAlpha = 1 - sinAlpha * sinAlpha;
    double uSq = cosSqAlpha * (a * a - b * b) / (b * b);
    double A = 1 + uSq / 16384
        * (4096 + uSq * (-768 + uSq * (320 - 175 * uSq)));
    double B = uSq / 1024 * (256 + uSq * (-128 + uSq * (74 - 47 * uSq)));
    double sinSigma = 0, cosSigma = 0, deltaSigma = 0, cos2SigmaM = 0;
    double sigma = s / (b * A), sigmaP = 2 * Math.PI;

    while (Math.abs(sigma - sigmaP) > 1e-12) {
      cos2SigmaM = Math.cos(2 * sigma1 + sigma);
      sinSigma = Math.sin(sigma);
      cosSigma = Math.cos(sigma);
      deltaSigma = B
          * sinSigma
          * (cos2SigmaM + B
              / 4
              * (cosSigma * (-1 + 2 * cos2SigmaM * cos2SigmaM) - B / 6
                  * cos2SigmaM * (-3 + 4 * sinSigma * sinSigma)
                  * (-3 + 4 * cos2SigmaM * cos2SigmaM)));
      sigmaP = sigma;
      sigma = s / (b * A) + deltaSigma;
    }

    double tmp = sinU1 * sinSigma - cosU1 * cosSigma * cosAlpha1;
    double lat2 = Math.atan2(sinU1 * cosSigma + cosU1 * sinSigma * cosAlpha1,
        (1 - f) * Math.sqrt(sinAlpha * sinAlpha + tmp * tmp));
    double lambda = Math.atan2(sinSigma * sinAlpha1, cosU1 * cosSigma - sinU1
        * sinSigma * cosAlpha1);
    double C = f / 16 * cosSqAlpha * (4 + f * (4 - 3 * cosSqAlpha));
    double L = lambda
        - (1 - C)
        * f
        * sinAlpha
        * (sigma + C * sinSigma
            * (cos2SigmaM + C * cosSigma * (-1 + 2 * cos2SigmaM * cos2SigmaM)));
    double lon2 = (toRad(lon1) + L + 3 * Math.PI) % (2 * Math.PI) - Math.PI; // normalise
                                                                             // to
                                                                             // -180...+180

    double revAz = Math.atan2(sinAlpha, -tmp); // final bearing, if required

    results[0] = toDegrees(lat2);
    results[1] = toDegrees(lon2);
    results[2] = toDegrees(revAz);

  }

  private static double toRad(double angle) {
    return angle * Math.PI / 180;
  }

  private static double toDegrees(double radians) {
    return radians * 180 / Math.PI;
  }

}
请注意:Vincenty 的公式在使用的椭球上精确到 0.5 毫米或 0.000015 英寸 (!) 以内。基于球形模型的计算,例如(更简单的)Haversine,精确到 0.3% 左右。所以之前的javascript解决方案可能是大多数人需要的。
2021-03-25 12:55:24

这个问题的答案和更多信息可以在这里找到:http : //www.edwilliams.org/avform.htm

我确实在该页面上看到了很多 acos、sin、tan 等。但我会坚持丹尼尔的回答。谢谢你的帮助。
2021-03-23 12:55:24
现在链接已断开。
2021-03-28 12:55:24
@ReinhardMänner 已修复。
2021-04-12 12:55:24

用于许多测地线计算(正反问题、面积计算等)的 Javascript。可在

http://geographiclib.sourceforge.net/scripts/geographiclib.js

示例用法显示在

http://geographiclib.sourceforge.net/scripts/geod-calc.html

谷歌地图的界面提供在

http://geographiclib.sourceforge.net/scripts/geod-google.html

这包括绘制测地线(蓝色)、测地线圆(绿色)和测地线包络线(红色)。

样本测地线从新西兰开始,绕世界 1.5 圈