【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)<