經緯度和球體三維坐标換算

引言

在球面坐标系上，坐标是用經緯度标示，緯線和經線一樣是為度量友善而假設出來的輔助線，定義為地球表面某點随地球自轉所形成的軌迹。

任何一根緯線都是圓形而且兩兩平行。緯線的長度是赤道的周長乘以緯線的緯度的餘弦，是以赤道最長，離赤道越遠的緯線，周長越短，到了兩極就縮為0。從赤道向北和向南，各分90°，稱為北緯和南緯，分别用“N”和“S”表示。

在實際的計算中，經緯度一般是不能拿來直接計算的，都需要進行一定的轉換，比方說計算距離等，假設把地球是一個正圓的，放在一個很大三維坐标系中，原心為地球的球心，坐标換算成三維的（笛卡爾）坐标數字，就友善計算了。

C#算法代碼

//設定地球半徑
 private  double WGS_84_RADIUS_EQUATOR = 6378137.0; 
//經緯度轉立體坐标
     public CVector WGS84ToWorld(CVector _wgscood)
        {
            if (_wgscood.VZ == 0 || Double.IsNaN(_wgscood.VZ))
            {
                _wgscood.VZ = 1;
            }
            _wgscood.VZ += WGS_84_RADIUS_EQUATOR;
            double tempx = 0;
            double tempy = 0;
            double tempz = 0;
            tempy = _wgscood.VZ * Math.Sin(_wgscood.VY / 180 * Math.PI);
            tempx = _wgscood.VZ * Math.Cos(_wgscood.VY / 180 * Math.PI) * Math.Sin(_wgscood.VX / 180 * Math.PI);
            tempz = _wgscood.VZ * Math.Cos(_wgscood.VY / 180 * Math.PI) * Math.Cos(_wgscood.VX / 180 * Math.PI);
            CVector rexyz = new CVector { VX = tempx, VY = tempy, VZ = tempz };
            return rexyz;
        }
     //立體坐标轉回經緯度
        public CVector WorldToWGS84(CVector _worldcood)
        {
            double tempx = 0;
            double tempy = 0;
            double tempz = 0;
            if (_worldcood.VX >= 0)
            {
                tempx = Math.Acos(_worldcood.VZ / Math.Sqrt(_worldcood.VX * _worldcood.VX + _worldcood.VZ * _worldcood.VZ)) * 180 / Math.PI;
                tempx = tempx < 0 ? tempx + 180 : tempx;
            }
            else
            {
                tempx = Math.Acos(_worldcood.VZ / Math.Sqrt(_worldcood.VX * _worldcood.VX + _worldcood.VZ * _worldcood.VZ)) * 180 / Math.PI;
                tempx = tempx < 0 ? -(tempx + 180) : -tempx;
            }
            tempy = Math.Asin(_worldcood.VY / Math.Sqrt(_worldcood.VX * _worldcood.VX + _worldcood.VY * _worldcood.VY + _worldcood.VZ * _worldcood.VZ)) * 180 / Math.PI;
            tempy = tempy > 90 ? 90 - tempy : tempy;

            tempz = (Math.Sqrt(_worldcood.VX * _worldcood.VX + _worldcood.VY * _worldcood.VY + _worldcood.VZ * _worldcood.VZ)) - WGS_84_RADIUS_EQUATOR;
            CVector relonlat = new CVector { VX = tempx, VY = tempy, VZ = tempz };
            return relonlat;
        }