四参数函数及应用-转点及全图坐标系转换
本帖最后由 wzg356 于 2020-11-12 11:27 编辑问题的提出，老老帖了，最近觉得这类资料少，算是补充
http://bbs.mjtd.com/forum.php?mo ... =%CB%C4%B2%CE%CA%FD
已知这四参数
纵坐标平移参数X=59.279033
横坐标平移参数Y=59.191856
旋转角度T=0.0001852315256（度）
伸缩尺度K=0.999985766528
可以写一个这样的程序来进行转换吗？因为转出来后几公分偏差，不知道写是不是正确的
(defun c:xzb()
(setq sicanshu (ssget) )
(command"move" sicanshu "" "0,0" "59.279033,59.191856" "")
(command"_scale" sicanshu "" "0,0" "0.999985766528" "")
(command"_rotate" sicanshu "" "0,0" "0.0001852315256" "")
)
;写个函数验证一下
;四参数获取，公共点对数量需对应
(defun pts24cs (oldps newps / ps1 ps2 p01 p02 ls1 ls2 ag1 ag2 k rtdx dy)
(setq ps1 (mapcar '(lambda (a) (mapcar '+ (list 0.0 0.0)a)) oldps)
ps2 (mapcar '(lambda (a) (mapcar '+ (list 0.0 0.0)a)) newps)
);转二维
(while(and(setq p01(car ps1)ps1(cdr ps1))
(setq p02(car ps2)ps2(cdr ps2))
)
(setq ls1(cons (mapcar '(lambda(a)(distance p01 a)) ps1) ls1))
(setq ag1(cons (mapcar '(lambda(a)(angle p01 a)) ps1) ag1))
(setq ls2(cons (mapcar '(lambda(a)(distance p02 a)) ps2)ls2))
(setq ag2(cons (mapcar '(lambda(a)(angle p02 a)) ps2)ag2))
);所有可能的两点组合的距离及方位角
(setq ls1 (apply 'append ls1) ag1(apply 'append ag1))
(setq ls2(apply 'append ls2) ag2(apply 'append ag2))
(setq k(/ (apply '+ (mapcar '/ ls2 ls1)) (length ls1)));平均k
(setq rt(/ (apply '+(mapcar '- ag2 ag1)) (length ag1)));平均旋转弧度
(setq oldps (mapcar '(lambda (a)(Rot2D a rt)) oldps));老点旋转后
(setq oldps (mapcar '(lambda (a)(mapcar '* (list k k)a)) oldps));乘k
(setq dxys(mapcar '(lambda(a b)(mapcar '- a b)) newps oldps))
(setq dx (/ (apply '+(mapcar 'car dxys)) (length (car dxys))));平均dx
(setq dy (/ (apply '+(mapcar 'cadr dxys)) (length (cadr dxys))));平均dy
(list dx dy rt k)
)
;;; 旋转向量到指定角度 ;;; 输入: 一个向量和指定的角度
;;; 输出: 被旋转后的向量
(defun Rot2D (v a / c s x y)
(setq c (cos a) s (sin a))
(setq x (car v) y (cadr v))
(list (- (* x c) (* y s)) (+ (* x s) (* y c)))
)
;应用及验证
;曾用过的点对-公共已知点 ， 最少2个点对
(setq oldps (list(list 5087067.129 532083.741)(list 5087050.578 532228.286));老坐标系
newps (list(list 5087679.709 609542.699)(list 5087664.971 609687.453));新坐标系
)
;现场实测第三点对-公共已知点
(setq oldp (list 5086983.205 532095.119) newp (list 5087595.926 609555.129))
;以上三个点对因保密原因，关键地理信息去掉了
(setq 4cs(pts24cs oldps newps));取得四参数
(setq newp(mapcar '+ (list (car 4cs)(cadr 4cs))(mapcar '* (list (cadddr 4cs)(cadddr 4cs)) (Rot2D oldp (caddr 4cs)))))
;转换第三个点
(list (rtos (car newp)2 5)(rtos (cadr newp)2 5))=>("5087595.92692" "609555.12984")与现场一致
;与其他软件计算结果也基本一致
;验证一下图面操作
(defun c:yz()
(setq oldps (list(list 5087067.129 532083.741)(list 5087050.578 532228.286))
newps (list(list 5087679.709 609542.699)(list 5087664.971 609687.453))
);坐标点对-公共点，老-新坐标系数据
(setq 4cs(pts24cs oldps newps));;取得四参数(list dx dy rt k)
(foreach cs 4cs (princ (strcat "\n"(rtos cs 2 10))))
(setq en(entmakex(list '(0 . "line") (cons 10 (car oldps)) (cons 11 (cadr oldps)))))
;老坐标画的线，黑
(command"_rotate" en "" "0,0" (* (/ (caddr 4cs) pi)180));弧度转为度
(command"_scale" en "" "0,0" (cadddr 4cs))
(command"move" en "" "0,0" (list (car 4cs)(cadr 4cs)))
(entmakex(list '(0 . "line") '(62 . 1)(cons 10 (car newps)) (cons 11 (cadr newps))))
;新坐标画的线，红
;两线重合
)
;说明依据四参数，可以整图转换的
；需要注意的是，转换图形的时候，command “move”是最后一步
很棒，感谢分享，学习一下
页:
[1]