从GPS离散点统计活动圈子

背景

数仓里面有一些IP和GPS共现的数据，希望从这些数据中汇总出IP的大致活动范围

方案

将IP的离散分布点聚合后，排除离群点，再对非离群点求得最小圆覆盖，考虑到同一个IP可能会在多个地方都有GPS分布，因此可能会有多个聚集点，在求最小圆覆盖前需要多GPS点分类，再对不同的类求最小覆盖圆。 (也可以不分类，直接求一个圆)

此外在计算最小覆盖圆时需要将GPS坐标投影到平面坐标上去，来计算距离（python中可以用pyproj来进行坐标转换，但是我在spark中使用时，出现了一些问题：pyspark中使用的pyproj是1.9,导致无法从wgs84 转到Web 墨卡托坐标，因此代码中我写了手动转的方法，不过会有一点点误差）。

如果需要聚成多个圈子的话，在对GPS分类前，需要计算分为几类，可以先用层次聚类求得分为几个类，再用Kmeans来进行聚类。

最小覆盖圆的算法网上可以搜索到一些，本文参考了
https://github.com/platinhom/platinhom.github.com/blob/master/_posts/2015-11-11-MinCircleforPoints.md

代码

# -*- coding: utf-8 -*-
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import scipy.cluster.hierarchy as sch  # 用于进行层次聚类，画层次聚类图的工具包
import scipy.spatial.distance as ssd
from scipy.cluster.vq import vq, kmeans, whiten
import numpy as np
from matplotlib.patches import Ellipse, Circle
import math

# 长半轴
a = 6378245.0;
# 扁率
ee = 0.00669342162296594323;


# 纬度转换
def transformlat(lng, lat):
    ret = -100.0 + 2.0 * lng + 3.0 * lat + 0.2 * lat * lat + 0.1 * lng * lat + 0.2 * math.sqrt(abs(lng))
    ret += (20.0 * math.sin(6.0 * lng * math.pi) + 20.0 * math.sin(2.0 * lng * math.pi)) * 2.0 / 3.0
    ret += (20.0 * math.sin(lat * math.pi) + 40.0 * math.sin(lat / 3.0 * math.pi)) * 2.0 / 3.0
    ret += (160.0 * math.sin(lat / 12.0 * math.pi) + 320 * math.sin(lat * math.pi / 30.0)) * 2.0 / 3.0
    return ret;


def transformlng(lng, lat):
    ret = 300.0 + lng + 2.0 * lat + 0.1 * lng * lng + 0.1 * lng * lat + 0.1 * math.sqrt(abs(lng));
    ret += (20.0 * math.sin(6.0 * lng * math.pi) + 20.0 * math.sin(2.0 * lng * math.pi)) * 2.0 / 3.0;
    ret += (20.0 * math.sin(lng * math.pi) + 40.0 * math.sin(lng / 3.0 * math.pi)) * 2.0 / 3.0;
    ret += (150.0 * math.sin(lng / 12.0 * math.pi) + 300.0 * math.sin(lng / 30.0 * math.pi)) * 2.0 / 3.0;
    return ret;


def wgs84togcj02(lng, lat):
    dlat = transformlat(lng - 105.0, lat - 35.0)
    dlng = transformlng(lng - 105.0, lat - 35.0)
    radlat = lat / 180.0 * math.pi;
    magic = math.sin(radlat);
    magic = 1 - ee * magic * magic;
    sqrtmagic = math.sqrt(magic);
    dlat = (dlat * 180.0) / ((a * (1 - ee)) / (magic * sqrtmagic) * math.pi);
    dlng = (dlng * 180.0) / (a / sqrtmagic * math.cos(radlat) * math.pi);
    mglat = lat + dlat;
    mglng = lng + dlng;
    return mglng, mglat;


def gcj02towgs84(lng, lat):
    dlat = transformlat(lng - 105.0, lat - 35.0);
    dlng = transformlng(lng - 105.0, lat - 35.0);
    radlat = lat / 180.0 * math.pi;
    magic = math.sin(radlat);
    magic = 1 - ee * magic * magic;
    sqrtmagic = math.sqrt(magic);
    dlat = (dlat * 180.0) / ((a * (1 - ee)) / (magic * sqrtmagic) * math.pi);
    dlng = (dlng * 180.0) / (a / sqrtmagic * math.cos(radlat) * math.pi);
    mglat = lat + dlat;
    mglng = lng + dlng;
    return lng * 2 - mglng, lat * 2 - mglat;


# 经纬度转Wev墨卡托
def lonLat2WebMercator(lng, lat):
    lng, lat = gcj02towgs84(lng, lat)
    x = lng * 20037508.34 / 180;
    y = math.log(math.tan((90 + lat) * math.pi / 360)) / (math.pi / 180);
    y = y * 20037508.34 / 180
    return [x, y]


# Web墨卡托转经纬度
def WebMercator2lonLat(x, y):
    lng = x / 20037508.34 * 180;
    lat = y / 20037508.34 * 180;
    lat = 180 / math.pi * (2 * math.atan(math.exp(lat * math.pi / 180)) - math.pi / 2);
    lng, lat = wgs84togcj02(lng, lat)
    return [lng, lat]


def lonLatArray2WebMercator(lng_lat_array):
    array = []
    for (lng, lat) in lng_lat_array:
        array.append(lonLat2WebMercator(lng, lat))
    return array


# 离群点删除
def detect_outliers(data, threshold=3):
    mean_d = np.mean(data[:, 0])
    std_d = np.std(data[:, 0])

    mean_d_1 = np.mean(data[:, 1])
    std_d_1 = np.std(data[:, 1])

    data_filter = []
    for y in data:
        if std_d == 0 and std_d_1 == 0:
            data_filter.append(True)
        else:
            z_score = (y[0] - mean_d) / std_d
            z_score_1 = (y[1] - mean_d_1) / std_d_1
            if (std_d == 0 or abs(z_score) <= threshold) and (abs(z_score_1) <= threshold or std_d_1 == 0):
                data_