UTM投影正算(求解)

gzxl

根据书上公式计算出的结果与笑脸软件有差别，不知为何，求高人解答

2. ;;;UTM投影正算
   
3. (defun c:tt ( / B L L0 TYA TYB XY)
   
4.   (setq TYa 6378137              ;WGS84 长半轴
   
5.         TYb 6356752.3142451      ;WGS84 短半轴
   
6.         B   18.0700              ;经度
   
7.         L   110.1900             ;纬度
   
8.         L0  111                  ;中央子午线
   
9.   )
   
10.   (setq B (BLTK_Du_Hu B))       ;以弧度为单位
    
11.   ;UTM投影自西经180°起每隔经差6度自西向东分带，第1带的中央经度为-177°，因此高斯-克吕格投影的第1带是UTM的第31带
    
12.   ;确保东经位于-180.00----179.9之间
    
13.   ;(setq L (- (+ L 180) (fix (* 360.0 (/ (+ L 180) 360.0))) 180))
    
14.   ;(setq L (BLTK_Du_Hu L))       ;以弧度为单位
    
15.   (setq L (* L (/ pi 180.0)))
    
16.   (setq L0 (* L0 (/ pi 180.0)))
    
17.   (setq XY (UTMTY_BL_XY TYa TYb B L L0))
    
18.   (princ XY)
    
19.   (princ)
    
20. )
    
21. (defun UTMTY_BL_XY (TYa TYb B L L0 / A a1 a2 a3 a4 b1 b2 b3 b4 b5 C e1 e2 FE FN M Nn Tn X Y k0)
    
22.   (setq FN 0)        ;北纬偏移
    
23.   (setq FE 500000)   ;东纬偏移
    
24.   (setq k0 0.9996)
    
25.   ;计算椭圆元素
    
26.   (setq Tn (expt (/ (sin B) (cos B)) 2.0)) ;0.107041764
    
27.   (setq e1 (sqrt (- 1.0 (expt (/ TYb TYa) 2.0))))  ;第一偏心率 0.081819197
    
28.   (setq e2 (sqrt (- (expt (/ TYa TYb) 2.0) 1.0)))  ;第二偏心率 0.082094440
    
29.   (setq C (* (expt e2 2.0) (expt (cos B) 2.0))) ;0.006087843
    
30.   ;计算正算系数
    
31.   (setq A (* (- L L0) (cos B)))
    
32.   (setq a1 (- 1.0 (/ (expt e1 2.0) 4.0) (/ (* 3.0 (expt e1 4.0)) 64.0) (/ (* 5.0 (expt e1 6.0)) 256.0))
    
33.         a2 (+ (/ (* 3.0 (expt e1 2.0)) 8.0) (/ (* 3.0 (expt e1 4.0)) 32.0) (/ (* 45.0 (expt e1 6.0)) 1024.0))
    
34.         a4 (/ (* 35.0 (expt e1 6.0)) 3072.0)
    
35.   )
    
36.   (setq M (* TYa (- (+ (* a1 B) (* a3 (sin (* 4.0 B)))) (* a2 (sin (* 2.0 B))) (* a4 (sin (* 6.0 B))))))
    
37.   (setq Nn (/ (/ (expt TYa 2.0) TYb) (sqrt (+ 1.0 (* (expt e2 2.0) (expt (cos B) 2.0)))))) ;卯酉圈半径
    
38.   (setq b1 (* Nn (/ (sin B) (cos B)))
    
39.         b2 (+ (- 5.0 Tn) (* 9.0 C) (* 4.0 (expt C 2.0)))
    
40.         b3 (- (+ 61.0 (expt Tn 2.0) (* 600.0 C)) (* 58.0 Tn) (* 330.0 (expt e2 2.0)))
    
41.         b4 (+ (- 1.0 Tn) C)
    
42.         b5 (- (+ 5.0 (expt Tn 2.0) (* 72.0 C)) (* 18.0 Tn) (* 58.0 (expt e2 2.0)))
    
43.   )
    
44.   (setq X (+ FN (* k0 (+ M (* b1 (+ (/ (expt A 2.0) 2.0) (* b2 (/ (expt A 4.0) 24.0)))) (* b3 (/ (expt A 6.0) 720.0))))))
    
45.   (setq Y (+ FE (* k0 Nn (+ A (* b4 (/ (expt A 3.0) 6.0)) (* b5 (/ (expt A 5.0) 120.0))))))
    
46.   (list (rtos X 2 6) (rtos Y 2 6))
    
47. )
    
48. ;;;将度.分秒形式化为弧度：输入为度.分秒形式，输出为弧度
    
49. (defun BLTK_Du_Hu (deg / du fe mi r)
    
50.   (setq du  (fix deg)
    
51.         deg (* (- deg du) 100.0)
    
52.         fe  (fix deg)
    
53.         mi  (* (- deg fe) 100.0)
    
54.         r   (+ du (/ fe 60.0) (/ mi 3600.0))
    
55.         r   (* r (/ pi 180.0))
    
56.   )
    
57. )

;;;UTM投影正算<br /> (defun c:tt ( / B L L0 TYA TYB XY)<br />   (setq TYa 6378137              ;WGS84 长半轴<br />         TYb 6356752.3142451      ;WGS84 短半轴<br />         B   18.0700              ;经度<br />         L   110.1900             ;纬度<br />         L0  111                  ;中央子午线<br />   )<br />   (setq B (BLTK_Du_Hu B))       ;以弧度为单位<br />   ;UTM投影自西经180°起每隔经差6度自西向东分带，第1带的中央经度为-177°，因此高斯-克吕格投影的第1带是UTM的第31带<br />   ;确保东经位于-180.00----179.9之间<br />   ;(setq L (- (+ L 180) (fix (* 360.0 (/ (+ L 180) 360.0))) 180))<br />   ;(setq L (BLTK_Du_Hu L))       ;以弧度为单位<br />   (setq L (* L (/ pi 180.0)))<br />   (setq L0 (* L0 (/ pi 180.0)))<br />   (setq XY (UTMTY_BL_XY TYa TYb B L L0))<br />   (princ XY)<br />   (princ)<br /> )<br /> (defun UTMTY_BL_XY (TYa TYb B L L0 / A a1 a2 a3 a4 b1 b2 b3 b4 b5 C e1 e2 FE FN M Nn Tn X Y k0)<br />   (setq FN 0)        ;北纬偏移<br />   (setq FE 500000)   ;东纬偏移<br />   (setq k0 0.9996)<br />   ;计算椭圆元素<br />   (setq Tn (expt (/ (sin B) (cos B)) 2.0)) ;0.107041764<br />   (setq e1 (sqrt (- 1.0 (expt (/ TYb TYa) 2.0))))  ;第一偏心率 0.081819197<br />   (setq e2 (sqrt (- (expt (/ TYa TYb) 2.0) 1.0)))  ;第二偏心率 0.082094440<br />   (setq C (* (expt e2 2.0) (expt (cos B) 2.0))) ;0.006087843<br />   ;计算正算系数<br />   (setq A (* (- L L0) (cos B)))<br />   (setq a1 (- 1.0 (/ (expt e1 2.0) 4.0) (/ (* 3.0 (expt e1 4.0)) 64.0) (/ (* 5.0 (expt e1 6.0)) 256.0))<br />         a2 (+ (/ (* 3.0 (expt e1 2.0)) 8.0) (/ (* 3.0 (expt e1 4.0)) 32.0) (/ (* 45.0 (expt e1 6.0)) 1024.0))<br />         a4 (/ (* 35.0 (expt e1 6.0)) 3072.0)<br />   )<br />   (setq M (* TYa (- (+ (* a1 B) (* a3 (sin (* 4.0 B)))) (* a2 (sin (* 2.0 B))) (* a4 (sin (* 6.0 B))))))<br />   (setq Nn (/ (/ (expt TYa 2.0) TYb) (sqrt (+ 1.0 (* (expt e2 2.0) (expt (cos B) 2.0)))))) ;卯酉圈半径<br />   (setq b1 (* Nn (/ (sin B) (cos B)))<br />         b2 (+ (- 5.0 Tn) (* 9.0 C) (* 4.0 (expt C 2.0)))<br />         b3 (- (+ 61.0 (expt Tn 2.0) (* 600.0 C)) (* 58.0 Tn) (* 330.0 (expt e2 2.0)))<br />         b4 (+ (- 1.0 Tn) C)<br />         b5 (- (+ 5.0 (expt Tn 2.0) (* 72.0 C)) (* 18.0 Tn) (* 58.0 (expt e2 2.0)))<br />   )<br />   (setq X (+ FN (* k0 (+ M (* b1 (+ (/ (expt A 2.0) 2.0) (* b2 (/ (expt A 4.0) 24.0)))) (* b3 (/ (expt A 6.0) 720.0))))))<br />   (setq Y (+ FE (* k0 Nn (+ A (* b4 (/ (expt A 3.0) 6.0)) (* b5 (/ (expt A 5.0) 120.0))))))<br />   (list (rtos X 2 6) (rtos Y 2 6))<br /> )<br /> ;;;将度.分秒形式化为弧度：输入为度.分秒形式，输出为弧度<br /> (defun BLTK_Du_Hu (deg / du fe mi r)<br />   (setq du  (fix deg)<br />         deg (* (- deg du) 100.0)<br />         fe  (fix deg)<br />         mi  (* (- deg fe) 100.0)<br />         r   (+ du (/ fe 60.0) (/ mi 3600.0))<br />         r   (* r (/ pi 180.0))<br />   )<br /> )

wangph

呵呵，我试过，结果是一样的，也是用你的程序

wzg356

不用专门写utm正反算了
就用你写的高斯正反算直接转
高斯正算结果bl-->xy  即(list x y)转换就行
即
近似采用 Xutm=0.9996 * X高斯，Yutm=0.9996 * Y高斯   进行坐标转换
这个就是utm正算结果(list (*(car xy)0.9996) (+ (* (- (cadr xy)500000)0.9996)500000))
已与国外的软件global mapper结果比对是正确的