从高德地图获取最新行政区划边界数据

August 2018 · 2 minute read

行政区划边界数据在很多地方都会用到。然而，国内公开的数据一般跟不上最新的变更，国外公开的数据因为问题很多（国界线、南海、台湾等等）所以最好别用。

另外，一般人可能没有注意到，很多软件自带的中国国界数据也是有问题的，最好别用，即便仅仅是作为 GIS 底图。

不少文献中的底图，中国的行政边界都有问题。

有没有办法获取 正确 且 精细 且 更新及时 的行政区划边界数据呢？几经搜寻，发现高德地图满足要求。

现在我们需要从高德的 API 里面查出 country、province、city、district，也就是国、省、市、区县四个级别的行政区划边界。 当然，还有乡镇一级的行政区划边界，有需要的可以参照以下的方法自行处理。

amap_CN.png

用浏览器在找到的连接处点击鼠标右键，在弹出的菜单中找到 Copy 下的 Copy as cURL(bash)，然后粘贴进 shell 命令行回车，我们就看到了原始数据的样子： 返回的是 jsonp，里面包含了一段 json 数据，看到 level: "country", name: "中华人民共和国", polyline: "120.823872,40.530257;..."， 显然，polyline 就是边界数据，再仔细看，; 分隔的是经纬度点，| 分割的是不同的多边形。

找到了？是的，我们找到了行政区划边界数据的 API，下面就好办了。

从国家到区县，很多行政边界并不仅仅是一个多边形，而是由多个多边形组成的，这一点要特别注意。

接下来我们可以编写一段程序，自动查询任意行政边界，并利用 pyshp 转换成 shp 格式保存到本地：

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import os
import re
import json
import requests
import shapefile
from xpinyin import Pinyin  # 便于转换中文名称为拼音


BASE_URL = u'http://restapi.amap.com/v3/config/district?subdistrict=1&extensions=all&level=district&key=608d75903d29ad471362f8c58c550daf&s=rsv3&output=json&keywords={name}&callback=jsonp_923427_&platform=JS&logversion=2.0&sdkversion=1.3&appname=http%3A%2F%2Flbs.amap.com%2Fapi%2Fjavascript-api%2Fexample%2Fdistrict-search%2Fdraw-district-boundaries%2F&csid=7EB50235-A6C0-45C3-B08C-2F77FE1837AC'  # .format(name="中国")

HEADERS = {
    'Accept-Encoding': 'gzip, deflate, sdch',
    'Accept-Language': 'zh-CN,zh;q=0.8',
    'User-Agent': 'Mozilla/5.0 (Windows NT 10.0; WOW64) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/52.0.2743.116 Safari/537.36',
    'Accept': '*/*',
    'Referer': 'http://lbs.amap.com/api/javascript-api/example/district-search/draw-district-boundaries/',
    'Cookie': 'guid=70c9-87fe-5807-3f17; cna=YKn+D5QKFQ0CAd3vO4CwxB9s; l=Ara234AMAU3gNSreb6UFtrzXhua4gPoX; isg=AqmphLSq2mevHubTbJlaWcpIuFVMnZ2o89H-lkuechDJEskkk8ateJdKogHf; key=608d75903d29ad471362f8c58c550daf',
    'Connection': 'keep-alive',
}


def get_bdr_by_name_allpoly(name, out_name='', out_dir='.'):
    if not out_name:
        out_name = Pinyin().get_pinyin(name, '')
    url = BASE_URL.format(name=name)
    r = requests.get(url=url, headers=HEADERS)
    if r.status_code == requests.codes.ok:
        js = r.text
    else:
        print('Fail with status_code = {}'.format(r.status_code))
    with open(os.path.join(out_dir, out_name), 'w') as o:
        o.write(js.encode("utf-8"))
    m = re.findall(r"jsonp_\d{1,10}_\((?P<json>.*)\)", js)
    js = m[0]
    js = json.loads(js)
    polylines = js.get('districts')[0].get('polyline')
    polylines = polylines.split('|')
    poly_list = []
    w = shapefile.Writer(shapefile.POLYGON)
    w.field(b'NAME', b'C')
    for polyline in polylines:
        polyline = polyline.split(';')
        x, y = [[float(p.split(',')[i]) for p in polyline] for i in (0, 1)]
        polygon = [[x[j], y[j]] for j in range(len(x))]
        w.poly(parts=[polygon])
        w.record(b'part')
        poly_list.append(polygon)
    w.save(os.path.join(out_dir, b'{out_name}'.format(out_name=out_name)))
    return poly_list


if __name__ == '__main__':
    get_bdr_by_name_allpoly(name=u'中国')
    get_bdr_by_name_allpoly(name=u'河北省')
    get_bdr_by_name_allpoly(name=u'武汉市')
    get_bdr_by_name_allpoly(name=u'若羌县')

    # 对于同名的地区，可以利用行政编码查询，避免混淆
    get_bdr_by_name_allpoly(name='210102')  # 沈阳市和平区行政编码
    get_bdr_by_name_allpoly(name='120101')  # 天津市和平区行政编码

我们用 ArcGIS Earth 打开看一下（别的 GIS 软件一样的）：

amap_shp_view.png

当然，还可以另外写为 geojson 或 kml 等格式，不难。

最后，需要 特别说明 的是，由于高德地图在国内采用的是 GCJ-02 坐标系（是在 WGS84 坐标上进行的随机加密，有几十到几百米的偏移），因此，如果对精度要求非常高，就必须经过相应的坐标转换（可自行从网上搜寻坐标转换方法），如果不在意，可以忽略。