[ArcPy] kml转shp文件 | Arcpy与shp文件 | ArcGIS脚本工具


涉及内容:ArcPy操作shape文件;python2解析xml文档

实验课中需要kml文件转shp文件
ArcGIS Desktop工具箱中有 “kml转图层” 与 “要素类转shp”两个工具,我们可以利用这两个工具制作一个模型,这个比较简单
因为小生最近接触ArcPy,所以想用Python在文件转换的角度进行转换

相关文件下载: http://download.csdn.net/download/summer_dew/9992151
若要运行此示例需要在Python2安装BeautifulSoup模块:https://blog.csdn.net/summer_dew/article/details/78045134

ArcPy操作shp文件

# 输出地址
outPath = r"F:\study\school\analyze\lab1\workspace\polygon.shp"
# 获取shape文件保存的地址
outWorkspace = os.path.split(outPath)[0]
# shape文件保存的名字
shpName = os.path.split(outPath)[-1]
# 坐标系
spat_ref = "4326"
# 创建shape文件
gp.CreateFeatureClass_management(outWorkspace,shpName,featureType,"","","",spat_ref)
# 创建字段
gp.addfield(outPath, "NAME", "text", "", "", "100")
gp.addfield(outPath, "DESCR", "text", "", "", "250")
gp.addfield(outPath, "FOLDER", "text", "", "", "100")
# 获取插入游标
cur = gp.InsertCursor(outPath)
# 新建一行
newRow = cur.newRow()
# 遍历数据
for feature in kmlDataList:
    name = feature[0] #名称字段
    desc = feature[1] #描述字段
    XYList = feature[2] #XY坐标
    folderName = feature[3] #文件夹名称字段

    if len(XYList)==0: #如果没有XY坐标
        print "no find coordinate -- feature's name :%s  " % name
        continue

    # 创建点
    pnt = gp.CreateObject("point")
    # 对shp文件中要素赋予空间坐标
    if featureType=="POLYGON" or featureType == "POLYLINE":
        array = gp.CreateObject("array")
        for XY in XYList:
            pnt.X,pnt.Y = XY[0],XY[1]
            array.add(pnt)
        newRow.Shape = array
    else:
        pnt.X,pnt.Y = XYList[0][0],XYList[0][1]
        newRow.Shape = pnt
    # 对该行的字段赋值
    newRow.NAME = name
    newRow.DESCR = desc
    newRow.FOLDER = folderName
    #插入数据,newRow自动下移
    cur.InsertRow(newRow)

完整代码

若该脚本要创建为工具箱,内容一定不能出现中文!连注释都不行!(ArcGIS10.2)
其他环境下可以出现中文

# -*- coding:utf-8 -*-
# Author: PasserQi
# time: 2017/9/23
# Dependence: BeautifulSoup
# function: kml to shp

import arcpy,arcgisscripting,os,sys
from bs4 import BeautifulSoup

# create gp
gp = arcgisscripting.create()
gp.OverwriteOutput = 1

# -------------------------------------------------------------------
# init param

gp.AddMessage("get init param....")
print "get init param"

# get param
kmlFilePath = sys.argv[1]
featureType = sys.argv[2]
outPath = sys.argv[3]
# test param
# kmlFilePath = r"F:\study\school\analyze\lab1\songbai.kml"
# featureType = "POLYGON"
# outPath = r"F:\study\school\analyze\lab1\workspace\polygon.shp"

# UI
gp.AddMessage("kml file path %s \n feature type %s \n out path %s" % (kmlFilePath,featureType,outPath) )
print "kml file path %s \n feature type %s \n out path %s" % (kmlFilePath,featureType,outPath)

# ------------------------------------------------
# read kml and analysis

openKmlFile = open(kmlFilePath,"r")
kmlDom = openKmlFile.read()
bsObj = BeautifulSoup(kmlDom,"html.parser")
openKmlFile.close()

# ---------------------------------------------------
# get kml data

kmlDataList = [] # list kml feature data
# search feature type
if featureType=="POLYGON": type = "polygon"
elif featureType=="POLYLINE" : type = "linestring"
elif featureType=="POINT" : type = "point"

folders = bsObj.findAll("folder")
for folder in folders:
    folderName = folder.find("name").get_text() # folder name
    print folderName
    features = bsObj.findAll(type) #Search for specified type
    for feature in features:
        name = description = ""
        XYList = []

        coordinateStr = feature.find("coordinates").get_text() #get coordinates
        coordinateStr = coordinateStr.replace('\t','').replace('\n','') #remove \n \t
        coordinateList = coordinateStr.split(' ')
        if coordinateList:
            for XYZstr in coordinateList:
                if XYZstr=='': #coordinateList:the last one is ''
                    continue
                XYZ = XYZstr.split(',')
                x = float(XYZ[0])
                y = float(XYZ[1])
                XYList.append([x,y])  #coordinate string to float
        else:
            XYList = []

        if (feature.parent.find("name")):
            name = feature.parent.find("name").get_text()
            print "name %s " % name
        if (feature.parent.find("description")):
            description = feature.parent.find("description").get_text()
            print "description %s" %description

        folderPath = '/'.join(folderName)
        kmlDataList.append([name,description,XYList,folderPath])


# -----------------------------------------------------------------------
# create shapefile
if kmlDataList:

    # create shp
    outWorkspace = os.path.split(outPath)[0]
    shpName = os.path.split(outPath)[-1]
    spat_ref = "4326"
    gp.CreateFeatureClass_management(outWorkspace,shpName,featureType,"","","",spat_ref)
    # create field
    gp.addfield(outPath, "NAME", "text", "", "", "100")
    gp.addfield(outPath, "DESCR", "text", "", "", "250")
    gp.addfield(outPath, "FOLDER", "text", "", "", "100")
    # get cursor
    cur = gp.InsertCursor(outPath)
    newRow = cur.newRow()
    for feature in kmlDataList:
        name = feature[0]
        desc = feature[1]
        XYList = feature[2]
        folderName = feature[3]

        if len(XYList)==0:
            print "no find coordinate -- feature's name :%s  " % name
            continue

        pnt = gp.CreateObject("point")
        if featureType=="POLYGON" or featureType == "POLYLINE":
            array = gp.CreateObject("array")
            for XY in XYList:
                pnt.X,pnt.Y = XY[0],XY[1]
                array.add(pnt)
            newRow.Shape = array
        else:
            pnt.X,pnt.Y = XYList[0][0],XYList[0][1]
            newRow.Shape = pnt

        newRow.NAME = name
        newRow.DESCR = desc
        newRow.FOLDER = folderName

        cur.InsertRow(newRow)
    print "finish"
    del cur,newRow

详细 为脚本创建工具箱

  1. 查看.kml文件
    在谷歌地球上随意画两个面,然后导出成KML格式
    这里写图片描述
    使用记事本查看
    这里写图片描述
    结论:其实kml文件是自定义的xml文件,观察到我们所需要的信息在<Folder>标签中,查看数据特点,便于读取

  2. 查看ArcPy如何操作shp文件
    ArcPy帮助文档查找shape文件操作,找了好久。原来不在ArcPy下
    结论:创建shape文件,详细定义请看
    http://resources.arcgis.com/zh-cn/help/main/10.2/index.html#/na/00170000002p000000/

    arcpy.CreateFeatureclass_management (输出的位置, 输出的文件名, {要素类型}, {要素类属性方案}, {素类是否包含线性测量值(m 值)}, {要素类是否包含高程值}, {空间参考方案}, ...)
    # 例子
    arcpy.CreateFeatureclass_management("C:/output", "habitatareas.shp", "POLYGON", "study_quads.shp", "DISABLED", "DISABLED", "C:/workspace/landuse.shp")
    
  3. 代码
    因kml为自定义的xml文件,需要解析它,这里使用BeautifulSoup包,ArcPy下(Python2版本)安装BeautifulSoup
    需要对文件进行写入,使用到ArcPy 游标 与 插入游标:http://resources.arcgis.com/zh-cn/help/main/10.2/index.html#/na/03q300000044000000/

  4. 创建新的工具箱
    在ArcCatalog新建一个工具箱–>添加一个脚本
    这里写图片描述
    脚本内容
    这里写图片描述
    添加脚本路径
    这里写图片描述
    设置三个参数
    这里写图片描述
    这里写图片描述
    这里写图片描述

  5. 测试
    这里写图片描述
  6. 结果
    这里写图片描述

注意

  1. 工具运行异常
    Python脚本如果要添加成工具,脚本内一定不能出现中文!会报异常错误!
  2. ValueError: could not convert string to float

    • 问题一
      kml文件标签<coordinates></coordinates>中存放坐标集,坐标集内一个坐标的格式是:X,Y,Z(Z后面有一个空格字符)
      所以split字符串时会多解析一个坐标,是一个空的字符串
      所以会出现以上string转float的报错
      解决方法:

      coordinateStr = feature.find("coordinates").get_text() #得到坐标集
      coordinateStr = coordinateStr.replace('\t','').replace('\n','') #移除\n与\t
      coordinateList = coordinateStr.split(' ') #以空格分隔字符串,提取坐标
      if coordinateList: #坐标不为空
          for XYZstr in coordinateList: #遍历坐标集
              if XYZstr=='': #遇到""退出该次循环
                  continue
              XYZ = XYZstr.split(',')
              x = float(XYZ[0])
              y = float(XYZ[1])
              XYList.append([x,y])  #coordinate string to float
      else:
          XYList = []
      
    • 问题2
      kml文件中<coordinates></coordinates>中值会有\n与\t不必要字符
      解决方法:

      coordinateStr = coordinateStr.replace('\t','').replace('\n','') #移除\n与\t
      

转载自:https://blog.csdn.net/summer_dew/article/details/78072870

You may also like...