UTM投影正算(求解)
根据书上公式计算出的结果与笑脸软件有差别，不知为何，求高人解答;;;UTM投影正算
(defun c:tt ( / B L L0 TYA TYB XY)
(setq TYa 6378137 ;WGS84 长半轴
TYb 6356752.3142451 ;WGS84 短半轴
B 18.0700 ;经度
L 110.1900 ;纬度
L0111 ;中央子午线
)
(setq B (BLTK_Du_Hu B)) ;以弧度为单位
;UTM投影自西经180°起每隔经差6度自西向东分带，第1带的中央经度为-177°，因此高斯-克吕格投影的第1带是UTM的第31带
;确保东经位于-180.00----179.9之间
;(setq L (- (+ L 180) (fix (* 360.0 (/ (+ L 180) 360.0))) 180))
;(setq L (BLTK_Du_Hu L)) ;以弧度为单位
(setq L (* L (/ pi 180.0)))
(setq L0 (* L0 (/ pi 180.0)))
(setq XY (UTMTY_BL_XY TYa TYb B L L0))
(princ XY)
(princ)
)
(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)
(setq FN 0) ;北纬偏移
(setq FE 500000) ;东纬偏移
(setq k0 0.9996)
;计算椭圆元素
(setq Tn (expt (/ (sin B) (cos B)) 2.0)) ;0.107041764
(setq e1 (sqrt (- 1.0 (expt (/ TYb TYa) 2.0))));第一偏心率 0.081819197
(setq e2 (sqrt (- (expt (/ TYa TYb) 2.0) 1.0)));第二偏心率 0.082094440
(setq C (* (expt e2 2.0) (expt (cos B) 2.0))) ;0.006087843
;计算正算系数
(setq A (* (- L L0) (cos B)))
(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))
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))
a4 (/ (* 35.0 (expt e1 6.0)) 3072.0)
)
(setq M (* TYa (- (+ (* a1 B) (* a3 (sin (* 4.0 B)))) (* a2 (sin (* 2.0 B))) (* a4 (sin (* 6.0 B))))))
(setq Nn (/ (/ (expt TYa 2.0) TYb) (sqrt (+ 1.0 (* (expt e2 2.0) (expt (cos B) 2.0)))))) ;卯酉圈半径
(setq b1 (* Nn (/ (sin B) (cos B)))
b2 (+ (- 5.0 Tn) (* 9.0 C) (* 4.0 (expt C 2.0)))
b3 (- (+ 61.0 (expt Tn 2.0) (* 600.0 C)) (* 58.0 Tn) (* 330.0 (expt e2 2.0)))
b4 (+ (- 1.0 Tn) C)
b5 (- (+ 5.0 (expt Tn 2.0) (* 72.0 C)) (* 18.0 Tn) (* 58.0 (expt e2 2.0)))
)
(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))))))
(setq Y (+ FE (* k0 Nn (+ A (* b4 (/ (expt A 3.0) 6.0)) (* b5 (/ (expt A 5.0) 120.0))))))
(list (rtos X 2 6) (rtos Y 2 6))
)
;;;将度.分秒形式化为弧度：输入为度.分秒形式，输出为弧度
(defun BLTK_Du_Hu (deg / du fe mi r)
(setq du(fix deg)
deg (* (- deg du) 100.0)
fe(fix deg)
mi(* (- deg fe) 100.0)
r (+ du (/ fe 60.0) (/ mi 3600.0))
r (* r (/ pi 180.0))
)
) 呵呵，我试过，结果是一样的，也是用你的程序 不用专门写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结果比对是正确的
页:
[1]