社区 › AutoCAD与GIS测绘 › 利用大神函数经纬度展点，多段线到十进制经纬度

树櫴希德 发表于 2021-6-21 18:18:57

利用大神函数经纬度展点，多段线到十进制经纬度

利用大神函数经纬度展点，多段线到十进制经纬度
谷歌地图KML文件经纬度是十进制的 奥维也是
以下是KML文件样表红色子可以改动
<?xml version="1.0" encoding="UTF-8"?>
<kml xmlns="http://www.opengis.net/kml/2.2" xmlns:gx="http://www.google.com/kml/ext/2.2" xmlns:kml="http://www.opengis.net/kml/2.2" xmlns:atom="http://www.w3.org/2005/Atom">
<Document>
        <name>CCCC巴中停车场.kml</name>
        <Style id="s_ylw-pushpin_hl21">
                <IconStyle>
                        <scale>1.3</scale>
                        <Icon>
                                <href>http://maps.google.com/mapfiles/kml/pushpin/ylw-pushpin.png</href>
                        </Icon>
                        <hotSpot x="20" y="2" xunits="pixels" yunits="pixels"/>
                </IconStyle>
                <LineStyle>
                        <color>ff0000ff</color>
                        <width>3</width>
                </LineStyle>
        </Style>
        <Style id="s_ylw-pushpin11">
                <IconStyle>
                        <scale>1.1</scale>
                        <Icon>
                                <href>http://maps.google.com/mapfiles/kml/pushpin/ylw-pushpin.png</href>
                        </Icon>
                        <hotSpot x="20" y="2" xunits="pixels" yunits="pixels"/>
                </IconStyle>
                <LineStyle>
                        <color>ff0000ff</color>
                        <width>3</width>
                </LineStyle>
        </Style>
        <StyleMap id="m_ylw-pushpin21">
                <Pair>
                        <key>normal</key>
                        <styleUrl>#s_ylw-pushpin11</styleUrl>
                </Pair>
                <Pair>
                        <key>highlight</key>
                        <styleUrl>#s_ylw-pushpin_hl21</styleUrl>
                </Pair>
        </StyleMap>
        <Placemark>
                <name>未命CCCCBB名路径</name>
                <styleUrl>#m_ylw-pushpin21</styleUrl>
                <LineString>
                        <tessellate>1</tessellate>
                        <coordinates>
                                113.235146862,22.957426436
113.235051235,22.957391803
113.235063160,22.957363577
113.235074058,22.957367525
113.235084711,22.957342288
113.235071549,22.957337522
                </coordinates>
                </LineString>
        </Placemark>
</Document>
</kml>


;;;功能：投影换带
;;作者: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
)
;;;十进制角度转DD.MMSS

(DEFUN 10-DDMMSS ( A / B C D E )
(SETQ B (Rad2Dms (Degree2Rad A)) )
(setq c (fix b)) (setq d (fix(*(- b c) 100)))
(setq e (- (* 10000 b) (* c 10000)(* d 100) ) )
(strcat (rtos c 2 0) "度" (rtos d 2 0) "分"(rtos e 2 6) "秒")
)

;;;弧度转十进制角度 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))
)
;;;
;;;cadx,cady转BL
;;;输入参数
;;;;pt为点位cad点坐标,
;;;L0为中央经线，单位为弧度
;;;prjType取值为0:北京54坐标系,1:西安80坐标系，2：CGCS2000坐标系
;;;输出经纬度B,L表，单位为十进制
(defun PrjPoint:xy2BL2(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 (Rad2Degree L) (Rad2Degree B) )
)

;;;;;;;;;;;;;;;;;;;;;;;;;;;


(defun c:jw-2000 (   / l0 bb ll vvv str)

(setq l0 (getreal "\n请输入中央子午线经度dd.mmss:"))
(setq bb (getreal "\n请输入待求点纬度dd.mmss:"))
(setq ll (getreal "\n请输入待求点经度dd.mmss:"))

(setq vvv(append(PrjPoint:BL2xy (dms2rad bb) (dms2rad ll) (dms2rad l0)2)'(0)))
(setq str(strcat (rtos bb 2 12 )","(rtos ll 2 12 )))
(entmake (list '(0 . "TEXT")'(8 . "jwd") (cons 1 str) (cons 10 vvv) (cons 40 2)))
(entmake (list '(0 . "POINT")'(8 . "jwd") (cons 10 vvv)))

)

(defun vxs(e / p a b n ob q et d d1 en et)
    (setq a(entget e)ob(vlax-ename->vla-object e)et(cdr(assoc 0 a))n 0 p nil d nil)
    (cond((="LWPOLYLINE"et)
    (repeat(length a)(setq b (nth n a) n (+ n 1))
      (if (= 10 (car b))(progn
      (setq q(list (cadr b) (caddr b))d1(vlax-curve-getDistAtPoint ob q))
      (if p (if (not(member d1 d)) (setq p (append p (list q))d (append d (list d1))))
          (setq p (list q)))))))
   ((="POLYLINE"et)
    (SETQ EN (ENTGET (SETQ E (ENTNEXT E))))
    (WHILE (/= (CDR (ASSOC 0 EN)) "SEQEND")
      (SETQ q (CDR (ASSOC 10 EN))d1(vlax-curve-getDistAtPoint ob q)q(reverse(cdr(reverse q))))
      (if p(if (not(member d1 d)) (setq p (append p (list q))d (append d (list d1))))(setq p(list q)))
      (SETQ EN (ENTGET (SETQ E (ENTNEXT E)))))
    (setq p(reverse p))))P)

(defun c:lw-bl2000 ( / p l0 zb jwd );;;2000坐标多段线TXT转经纬度文件，经纬度为10进制

(vl-load-com)
(setq l0 (getreal "\n请输入中央子午线经度dd.mmss:"))
(setq p (car(entsel "\n请选择2000或者wgs84坐标多段线：")))

(setq zb (vxs p))
(setq ffn (getfiled "选取/建立数据导出10进制经纬度文件" "" "txt" 1))
(setq ff (open ffn "w"))
(foreach x zb
(setq jwd (PrjPoint:xy2BL2 x (Dms2Rad l0) 2))

(write-line (strcat (rtos (car jwd)2 9)","(rtos (cadr jwd)2 9)) ff)

)
(close ff)
       )

;;;;;;;;;;;;;;;;;;;;;;;;;
(defun c:bl-2000 (/l0 fn f ll k lll bbb lllbbb str);;;经纬度文件TXT转2000坐标，图上经纬度为度分秒
(vl-load-com)
(setq l0 (getreal "\n请输入中央子午线经度dd.mmss:"))
(setq fn (getfiled "Select Log file请打开10进制经纬度数据文件：格式为：经度，纬度" "" "txt" 8))

(setq f (open (findfile fn) "r"))

(while (setq k (read-line f))
    (setq ll (cons kll ) )
)
(close f)
;(String:Split (car ll) ",")
(foreach x ll
(setq lll (atof (car(String:Split x ","))))
(setq bbb (atof (cadr(String:Split x ","))))

(setq lllbbb( PrjPoint:BL2xy(Degree2Rad bbb) (Degree2Rad lll) (dms2rad l0) 2))
(setq str(strcat (10-DDMMSS bbb )","(10-DDMMSS lll )))
(entmake (list '(0 . "TEXT")'(8 . "jwd") (cons 1 str) (cons 10 lllbbb) (cons 40 2)))
(entmake (list '(0 . "POINT")'(8 . "jwd") (cons 10 lllbbb)))
)
(princ)
)

(defun String:Split (str delimiter / post strlst stl)
    (if      str
      (progn
      (setq stl (strlen delimiter))
      (while (vl-string-search delimiter str)
          (setq      post   (vl-string-search delimiter str)
                strlst (cons (substr str 1 post) strlst)
                str    (substr str (+ 1 post stl))
          )
      )
      (reverse (vl-remove "" (cons str strlst)))
      )
    )
)

烟盒迷唇 发表于 2021-6-22 08:40:33

谢谢分享啊

664571221 发表于 2021-6-22 10:29:08

兄弟你qq多少可以加你一下吗

树櫴希德 发表于 2021-6-22 14:47:47

根据4参数算坐标？


;;;角度转弧度 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
)
;;;
(setq pp'(40 50))

(setq aa(Dms2Rad(getreal "\n请输入旋转角度dd.mmss:")))
(setq kk(getreal "\n请输入尺度1.00000:"))

;;;;;;;;;;;;;;;;;;;
(defun zbzh ( pp aa kk /   xx yy x0 y0 )


(setq xx (cadr pp) yy(car pp) )

(setq x0   (* (cos aa)(- (* xx kk)    (* (/ (sin aa) (cos aa)) (* kkyy) )            ))   )

;

(setq y0 (+(/ (* kkyy) (cos aa))   (* (/ (sin aa) (cos aa)) x0 )      ))

(list y0 x0)

)

树櫴希德 发表于 2021-6-22 17:16:19

地方坐标系或者自定义坐标根据4参数多段线转经纬度文件TXT，经纬度为10进制，可以复制到KML文件用\n""可以将手薄里面参与求控制点的WGS84经纬度换算成坐标\n" "然后用COORD笑脸软件求地方坐标--》WGS84坐标的4参数抄下来用113.47228431,22.986922997
113.472186361,22.986816298
113.472054915,22.986774606
113.4719921,22.986754673
113.471037673,22.984983334
113.471092087,22.984715313
113.470541912,22.985446895
113.470968935,22.985462633
113.471749262,22.986910865
113.471755997,22.986923274
113.471762811,22.986935665
113.471769664,22.98694801
113.471776604,22.986960328
113.471783612,22.98697261
113.471790689,22.986984865
113.471797844,22.986997093
113.471805047,22.987009266
113.471812319,22.987021422
113.471819678,22.987033532
113.471827087,22.987045624
113.471834553,22.987057671
113.471842108,22.987069681
113.47184972,22.987081638
113.471857401,22.987093585
113.471865131,22.987105478
113.471872949,22.987117353
113.471880815,22.987129182
113.471888769,22.987140957
113.471896772,22.987152715
113.471904843,22.987164427
113.471912982,22.987176102
113.47192117,22.987187742
113.471929445,22.987199327
113.471937769,22.987210886
113.471946172,22.987222408
113.471954613,22.987233894
113.471963142,22.987245325
113.471971719,22.987256721


;;;功能：投影换带
;;作者: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
)
;;;十进制角度转DD.MMSS

(DEFUN 10-DDMMSS ( A / B C D E )
(SETQ B (Rad2Dms (Degree2Rad A)) )
(setq c (fix b)) (setq d (fix(*(- b c) 100)))
(setq e (- (* 10000 b) (* c 10000)(* d 100) ) )
(strcat (rtos c 2 0) "度" (rtos d 2 0) "分"(rtos e 2 6) "秒")
)

;;;弧度转十进制角度 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))
)
;;;
;;;cadx,cady转BL
;;;输入参数
;;;;pt为点位cad点坐标,
;;;L0为中央经线，单位为弧度
;;;prjType取值为0:北京54坐标系,1:西安80坐标系，2：CGCS2000坐标系
;;;输出经纬度B,L表，单位为十进制
(defun PrjPoint:xy2BL2(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 (Rad2Degree L) (Rad2Degree B) )
)

;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;

(defun vxs(e / p a b n ob q et d d1 en et)
    (setq a(entget e)ob(vlax-ename->vla-object e)et(cdr(assoc 0 a))n 0 p nil d nil)
    (cond((="LWPOLYLINE"et)
    (repeat(length a)(setq b (nth n a) n (+ n 1))
      (if (= 10 (car b))(progn
      (setq q(list (cadr b) (caddr b))d1(vlax-curve-getDistAtPoint ob q))
      (if p (if (not(member d1 d)) (setq p (append p (list q))d (append d (list d1))))
          (setq p (list q)))))))
   ((="POLYLINE"et)
    (SETQ EN (ENTGET (SETQ E (ENTNEXT E))))
    (WHILE (/= (CDR (ASSOC 0 EN)) "SEQEND")
      (SETQ q (CDR (ASSOC 10 EN))d1(vlax-curve-getDistAtPoint ob q)q(reverse(cdr(reverse q))))
      (if p(if (not(member d1 d)) (setq p (append p (list q))d (append d (list d1))))(setq p(list q)))
      (SETQ EN (ENTGET (SETQ E (ENTNEXT E)))))
    (setq p(reverse p))))P)


;;;;;;;;;;;;;;;;;;;
(defun zbzh ( pp aa kk /   xx yy x0 y0 )


(setq xx (cadr pp) yy(car pp) )

(setq x0   (* (cos aa)(- (* xx kk)    (* (/ (sin aa) (cos aa)) (* kkyy) )            ))   )

;

(setq y0 (+(/ (* kkyy) (cos aa))   (* (/ (sin aa) (cos aa)) x0 )      ))

(list (+ y0 (car syd))(+ x0 (cadr syd)))

)
;3、点表生成多段线
(defun makepl (lst / pt)
(entmake (append    (list '(0 . "LWPOLYLINE") '(100 . "AcDbEntity") '(100 . "AcDbPolyline") '(8 . "ctpzx") (cons 90 (length lst)) (cons 70 129))
      (mapcar '(lambda (pt)(cons 10 pt)) lst ))
) )
;;;;;;

(defun c:dflw-bl2000 ( / p l0 zb jwd syd aa kk xzblll);;;地方坐标系或者自定义坐标根据4参数多段线转经纬度文件TXT，经纬度为10进制，可以复制到KML文件用
(alert (strcat "地方坐标系或者自定义坐标根据4参数多段线转经纬度文件TXT，经纬度为10进制，可以复制到KML文件用\n""可以将手薄里面参与求控制点的WGS84经纬度换算成坐标\n" "然后用COORD笑脸软件求地方坐标--》WGS84坐标的4参数抄下来用"))
(vl-load-com)
(setq l0 (getreal "\n请输入中央子午线经度dd.mmss:"))
(setq p (car(entsel "\n请选择地方坐标系或者自定义坐标多段线：")))

;(setq pp '(62295.5790 17258.584))

(setq syd'(454891.370478 2529629.687669 )) ;4参数x y位移人工输入

(setq aa (Dms2Rad 0.0047650491))   ;4参数旋转参数人工输入dd.mmss
(setq kk 1.00000515799711)      ;4参数缩放比列因子人工输入
    (setq zb (vxs p))
(setq lll '())
(foreach y zb
    (setq lll (zbzh y aa kk))
    (setq xzb (cons lll xzb))
)
(makepl xzb)
(setq ffn (getfiled "选取/建立数据导出10进制经纬度文件" "" "txt" 1))
(setq ff (open ffn "w"))
(foreach x xzb
(setq jwd (PrjPoint:xy2BL2 x (Dms2Rad l0) 2))

(write-line (strcat (rtos (car jwd)2 9)","(rtos (cadr jwd)2 9)) ff)

)
(close ff)
       )

流氓兔 发表于 2021-6-22 21:32:33

牛逼的哥们，每次都是好代码

yshf 发表于 2021-6-22 23:49:16

谢谢分享

w379106181 发表于 2021-8-24 08:07:31

谢谢分享啊
谢谢分享啊

xujinhua 发表于 2024-3-31 20:58:20

厉害，正需要转经纬度

查看完整版本: 利用大神函数经纬度展点，多段线到十进制经纬度