基于dbscan和k-means算法的北京某品牌门店聚类实战

项目背景：
某个品牌全国经营着将近3000家门店，无法及时发现哪个片区的门店问题比较多，希望能对门店进行分级，按照商圈的维度去看问题，那么如何把门店进行商圈分类，就成了比较棘手的问题。设想的是有以下几种方法：
1.用门店地址经纬度返回其所在的business area
2.获取商圈中心坐标经纬度，再计算各门店到各商圈中心点的距离，按门店GROUPBY取最小距离的前3名进行观察；
3.基于门店经纬度的聚类分析，将门店分成簇，再通过爬取大众点评商圈类型（例如北京的王府井商圈，西单商圈，上海地陆家嘴商圈等)，把聚类成的门店簇对应到爬回来的商圈，可视化观察；
第三种方法又采取了k-means和dbscan两种方法以及其混合算法进行尝试；

因为北京的门店最多，我们先尝试用北京的282家门店进行聚类；

1.用门店地址经纬度返回其所在的business area

#方法1：用门店地址经纬度返回其所在的business area

import numpy as np
import pandas as pd
import requests

def regeocode(location):
    parameters1 = {'location': location, 'key': 'yourkey','radius': '1000', 'extensions' :'all','batch':'false','roadlevel':'0'}
    base1 = 'http://restapi.amap.com/v3/geocode/regeo'
    response1 = requests.get(base1, parameters1)
    answer1 = response1.json()
    #print (answer1)
    try:
        if answer1['regeocode']['addressComponent']['businessAreas'][0]['name'] ==[]:
            businessAreas=''
        else:
            businessAreas=answer1['regeocode']['addressComponent']['businessAreas'][0]['name']
    except:
            businessAreas = ''
    try:
        if answer1['regeocode']['pois'][0]['name'] ==[]:
            pois  =''
        else:
            pois = answer1['regeocode']['pois'][0]['name']  # 兴趣点
    except:
        pois =  ''
    try:
        if answer1['regeocode']['pois'][0]['businessarea'] ==[]: #兴趣点商圈
            poi_businessarea=''
        else:
            poi_businessarea=answer1['regeocode']['pois'][0]['businessarea']
    except:
        poi_businessarea = ''

    return businessAreas,pois,poi_businessarea

if __name__=='__main__':
    path=r'D:\商圈门店\北京商圈门店坐标运营指标(1).xlsx'
    df1=pd.read_excel(path,sheet_name='门店list及经纬度')
    result=[]
    for i in range(len(df1)):
        dic1={}
        jingweidu=df1['经纬度'].iloc[i]
        shop_name=df1['StoreName_chn'].iloc[i]
        businessAreas,pois,poi_businessarea=regeocode(jingweidu)
        dic1['jingweidu']=jingweidu
        dic1['shp_name']=shop_name
        dic1['businessAreas']=businessAreas #返回商圈名称，具体可以查看高德api的《逆地理编码》的文档介绍
        dic1['pois']=pois #兴趣点
        dic1['poi_businessarea']=poi_businessarea
        result.append(dic1)
        print ('第%s条执行完毕'%i)
    result_df=pd.DataFrame(result)
    writer=pd.ExcelWriter(r'D:\商圈门店\门店商圈类型方案1.xlsx')
    result_df.to_excel(writer,index=False)
    writer.save()
    #print (businessAreas)
    #print (pois)
    #print (poi_businessarea)
    print ('success')

2.获取商圈中心坐标经纬度，再计算各门店到各商圈中心点的距离；
通过大众点评划分好的商圈，通过调取高德API返回该商圈关键地标的经纬度（例如天坛商圈，获取的就是天坛建筑的经纬度），
再分别计算每个门店的经纬度距离所有商圈的距离，将门店划分到距离最近的商圈中；

import numpy as np
import pandas as pd
import requests
import math
#获取商圈中心坐标经纬度
def geocode(address):
    parameters = {'address': address, 'key': 'yourkey','city':'北京'}
    base = 'http://restapi.amap.com/v3/geocode/geo'
    response = requests.get(base, parameters)
    answer = response.json()
    try:
        location =  answer['geocodes'][0]['location']    #  经纬度
        jingdu,weidu = location.split(',')
    except:
        location = ''
        jingdu =''
        weidu  =''
    try:
        district=answer['geocodes'][0]['district']
    except:
        district=''
    try:
        neighbor_name=answer['geocodes'][0]['neighborhood']['name']
    except:
        neighbor_name=''
    try:
        neighbor_type=answer['geocodes'][0]['neighborhood']['type']
    except:
        neighbor_type=''
    try:
        building_name=answer['geocodes'][0]['building']['name']
    except:
        building_name=''
    try:
        building_type=answer['geocodes'][0]['building']['type']
    except:
        building_type=''
    return  address,jingdu,weidu,district

#计算门店到各商圈中心点的距离，定义距离函数，（x1,y1）和（x2,y2）分别为商圈和门店的经纬度坐标
def fun1(x1,y1,x2,y2):
    distance=6378137.0*2.0*math.asin(math.sqrt(math.pow(math.sin((y1-y2)*math.acos(-1)/360),2)+math.cos(y1*math.acos(-1)/180)*math.cos(y2*math.acos(-1)/180)*math.pow(math.sin((x1-x2)*math.acos(-1)/360),2)))
    return distance
	
#按门店GROUPBY取最小距离的前3名进行观察

if __name__=='__main__':
    path1=r'D:\商圈门店\北京商圈.xlsx' #读取商圈经纬度
    path2=r'D:\商圈门店\北京商圈门店坐标运营指标(1).xlsx' #读取门店经纬度
    df1=pd.read_excel(path1,sheet_name='Sheet1')
    df2=pd.read_excel(path2,sheet_name='门店list及经纬度')
    result=[]
    for i in range(len(df1)):
        for j in range(len(df2)):
            dic1={}
            sq_name=df1['商圈'].iloc[i]
            address,jingdu,weidu,district=geocode(sq_name)
            if jingdu != '':
                x1=float(jingdu)
                y1=float(weidu)
                jingweidu=df2['经纬度'].iloc[j]
                shop_name=df2['StoreName_chn'].iloc[j]
                x2 = float(jingweidu.split(',')[0])
                y2= float(jingweidu.split(',')[1])
                distance=fun1(x1,y1,x2,y2)
                #为了防止有些偏远门店是孤立点，与他最近的商圈可能较远，被错误的划分到该商圈，但实际我们认为他就是一个单独的孤立门店，因此我们这里将距离超过5公里以外的都不算做该商圈；
                if distance<5000:
                    dic1['sq_name']=sq_name
                    dic1['sq_jingdu']=jingdu
                    dic1['sq_weidu']=weidu
                    dic1['district']=district
                    dic1['shop_name']=shop_name
                    dic1['shop_jingweidu']=jingweidu
                    dic1['distance']=distance
                    result.append(dic1)
            else:
                continue
        print ('第%s个商圈执行完成'%i)
    print ('step1:sucess')
    result_df=pd.DataFrame(result)
    writer=pd.ExcelWriter(r'D:\商圈门店\门店商圈类型方案2.xlsx')
    result_df.to_excel(writer,index=False)
    writer.save()
    print ('end')

最后，通过对2种方案的结果进行对比，并抽查了几个不一致门店商圈类型，进行了门店地理位置搜索，最终确定选择方案2的划分方式；
与品牌方爸爸业务会议后，品牌方认为最后划分的商圈仍然过去细，比如世贸天阶和建外大街这种其实是非常近的，可以把他认为是一个，这就让我想到了通过密度去聚类，因此有了下面的第三种方法；

3.门店聚类
3.1 dbscan聚类

#DBSCAN算法的重点是选取的聚合半径参数和聚合所需指定的MinPts数目。
#在此使用球面距离来衡量地理位置的距离，来作为聚合的半径参数。
#如下实验，选取2公里作为密度聚合的半径参数，MinPts个数为5.
import numpy as np
import pandas as pd
from sklearn.cluster import DBSCAN
from shapely.geometry  import MultiPoint
import shapefile
import math
from scipy.spatial.distance import pdist,squareform
from scipy.cluster.vq import vq,kmeans,whiten,kmeans2
from matplotlib import pyplot as plt

def my_get_address_text_by_location(location):
    parameters1 = {'location': location, 'key': '329937aa1c76185568ab0fcda2cf4680','radius': '1000', 'extensions' :'all','batch':'false','roadlevel':'0'}
    base1 = 'http://restapi.amap.com/v3/geocode/regeo'
    response1 = requests.get(base1, parameters1)
    answer1 = response1.json()
    try:
        if answer1['regeocode']['formatted_address'] ==[]:
            formatted_address=''
        else:
            formatted_address=answer1['regeocode']['formatted_address']
    except:
        formatted_address = ''
    return formatted_address

def haversine(lonlat1, lonlat2):
    lat1, lon1 = lonlat1
    lat2, lon2 = lonlat2
    lon1, lat1, lon2, lat2 = map(math.radians, [lon1, lat1, lon2, lat2])
    dlon = lon2 - lon1
    dlat = lat2 - lat1
    a = math.sin(dlat / 2) ** 2 + math.cos(lat1) * math.cos(lat2) * math.sin(dlon / 2) ** 2
    c = 2 * math.asin(math.sqrt(a))
    r = 6371  # Radius of earth in kilometers. Use 3956 for miles
    return c * r

def clustering_by_dbscan(X):
    distance_matrix = squareform(pdist(X, (lambda u, v: haversine(u, v))))
    db = DBSCAN(eps=1, min_samples=2, metric='precomputed')  
    y_db = db.fit_predict(distance_matrix)
    X['cluster'] = y_db
    plt.scatter(X['lat'], X['lng'], c=X['cluster'])
    plt.show()
    return X

3.2基于dbscan和k-means混合算法聚类


其中，望京，王府井等地的门店数是最多的。如下是各门店分布地图，同一个颜色为相同的商圈，从颜色分布来看划分效果不错；

具体代码实现：

'''
.....
同上脚本复制过来
'''
def clustering_by_dbscan_and_kmeans2(X):
    distance_matrix = squareform(pdist(X, (lambda u, v: haversine(u, v))))
    db = DBSCAN(eps=1, min_samples=2, metric='precomputed')  
    y_db = db.fit_predict(distance_matrix)
    X['cluster'] = y_db
    results = {}
    #将X从dataframe转成np.array x.values
    for i in X.values:
        if i[2] not in results.keys():
            results[i[2]] = [[i[1], i[0]]]
        else:
            if results[i[2]]:
                results[i[2]].append([i[1], i[0]])
            else:
                results[i[2]] = [[i[1], i[0]]]
    print ("DBSCAN output: ", len(results), results.keys())
    #print (results['0'])
    print ("KMeans calc center as below: ")
    for k in results.keys():
        xy = np.array(results[k])
        #print (xy)
        z = np.sin(xy[:, 1] - 0.2 * xy[:, 1])
        z = whiten(z)
        tt=np.array(zip(xy[:, 0], xy[:, 1], z))
        print (type(tt))
        print (tt)
        res, idx = kmeans2(tt, 1, iter=20, minit='points')
        location=res[0][1]+','+res[0][0]
        print (location)
        address_text = my_get_address_text_by_location(location)
        print (res, address_text)
        #return results

if __name__=='__main__':
    path=r'D:\商圈门店\北京商圈门店坐标运营指标(1).xlsx'
    df1=pd.read_excel(path,sheet_name='门店list及经纬度')
    x=df1[['lat','lng']]
    clustering_by_dbscan_and_kmeans2(x)
    #print ('execute end')
    #result=clustering_by_dbscan(x)
    #writer=pd.ExcelWriter(r'D:\商圈门店\dbscanresult.xlsx')
    #result.to_excel(writer,index=False)
    #writer.save()

该种方法分类结果受到了品牌方的认可，后面在此商圈标签基础上和业务一起做了一些探索性分析；