从GPS模块中获得的坐标是属于球面坐标，并不是我们平常接触到的平面坐标，所以在计算两点间距离的时候不能简单的利用勾股定理来计算，而是要将GPS中的坐标转换为常用的平面坐标。由于地球是一个不规则的椭圆，精确计算方法很复杂，计算量也非常大。由于我们平常使用中不需要非常精确的转换，这里就利用一个近似的计算方法来进行转换，经过测试对比，比百度的计算方式更精确。（关于坐标转换的更多知识可以百度上查找，本计算方法也是在别人的基础上修改而来。）
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 |
type // 笛卡尔坐标系 tagCRDCARTESIAN = record X, Y, z: double; end; CRDCARTESIAN = tagCRDCARTESIAN; // 大地坐标系 tagCRDGEODETIC = record longitude: double; // 经度 latitude: double; // 纬度 height: double; // 大地高,可设为0 end; CRDGEODETIC = tagCRDGEODETIC; function GeodeticToCartesian(pcg: CRDGEODETIC): CRDCARTESIAN; var b: double; // 纬度度数 L: double; // 经度度数 L0: double; // 中央经线度数 L1: double; // L-L0 t: double; // tanB m: double; // ltanB N: double; // 卯酉圈曲率半径 q2: double; X: double; // 高斯平面纵坐标 Y: double; // 高斯平面横坐标 s: double; // 赤道至纬度B的经线弧长 f: double; // 参考椭球体扁率 e2: double; // 椭球第一偏心率 a: double; // 参考椭球体长半轴 // b:double; //参考椭球体短半轴 a1: double; a2: double; a3: double; a4: double; b1: double; b2: double; b3: double; b4: double; c0: double; c1: double; c2: double; c3: double; Datum, prjno, zonewide: Integer; IPI: double; begin Datum := 84; // 投影基准面类型：北京54基准面为54，西安80基准面为80，WGS84基准面为84 prjno := 0; // 投影带号 zonewide := 3; IPI := 0.0174532925199433333333; // 3.1415926535898/180.0 b := pcg.latitude; // 纬度 L := pcg.longitude; // 经度 if (zonewide = 6) then begin prjno := trunc(L / zonewide) + 1; L0 := prjno * zonewide - 3; end else begin prjno := trunc((L - 1.5) / 3) + 1; L0 := prjno * 3; end; if (Datum = 54) then begin a := 6378245; f := 1 / 298.3; end else if (Datum = 84) then begin a := 6378137; f := 1 / 298.257223563; end; L0 := L0 * IPI; L := L * IPI; b := b * IPI; e2 := 2 * f - f * f; // (a*a-b*b)/(a*a); L1 := L - L0; t := Tan(b); m := L1 * cos(b); N := a / sqrt(1 - e2 * sin(b) * sin(b)); q2 := e2 / (1 - e2) * cos(b) * cos(b); a1 := 1 + 3 / 4 * e2 + 45 / 64 * e2 * e2 + 175 / 256 * e2 * e2 * e2 + 11025 / 16384 * e2 * e2 * e2 * e2 + 43659 / 65536 * e2 * e2 * e2 * e2 * e2; a2 := 3 / 4 * e2 + 15 / 16 * e2 * e2 + 525 / 512 * e2 * e2 * e2 + 2205 / 2048 * e2 * e2 * e2 * e2 + 72765 / 65536 * e2 * e2 * e2 * e2 * e2; a3 := 15 / 64 * e2 * e2 + 105 / 256 * e2 * e2 * e2 + 2205 / 4096 * e2 * e2 * e2 * e2 + 10359 / 16384 * e2 * e2 * e2 * e2 * e2; a4 := 35 / 512 * e2 * e2 * e2 + 315 / 2048 * e2 * e2 * e2 * e2 + 31185 / 13072 * e2 * e2 * e2 * e2 * e2; b1 := a1 * a * (1 - e2); b2 := -1 / 2 * a2 * a * (1 - e2); b3 := 1 / 4 * a3 * a * (1 - e2); b4 := -1 / 6 * a4 * a * (1 - e2); c0 := b1; c1 := 2 * b2 + 4 * b3 + 6 * b4; c2 := -(8 * b3 + 32 * b4); c3 := 32 * b4; s := c0 * b + cos(b) * (c1 * sin(b) + c2 * sin(b) * sin(b) * sin(b) + c3 * sin (b) * sin(b) * sin(b) * sin(b) * sin(b)); X := s + 1 / 2 * N * t * m * m + 1 / 24 * (5 - t * t + 9 * q2 + 4 * q2 * q2) * N * t * m * m * m * m + 1 / 720 * (61 - 58 * t * t + t * t * t * t) * N * t * m * m * m * m * m * m; Y := N * m + 1 / 6 * (1 - t * t + q2) * N * m * m * m + 1 / 120 * (5 - 18 * t * t + t * t * t * t - 14 * q2 - 58 * q2 * t * t) * N * m * m * m * m * m; Y := Y + 1000000 * prjno + 500000; Result.X := X; Result.Y := Y - 38000000; Result.z := 0; end; |
