qgis自动提取百度地图aoi插件

• 背景

作为一位gis工作者，有时候会遇到这样一个难题就是我们想在互联网地图查找某个位置想把这个位置的poi点数据或者aoi面数据绘制到gis软件中，往往是很困慢的，经过我不断考察验证决定在qgis上实现这一功能，于是我开发了一个qgis插件轻松绘制地图查找到的gis数据（由于百度地图aoi面数据比高德、腾讯地图多一些就先从百度入手）。

• 相关工具

百度地图ak、sk

Python（3.9+）

PyCharm

Qgis（3.20+）

Qt（qgis安装自带）

• 申请百度地图ak

进入百度地图开发者平台官网：https://lbsyun.baidu.com，点击控制台，先经过一番操作获得百度地图开放者账号；

然后在控制台》应用管理》我的应用》创建应用；

最终你就会得到你要的百度地图ak和sk

• 获取aoi数据代码

import requests
import hashlib
import urllib
import json


#bd墨卡托转BD-09
import math
pi = 3.1415926535897932384626
def Yr(lnglat,b):
    if b!='':
        c=b[0]+b[1]*abs(lnglat[0])
        d=abs(lnglat[1]/b[9])
        d=b[2]+b[3]*d+b[4]*d*d+b[5]*d*d*d+b[6]*d*d*d*d+b[7]*d*d*d*d*d+b[8]*d*d*d*d*d*d
        if 0>lnglat[0]:
            bd=-1*c
        else:
            bd=c
        lnglat[0]=bd
        if 0 > lnglat[0]:
            bd2 = -1 * d
        else:
            bd2 = d
        lnglat[1] = bd2
        return lnglat
    return
def Mecator2BD09(lng,lat):
    lnglat=[0,0]
    Au=[[1.410526172116255E-8, 8.98305509648872E-6, -1.9939833816331, 200.9824383106796, -187.2403703815547,
          91.6087516669843, -23.38765649603339, 2.57121317296198, -0.03801003308653, 1.73379812E7],
         [- 7.435856389565537E-9, 8.983055097726239E-6, -0.78625201886289, 96.32687599759846, -1.85204757529826,
          -59.36935905485877, 47.40033549296737, -16.50741931063887, 2.28786674699375, 1.026014486E7],
         [- 3.030883460898826E-8, 8.98305509983578E-6, 0.30071316287616, 59.74293618442277, 7.357984074871,
          -25.38371002664745, 13.45380521110908, -3.29883767235584, 0.32710905363475, 6856817.37],
         [- 1.981981304930552E-8, 8.983055099779535E-6, 0.03278182852591, 40.31678527705744, 0.65659298677277,
          -4.44255534477492, 0.85341911805263, 0.12923347998204, -0.04625736007561, 4482777.06],
         [3.09191371068437E-9, 8.983055096812155E-6, 6.995724062E-5, 23.10934304144901, -2.3663490511E-4,
          -0.6321817810242, -0.00663494467273, 0.03430082397953, -0.00466043876332, 2555164.4],
         [2.890871144776878E-9, 8.983055095805407E-6, -3.068298E-8, 7.47137025468032, -3.53937994E-6, -0.02145144861037,
          -1.234426596E-5, 1.0322952773E-4, -3.23890364E-6, 826088.5]]
    Sp=[1.289059486E7, 8362377.87, 5591021, 3481989.83, 1678043.12, 0 ]
    lnglat[0]=math.fabs(lng)
    lnglat[1] =abs(lat)
    for d in range(0,6):
        if lnglat[1]>=Sp[d]:
            c=Au[d]
            break
    lnglat=Yr(lnglat,c)
    return lnglat
def BD092WGS84(lnglat):

    x_pi = 3.14159265358979324 * 3000.0 / 180.0
    pi = 3.1415926535897932384626  # π
    a = 6378245.0  # 长半轴
    ee = 0.00669342162296594323  # 扁率
    x = lnglat[0] - 0.0065
    y = lnglat[1] - 0.006
    z = math.sqrt(x * x + y * y) - 0.00002 * math.sin(y * x_pi)
    theta = math.atan2(y, x) - 0.000003 * math.cos(x * x_pi)
    lnglat[0] = z * math.cos(theta)
    lnglat[1] = z * math.sin(theta)

    dlat = tranlat1(lnglat[0] - 105.0, lnglat[1] - 35.0)
    dlng = tranlng1(lnglat[0] - 105.0, lnglat[1] - 35.0)
    radlat = lnglat[1] / 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 = lnglat[1] + dlat
    mglng = lnglat[0] + dlng
    return [lnglat[0]* 2 - mglng, lnglat[1] * 2 - mglat]
def tranlat1(lng, lat):
    ret = -100.0 + 2.0 * lng + 3.0 * lat + 0.2 * lat * lat + 0.1 * lng * lat + 0.2 * math.sqrt(math.fabs(lng))
    ret += (20.0 * math.sin(6.0 * lng * pi) + 20.0 *
            math.sin(2.0 * lng * pi)) * 2.0 / 3.0
    ret += (20.0 * math.sin(lat * pi) + 40.0 *
            math.sin(lat / 3.0 * pi)) * 2.0 / 3.0
    ret += (160.0 * math.sin(lat / 12.0 * pi) + 320 *
            math.sin(lat * pi / 30.0)) * 2.0 / 3.0
    return ret
def tranlng1(lng, lat):
    ret = 300.0 + lng + 2.0 * lat + 0.1 * lng * lng + \
          0.1 * lng * lat + 0.1 * math.sqrt(math.fabs(lng))
    ret += (20.0 * math.sin(6.0 * lng * pi) + 20.0 *
            math.sin(2.0 * lng * pi)) * 2.0 / 3.0
    ret += (20.0 * math.sin(lng * pi) + 40.0 *
            math.sin(lng / 3.0 * pi)) * 2.0 / 3.0
    ret += (150.0 * math.sin(lng / 12.0 * pi) + 300.0 *
            math.sin(lng / 30.0 * pi)) * 2.0 / 3.0
    return ret

#geojson模版
mb = {
    "type": "FeatureCollection",
    "name": "cs",
    "crs": {
        "type": "name",
        "properties": {
            "name": "urn:ogc:def:crs:EPSG::4326"
        }
    },
    "features": [{
            "type": "Feature",
            "properties": {
                "ID": 0
            },
            "geometry": {
                "type": "MultiPolygon",
                "coordinates": [[[]]]
            }
        }
    ]
}

#提取函数
def getbdwzi(dz,ct,aky,sky):
    datab = mb
    # 参数拼接sign电子签名验证
    baiduurl = "https://api.map.baidu.com"
    queryStr = "/place/v2/suggestion?query=" + dz + "&region=" + ct + "&city_limit=true&coord_type=1&output=json&ak=" + aky
    encodedStr = urllib.parse.quote(queryStr, safe="/:=&?#+!$,;'@()*[]")
    rawStr = encodedStr + sky
    sign = hashlib.md5(urllib.parse.quote_plus(rawStr).encode("utf8")).hexdigest()
    qqurl = baiduurl + queryStr + '&sn=' + sign
    # 获取请求数据
    data = requests.request('GET', qqurl).content.decode('utf-8')
    dataa = json.loads(data)
    # 判断请求结果
    if dataa['status'] == 0:
        # aoi参数拼接获取aoi
        AOI_id = dataa['result'][0]['uid']
        uel_AOI = 'https://map.baidu.com/?newmap=1&qt=ext&uid=' + AOI_id + '&ext_ver=new&ie=utf-8&l=11'
        r_AOI = requests.request('GET', uel_AOI).content.decode('utf-8')
        data_AOI = json.loads(r_AOI)
        point_transform = []
        # poi点位坐标转换
        xy = dataa['result'][0]['location']
        nxy = BD092WGS84([float(xy['lng']), float(xy['lat'])])
        dataa['result'][0]['location'] = nxy
        # 判断是否有aoi
        if 'geo' in data_AOI['content']:
            data_AOI['content']['geo']
            geo_AOI = data_AOI['content']['geo']
            geo_AOI = geo_AOI.split('|')
            point = geo_AOI[2].split(",")
            #aoi坐标转换
            for i in range(int(len(point) / 2)):  # 全部点的坐标，分别是x,y,的形式
                if i == 0:  # 第一个点的x坐标删除‘1-’
                    point[2 * i] = point[2 * i][2:]
                if i == int((len(point) / 2) - 1):  # 最后的点的y坐标删除‘;’
                    point[2 * i + 1] = point[2 * i + 1][:-1]
                point_Mecator2BD09 = Mecator2BD09(float(point[2 * i]), float(point[2 * i + 1]))
                point_BD092WGS84 = BD092WGS84(point_Mecator2BD09)
                point_transform.append(point_BD092WGS84)

            datab['features'][0]['geometry']['coordinates'][0][0] = point_transform
        else: # 若没有aoi则返回poi点数据
            datab['features'][0]['geometry']['type'] = "Point"
            datab['features'][0]['geometry']['coordinates'] = nxy

        datab['name'] = dz
        datab['features'][0]['properties'] = dataa['result'][0]
        dataz = json.dumps(datab, ensure_ascii=False)

        # 数据写的geojson
        with open(dz + ".geojson", 'w', encoding='utf-8') as geo:
            geo.write(dataz)

    else:
        print('未找到')


if __name__ == '__main__':
    dz = '深圳北站'
    ct ='深圳市'
    aky = '4wayPmjAKZzmzzr1YXSjqVpQ70V9Ct5U'
    sky = 'dciEGgSKzMiTe5myW9Y62jrFGZL0zhTS'
    getbdwzi(dz,ct,aky,sky)

• Qgis插件制作

首先Plugin Builder3，这个制作插件的插件

剩下步骤参考这位大佬就不一一详述了https://www.cnblogs.com/wsh233/p/16976884.html

• 插件安装使用

到我的git网站下载插件zip包https://github.com/srmaper/baiduAOI

安装插件：插件》安装并管理插件》从zip文件安装

使用插件：矢量》bdaoi》aoi

输入查找地点、城市、百度ak、百度sk点击确定自动定位到要找的位置并且生成矢量图层

大功告成！！！（技术有限仅供参考）

*********************************************************************************

承接各种gis处理制图分析爬取任务及小工具开发