【Matplotlib 使用】k-d tree 实现用户划分（聚类）；结合 Basemap 绘制划分效果图

目的

最近要绘制一张地图的分割效果图。给定所有点的经纬度坐标，使用 k-d tree 、DBSCAN 之类的聚类算法（分割算法）将所有点划分至不同的区域，并将区域的编号作为该点的类别标签。为了观察特定分割算法的分割效果，因此需要将分割后的各个区域的轮廓画出来。

先前的方法是将所有的轮廓看作四边形，在区域内部的所有点中直接去找矩形的4个“顶点”（即左上角、左下角、右上角、右下角），然后以地图为背景，在上面绘制一个个四边形得到地图分割效果图。这种做法最不好的一点就不精确，分割后的区域不一定都是四边形呀。

为了将多边形的轮廓也画出来，关键在于找到划分后的每个区域的轮廓，这便是 凸包问题 （指在 n 个点中，寻找一个凸多边形，使所有的点在凸多边形的边界或者内部），然后利用 Basemap 就可以将一个个多边形画出来了。

效果

• 原始的四边形轮廓绘制
  
• 寻找经纬度的凸包并在地图上绘制
  

实现

Graham 扫描法

凸包问题属于计算机图形学问题，一种常用的解法是 Graham 扫描法，运行时间为 $\mathcal{O}(n\log (n))$。

算法流程

1. 寻找 y 轴最小的点，如果 y 轴位置是相同的，那个找 x 轴位置最小的，称之为基准点。
2. 计算 1 中找到基准点与其他点的极角（即过此2点的直线与x轴正方向的夹角，代码中以弧度表示），将这些点按极角的大小正序排列。
3. 进行基准点与 2 中点的连线迭代，对新连线的点计算其是否符合凸多边形的定义，如果不满足舍弃此点。判断的方法是计算三点组成线段的叉乘，值为正表示满足条件。

动态模拟

参考

• Python求凸包及多边形面积
• 凸包问题，Graham扫描法，python实现
• matplotlib高级教程之形状与路径——patches和path

代码

版本

• matplotlib 2.2.2
• basemap 1.2.1
• python 2.7

思路

1. 调用 k-d tree 对经纬度数据点进行分割，得到若干个蔟类及其对应的数据点，可以设置 bucket_size 控制每个类别中数据点个数的最大值。

   clusterer = kdtree.KDTreeClustering(bucket_size=bucket_size)
   

2. 对每个蔟类调用 graham 算法得到凸边形的所有顶点位置。convex_hull 的具体实现在 graham_morton.py 中。

   hull_point = convex_hull(points)
   

3. 使用 mpath 将凸边形顶点按顺序连接起来形成一条首尾相连的路径。其中，MOVETO、LINETO、CURVE3、CURVE4、CLOSEPOLY 用来控制路径的线型。
4. 对路径对应的凸边形上色，并在 Basemap 预设的背景图之上将其绘制出来。

完整代码

• 主文件：需要准备好的经纬度数据。上述示例图中的数据文件 cmu_train_loc.pkl 的下载地址 点击下载

  # -*- coding: UTF-8 -*-
  import numpy as np
  import random
  import kdtree
  import pickle
  import gzip
  
  
  def load_obj(filename):
      with gzip.open(filename, 'rb') as fin:
          obj = pickle.load(fin)
      return obj
  
  
  def draw_kd_clusters_in_pathpatch(locs_data="None", bucket_size=50, savefile="./pic.png"):
      from collections import defaultdict as dd
      from matplotlib.collections import PatchCollection
      import matplotlib as mpl
      mpl.use('Agg')
      import matplotlib.pyplot as plt
      from mpl_toolkits.basemap import Basemap
  
      lllat = 24.396308
      lllon = -124.848974
      urlat = 49.384358
      urlon = -66.885444
  
      height = 25
      weight = 58
  
      fig = plt.figure(figsize=(weight, height), edgecolor='r')
      print fig.get_size_inches()
  
      m = Basemap(llcrnrlat=lllat, urcrnrlat=urlat, llcrnrlon=lllon,  urcrnrlon=urlon,
                  resolution='c', projection='cyl')
      m.drawmapboundary(fill_color='white')   # 背景
      m.drawcoastlines(linewidth=2)           # 海岸线
      m.drawcountries(linewidth=2)            # 国界
      m.drawstates(linewidth=0.5)             # 州界
      m.drawlsmask(land_color="0.8", ocean_color="w", lakes=True)  # 陆地和海洋的颜色
  
      '''split the point into multi-cluster using kd-tree.'''
      clusterer = kdtree.KDTreeClustering(bucket_size=bucket_size)
      train_locs = locs_data
      mlon, mlat = m(*(train_locs[:, 1], train_locs[:, 0]))
      train_locs = np.transpose(np.vstack((mlat, mlon)))
      clusterer.fit(train_locs)
      clusters = clusterer.get_clusters()
      cluster_points = dd(list)
      for i, cluster in enumerate(clusters):
          '''filter right-center points for cmu.'''
          if train_locs[i][1] > urlon:
              print("filter the out range sample:", train_locs[i])
              continue
  
          cluster_points[cluster].append(train_locs[i])
  
      '''get hull point and put patches into collection with specific colors.'''
      hull_point_list, number_of_points = get_convex_hull_of_discrete_points(cluster_with_point=cluster_points)
      hull_patches_list = get_path_of_hull_point(hull_point_list=hull_point_list)
      colors = [random.random() for i in range(0, len(hull_patches_list))]
      collection = PatchCollection(hull_patches_list, cmap=plt.cm.hsv, alpha=0.3)
      collection.set_array(np.array(colors))
  
      '''plot the collections.'''
      ax = plt.gca()
      ax.add_collection(collection)
      ax.set_xlim([lllon, urlon])
      ax.set_ylim([lllat, urlat])
  
      ax.spines['left'].set_linewidth(5)
      ax.spines['right'].set_linewidth(5)
      ax.spines['top'].set_linewidth(5)
      ax.spines['bottom'].set_linewidth(5)
  
      plt.tight_layout()
      plt.subplots_adjust(top=1, bottom=0, left=0, right=1, hspace=0, wspace=0)
      plt.savefig(savefile)
      print "the plot saved in " + savefile
  
  
  def get_path_of_hull_point(hull_point_list=list()):
      import matplotlib.path as mpath
      import matplotlib.patches as mpatches
  
      hull_patches_list = list()
      Path = mpath.Path
  
      for point_list in hull_point_list:
          path_data = list()
          num = len(point_list) - 1
          for id, point in enumerate(point_list):
              if id == 0:
                  path_data.append((Path.MOVETO, [point[0], point[1]]))
              elif id == num:
                  path_data.append((Path.CLOSEPOLY, [point[0], point[1]]))
              else:
                  path_data.append((Path.LINETO, [point[0], point[1]]))
          codes, verts = zip(*path_data)
          path = mpath.Path(verts, codes)
          patch = mpatches.PathPatch(path, linewidth=3)
          hull_patches_list.append(patch)
  
      # Path.MOVETO :拿起钢笔, 移动到给定的顶点。一般指的是 “起始点”
      # Path.LINETO :从当前位置绘制直线到给定顶点。
      # Path.CURVE3 :从当前位置 (用给定控制点) 绘制一个二次贝塞尔曲线到给定端点。
      # Path.CURVE4 :从当前位置 (与给定控制点) 绘制三次贝塞尔曲线到给定端点。
      # Path.CLOSEPOLY :将线段绘制到当前折线的起始点。
      # Path.STOP :整个路径末尾的标记
  
      print("hull_patches_list len is:", len(hull_patches_list))
      return hull_patches_list
  
  
  def get_convex_hull_of_discrete_points(cluster_with_point=list()):
      from graham_morton import convex_hull
  
      def list_to_tuple(temp):
          return (temp[1], temp[0])
      # test_data = [(220, -100), (0, 0), (-40, -170), (240, 50), (-160, 150), (-210, -150)]
      hull_point_list = list()
      c_number = len(cluster_with_point)
      number_of_points = list()
  
      for i in range(0, c_number):
          points = np.vstack(cluster_with_point[i]).tolist()
          points = list(map(list_to_tuple, points))
          number_of_points.append(len(points))
          hull_point = convex_hull(points)
          if len(hull_point) != 0:
              hull_point_list.append(hull_point)
      print("hull_point_list len is:", len(hull_point_list))
      return hull_point_list, number_of_points
  
  
  if __name__ == '__main__':
  
      train_loc = load_obj(filename="./data/cmu/cmu_train_loc.pkl")
      draw_kd_clusters_in_pathpatch(locs_data=train_loc, bucket_size=50, savefile="./pic_morton/kd-cmu-qxx.png")
  

• graham_morton.py 使用 Graham 算法寻找给定点集的凸包 （代码来自：Python求凸包及多边形面积）

  # -*- coding: UTF-8 -*-
  import math
  
  
  # 获取基准点的下标,基准点是p[k]
  def get_leftbottompoint(p):
      k = 0
      for i in range(1, len(p)):
          if p[i][1] < p[k][1] or (p[i][1] == p[k][1] and p[i][0] < p[k][0]):
              k = i
      return k
  
  
  # 叉乘计算方法
  def multiply(p1, p2, p0):
      return (p1[0] - p0[0]) * (p2[1] - p0[1]) - (p2[0] - p0[0]) * (p1[1] - p0[1])
  
  
  # 获取极角，通过求反正切得出，考虑pi/2的情况
  def get_arc(p1, p0):
      # 兼容sort_points_tan的考虑
      if (p1[0] - p0[0]) == 0:
          if (p1[1] - p0[1]) == 0:
              return -1
          else:
              return math.pi / 2
      tan = float((p1[1] - p0[1])) / float((p1[0] - p0[0]))
      arc = math.atan(tan)
      if arc >= 0:
          return arc
      else:
          return math.pi + arc
  
  
  # 对极角进行排序,排序结果list不包含基准点
  def sort_points_tan(p, pk):
      p2 = []
      for i in range(0, len(p)):
          p2.append({"index": i, "arc": get_arc(p[i], pk)})
      # print('排序前:',p2)
      p2.sort(key=lambda k: (k.get('arc')))
      # print('排序后:',p2)
      p_out = []
      for i in range(0, len(p2)):
          p_out.append(p[p2[i]["index"]])
      return p_out
  
  
  def convex_hull(p):
      p = list(set(p))
      # print('全部点:',p)
      k = get_leftbottompoint(p)
      pk = p[k]
      p.remove(p[k])
      # print('排序前去除基准点的所有点:',p,'基准点:',pk)
      p_sort = sort_points_tan(p, pk)   # 按与基准点连线和x轴正向的夹角排序后的点坐标
      # print('其余点与基准点夹角排序:',p_sort)
      if len(p_sort) != 0:
          p_result = [pk, p_sort[0]]
      else:
          p_result = [pk]
  
      for i in range(1, len(p_sort)):
          # 叉乘为正,向前递归删点;叉乘为负,序列追加新点
          while(multiply(p_result[-2], p_sort[i], p_result[-1]) > 0):
              p_result.pop()
          p_result.append(p_sort[i])
      return p_result
  
  
  if __name__ == '__main__':
      test_data = [(220, -100), (0, 0), (-40, -170), (240, 50), (-160, 150), (-210, -150)]
      print(test_data)
  
      result = convex_hull(test_data)
      print(result)
      x1, y1 = [], []
      for i in range(len(test_data)):
          ri = test_data[i]
          x1.append(ri[0])
          y1.append(ri[1])
  
      import matplotlib.pyplot as plt
      plt.plot(x1, y1, linestyle=' ', marker='.')
      xx, yy = [], []
      for i in range(len(result)):
          ri = result[i]
          xx.append(ri[0])
          yy.append(ri[1])
  
      plt.plot(xx, yy, linestyle=' ', marker='*')