cad（2000坐标系）转kml

拼凑来的，找不到原文章了

import win32com.client
import pythoncom
import simplekml
from tkinter import messagebox,filedialog
import math

def XY2LatLon(X, Y, L0):

  iPI = 0.0174532925199433
  a = 6378137.0
  f= 0.00335281006247
  ZoneWide = 3 #按3度带进行投影

  ProjNo = int(X / 1000000)
  L0 = L0 * iPI
  X0 = ProjNo * 1000000 + 500000
  Y0 = 0
  xval = X - X0
  yval = Y - Y0

  e2 = 2 * f - f * f #第一偏心率平方
  e1 = (1.0 - math.sqrt(1 - e2)) / (1.0 + math.sqrt(1 - e2))
  ee = e2 / (1 - e2) #第二偏心率平方

  M = yval
  u = M / (a * (1 - e2 / 4 - 3 * e2 * e2 / 64 - 5 * e2 * e2 * e2 / 256))

  fai = u + (3 * e1 / 2 - 27 * e1 * e1 * e1 / 32) * math.sin(2 * u) + (21 * e1 * e1 / 16 - 55 * e1 * e1 * e1 * e1 / 32) * math.sin(4 * u) + (151 * e1 * e1 * e1 / 96) * math.sin(6 * u)+ (1097 * e1 * e1 * e1 * e1 / 512) * math.sin(8 * u)
  C = ee * math.cos(fai) * math.cos(fai)
  T = math.tan(fai) * math.tan(fai)
  NN = a / math.sqrt(1.0 - e2 * math.sin(fai) * math.sin(fai))
  R = a * (1 - e2) / math.sqrt(
    (1 - e2 * math.sin(fai) * math.sin(fai)) * (1 - e2 * math.sin(fai) * math.sin(fai)) * (1 - e2 * math.sin(fai) * math.sin(fai)))
  D = xval / NN

  #计算经纬度（弧度单位的经纬度）
  longitude1 = L0 + (D - (1 + 2 * T + C) * D * D * D / 6 + (
  5 - 2 * C + 28 * T - 3 * C * C + 8 * ee + 24 * T * T) * D * D * D * D * D / 120) / math.cos(fai)
  latitude1 = fai - (NN * math.tan(fai) / R) * (
  D * D / 2 - (5 + 3 * T + 10 * C - 4 * C * C - 9 * ee) * D * D * D * D / 24 + (
  61 + 90 * T + 298 * C + 45 * T * T - 256 * ee - 3 * C * C) * D * D * D * D * D * D / 720)

  #换换为deg
  longitude = longitude1 / iPI
  latitude = latitude1 / iPI
  return  longitude,latitude
acad = win32com.client.Dispatch("AutoCAD.Application")
doc = acad.ActiveDocument
doc.Utility.Prompt("\n醉后不知天在水\n满船清梦压星河\n")
mp = doc.ModelSpace  # 模型空间
messagebox.showinfo('提示', "请在屏幕拾取图元，以Enter键结束")
try:
    doc.SelectionSets.Item("SS1").Delete()
except:
    messagebox.showinfo('警告', "Delete selection failed")
slt = doc.SelectionSets.Add("SS1")
slt.SelectOnScreen()
kml = simplekml.Kml()

for x in slt:
    if x.ObjectName == "AcDbPolyline":
        pts = x.Coordinates
        points = [(XY2LatLon(pts[x],pts[x+1],120)) for x in range(0,len(pts),2)]
        pol = kml.newpolygon(name="czh", outerboundaryis=points)
        pol.style.polystyle.color = '550000FF' #aabbggrr
        pol.style.polystyle.outline = 1
        pnt = kml.newpoint(name="czh", coords=[points[0]])
        pnt.description = "HPPs in Nepal"
        pnt.style.labelstyle.color = 'FFFFFF'  # html color code
savepath = filedialog.asksaveasfilename(filetypes=[("kml", "kml")], title="保存kml")+".kml"
kml.save(savepath)