KML文件读取多个点，分析密度

#geoogel earth 中点位置存储在一个文件夹，直接读取.kml文件
# _*_ coding: cp936 _*_
import re,arcpy,os,codecs
from arcpy import env
txtfile=''
workspace=os.path.dirname(txtfile)
env.workspace=workspace
tf=os.path.basename(txtfile)
reader=codecs.open(txtfile,'r','utf-8')
inputname=tf.split('.')[0]+'.shp'
prjFile='WGS 1984.prj'
patA=re.compile('<name>.*</name>')
patB=re.compile('<coordinates>.*</coordinates>')
ncdic={}
while True:
    line=reader.readline()
    if len(line)==0:
        break
    line=line.strip()
    if line=='<Placemark>':
        while True:
            line=reader.readline()
            line=line.strip()
            if line=='</Placemark>':
                break
            m=patA.match(line)
            if m:
                names=re.sub('<name>(.*?)</name>',r'\1',line)
            n=patB.match(line)
            if n:
                coordi=re.sub('<coordinates>(.*?)</coordinates>',r'\1',line)
        ncdic[names]=coordi
reader.close()
arcpy.CreateFeatureclass_management(workspace,inputname,'Point','','','',prjFile)
fieldName1='Name'
arcpy.AddField_management(inputname,fieldName1,'TEXT')
keys=ncdic.keys()
newdic={}
for key in keys:
    newdic[key]=[float(a) for a in ncdic[key].split(',')]
cursor=arcpy.da.InsertCursor(inputname,['SHAPE@','Name'])
point=arcpy.Point()
for key in keys:
    point.X=newdic[key][0]
    point.Y=newdic[key][1]
    cursor.insertRow([point,key])
del cursor
    


#定义函数
txtfile=arcpy.GetParameterAsText(0)
prjFile=arcpy.GetParameterAsText(1)
def KMLLoad(txtfile,prjFile):
    import re,arcpy,os,codecs
from arcpy import env

workspace=os.path.dirname(txtfile)
env.workspace=workspace
tf=os.path.basename(txtfile)
reader=codecs.open(txtfile,'r','utf-8')
inputname=tf.split('.')[0]+'.shp'

patA=re.compile('<name>.*</name>')
patB=re.compile('<coordinates>.*</coordinates>')
ncdic={}
while True:
    line=reader.readline()
    if len(line)==0:
        break
    line=line.strip()
    if line=='<Placemark>':
        while True:
            line=reader.readline()
            line=line.strip()
            if line=='</Placemark>':
                break
            m=patA.match(line)
            if m:
                names=re.sub('<name>(.*?)</name>',r'\1',line)
            n=patB.match(line)
            if n:
                coordi=re.sub('<coordinates>(.*?)</coordinates>',r'\1',line)
        ncdic[names]=coordi
reader.close()
arcpy.CreateFeatureclass_management(workspace,inputname,'Point','','','',prjFile)
fieldName1='Name'
arcpy.AddField_management(inputname,fieldName1,'TEXT')
keys=ncdic.keys()
newdic={}
for key in keys:
    newdic[key]=[float(a) for a in ncdic[key].split(',')]
cursor=arcpy.da.InsertCursor(inputname,['SHAPE@','Name'])
point=arcpy.Point()
for key in keys:
    point.X=newdic[key][0]
    point.Y=newdic[key][1]
    cursor.insertRow([point,key])
del cursor
return inputname,workspace
pg,nd=KMLLoad(txtfile,prjFile)
arcpy.SetParameter(2,nd+'\\'+pg)

#计算点密度
txtfile=''#需要处理的,shp文件路径
prj=''#投影文件，可以是具有投影信息的文件
def prjdensity(txtfile,prj,pg):#转化投影，计算点密度
    import arcpy,os
    from arcpy import env
    from arcpy.sa import *
    workspace=os.path.dirname(txtfile)
    env.workspace=workspace
    prjfc=workspace+'\\'+pg+'_prj'+'.shp'#具有投影信息的文件名
    arcpy.Project_management(pg,prjfc,prj)
    densityR=PointDensity(prjfc,'NONE',5,NbrCircle(3500,'MAP'))
    densityR.save(workspace+'\\'+'densityR')
    return densityR