如何計算以GPS座標爲中心的地球上的一個圓上的點？

Question

在KML中畫一個圓

如何獲取地球上某點的GPS座標（比如說十進制格式）並生成一個多邊形的座標，該座標近似以該點爲中心的圓？

具有20+個數據點的多邊形看起來像一個圓圈。數據點越多 - 看起來越圓。

我正在寫一個程序，它將生成KML，並且不知道如何計算多邊形頂點的座標。

的數據輸入示例：

經緯度，圓弧半徑（以英尺爲單位），NumberOfDataPoints

26.128477，-80.105149，500個，20

Answer 1

我不知道這是否是最簡單的解決方案，它假定世界是一個球體。

定義：

R是球體（即地球）的半徑。

r是圓的半徑（以相同的單位）。

t是球體中心長度爲r的大圓弧所對應的角度，因此t = r/R弧度。

現在假設球體的半徑爲1並且以原點爲中心。

C是代表圓的中心的單位向量。

想象一下圍繞北極的一個圓，並考慮圓的平面與從地球中心到北極的直線相交的點。顯然這一點將在北極以下的某個地方。

K是對應點「下方」 C（即，其中的圓的平面相交C）所以K = COS（T）C

s是在三維空間中測量的圓的半徑（即不在球體上），所以S = SIN（t）的

現在，我們希望在3D空間中的圓圈中心K，半徑S和點躺在通過並垂直於K.

This answer平面（忽略旋轉的東西）解釋瞭如何找到平面的基本矢量（即與正常的K或C正交的矢量）。使用交叉產品找到第二個。

稱這些基向量U和V

// Pseudo-code to calculate 20 points on the circle 
for (a = 0; a != 360; a += 18) 
{ 
    // A point on the circle and the unit sphere 
    P = K + s * (U * sin(a) + V * cos(a)) 
} 

轉換每個點球面座標和你做。

很無聊，我用C＃編碼了這個。結果是合理的：他們在一個圈子裏，躺在球體上。大多數代碼實現代表向量的struct。實際計算非常簡單。

using System; 

namespace gpsCircle 
{ 
    struct Gps 
    { 
     // In degrees 
     public readonly double Latitude; 
     public readonly double Longtitude; 

     public Gps(double latitude, double longtitude) 
     { 
      Latitude = latitude; 
      Longtitude = longtitude; 
     } 

     public override string ToString() 
     { 
      return string.Format("({0},{1})", Latitude, Longtitude); 
     } 

     public Vector ToUnitVector() 
     { 
      double lat = Latitude/180 * Math.PI; 
      double lng = Longtitude/180 * Math.PI; 

      // Z is North 
      // X points at the Greenwich meridian 
      return new Vector(Math.Cos(lng) * Math.Cos(lat), Math.Sin(lng) * Math.Cos(lat), Math.Sin(lat)); 
     } 
    } 

    struct Vector 
    { 
     public readonly double X; 
     public readonly double Y; 
     public readonly double Z; 

     public Vector(double x, double y, double z) 
     { 
      X = x; 
      Y = y; 
      Z = z; 
     } 

     public double MagnitudeSquared() 
     { 
      return X * X + Y * Y + Z * Z; 
     } 

     public double Magnitude() 
     { 
      return Math.Sqrt(MagnitudeSquared()); 
     } 

     public Vector ToUnit() 
     { 
      double m = Magnitude(); 

      return new Vector(X/m, Y/m, Z/m); 
     } 

     public Gps ToGps() 
     { 
      Vector unit = ToUnit(); 
      // Rounding errors 
      double z = unit.Z; 
      if (z > 1) 
       z = 1; 

      double lat = Math.Asin(z); 

      double lng = Math.Atan2(unit.Y, unit.X); 

      return new Gps(lat * 180/Math.PI, lng * 180/Math.PI); 
     } 

     public static Vector operator*(double m, Vector v) 
     { 
      return new Vector(m * v.X, m * v.Y, m * v.Z); 
     } 

     public static Vector operator-(Vector a, Vector b) 
     { 
      return new Vector(a.X - b.X, a.Y - b.Y, a.Z - b.Z); 
     } 

     public static Vector operator+(Vector a, Vector b) 
     { 
      return new Vector(a.X + b.X, a.Y + b.Y, a.Z + b.Z); 
     } 

     public override string ToString() 
     { 
      return string.Format("({0},{1},{2})", X, Y, Z); 
     } 

     public double Dot(Vector that) 
     { 
      return X * that.X + Y * that.Y + Z * that.Z; 
     } 

     public Vector Cross(Vector that) 
     { 
      return new Vector(Y * that.Z - Z * that.Y, Z * that.X - X * that.Z, X * that.Y - Y * that.X); 
     } 

     // Pick a random orthogonal vector 
     public Vector Orthogonal() 
     { 
      double minNormal = Math.Abs(X); 
      int minIndex = 0; 
      if (Math.Abs(Y) < minNormal) 
      { 
       minNormal = Math.Abs(Y); 
       minIndex = 1; 
      } 
      if (Math.Abs(Z) < minNormal) 
      { 
       minNormal = Math.Abs(Z); 
       minIndex = 2; 
      } 

      Vector B; 
      switch (minIndex) 
      { 
       case 0: 
        B = new Vector(1, 0, 0); 
        break; 
       case 1: 
        B = new Vector(0, 1, 0); 
        break; 
       default: 
        B = new Vector(0, 0, 1); 
        break; 
      } 

      return (B - minNormal * this).ToUnit(); 
     } 
    } 

    class Program 
    { 
     static void Main(string[] args) 
     { 
      // Phnom Penh 
      Gps centre = new Gps(11.55, 104.916667); 

      // In metres 
      double worldRadius = 6371000; 
      // In metres 
      double circleRadius = 1000; 

      // Points representing circle of radius circleRadius round centre. 
      Gps[] points = new Gps[20]; 

      CirclePoints(points, centre, worldRadius, circleRadius); 
     } 

     static void CirclePoints(Gps[] points, Gps centre, double R, double r) 
     { 
      int count = points.Length; 

      Vector C = centre.ToUnitVector(); 
      double t = r/R; 
      Vector K = Math.Cos(t) * C; 
      double s = Math.Sin(t); 

      Vector U = K.Orthogonal(); 
      Vector V = K.Cross(U); 
      // Improve orthogonality 
      U = K.Cross(V); 

      for (int point = 0; point != count; ++point) 
      { 
       double a = 2 * Math.PI * point/count; 
       Vector P = K + s * (Math.Sin(a) * U + Math.Cos(a) * V); 
       points[point] = P.ToGps(); 
      } 
     } 
    } 
} 

Answer 2

我寫了Polycircles，它是Python中的一個小型開源軟件包。它使用geographiclib進行地理空間計算。