使用高德地图API获取南京市内的医院POI（python，含代码）

前言

      python在获取数据方面具有较为强大的功能，高德地图API则提供了相当充足的数据。通过编译python代码，使用多边形搜索对南京市范围内的医院POI点进行获取，再使用地理方法对其进行可视化展现和核密度分析，可以对南京市的医院种类、分布进行较为直观的了解、分析。

（自己学习总结的，包能跑的捏）

技术路线

      主要利用python、excel和ArcGIS三个平台，依次进行代码编写和运行、数据处理以及出图三个步骤。在代码实现中，通过高德地图api所提供的服务和参数，确认爬虫范围，并在外接矩形中分化四叉树，在分化的子区域中进行POI信息提取，并根据POI点密度决定是否再分化，将所获取的数据存储在csv文件中。

      在数据处理中，由于编码格式问题，csv文件直接显示存在乱码问题，所以将其转化为GB18030格式，并作适当修改使其符合xy数据格式；在ArcGIS中加载南京市行政区域底图和csv数据，csv数据转化为点格式的shp文件并导出保存，最后进行核密度分析和出图处理。

方法与过程

    1.代码实现

    (1) GetPoi_keywords文件   

      主要包括所需功能函数：gcj02towgs84用于将高德地图坐标系转化为WGS84坐标系，Get_poi_polygon用于在单个格网中获取poi数据，Get_times_polygon用于在多个格网中申请循环获取，writecsv将读取的数据转化为csv文件，get_city_scope获取南京市行政区范围。

import requests
import json
import csv
import re
def gcj02towgs84(localStr): # 高德坐标系转换为 WGS84 
    lng = float(localStr.split(',')[0])
    lat = float(localStr.split(',')[1])
    PI = 3.1415926535897932384626
    ee = 0.00669342162296594323
    a = 6378245.0
    dlat = transformlat(lng - 105.0, lat - 35.0)
    dlng = transformlng(lng - 105.0, lat - 35.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)
    dlng = (dlng * 180.0) / (a / sqrtmagic * math.cos(radlat) * PI)
    mglat = lat + dlat
    mglng = lng + dlng
    return [lng * 2 - mglng,lat * 2 - mglat]
#单个格网内申请数据
def Get_poi_polygon(key,polygon,keywords,page):
    #设置header
    header = {'User-Agent': "Mozilla/5.0 (Windows; U; Windows NT 6.1; en-us) AppleWebKit/534.50 (KHTML, like Gecko) Version/5.1 Safari/534.50"}
    #格式化输入矩形
    Polygonstr = str(polygon[0]) + ',' + str(polygon[1]) + '|' + str(polygon[2]) + ',' + str(polygon[3])
    #构建url
url = 'https://restapi.amap.com/v3/place/polygon?polygon={}&key={}&keywords={}&page={}'.format(Polygonstr, key, keywords, page)
    print(url)
    #用get函数请求数据
    r = requests.get(url, headers=header)
    #设置数据编码格式
    r.encoding = 'utf-8'
    data = r.text
    return data
#多个格网循环申请
def Get_times_polygon(key,polygon,keywords):
    page = 1
    # count为0的时候跳出循环
    while True:
        result = Get_poi_polygon(key,polygon, keywords, page)
        # json.loads解码json
        content = json.loads(result)
        pois = content['pois']
        count = content['count']
        print(count)
        #如果区域内poi的数量大于800，则认为超过上限，返回False请求对区域进行切割
        if int(count) > 800:
            return False
        else:
            for i in range(len(pois)):
                name = pois[i]['name']
                location = pois[i]['location']
                if 'address' not in pois[i].keys():
                    address = str(-1)
                else:
                    address = pois[i]['address']
                adname = pois[i]['adname']
                result = location 
                lng = result[0]
                lat = result[1]
                row = [name, address, adname, lng, lat]
                #保存数据为csv
                writecsv(row, keywords)
            if count == '0':
                break
            page = page + 1
def writecsv(poilist,keywords):
    # 保存为csv
    with open('{}.csv'.format(keywords),'a',newline='',encoding='utf-8') as csvfile:
        writer = csv.writer(csvfile)
        writer.writerow(poilist)
# 获取行政区范围
def get_city_scope(key, cityname):
    parameters = 'key={}&keywords={}&subdistrict={}&output=JSON&extensions=all'.format(key, cityname, 0)
    url = 'https://restapi.amap.com/v3/config/district?'
    # 设置header
    header = {'User-Agent': "Mozilla/5.0 (Windows; U; Windows NT 6.1; en-us) AppleWebKit/534.50 (KHTML, like Gecko) Version/5.1 Safari/534.50"}
    res = requests.get(url,params=parameters)
    jsonData = res.json()
    if jsonData['status'] == '1':
        district = jsonData['districts'][0]['polyline']
        district_list = re.split(';|\|',district)
        xlist, ylist = [], []
        for d in district_list:
            xlist.append(float(d.split(',')[0]))
            ylist.append(float(d.split(',')[1]))
        xmax = max(xlist)
        xmin = min(xlist)
        ymax = max(ylist)
        ymin = min(ylist)
        return [xmin, ymax, xmax, ymin]
    else:
        print ('fail to acquire: {}'.format(jsonData['info']))

(2) RectanSearch文件

import GetPoi_keywords as gp
#四叉树分化并执行
def Quadrangle(key,polygon,keywords):
    #空列表，存放切割后的子区域
    PolygonList = []
    for i in range(len(polygon)):
        currentMinlat = round(polygon[i][0],6)#最小经度  #最大经度、最大纬度、最小纬度略
        cerrnt_list = [currentMinlat, currentMaxlon, currentMaxlat, currentMinlon]
        #将多边形输入获取函数中，判断区域内poi的数量
        status = gp.Get_times_polygon(key,cerrnt_list,keywords)
        #如果数量大于800，四叉树进行再次细分,否则返回区域的坐标对
        if status != False:
            print('该区域poi数量小于800，正在写入数据')
        else:
            #左上矩形
            PolygonList.append([
                currentMinlat, #左经
                currentMaxlon, #上纬
                (currentMaxlat+currentMinlat)/2, #右经
                (currentMaxlon+currentMinlon)/2]) #下纬
            #右上、左下、右下矩形略
            print(len(PolygonList))
            if len(PolygonList) == 0:
                break
            else:
                Quadrangle(key,PolygonList,keywords)
key =' '#在高德地图开放平台上申请
keywords = '医院'
city = '南京'
#调用高德查询行政区的API接口来返回矩形坐标对
Retance = gp.get_city_scope(key,city)
#存储区域矩形的列表
input_polygon = []
input_polygon.append(Retance)
Quadrangle(key,input_polygon,keywords)

2.数据处理与出图

      将乱码csv文件转存为GB18030格式，并添加序列号和表头，制作为每个Poi点的属性表。 

      将poi点导入ArcGIS，并使其坐标系与地图坐标系相吻合。结合空间信息或属性信息，去除范围在外接矩形内，但是在行政区外的poi点。删除多余点后导出该shp图层。

       后续出图按照标准要求即可，亦可自定义进行插值