用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的圆心角成正比。

d_{AB}=R\cdot\theta,其中\theta=\angle{AOB}

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

x_P=R\cdot\cos{\varphi_P}\cdot\cos{\lambda_P}\\y_P=R\cdot\cos{\varphi_P}\cdot\sin{\lambda_P}\\z_P=R\cdot\sin{\varphi_P}

利用向量的内积

\vec{OA}\cdot\vec{OB}=|\vec{OA}|\:|\vec{OB}|\:\cos{\theta}

就可以求出cosθ的值,

\cos{\theta}=\cos\varphi_A\cos{\varphi}_B\cos{(\lambda_A-\lambda_B)}+\sin{\varphi_A}\sin{\varphi_B}

进而得到θ的值。

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

\cos{\theta}=\cos\varphi_A\cos{\varphi}_P\cos{(\lambda_A-\lambda_P)}+\sin{\varphi_A}\sin{\varphi_P}

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

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

\sin^2{\theta}\cos{\zeta}+\cos^2{\theta}=\cos{(\varphi_{A}+\theta)}\cos{\varphi_P}\cos{(\lambda_{A}-\lambda_P)}+\sin{(\varphi_{A}+\theta)}\sin{\varphi_P}

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

\sin{\varphi_P}=\frac{(\sin^2{\theta}\cos{\zeta}+\cos^2{\theta})\cos{\varphi_A}-\cos{\theta}\cos{(\varphi_A+\theta)}}{\sin{\theta}}

\cos{(\lambda_P-\lambda_A)}=\frac{\cos{\theta}-\sin{\varphi_A}\sin{\varphi_P}}{\cos{\varphi_A}\cos{\varphi_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>"

《用KML在Google Earth上画圆的算法推导和程序实现》有3个想法

  1. 我的数学不好,但是我也想出了一个画扇形的算式。代码如下(javascript):
    //画一个扇形 半径,初始角度,旋转角度,纬度,经度,颜色,姓名
    function DrawShanxing(r, initDegree, Degree, lat, lon, color, name) {
    //创建一个画板可以让你开始画画
    var polygonPlacemark = ge.createPlacemark(”);
    polygonPlacemark.setGeometry(ge.createPolygon(”));
    var outer = ge.createLinearRing(”);
    polygonPlacemark.getGeometry().setOuterBoundary(outer);
    //设置扇形的颜色
    var style = ge.createStyle(name);
    style.getPolyStyle().getColor().set(color);
    style.getLineStyle().getColor().set(’00ffffff’);
    polygonPlacemark.setStyleSelector(style);

    //为画板添加事件
    google.earth.addEventListener(polygonPlacemark, ‘click’, function(event) { try{GAjax.getPlotInfo(name,callback);}catch(err){alert(err.description )} });
    //将画板加入到地图中
    ge.getFeatures().appendChild(polygonPlacemark);

    var coords = outer.getCoordinates();
    var bLat, bLng;
    //画一条边
    bLng = Math.cos(initDegree * 3.1415926 / 180) * r;
    bLat = Math.sin(initDegree * 3.1415926 / 180) * r;

    coords.pushLatLngAlt(lat + bLat, lon + bLng, 0);

    //画顶点

    coords.pushLatLngAlt(lat, lon, 0);
    //画二条边
    bLng = Math.cos(initDegree * 3.1415926 / 180 + Degree * 3.1415926 / 180) * r;
    bLat = Math.sin(initDegree * 3.1415926 / 180 + Degree * 3.1415926 / 180) * r;
    coords.pushLatLngAlt(lat + bLat, lon + bLng, 0);

    for (var i = 15; i > 0; i–) {
    bLng = Math.cos(initDegree *Math.PI / 180 + Degree / 15 * Math.PI / 180 * i) * r;
    bLat = Math.sin(initDegree * Math.PI / 180 + Degree / 15 * Math.PI / 180 * i) * r;
    coords.pushLatLngAlt(lat + bLat, lon + bLng, 0);
    }

    }

    但是仍然有一个问题没有找到解决方法,就是由于一经度的距离不等于一纬度的距离 这导致画出来的都不是正规的圆形,我想请问这个有没有什么好的解决方法。

    谢谢,有回复的话请回tianma8778@126.com

  2. 看了下你的代码,你在计算坐标的时候用的是平面直角坐标系下圆上坐标的计算方法,而带入的是球面坐标系下的经纬度坐标,所以会出问题。
    文章中已经给出了通过经纬度坐标求球面上两点间距离的公式以及已知一点坐标、方位角、距离求另一点坐标的计算公式,用来画扇也是完全可以的。

发表评论

电子邮件地址不会被公开。 必填项已用*标注

Time limit is exhausted. Please reload CAPTCHA.