获取行政区边界和利用GDAL生成shp文件

# coding=utf-8
try:
    from osgeo import gdal
    from osgeo import ogr
    import urllib2 as ulb
    import json
except ImportError:
    import gdal
    import ogr


# 调用高德地图API获取边界信息
def GetBoundary(keyword):
    # API的KEY
    key = "bb4c4785c81b07d286fc561879c6037b"

    # 拼接请求url
    url = "https://restapi.amap.com/v3/config/district?key=" + key + "&keywords=" + keyword + "&subdistrict=0&extensions=all"

    # 接收返回的结果
    response = ulb.urlopen(url)
    response = response.read()
    s = json.loads(response)
    polyline = s['districts'][0]['polyline']

    # 第一次拆分，将多段线拆分开
    first_separate = polyline.split("|")
    polylines = []
    str_polylines = []

    # 第二次拆分，将各多段线内的坐标拆开
    for item in first_separate:
        line = item.split(";")
        polylines.append(line)

    # 依次遍历每一条多段线，拼接对应格式的字符串
    for line in polylines:
        list_lat = []
        list_lon = []
        for item in line:
            res = item.split(",")
            list_lon.append(res[0])
            list_lat.append(res[1])
        result_str = ""
        print list_lat.__len__(), "points are read."
        for i in range(list_lat.__len__()):
            result_str += list_lon[i] + " " + list_lat[i] + ","
        result_str = "POLYGON ((" + result_str[:-1] + "))"
        str_polylines.append(result_str)
    return str_polylines


# 写形文件
def WriteVectorFile(filename, content):
    # 为了支持中文路径，添加下面这句代码
    gdal.SetConfigOption("GDAL_FILENAME_IS_UTF8", "NO")
    # 为了使属性表字段支持中文，添加下面这句
    gdal.SetConfigOption("SHAPE_ENCODING", "")

    # 注册所有的驱动
    ogr.RegisterAll()

    # 创建数据
    strDriverName = "ESRI Shapefile"
    oDriver = ogr.GetDriverByName(strDriverName)
    if oDriver is None:
        print strDriverName, "驱动不可用"
        return

    # 创建数据源
    oDS = oDriver.CreateDataSource(filename)
    if oDS == None:
        print filename, "创建文件失败！"
        return

    # 创建图层，创建一个多边形图层，这里没有指定空间参考，如果需要的话，需要在这里进行指定
    papszLCO = []
    oLayer = oDS.CreateLayer("Polygon", None, ogr.wkbPolygon, papszLCO)
    if oLayer == None:
        print("图层创建失败！\n")
        return

    # 下面创建属性表
    # 先创建一个叫FieldID的整型属性
    oFieldID = ogr.FieldDefn("FieldID", ogr.OFTInteger)
    oLayer.CreateField(oFieldID, 1)

    # 再创建一个叫FeatureName的字符型属性，字符长度为50
    oFieldName = ogr.FieldDefn("FieldName", ogr.OFTString)
    oFieldName.SetWidth(100)
    oLayer.CreateField(oFieldName, 1)

    oDefn = oLayer.GetLayerDefn()

    for item in content:
        # 创建要素
        oFeature = ogr.Feature(oDefn)
        oFeature.SetField(0, 0)
        oFeature.SetField(1, "Border")
        geom = ogr.CreateGeometryFromWkt(item)
        oFeature.SetGeometry(geom)
        oLayer.CreateFeature(oFeature)

    oDS.Destroy()
    print("形文件创建完成！\n")


# path为文本文件路径
# spl_str为坐标分隔符，如空格、逗号、制表符等
# format为坐标样式，int类型变量
# 1：latitude longitude
# 2：longitude latitude
def ReadContent(path, spl_str, format):
    str_polylines = []
    list_lat = []
    list_lon = []
    # 判断坐标顺序
    if format == 1:
        for line in open(path):
            res = line.split(spl_str)
            list_lat.append(res[0])
            list_lon.append(res[1].replace("\n", ""))
    elif format == 2:
        for line in open(path):
            res = line.split(spl_str)
            list_lon.append(res[0])
            list_lat.append(res[1].replace("\n", ""))

    result_str = ""
    print list_lat.__len__(), "points are read."
    for i in range(list_lat.__len__()):
        result_str += list_lon[i] + " " + list_lat[i] + ","
    result_str = "POLYGON ((" + result_str[:-1] + "))"
    str_polylines.append(result_str)
    return str_polylines


mode = input("Select mode:\n1:Load points from file.\n2:Get points from Internet.\n")
if mode == 1:
    filename = raw_input("Input file path(e.g. E:\input.txt):\n")

    separator = raw_input("Input separator of text(for \"\\t\" input \"t\"):\n")
    if separator == "t":
        separator = "\t"

    format = input("Select format of text:\n1:latitude longitude\n2:longitude latitude\n")
    str_points = ReadContent(filename, separator, format)

    out = raw_input("Input output filename(e.g. F:\out.shp):\n")
    WriteVectorFile(out, str_points)
elif mode == 2:
    keyword = raw_input("Input name of strict:\n")
    str_points = GetBoundary(keyword)
    out = raw_input("Input output filename(e.g. F:\out.shp):\n")
    WriteVectorFile(out, str_points)