具体需求

基于电信的位置融合数据，分析用户是否通过火车出行，以及火车出行的相关信息，如乘车车次、上车站、下车站等。

数据描述

1、用户定位数据手机号码业务开始时间城市编码经纬度基站标识网格号业务类型事件类型数据源天分区2、火车停靠数据车次车站序号车站名称行驶时间到站时间行驶里程经纬度

Python代码如下：

输入经纬度 和编码级别计算对应的geohash编码

from math import log10__base32 = '0123456789bcdefghjkmnpqrstuvwxyz'__decodemap = { }for i in range(len(__base32)):    __decodemap[__base32[i]] = idel idef decode_exactly(geohash):    lat_interval, lon_interval = (-90.0, 90.0), (-180.0, 180.0)    lat_err, lon_err = 90.0, 180.0    is_even = True    for c in geohash:        cd = __decodemap[c]        for mask in [16, 8, 4, 2, 1]:            if is_even: # adds longitude info                lon_err /= 2                if cd & mask:                    lon_interval = ((lon_interval[0]+lon_interval[1])/2, lon_interval[1])                else:                    lon_interval = (lon_interval[0], (lon_interval[0]+lon_interval[1])/2)            else:      # adds latitude info                lat_err /= 2                if cd & mask:                    lat_interval = ((lat_interval[0]+lat_interval[1])/2, lat_interval[1])                else:                    lat_interval = (lat_interval[0], (lat_interval[0]+lat_interval[1])/2)            is_even = not is_even    lat = (lat_interval[0] + lat_interval[1]) / 2    lon = (lon_interval[0] + lon_interval[1]) / 2    return lat, lon, lat_err, lon_errdef decode(geohash):    lat, lon, lat_err, lon_err = decode_exactly(geohash)    # Format to the number of decimals that are known    lats = "%.*f" % (max(1, int(round(-log10(lat_err)))) - 1, lat)    lons = "%.*f" % (max(1, int(round(-log10(lon_err)))) - 1, lon)    if '.' in lats: lats = lats.rstrip('0')    if '.' in lons: lons = lons.rstrip('0')    return lats, lonsdef encode(latitude, longitude, precision=12):    """    Encode a position given in float arguments latitude, longitude to    a geohash which will have the character count precision.    """    lat_interval, lon_interval = (-90.0, 90.0), (-180.0, 180.0)    geohash = []    bits = [ 16, 8, 4, 2, 1 ]    bit = 0    ch = 0    even = True    while len(geohash) < precision:        if even:            mid = (lon_interval[0] + lon_interval[1]) / 2            if longitude > mid:                ch |= bits[bit]                lon_interval = (mid, lon_interval[1])            else:                lon_interval = (lon_interval[0], mid)        else:            mid = (lat_interval[0] + lat_interval[1]) / 2            if latitude > mid:                ch |= bits[bit]                lat_interval = (mid, lat_interval[1])            else:                lat_interval = (lat_interval[0], mid)        even = not even        if bit < 4:            bit += 1        else:            geohash += __base32[ch]            bit = 0            ch = 0    return ''.join(geohash)

上面的decode是解码算法 encode是编码算法

Google S2地图索引算法

编码过程和解码过程如下：

hilbert_map = {    'a': {(0, 0): (0, 'd'), (0, 1): (1, 'a'), (1, 0): (3, 'b'), (1, 1): (2, 'a')},    'b': {(0, 0): (2, 'b'), (0, 1): (1, 'b'), (1, 0): (3, 'a'), (1, 1): (0, 'c')},    'c': {(0, 0): (2, 'c'), (0, 1): (3, 'd'), (1, 0): (1, 'c'), (1, 1): (0, 'b')},    'd': {(0, 0): (0, 'a'), (0, 1): (3, 'c'), (1, 0): (1, 'd'), (1, 1): (2, 'd')},}un_hilbert_map = {    'a': { 0: (0, 0,'d'), 1: (0, 1,'a'), 3: (1, 0,'b'),  2: (1, 1,'a')},    'b': { 2: (0, 0,'b'), 1: (0, 1,'b'), 3: (1, 0,'a'),  0: (1, 1,'c')},    'c': { 2: (0, 0,'c'), 3: (0, 1,'d'), 1: (1, 0,'c'),  0: (1, 1,'b')},    'd': { 0: (0, 0,'a'), 3: (0, 1,'c'), 1: (1, 0,'d'),  2: (1, 1,'d')}}#编码def point_to_hilbert(lng,lat, order=16):    print ('hilbert')    lng_range = [-180.0, 180.0]    lat_range = [-90.0, 90.0]    current_square = 'a'    position = 0    for i in range(order - 1, -1, -1):        position <<= 2        lng_mid = (lng_range[0]+lng_range[1])/2        lat_mid = (lat_range[0]+lat_range[1])/2        if lng >= lng_mid :            quad_x = 1            lng_range[0] = lng_mid        else:            quad_x = 0            lng_range[1] = lng_mid        if lat >= lat_mid :            quad_y = 1            lat_range[0] = lat_mid        else:            quad_y = 0            lat_range[1] = lat_mid        quad_position,current_square = hilbert_map[current_square][(quad_x, quad_y)]        position |= quad_position    return position#解码def hilbert_to_point( d , order=16):    print ('hilbert')    lng_range = [-180.0, 180.0]    lat_range = [-90.0, 90.0]    current_square = 'a'    lng=lat=lng_mid=lat_mid=0    for i in range(order - 1, -1, -1):        lng_mid = ( lng_range[0] + lng_range[1] ) / 2        lat_mid = ( lat_range[0] + lat_range[1] ) / 2        mask = 3 << (2*i)        quad_position = (d & mask) >> (2*i)        quad_x, quad_y, current_square= un_hilbert_map[current_square][quad_position]        if quad_x:            lng_range[0] = lng_mid        else:            lng_range[1] = lng_mid        if quad_y:            lat_range[0] = lat_mid        else:            lat_range[1] = lat_mid    lat = lat_range[0]    lng = lng_range[0]    return lng,latif __name__ == '__main__':    d = point_to_hilbert(-50.555443,77.776655,36)    print (d)    lng,lat = hilbert_to_point(d,36)    print (lng,lat)

获取百度路径规划算法结果的代码如下：

对应的百度地图API参考代码中的链接

import loggingimport sysimport tracebackfrom time import sleepimport jsonimport requestsfrom bs4 import BeautifulSoupimport pymysqlimport os.pathimport timedef getroadinfo(par):    linearray = par.split('#&;')    startlon = linearray[0]    startlat = linearray[1]    endlon = linearray[2]    endlat = linearray[3]    strcount = linearray[4]    #http: // api.map.baidu.com / direction / v1?mode = driving & origin = 上地五街 & destination = 北京大学 & origin_region = 北京 & destination_region = 北京 & output = json & ak = 您的ak    origin = str(startlat) + "," + str(startlon)    dest = str(endlat) + "," + str(endlon)    url = "http://api.map.baidu.com/direction/v1";    Param = "mode=driving&origin=" + origin + "&destination=" + dest + "&origin_region=武汉&destination_region=武汉&output=json&ak=" + "";    strurl = url + "?" + Param;    #请求百度接口返回数据    try:        web = requests.get(strurl, timeout=30)        data = json.loads(web.text)        if 'result' in data:            if 'routes' in data['result']:                routes = data['result']['routes']                listtext = []                writestr=""                for route in routes:                    if 'steps' in route:                        steps = route['steps']                        for step in steps:                            # area = 0                            # if 'area' in step:                            #     area = step['area']                            # direction = 0                            # if 'direction' in step:                            #     direction = step['direction']                            # distance = 0                            # if 'distance' in step:                            #     distance = step['distance']                            # duration = 0                            # if 'duration' in step:                            #     duration = step['duration']                            # instructions = ""                            # if 'instructions' in step:                            #     instructions = step['instructions']                            path = ""                            if 'path' in step:                                path = step['path']                            # turn = 0                            # if 'turn' in step:                            #     turn = step['turn']                            # type = 0                            # if 'type' in step:                            #     type = step['type']                            # stepOriginInstruction = ""                            # if 'stepOriginInstruction' in step:                            #     stepOriginInstruction = step['stepOriginInstruction']                            # stepDestinationInstruction = ""                            # if 'stepDestinationInstruction' in step:                            #     stepDestinationInstruction = step['stepDestinationInstruction']                            # traffic_condition = 0                            # if 'traffic_condition' in step:                            #     traffic_condition = step['traffic_condition']                            # currenttime = time.strftime('%Y-%m-%d %H:%M:%S',time.localtime(time.time()))                            # writestr = origin + "#&;" + dest + "#&;" + str(area) +"#&;"+ str(direction) +"#&;"+ str(distance) +"#&;"+ str(duration) +"#&;"+ instructions +"#&;"+ path +"#&;"+ str(turn) +"#&;"+ str(type) +"#&;"+ stepOriginInstruction +"#&;"+ stepDestinationInstruction +"#&;"+ str(traffic_condition) +"#&;"+ currenttime + "\n"                            writestr = writestr + path + ";"                print(time.strftime('%Y-%m-%d %H:%M:%S', time.localtime(time.time())))                print(writestr)                writestr = writestr + "\n"                listtext.append(writestr)                yushu = int(strcount)%20                text = "E:/roadinfo/roadinfo"+str(yushu)+".txt"                file_object = open(text, 'a')                lock(file_object, LOCK_EX)                file_object.writelines(listtext)                unlock(file_object)                file_object.close()    except:        logging.error('解析线路信息失败')        traceback.print_exc()

计算步骤

为了计算所有用户可能乘坐的火车班次和上下车车站具体的火车站，具体计算过程如下：方案1

1、通过火车班次数据中的第三列和第七列

计算出所有的火车运行时速

2、对所有的用户定位经纬度进行geohash编码或者S2编码存成另外一个列

3、参考geohash编码级别对应的定位精度表如下：Geohash 字符串长度 纬度 经度 纬度误差 经度误差 km误差1 2 3 ±23 ±23 ±25002 5 5 ±2.8 ±5.6 ±6303 7 8 ±0.70 ±0.70 ±784 10 10 ±0.087 ±0.18 ±205 12 13 ±0.022 ±0.022 ±2.46 15 15 ±0.0027 ±0.0055 ±0.617 17 18 ±0.00068 ±0.00068 ±0.0768 20 20 ±0.000085 ±0.00017 ±0.019

如果用户在某一个小时的时间间隔内 geohash的前三位有变动，我们可以认为该用户在高速移动中

4、筛选出所有geohash在一个小时内前三位有变动的点的用户数据，5、对所有的火车站进行geohash第五级编码6、如果某个用户被过滤后的数据的geohash编码和某个车站的geohash编码的前五位相同，则判断该用户的定位时间是否在该geohash的车站的到站时间前后半个小时之内，如果在半个小时之内，则判断该用户前后的车站是否在某一个火车形成序列中，如果是则输出。

方案2

1、因为考虑到用户的数据并不是以一定的频率匀速上传，中间有间断的可能，所以只能采用路网匹配的算法，即判断某个用户某个时间的数据点是否在道路网上：2、参考百度地图API，获取百度地图的铁路网路径规划信息，具体爬虫程序参考我上面的获取路径规划的程序。以火车时刻表的每个班车每个时间段的起点和终点（车站）作为输入获取道路网数据。3、对所有铁路道路网点进行geohash第五级编码 也就是定位精度为2.4公里4、对所有用户数据的经纬度进行geohash第五级编码5、对用户时间进行格式转化为小时和分钟两个列6、对道路网数据进行时间补充操作例如：

A B C

如上图A点经过B点到达C点 总共时间一个小时，A到B路程为四分之一，已知A，C时间 则B的时间点为A时间加上一刻钟。

之后将所有道路网时间转化为 小时和分钟两个列7、对用户数据表和道路网表的时间数据列和 geohash列添加索引注意是前缀索引，可以另外加两个列是geohash的前缀，通过前缀匹配会加速。

8、使用铁路网数据的geohash列对用户定位数据进行过滤，去除所有不在火车道路网附近的用户数据

9、对每一组用户数据，和道路网的数据进行 geohash列的join 10、判断时间是否和对应的时间点吻合，吻合即输出