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 seconddms)
(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)
(setqsign(<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)
(setqsign(<degree 0)
         degree   (abs degree )      
       rad (/ (* degreePI) 180.0)
)
(if sign (setq rad(- 0 rad)))
rad
)
;;;
;;;弧度转十进制角度 rad->dd
(defun Rad2Degree(rad / sign degree)
(setqsign(<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 )
(setqa 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
      A316.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)
(setqa    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)
(setqa    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))
)
(setqX       (+ (* 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)
      eta1
)
(while (> eta0.000000000048)
    (setqpreB0B0
      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

没搞懂怎么用
页: [1]
查看完整版本: 高斯投影正反算