高斯投影正反算

hothua2020 · 发表于 2020-8-20 11:00:22

本帖最后由 hothua2020 于 2020-9-6 15:28 编辑

;;;功能：投影换带
;;作者:HUA
;;;创建日期：2020-8-18
;;;最后编辑日期:2020-8-20
;;;
;;;秒转度分秒 s->dd.mmss
(defun Second2Dms(second / sign degree miniute second dms)
(setq sign (< second 0)
second (/ (abs second) 3600.0)
degree (fix second)
second (* (- second degree) 60.0)
miniute (fix second)
second (* (- second miniute) 60.0)
dms (+ degree (/ miniute 100.0) (/ second 10000.0))
)
(if sign (setq dms (- 0 dms)))
dms
)
;;;
;;;角度转弧度 dd.mmss->rad
(defun Dms2Rad(dms / sign degree miniute second rad)
(setq sign (< dms 0)
dms (abs dms)
degree (fix dms)
dms (* (- dms degree) 100)
miniute (fix dms)
second (* (- dms miniute) 100)
rad (/ (* (+ degree (/ miniute 60.0) (/ second 3600.0)) PI) 180.0)
)
(if sign (setq rad (- 0 rad)))
rad
)
;;;
;;;弧度转角度 rad->dd.mmss
(defun Rad2Dms(rad / sign degree miniute second dms)
(setq sign (< rad 0)
rad (abs rad)
rad (/ (* rad 180.0) PI)
degree (fix rad)
rad (* (- rad degree) 60.0)
miniute (fix rad)
second (* (- rad miniute) 60.0)
dms (+ degree (/ miniute 100.0) (/ second 10000.0))
)
(if sign (setq dms (- 0 dms)))
dms
)
;;;
;;;十进制角度转弧度 dd->rad
(defun Degree2Rad(degree / sign rad)
(setq sign (< degree 0)
degree (abs degree )
rad (/ (* degree PI) 180.0)
)
(if sign (setq rad (- 0 rad)))
rad
)
;;;
;;;弧度转十进制角度 rad->dd
(defun Rad2Degree(rad / sign degree)
(setq sign (< rad 0)
rad (abs rad)
degree (/ (* rad 180.0) PI)
)
(if sign (setq degree (- 0 degree)))
degree
)
;;;
;;;54坐标系参数
(defun PrjPoint_Krasovsky( / a f e2 e12 A1 A2 A3 A4 A5 A6 A7 tf )
(setq a 6378245.0
f 298.3
tf (/ (- f 1) f)
e2 (- 1 (* tf tf))
tf (/ f (- f 1))
e12 (- (* tf tf) 1)
A1 6367558.4968749800000000
A2 -16036.4802694116000000
A3 16.8280668841393000
A4 -0.0219753089548841
A5 0.0000311310026601
A6 -0.0000000459827472
A7 0.0000000000679857
)
(list a f e2 e12 A1 A2 A3 A4 A5 A6 A7)
)
;;;
;;;80坐标系参数
(defun PrjPoint_IUGG1975(/ a f e2 e12 A1 A2 A3 A4 A5 A6 A7 tf)
(setq a 6378140.0
f 298.257
tf (/ (- f 1) f)
e2 (- 1 (* tf tf))
tf (/ f (- f 1))
e12 (- (* tf tf) 1)
A1 6367452.1327884300000000
A2 -16038.5282286966000000
A3 16.8326464353061000
A4 -0.0219844636486152
A5 0.0000311484690908
A6 -0.0000000460151868
A7 0.0000000000680433
)
(list a f e2 e12 A1 A2 A3 A4 A5 A6 A7)
)
;;;
;;;2000坐标系参数
(defun PrjPoint_2000(/ a f e2 e12 A1 A2 A3 A4 A5 A6 A7 tf)
(setq a 6378137.0
f 298.257222101
tf (/ (- f 1) f)
e2 (- 1 (* tf tf))
tf (/ f (- f 1))
e12 (- (* tf tf) 1)
A1 6367449.1457710400000000
A2 -16038.5087415914000000
A3 16.8326134276660000
A4 -0.0219844041401628
A5 0.0000311483615430
A6 -0.0000000460149936
A7 0.0000000000680429
)
(list a f e2 e12 A1 A2 A3 A4 A5 A6 A7)
)
;;;
;;;BL转cadx,cady
;;;输入参数
;;;B,L为点位纬度和经度，单位为弧度
;;;L0为中央经线，单位为弧度
;;;prjType取值为0:北京54坐标系,1:西安80坐标系，2：CGCS2000坐标系
;;;输出CAD坐标x,y
(defun PrjPoint:BL2xy(B L L0 prjType / a f e2 e12 A1 A2 A3 A4 A5 A6 A7 X N tanB tanB2 m m2 ng2 sinB cosB cadx cady)
(cond
((= 0 prjType) (mapcar 'set (list 'a 'f 'e2 'e12 'A1 'A2 'A3 'A4 'A5 'A6 'A7) (PrjPoint_Krasovsky)))
((= 1 prjType) (mapcar 'set (list 'a 'f 'e2 'e12 'A1 'A2 'A3 'A4 'A5 'A6 'A7) (PrjPoint_IUGG1975)))
((= 2 prjType) (mapcar 'set (list 'a 'f 'e2 'e12 'A1 'A2 'A3 'A4 'A5 'A6 'A7) (PrjPoint_2000)))
(T (alert "未知prjType参数值")(exit))
)
(setq X (+ (* A1 B) (* A2 (sin (* 2 B))) (* A3 (sin (* 4 B))) (* A4 (sin (* 6 B))) (* A5 (sin (* 8 B))) (* A6 (sin (* 10 B))) (* A7 (sin (* 12 B))))
sinB (sin B) cosB (cos B) tanB (/ (sin B) (cos B)) tanB2 (* tanB tanB)
N (/ a (sqrt (- 1 (* e2 sinB sinB))))
m (* cosB (- L L0)) m2 (* m m)
ng2 (/ (* cosB cosB e2) (- 1 e2))
cady (+ X (* N tanB (* (+ 0.5 (* (+ (/ (+ (- 5 tanB2) (* 9 ng2) (* 4 ng2 ng2)) 24.0) (/ (* (+ (- 61 (* 58 tanB2)) (* tanB2 tanB2)) m2) 72.0)) m2)) m2)))
cadx (+ 500000 (* N m (+ 1 (* m2 (+ (/ (+ (- 1 tanB2) ng2) 6.0) (/ (* m2 (- (+ 5 (* tanB2 tanB2) (* 14 ng2)) (* 18 tanB2) (* 58 ng2 tanB2))) 120.0))))))
)
(list cadx cady)
)
;;;
;;;cadx,cady转BL
;;;输入参数
;;;;pt为点位cad点坐标,
;;;L0为中央经线，单位为弧度
;;;prjType取值为0:北京54坐标系,1:西安80坐标系，2：CGCS2000坐标系
;;;输出经纬度B,L表，单位为dd.mmss
(defun PrjPoint:xy2BL(pt L0 prjType / a f e2 e12 A1 A2 A3 A4 A5 A6 A7 sinB cosB tanB tanB2 N ng2 V yN preB0 B0 eta x y )
(cond
((= 0 prjType) (mapcar 'set (list 'a 'f 'e2 'e12 'A1 'A2 'A3 'A4 'A5 'A6 'A7) (PrjPoint_Krasovsky)))
((= 1 prjType) (mapcar 'set (list 'a 'f 'e2 'e12 'A1 'A2 'A3 'A4 'A5 'A6 'A7) (PrjPoint_IUGG1975)))
((= 2 prjType) (mapcar 'set (list 'a 'f 'e2 'e12 'A1 'A2 'A3 'A4 'A5 'A6 'A7) (PrjPoint_2000)))
(T (alert "未知prjType参数值")(exit))
)
(setq x (cadr pt)
y (- (car pt) 500000)
B0 (/ x A1)
eta 1
)
(while (> eta 0.000000000048)
(setq preB0 B0
B0 (/ (- x (+ (* A2 (sin (* 2 B0))) (* A3 (sin (* 4 B0))) (* A4 (sin (* 6 B0))) (* A5 (sin (* 8 B0))) (* A6 (sin (* 10 B0))) (* A7 (sin (* 12 B0))))) A1)
eta (abs (- B0 preB0))
)
)
(setq sinB (sin B0)
cosB (cos B0)
tanB (/ (sin B0) (cos B0))
tanB2 (* tanB tanB)
N (/ a (sqrt (- 1 (* e2 sinB sinB))))
ng2 (/ (* cosB cosB e2) (- 1 e2))
V (sqrt (+ 1 ng2))
yN (/ y N)
B (- B0 (* (+ (- (* yN yN) (/ (* (- (+ 5 (* 3 tanB2) ng2) (* 9 ng2 tanB2)) yN yN yN yN) 12.0)) (/ (* (+ 61 (* 90 tanB2) (* 45 tanB2 tanB2)) yN yN yN yN yN yN) 360.0)) V V tanB 0.5))
L (+ L0 (/ (+ (- yN (/ (* (+ 1 (* 2 tanB2) ng2) yN yN yN) 6.0 )) (/ (* (+ 5 (* 28 tanB2) (* 24 tanB2 tanB2) (* 6 ng2) (* 8 ng2 tanB2)) yN yN yN yN yN) 120.0)) cosB))
)
(list (Rad2Dms B) (Rad2Dms L))
)

;;;功能：投影换带 
;;作者:HUA 
;;;创建日期：2020-8-18 
;;;最后编辑日期:2020-8-20 
;;; 
;;;秒转度分秒 s->dd.mmss 
(defun Second2Dms(second / sign degree miniute second  dms) 
  (setq   sign (< second 0) 
      second (/ (abs second) 3600.0) 
      degree (fix second) 
      second (* (- second degree) 60.0) 
      miniute (fix second) 
      second (* (- second miniute) 60.0) 
      dms (+ degree (/ miniute 100.0) (/ second 10000.0)) 
  ) 
  (if sign (setq dms (- 0 dms))) 
  dms 
) 
;;; 
;;;角度转弧度 dd.mmss->rad 
(defun Dms2Rad(dms / sign degree miniute second rad) 
  (setq  sign  (<  dms 0) 
         dms  (abs dms) 
       degree (fix dms) 
       dms (* (- dms degree) 100) 
       miniute (fix dms) 
       second (* (- dms miniute) 100)    
       rad (/ (* (+ degree (/ miniute 60.0) (/ second 3600.0)) PI) 180.0) 
  ) 
  (if sign (setq rad  (- 0 rad))) 
  rad 
) 
;;; 
;;;弧度转角度 rad->dd.mmss 
(defun Rad2Dms(rad / sign degree miniute second dms) 
  (setq   sign (< rad 0) 
          rad (abs rad) 
      rad (/ (* rad 180.0) PI) 
      degree (fix rad) 
      rad (* (- rad degree) 60.0) 
      miniute (fix rad) 
      second (* (- rad miniute) 60.0) 
      dms (+ degree (/ miniute 100.0) (/ second 10000.0)) 
  ) 
  (if sign (setq dms (- 0 dms))) 
  dms 
) 
;;; 
;;;十进制角度转弧度 dd->rad 
(defun Degree2Rad(degree / sign rad) 
  (setq  sign  (<  degree 0) 
         degree   (abs degree )       
       rad (/ (* degree  PI) 180.0) 
  ) 
  (if sign (setq rad  (- 0 rad))) 
  rad 
) 
;;; 
;;;弧度转十进制角度 rad->dd 
(defun Rad2Degree(rad / sign degree) 
  (setq  sign  (<  rad 0) 
         rad (abs rad) 
         degree (/ (* rad 180.0) PI) 
  ) 
  (if sign (setq degree  (- 0 degree))) 
  degree 
) 
;;; 
;;;54坐标系参数 
(defun PrjPoint_Krasovsky( / a f e2 e12 A1 A2 A3 A4 A5 A6 A7 tf ) 
  (setq  a 6378245.0 
      f 298.3 
      tf (/ (- f 1) f) 
      e2  (- 1 (* tf tf)) 
      tf (/ f (- f 1)) 
      e12 (- (* tf tf) 1) 
      A1 6367558.4968749800000000 
      A2  -16036.4802694116000000 
      A3  16.8280668841393000 
      A4  -0.0219753089548841 
      A5 0.0000311310026601 
      A6 -0.0000000459827472 
      A7 0.0000000000679857       
  ) 
  (list a f e2 e12 A1 A2 A3 A4 A5 A6 A7)   
) 
;;; 
;;;80坐标系参数 
(defun PrjPoint_IUGG1975(/ a f e2 e12 A1 A2 A3 A4 A5 A6 A7 tf) 
  (setq  a    6378140.0 
      f    298.257 
      tf (/ (- f 1) f) 
      e2  (- 1 (* tf tf)) 
      tf (/ f (- f 1)) 
      e12 (- (* tf tf) 1) 
      A1    6367452.1327884300000000 
      A2    -16038.5282286966000000 
      A3    16.8326464353061000 
      A4    -0.0219844636486152 
      A5    0.0000311484690908 
      A6    -0.0000000460151868 
      A7    0.0000000000680433       
  ) 
  (list a f e2 e12 A1 A2 A3 A4 A5 A6 A7)   
) 
;;; 
;;;2000坐标系参数 
(defun PrjPoint_2000(/ a f e2 e12 A1 A2 A3 A4 A5 A6 A7 tf) 
  (setq  a    6378137.0 
      f    298.257222101 
      tf (/ (- f 1) f) 
      e2  (- 1 (* tf tf)) 
      tf (/ f (- f 1)) 
      e12 (- (* tf tf) 1) 
      A1    6367449.1457710400000000 
      A2    -16038.5087415914000000 
      A3    16.8326134276660000 
      A4    -0.0219844041401628 
      A5    0.0000311483615430 
      A6    -0.0000000460149936 
      A7    0.0000000000680429       
  ) 
   (list a f e2 e12 A1 A2 A3 A4 A5 A6 A7)   
) 
;;; 
;;;BL转cadx,cady 
;;;输入参数 
;;;B,L为点位纬度和经度，单位为弧度 
;;;L0为中央经线，单位为弧度 
;;;prjType取值为0:北京54坐标系,1:西安80坐标系，2：CGCS2000坐标系 
;;;输出CAD坐标x,y 
(defun PrjPoint:BL2xy(B L L0 prjType / a f e2 e12 A1 A2 A3 A4 A5 A6 A7 X N tanB tanB2 m m2 ng2 sinB cosB cadx cady)    
    (cond 
    ((= 0 prjType) (mapcar 'set (list 'a 'f 'e2 'e12 'A1 'A2 'A3 'A4 'A5 'A6 'A7) (PrjPoint_Krasovsky))) 
    ((= 1 prjType) (mapcar 'set (list 'a 'f 'e2 'e12 'A1 'A2 'A3 'A4 'A5 'A6 'A7) (PrjPoint_IUGG1975))) 
    ((= 2 prjType) (mapcar 'set (list 'a 'f 'e2 'e12 'A1 'A2 'A3 'A4 'A5 'A6 'A7) (PrjPoint_2000))) 
    (T (alert "未知prjType参数值")(exit)) 
  ) 
  (setq  X       (+ (* A1 B) (* A2 (sin (* 2 B))) (*  A3 (sin (* 4 B))) (* A4 (sin (* 6 B))) (* A5  (sin (* 8 B))) (* A6 (sin (* 10 B))) (* A7 (sin (* 12 B)))) 
      sinB   (sin B)  cosB   (cos B) tanB   (/ (sin B) (cos B))   tanB2   (* tanB tanB) 
      N      (/ a (sqrt (- 1 (* e2 sinB sinB)))) 
      m      (* cosB (- L L0))  m2 (* m m) 
      ng2      (/ (* cosB cosB e2) (- 1 e2)) 
      cady    (+ X (* N tanB (* (+ 0.5 (* (+ (/ (+ (- 5 tanB2) (* 9 ng2) (* 4 ng2 ng2)) 24.0) (/ (* (+ (- 61 (* 58 tanB2)) (* tanB2 tanB2)) m2) 72.0)) m2)) m2))) 
      cadx    (+ 500000 (* N m  (+ 1 (* m2 (+ (/ (+ (- 1 tanB2) ng2) 6.0) (/ (* m2 (-  (+ 5 (* tanB2 tanB2) (* 14 ng2)) (* 18 tanB2) (* 58 ng2 tanB2))) 120.0))))))    
  ) 
  (list cadx cady)   
) 
;;; 
;;;cadx,cady转BL 
;;;输入参数 
;;;;pt为点位cad点坐标, 
;;;L0为中央经线，单位为弧度 
;;;prjType取值为0:北京54坐标系,1:西安80坐标系，2：CGCS2000坐标系 
;;;输出经纬度B,L表，单位为dd.mmss 
(defun PrjPoint:xy2BL(pt L0 prjType / a f e2 e12 A1 A2 A3 A4 A5 A6 A7 sinB cosB tanB tanB2 N ng2 V yN preB0 B0 eta x y ) 
  (cond 
    ((= 0 prjType) (mapcar 'set (list 'a 'f 'e2 'e12 'A1 'A2 'A3 'A4 'A5 'A6 'A7) (PrjPoint_Krasovsky))) 
    ((= 1 prjType) (mapcar 'set (list 'a 'f 'e2 'e12 'A1 'A2 'A3 'A4 'A5 'A6 'A7) (PrjPoint_IUGG1975))) 
    ((= 2 prjType) (mapcar 'set (list 'a 'f 'e2 'e12 'A1 'A2 'A3 'A4 'A5 'A6 'A7) (PrjPoint_2000))) 
    (T (alert "未知prjType参数值")(exit)) 
  ) 
  (setq   x    (cadr pt) 
      y    (- (car pt) 500000) 
      B0   (/ x A1) 
      eta  1 
  )   
  (while (> eta  0.000000000048) 
    (setq  preB0  B0 
        B0   (/ (- x (+ (* A2 (sin (* 2 B0))) (* A3 (sin (* 4 B0))) (* A4 (sin (* 6 B0))) (* A5 (sin (* 8 B0))) (* A6 (sin (* 10 B0))) (* A7 (sin (* 12 B0))))) A1) 
        eta (abs (- B0 preB0)) 
    )   
  )   
  (setq   sinB  (sin B0) 
      cosB   (cos B0) 
      tanB    (/ (sin B0) (cos B0)) 
      tanB2   (* tanB tanB) 
      N      (/ a (sqrt (- 1 (* e2 sinB sinB)))) 
      ng2      (/ (* cosB cosB e2) (- 1 e2)) 
      V    (sqrt (+ 1 ng2)) 
      yN    (/ y N)    
      B    (- B0 (* (+ (- (* yN yN) (/ (* (- (+ 5 (* 3 tanB2) ng2) (* 9 ng2 tanB2)) yN yN yN yN) 12.0)) (/ (* (+ 61 (* 90 tanB2) (* 45 tanB2 tanB2)) yN yN yN yN yN yN) 360.0)) V V tanB 0.5))       
      L    (+ L0 (/ (+ (- yN (/ (* (+ 1 (* 2 tanB2) ng2) yN yN yN) 6.0 )) (/ (* (+ 5 (* 28 tanB2) (* 24 tanB2 tanB2) (* 6 ng2) (* 8 ng2 tanB2)) yN yN yN yN yN) 120.0)) cosB))    
  ) 
  (list (Rad2Dms B) (Rad2Dms L)) 
)

434939575 · 发表于 2020-8-21 10:02:07

这个做什么用的啊。

yshf · 发表于 2020-8-21 17:15:28

请问：能否给出A1～A7的计算公式？

skymudy · 发表于 2020-8-23 09:24:34

牛！！！学习，收藏了！

hothua2020 · 发表于 2020-8-24 08:19:27

yshf 发表于 2020-8-21 17:15
请问：能否给出A1～A7的计算公式？

; //A1-A7计算公式
; //e4 = pow(e2, 2); // e2 * e2; e6 = pow(e2, 3); //e4 * e2; e8 = pow(e2, 4); //e6 * e2;e10 = pow(e2, 5); //e8 * e2;e_12 = pow(e2, 6); //e10 * e2;
; //A1 = a*(1-e2)*(1 + 3 * e2 / 4 + 45 * e4 / 64 + 175 * e6 / 256 + 11025 * e8 / 16384 + 43659 * e10 / 65536 + 693693 * e_12 / 1048576);
; //A2 = a*(1-e2)*(3 * e2 / 8 + 15 * e4 / 32 + 525 * e6 / 1024 + 2205 * e8 / 4096 + 72765 * e10 / 131072 + 297297 * e_12 / 524288);
; //A3 = a*(1-e2)*(15 * e4 / 256 + 105 * e6 / 1024 + 2205 * e8 / 16384 + 10395 * e10 / 65536 + 1486485 * e_12 / 8388608);
; //A4 = a*(1-e2)*(35 * e6 / 3072 + 105 * e8 / 4096 + 10395 * e10 / 262144 + 55055 * e_12 / 1048576);
; //A5 = a*(1-e2)*(315 * e8 / 131072 + 3465 * e10 / 524288 + 99099 * e_12 / 8388608);
; //A6 = a*(1-e2)*(693 * e10 / 1310720 + 9009 * e_12 / 5242880);
; //A7 = a*(1-e2)*(1001 * e_12 / 8388608);

yshf · 发表于 2020-8-24 14:32:38

谢谢hothua2020提供的计算式

wzg356 · 发表于 2020-8-25 18:25:43

这个清爽.................

yuanziyou · 发表于 2020-8-29 11:44:54

哈哈哈，这个世界太小了，我认识楼主

树櫴希德 · 发表于 2021-6-21 09:39:28

用笑脸软件COORD核算过果然不差

Siuzf · 发表于 2022-5-26 16:45:56

没搞懂怎么用

账号		自动登录	找回密码
密码			注册

高斯投影正反算

评分