基于时空轨迹的可达区域计算（论文待发表）

import psycopg2
import geopandas as gpd
import numpy as np
from fitter import Fitter
import matplotlib.pyplot as plt
from utils_cls import RRP, RPA, Trj, Road, CRoad

pgisCon = psycopg2.connect(database="beijing", user="jiangshan", password="jiangshan", host="localhost", port="5432")
pgisCursor = pgisCon.cursor()

# 时空查询点
point = (116.306251, 39.98070)# 新中关商业区
point_r = 2000.0
time_interal = ['2012-10-24 10:00:00', '2012-10-24 10:10:00']
EPS = 0.001

# 创建视图
sql = "CREATE VIEW st_trjs AS " \
    "with taxi as ("\
    "select tb.unit_id, tb.ts_s, tb.the_geom "\
    "from taxi2012_bj as tb "\
    "where ST_DWithin(ST_Transform(tb.the_geom, 2436),ST_Transform(st_setsrid(st_makePoint({0}, {1}),4326),2436), {2})) " \
"select taxi.unit_id as u_id, taxi.ts_s as tstamp, taxi.the_geom as the_geom "\
"from taxi "\
"WHERE \'{3}\' < taxi.ts_s and taxi.ts_s < \'{4}\'".format(point[0], point[1], point_r, time_interal[0], time_interal[1])

print(sql)
pgisCursor.execute(sql)
pgisCon.commit()
# 基于时空查询视图 目标轨迹
sql = 'SELECT st_trjs.u_id, st_trjs.tstamp,st_trjs.the_geom FROM st_trjs order by st_trjs.tstamp'
print(sql)
g_out2 = gpd.read_postgis(sql=sql, con=pgisCon, geom_col="the_geom")
trj_id_list = np.unique(g_out2['u_id'].values).tolist()
trj_list = []
road_id_list = []
road_list = []
RRP_list = []
for ui_d in trj_id_list:
    trj = Trj()
    rd = Road()
    rp_i = RRP()
    fugai_road = {}
    # trj_rp = {}
    # 基于时空查询视图 轨迹可达点（指定轨迹的几何中心）
    sql = 'select ST_Transform(ST_Centroid(ST_Collect(ST_Transform(st_trjs.the_geom,2436))), 4326) as the_ctgeom from st_trjs where st_trjs.u_id = \'{0}\''.format(ui_d)
    print(sql)
    g_out3 = gpd.read_postgis(sql=sql, con=pgisCon, geom_col="the_ctgeom")
    ct = g_out3['the_ctgeom'].values[0]
    ct_pointstr = '{}'.format(ct)# 'POINT (116.3012000856913 39.98411309689593)'
    lon_ct = ct.x
    lat_ct = ct.y

    # --查询一个点最近(500.5米)的一条道路 --- 可达路段(bj_ways_car表)
    sql = 'SELECT ways.gid, ways.the_geom, ' \
          'ST_Distance(ST_Transform(st_setsrid(st_makePoint({0}, {1}),4326),2436), ST_Transform(ways.the_geom,2436)) as dis, ' \
          'degrees(ST_azimuth(ST_Transform(st_setsrid(st_makePoint({2}, {3}), 4326), 2436),' \
          'ST_Transform(st_setsrid(ST_StartPoint(ways.the_geom), 4326), 2436))) as h_deg_s,' \
          'degrees(ST_azimuth(ST_Transform(st_setsrid(st_makePoint({2}, {3}), 4326), 2436),' \
          'ST_Transform(st_setsrid(ST_EndPoint(ways.the_geom), 4326), 2436))) as h_deg_e, ' \
          'ST_Distance(ST_Transform(st_setsrid(st_makePoint({2}, {3}),4326),2436), ST_Transform(ways.the_geom,2436)) as dis_q '\
          'FROM bj_ways_car as ways ' \
          'WHERE st_dwithin(ST_Transform(ways.the_geom, 2436), ST_Transform(st_setsrid(st_makePoint({0}, {1}),4326),2436),500.5)=true ' \
          'order by dis ' \
          'limit 1'.format(lon_ct, lat_ct, point[0], point[1])
    print(sql)
    g_out4 = gpd.read_postgis(sql=sql, con=pgisCon, geom_col="the_geom")
    # print("g_out4[\"gid\"].values: {0}".format(g_out4["gid"].values.shape))
    road_id = g_out4["gid"].values[0]
    road_dis = g_out4["dis_q"].values[0]
    road_h_deg_s = g_out4["h_deg_s"].values[0]
    road_h_deg_e = g_out4["h_deg_e"].values[0]
    road_line = g_out4["the_geom"].values[0]
    linestring = '{}'.format(road_line)  # 'LINESTRING (116.3012722 39.9841254, 116.3006655 39.9841259, 116.3005954 39.9841339, 116.3005336 39.9841491, 116.300455 39.9841788, 116.3003859 39.9842138, 116.3003198 39.9842554, 116.3002784 39.984309, 116.3002558 39.9843802, 116.3002544 39.9844514, 116.3002732 39.9845311, 116.3003132 39.9846039, 116.3003756 39.9846472, 116.3004554 39.9846892, 116.3005528 39.9847152, 116.3006465 39.9847316, 116.3007301 39.9847408, 116.3008259 39.9847431, 116.3017861 39.9847709, 116.3019589 39.9847396, 116.3020805 39.9846675, 116.3021492 39.9845716, 116.302164 39.9844806, 116.3021542 39.9841369, 116.3021107 39.9826122)'
    if road_id not in road_id_list:
        rd.road_id = road_id
        rd.h_deg_e = 360+90-road_h_deg_e
        rd.h_deg_s = 360+90-road_h_deg_s
        rd.dis = road_dis
        rd.linestring = linestring
        road_list.append(rd)
    # h_deg
    sql = 'select degrees(ST_azimuth(ST_Transform(st_setsrid(st_makePoint({0}, {1}),4326),2436), ' \
          'ST_Transform(st_setsrid(st_makePoint({2}, {3}),4326),2436))) AS h_deg, st_setsrid(st_makePoint({0}, {1}),4326) as geom'.format(point[0], point[1], lon_ct, lat_ct)
    print(sql)
    g_out5 = gpd.read_postgis(sql=sql, con=pgisCon, geom_col="geom")
    h_deg = g_out5["h_deg"].values[0]
    trj.ct_pointstr = ct_pointstr
    trj.h_deg = 360+90-h_deg
    trj.lat_ct = lat_ct
    trj.lon_ct = lon_ct
    trj.road_id = road_id
    trj.trj_num = 1.0
    trj.ui_d = ui_d
    trj_list.append(trj)
    road_id_list.append(road_id)
    rp_i.h_deg = 360+90-h_deg
    rp_i.num = 1.0
    rp_i.point = ct_pointstr
    RRP_list.append(rp_i)

# road table  道路段
road_list
# trj table  目标轨迹
trj_list
# road_list = road_list*10 # test
# trj_list = trj_list * 10 # test
# cr table  覆盖路段
CR = []# # 覆盖路段
for rd in road_list:
    cr = CRoad()
    cr.road_id = rd.road_id
    cr.linestring = rd.linestring
    cr.h_deg_s = rd.h_deg_s
    cr.h_deg_e = rd.h_deg_e
    cr.dis = rd.dis
    for trj in trj_list:
        if trj.road_id == rd.road_id:
            rp = RRP()
            rp.point = trj.ct_pointstr
            rp.num = trj.trj_num
            rp.ui_d = trj.ui_d
            rp.road_id = trj.road_id
            rp.lon_ct = trj.lon_ct
            rp.lat_ct = trj.lat_ct
            cr.rrp.append(rp)
            rpa = RPA()
            rpa.road_id = trj.road_id
            rpa.ui_d = trj.ui_d
            rpa.h_deg = trj.h_deg
            cr.rpa.append(rpa)
    cr.update_rrp_num()
    cr.updata_cra()
    CR.append(cr)

def get_crmax(cr_list):
    num_max = 1.0
    crm_ind = 0
    i = 0
    for cr in cr_list:
        if num_max < cr.rrp_num:
            num_max = cr.rrp_num
            crm_ind = i
        i += 1
    return crm_ind, num_max
crm_index, num_max = get_crmax(CR)

def softmax(x):
    x_row_max = x.max(axis=-1)
    x_row_max = x_row_max.reshape(list(x.shape)[:-1]+[1])
    x = x - x_row_max
    x_exp = np.exp(x)
    x_exp_row_sum = x_exp.sum(axis=-1).reshape(list(x.shape)[:-1]+[1])
    softmax = x_exp / x_exp_row_sum
    return softmax

dis_list = np.array([cr.dis for cr in CR])
print(dis_list)
rpnum_list = np.array([cr.rrp_num for cr in CR])
print(rpnum_list)

#########################################################################
# 利用fitter拟合数据样本的分布
f = Fitter(dis_list, distributions=['norm', 't', 'laplace', 'rayleigh'])
f.fit()
f.hist() #绘制组数=bins的标准化直方图
plt.show()
print("==="*10)
f = Fitter(rpnum_list, distributions=['norm', 't', 'laplace', 'rayleigh'])
f.fit()
f.hist() #绘制组数=bins的标准化直方图
plt.show()
#########################################################################

dis_soft = softmax(dis_list)
rpnum_soft = softmax(rpnum_list)
result_soft = dis_soft * rpnum_soft

cr_result = []# 最终结果
eps = EPS# 超参数
i_ind = 0
for rs in result_soft:
    if rs > eps:
        cr_result.append(CR[i_ind])
        print(CR[i_ind].road_id)
        print(CR[i_ind].rrp_num)
        print(CR[i_ind].linestring)
    i_ind += 1
print(cr_result)

# 之前有视图，删除视图
sql = 'drop view st_trjs'
pgisCursor.execute(sql)
pgisCon.commit()

pgisCursor.close()
pgisCon.close()