用KML在Google Earth上画圆的算法推导和程序实现

年初在做毕业设计时遇到过一个问题，就是给定一个点的经纬度坐标，给定半径的长度，在Google Earth上围绕这个点画一个圆。这个问题的难点在于，Google Earth上画几何对象要用到KML，而KML中没有圆这个对象，只能用LineString来画，这样问题就转换成了找到圆上点的经纬度坐标并将它们连成一个LineString，取得点越细密，LineString就越接近圆。接下来问题就变为如何取圆上的点、怎样求出圆上点的坐标。

首先看如何取圆上的点，取法当然是越平均越好了。圆上的点，当然就弧度平均分了，弧度对应着圆心角。圆心角有何圆心点处的方位角有关系。我们可以定义正北方向为0°方位角，沿顺时针方向方位角递增。这样，我们把就问题转化成了，已知一个点A的经纬度坐标，求一个距A点距离为r的方位角为ζ的点P的经纬度坐标。

首先，我们需要知道球面上两点间球面距离的定义。球面上两点A、B之间的距离为经过A、B的大圆上以A、B为端点的劣弧的弧长。AB间的球面距离与弧AB的圆心角成正比。

，其中

如何根据两个点的经纬度坐标来求出θ呢。这需要用将点的经纬度坐标(λ,φ)转换为直角坐标系下的坐标(x,y,z)。用到下面公式：

利用向量的内积

就可以求出cosθ的值，

进而得到θ的值。

下面回到求点P经纬度坐标的问题上来。通过A和P之间的距离的关系，我们可以列出一个等式：

但是要求出方程的解还需要利用方位角的关系再列一个方程。

这里我们要引入一个特殊的点P₀。P₀位于A点的正北方向且距A点的距离为θ，很容求出P点的坐标为(λ_A,φ_A+θ)。通过P₀和P点之间的距离可以列出一个等式来。但是PP₀的长度还是未知的。为了求出PP₀的值，我们要进行一下坐标转换，将A点看作球面坐标系的北极点，将AP₀所在的大圆弧为0度经线，则有P₀点坐标为(0,π/2-θ)，P点坐标为(ζ,π/2-θ)。通过两种PP₀长度的等式可以列出方程：

通过解方程可以求出P点的坐标：

知道了P点坐标计算公式，就可以通过编写程序来算出圆上点的坐标，生成KML文件。以下是用Python语言编写的代码。在实际的程序中，需要考虑一些极端情况，如0值判断、编程语言中三角函数的参数取值范围、弧度角度转化等情况。

import math

class PointOnEarth:
    tolerance=0.000000001
    R=6371300.0
    def __init__(self,lon=0.0,lat=0.0):
        self.lon=lon
        self.lat=lat
    def radToDeg(self,radian):
        return radian*180.0/math.pi
    def degToRad(self,degree):
        return math.pi*degree/180.0
    def disToRad(self,distance):
        return distance/self.R
    def radToDis(self,radian):
        return radian*self.R
    def distanceTo(self,p):
        return self.radToDis(
            math.acos(
                math.cos(
                    self.degToRad(self.lat)
                )*math.cos(
                    self.degToRad(p.lat)
                )*math.cos(
                    self.degToRad(self.lon-p.lon)
                )+math.sin(
                    self.degToRad(self.lat)
                )*math.sin(
                    self.degToRad(p.lat)
                )
            )
        )
    def getPointBydirection(self,direction,distance):
        direct=self.degToRad(direction)
        r=self.disToRad(distance)
        lon=self.degToRad(self.lon)
        lat=self.degToRad(self.lat)
        lonDiff = 0.0
        p= PointOnEarth()
        s = math.sin(r) * math.sin(r) * math.cos(direct) + math.cos(r)* math.cos(r)
        sinLatX = (s * math.cos(lat) - math.cos(r) * math.cos(lat + r))/ math.sin(r)
        if -sinLatX*sinLatX+1.0<self.tolerance*self.tolerance:
            if direct>(math.pi/2) and direct<(math.pi*3/2):
                p.lat=-math.pi/2
            else:
                p.lat=math.pi/2
            lonDiff = 0.0
        else:
            p.lat=math.asin(sinLatX)
            if (math.cos(lat)*math.cos(lat)<self.tolerance*self.tolerance):#Point p is polar.
                lonDiff=math.pi-direct
            else:# Point p is not polar.
                cosLonDiff=(math.cos(r)-math.sin(p.lat)*math.sin(lat))/(math.cos(p.lat)*math.cos(lat))
                if -cosLonDiff*cosLonDiff+1.0<self.tolerance*self.tolerance:
                    if cosLonDiff>0.0:
                        lonDiff = 0.0
                    else:
                        lonDiff=math.pi
                else:
                    lonDiff=math.acos(cosLonDiff)
                    if direct>math.pi:
                        lonDiff = -lonDiff
        p.lon=lon+lonDiff
        while p.lon>math.pi or p.lon<-math.pi:
            if p.lon>0.0:
                p.lon=p.lon-math.pi*2
            else:
                p.lon=p.lon+math.pi*2
        p.lon=self.radToDeg(p.lon)
        p.lat=self.radToDeg(p.lat)
        return p

if __name__=="__main__":
    p=PointOnEarth(45.0,45.0)
    print "<Folder><name>KML Circle</name>"
    print "<Placemark><name>circle</name>"
    print "<visibility>1</visibility>"
    print "<Style><geomColor>ff0000ff</geomColor>"
    print "<geomScale>1</geomScale></Style>"
    print "<LineString><coordinates>"
    for i in range(0,361,5):
        p2=p.getPointBydirection(i,p.R*math.pi/2.0)
        print "%.9f, %.9f" % (p2.lon,p2.lat)
    print "</coordinates></LineString></Placemark></Folder>"