python爬取高德地图乡镇区行政区划

适用于爬取一个省、市、县、区等各个行政级别的行政区划数据，数据精确到区

级别。

以湖南省为例：

# -*- coding:utf-8 -*-

import arcpy

import urllib2

import sys

import json

import math

# 申请：高德API的开发者KEY

# 版本要求：ArcGis 10.2 Python 2.7

# 输入：输入省的名称 存储路径 文件名称

# 输出：该省的行政区划的shp 字段 city citycode district adcode

# 获取行政区划范围内的乡镇区划列表 并插入数据

def get_ad(name):

    url = "http://restapi.amap.com/v3/config/district?extensions=base&subdistrict=3&key="+key+"&keywords="+name

    # extensions说明：all只返回当前查询的district边界坐标点

    # subdistrict说明：1表示下一级行政区（到市），2表示下两级行政区（到乡镇区县），3表示下三级行政区（到街道）

    try:

        req = urllib2.Request(url)

        res_data = urllib2.urlopen(req)

        res = res_data.read()

        data = json.loads(res)

        # 返回查询字段的行政区列表

        city_list = data['districts'][0]['districts']

        for city in city_list:

            district_list = city['districts']

            districts = {}  # 存储district列表

            for district in district_list:

                districts[district['adcode'].encode('utf-8')] = district['name'].encode('utf-8')

            # 添加数据

            set_polyline(districts, city['name'].encode('utf-8'), city['citycode'].encode('utf-8'))

    except urllib2.HTTPError, e:

        if hasattr(e, "code"):

            print e.code

        if hasattr(e, "reason"):

            print e.reason

def get_polyline(adcode):

    url = "http://restapi.amap.com/v3/config/district?extensions=all&subdistrict=0&key=" + key

    _url = url + "&keywords=" + adcode

    try:

        req = urllib2.Request(_url)

        res_data = urllib2.urlopen(req)

        res = res_data.read()

        data = json.loads(res)

        polyline = data['districts'][0]['polyline']

        return polyline

    except urllib2.HTTPError, e:

        if hasattr(e, "code"):

            print e.code

        if hasattr(e, "reason"):

            print e.reason

def prepare_shp():

    fields_long = ['city', 'district']

    fields_short = ['citycode', 'adcode']

    length_long = 20

    length_short = 10

    sr = arcpy.SpatialReference(4326)

    if arcpy.Exists(file):

        arcpy.Delete_management(file)

    polygon = arcpy.CreateFeatureclass_management(path, shp_name, 'POLYGON', '', '', '', sr)

    print 'create:' + shp_name

    for field in fields_long:

        arcpy.AddField_management(polygon, field, 'TEXT', '', '', length_long)

        print 'add filed:' + field

    for field in fields_short:

        arcpy.AddField_management(polygon, field, 'TEXT', '', '', length_short)

        print 'add filed:' + field

def set_polyline(districts_map, city, citycode):

    insert_cursor = arcpy.da.InsertCursor(file, ('SHAPE@', 'city', 'district', 'citycode', 'adcode'))

    for adcode in districts_map.keys():

        polyline = get_polyline(adcode)

        array = arcpy.Array()

        for j in polyline.split('|'):

            arr = arcpy.Array()

            for i in j.split(';'):

                pt = arcpy.Point()

                pt.X = gcj2wgs(i)[0]

                pt.Y = gcj2wgs(i)[1]

                arr.append(pt)

            array.append(arr)

        polygon = arcpy.Polygon(array)

        insert_cursor.insertRow((polygon, city, districts_map[adcode], citycode, adcode))

    print "data create:" + city

    del insert_cursor

# 高德火星坐标系转化为WGS84地理坐标系

def gcj2wgs(loc):

    lon = float(loc.split(',')[0])

    lat = float(loc.split(',')[1])

    a = 6378245.0  # 克拉索夫斯基椭球参数长半轴a

    ee = 0.00669342162296594323  # 克拉索夫斯基椭球参数第一偏心率平方

    PI = 3.14159265358979324  # 圆周率

    # 以下为转换公式

    x = lon - 105.0

    y = lat - 35.0

    # 经度

    dLon = 300.0 + x + 2.0 * y + 0.1 * x * x + 0.1 * x * y + 0.1 * math.sqrt(abs(x))

    dLon += (20.0 * math.sin(6.0 * x * PI) + 20.0 * math.sin(2.0 * x * PI)) * 2.0 / 3.0

    dLon += (20.0 * math.sin(x * PI) + 40.0 * math.sin(x / 3.0 * PI)) * 2.0 / 3.0

    dLon += (150.0 * math.sin(x / 12.0 * PI) + 300.0 * math.sin(x / 30.0 * PI)) * 2.0 / 3.0

    # 纬度

    dLat = -100.0 + 2.0 * x + 3.0 * y + 0.2 * y * y + 0.1 * x * y + 0.2 * math.sqrt(abs(x))

    dLat += (20.0 * math.sin(6.0 * x * PI) + 20.0 * math.sin(2.0 * x * PI)) * 2.0 / 3.0

    dLat += (20.0 * math.sin(y * PI) + 40.0 * math.sin(y / 3.0 * PI)) * 2.0 / 3.0

    dLat += (160.0 * math.sin(y / 12.0 * PI) + 320 * math.sin(y * PI / 30.0)) * 2.0 / 3.0

    radLat = lat / 180.0 * PI

    magic = math.sin(radLat)

    magic = 1 - ee * magic * magic

    sqrtMagic = math.sqrt(magic)

    dLat = (dLat * 180.0) / ((a * (1 - ee)) / (magic * sqrtMagic) * PI);

    dLon = (dLon * 180.0) / (a / sqrtMagic * math.cos(radLat) * PI);

    wgsLon = lon - dLon

    wgsLat = lat - dLat

    return wgsLon, wgsLat

reload(sys)

sys.setdefaultencoding('utf-8')

key = "你的KEY"

search_name = "湖南省"  # 输入要获取的行政区划名称

path = "C:/ArcgisData/data" //选择一个存储路径

shp_name = search_name + ".shp"

file = path + "/" + shp_name

prepare_shp()

get_ad(search_name)

print "create success"

————————————————

版权声明：本文为CSDN博主「K.I.O」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。

原文链接：https://blog.csdn.net/sinat_29675423/article/details/87434839